2

I am working with a large(ish) shapefile containing 30,000 features. I am trying to run a zonal histogram and have been running into many problems. Running this function with the entire set of features crashes the tool. I've tried rasterizing the feature by FID but this also gets hung up and never finishes. I've written some python code that runs a zonal histogram on each feature but this has not really been finishing and is taking a long time.

I would now like to try running the tool by taking groups of 100 features and running the zonal histogram on each of those groups, outputting a uniquely named DBF file. I can then process all the DBF files in R to create one final table.

Below is my code for running the tool feature by feature. I found writing a temporary shapefile is more efficient than creating a temporary lyr file so I'd like to keep this approach. How could I adapt the code below to run in groups of 100 features?

import arcpy
import os
import os.path
from arcpy import env
from arcpy.sa import *
import gc

rasterLayer = "D:/pathToMyRaster/rasterLayer.tif"
outputDirectory = "D:/myOutputPath/"

newlayer = arcpy.mapping.Layer('D:/pathToShapefile/parcels.shp')

arcpy.env.parallelProcessingFactor = "100%"

with arcpy.da.SearchCursor(newlayer,['OID@','statsZone']) as cursor:
    arcpy.env.addOutputsToMap = False
    for row in cursor:
        print(row[1])
        tempFileName = outputDirectory+row[1]+".dbf"        
        sql="""{0} = {1}""".format(arcpy.AddFieldDelimiters(newlayer, arcpy.Describe(
            newlayer).OIDFieldName),row[0])
        arcpy.Select_analysis(in_features=newlayer, out_feature_class=os.path.join(outputDirectory,'TempShapefile.shp'.format(row[0])),
                             where_clause=sql)
        try:
            arcpy.gp.ZonalHistogram_sa(os.path.join(outputDirectory,'TempShapefile.shp'.format(row[0])), 'statsZone', rasterLayer, outputDirectory+row[1]+".dbf", "")
            del row
            gc.collect()
        except:
            print(str(errors))
            del row
            gc.collect()
  • I would first make sure my geometries are valid by running RepairGeometry tool since it is weird that even raster conversion crashes. – fatih_dur Jun 9 '18 at 7:29
  • I think the issues with crashing come from the number of columns generated in the output. Using a raster of parcels with OID as the value would generate a dbf of 27,000 columns. ArcMap fails and cannot handle this as an output – jotamon Jun 11 '18 at 15:54
2

I'd create a new short field in original shapefile, e.g. GROUP and populate it using:

!FID!//100

and create temporary subset using:

for i in range(30):
    query='"GROUP" = %s' %i
    arcpy.Select_analysis("ORIGINAL", r"in_memory/extract", query)
2

You can use a generator to yield chunks of 100 items at a time from a list of all rows. Then use each chunk to build sql query using the IN operator:

...
def chunks(l, n):
    """Yield successive n-sized chunks from l."""
    for i in range(0, len(l), n):
        yield l[i:i + n]

all_rows = [i for i in arcpy.da.SearchCursor(newlayer,['OID@','statsZone'])]

for p in chunks(all_rows,100):
    sql="""{0} IN({1})""".format(arcpy.AddFieldDelimiters(newlayer, arcpy.Describe(
                    newlayer).OIDFieldName),','.join([str(i) for i in [j[0] for j in p]]))
    arcpy.Select_analysis(..., where_clause=sql)
...
0

With the extra information given in your comment and presuming you are OK with using rasterised parcels, I think what you are looking for is Combine tool followed by Summary Statistics.

Combine will replicate what intersect tool does with two rasters and summary statistics is to dissolve your parcel OIDs (case field) as COUNTing the non-parcel raster (values). This output can be exported to a spreadsheet format to draw the histogram.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.