Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have several projects with tree point features (SpecialStatusSpeciesPoint) that contain up to 2,000,000 records. I am trying to record tree height (Hgt) and base elevation (BaseElev) values, based on DEMs representing vegetation height (veg_HGT) and bare earth (BEm). I am using ArcMap 10.1 with Spatial Analyst. Features are in a File GDB. DEMs are initially IMGs, but moved into a scratch file gdb as a GRID.

I initially tried Sample but was left with inconsistently null values in the attribute table. The same was true for Extract Multi Values to Points as well. Ultimately, I found Zonal Statistics as Table to work, but at a painstaking slow speed. Sometimes it would crash or not complete all records. I thought if I broke the points up into chunks, it would ensure all records were copied and perhaps in less time.

I wrote the following code (with the aid of snippets from you guys) as part of a larger script. All necessary modules and variables are taken care of in the code. For the sake of brevity, I only included the following. The processing time for only 10,000 records is on average 55 mins, 83 mins for 100,000 records, and 88 mins for 250,000 records (the code below).

# Define workspaces
arcpy.env.workspace = sdsfie_db
arcpy.env.scratchWorkspace = scratchGDB
arcpy.env.overwriteOutput = True
print("Adding attributes to SpecialStatusSpeciesPoint...")

# Add Field Id
if (not FieldExist(SpecialStatusSpeciesPoint, "Id")):
    arcpy.AddField_management(SpecialStatusSpeciesPoint, "Id", "LONG", "", "", "", "Id", "NULLABLE", "NON_REQUIRED", "")

# Use updateCursor to populate Id field
print("      Copying Id field...")
rows = arcpy.UpdateCursor(SpecialStatusSpeciesPoint)
for row in rows:
    OID = row.getValue("OBJECTID")
    row.setValue("Id", OID)
    rows.updateRow(row)
del row, rows

# Count number of records
recCount = int(arcpy.GetCount_management(SpecialStatusSpeciesPoint).getOutput(0))

# Decide whether to break data up into chunks
if recCount > 250000:
    # Break data into chunks
    print("   Breaking data into chunks of 250,000...")
    inDesc = arcpy.Describe(SpecialStatusSpeciesPoint)
    oidName = arcpy.AddFieldDelimiters(SpecialStatusSpeciesPoint, inDesc.oidFieldName)
    sql = '%s = (select min(%s) from %s)' % (oidName,oidName,os.path.basename(SpecialStatusSpeciesPoint))
    cur = arcpy.da.SearchCursor(SpecialStatusSpeciesPoint,[inDesc.oidFieldName],sql)
    minOID = cur.next()[0]
    del cur, sql
    sql = '%s = (select max(%s) from %s)' % (oidName,oidName,os.path.basename(SpecialStatusSpeciesPoint))
    cur = arcpy.da.SearchCursor(SpecialStatusSpeciesPoint,[inDesc.oidFieldName],sql)
    maxOID = cur.next()[0]
    del cur, sql
    breaks = range(minOID,maxOID)[0:-1:250000] # makes slices of 250000
    breaks.append(maxOID+1)
    ranges = []
    for b in range(len(breaks)-1):
        ranges.append([breaks[b], breaks[b+1]-1])
    numChunks = str(len(ranges))
    print("      Broken into " + numChunks + " chunks...")

    # Begin timer
    timerStart = time.time()

    tempOut = scratchGDB + "\\tempOut"
    counter = 0

    # Iterate through chunks
    for i,j in ranges:
        counter += 1
        lapStart = time.time()
        print("   Processing chunk " + str(counter) + " of " + numChunks + ". ID " + str(i) + " - " + str(j) + "...")

        # Build sql where_clause independent of GDB type
        fieldDelimited = arcpy.AddFieldDelimiters(SpecialStatusSpeciesPoint, "OBJECTID")
        fieldType = arcpy.ListFields(SpecialStatusSpeciesPoint, "OBJECTID")[0].type
        if str(fieldType) == 'String':
            i = "'%s'" % i
            j = "'%s'" % j
        whereClause = "%s BETWEEN %s AND %s" % (fieldDelimited, i, j)

        arcpy.env.workspace = sdsfie_db
        arcpy.env.scratchWorkspace = scratchGDB
        arcpy.env.overwriteOutput = True

        # Create layer of chunk
        print("      Creating temp layer...")
        # Create temp layer slice
        arcpy.Select_analysis(SpecialStatusSpeciesPoint, tempOut, whereClause)
        # Count features in temp layer
        tempOutCount = arcpy.GetCount_management(tempOut).getOutput(0)
        # Zonal Statistics as Table
        print("      Creating zoneTable7...")
        arcpy.gp.ZonalStatisticsAsTable_sa(tempOut, "Id", veg_HGT, zoneTable7, "DATA", "MAXIMUM")
        # Join Field
        print("      Joining MAX field...")
        arcpy.JoinField_management(SpecialStatusSpeciesPoint, "Id", zoneTable7, "Id", "MAX")
        # Zonal Statistics as Table
        print("      Creating zoneTable8...")
        arcpy.gp.ZonalStatisticsAsTable_sa(tempOut, "Id", BEm, zoneTable8, "DATA", "MINIMUM")
        # Join Field
        print("      Joining MIN field...")
        arcpy.JoinField_management(SpecialStatusSpeciesPoint, "Id", zoneTable8, "Id", "MIN")
        # Use updateCursor instead of calculate field
        print("      Copying MAX and MIN fields...")
        rows = arcpy.UpdateCursor(SpecialStatusSpeciesPoint, whereClause)
        for row in rows:
            veg_HGT_max = row.getValue("MAX")
            BEm_min = row.getValue("MIN")
            row.setValue("Hgt", veg_HGT_max)
            row.setValue("BaseElev", BEm_min)
            rows.updateRow(row)
        del row, rows
        # Clean up for next iteration
        arcpy.DeleteField_management(SpecialStatusSpeciesPoint, "MAX;MIN")
        arcpy.Delete_management("in_memory")
        # Print lap time
        lapStop = time.time()
        lapTime = lapStop - lapStart
        print("   Chunk run time: " + str(lapTime/60) + " minutes...")
else:
    rows = arcpy.UpdateCursor(SpecialStatusSpeciesPoint)
    for row in rows:
        OID = row.getValue("OBJECTID")
        row.setValue("Id", OID)
        rows.updateRow(row)
    del row, rows
    # Process: Zonal Statistics as Table
    arcpy.gp.ZonalStatisticsAsTable_sa(SpecialStatusSpeciesPoint, "Id", veg_HGT, zoneTable7, "DATA", "MAXIMUM")
    # Process: Join Field
    arcpy.JoinField_management(SpecialStatusSpeciesPoint, "Id", zoneTable7, "Id", "MAX")
    # Process: Zonal Statistics as Table
    arcpy.gp.ZonalStatisticsAsTable_sa(SpecialStatusSpeciesPoint, "Id", BEm, zoneTable8, "DATA", "MINIMUM")
    # Process: Join Field
    arcpy.JoinField_management(SpecialStatusSpeciesPoint, "Id", zoneTable8, "Id", "MIN")
    rows = arcpy.UpdateCursor(SpecialStatusSpeciesPoint)
    for row in rows:
        veg_HGT_max = row.getValue("MAX")
        BEm_min = row.getValue("MIN")
        row.setValue("Hgt", veg_HGT_max)
        row.setValue("BaseElev", BEm_min)
        rows.updateRow(row)
    del row, rows

# Clean up
arcpy.DeleteField_management(SpecialStatusSpeciesPoint, "Id;MAX;MIN")
arcpy.Delete_management("in_memory")
arcpy.Delete_management(scratchGDB)

# Stop timer and print results
timerStop = time.time()
elapsed = timerStop - timerStart
hrs, remain = divmod(elapsed, 3600)
mins, secs = divmod(remain, 60)
print("   SpecialStatusSpeciesPoint processing time: %d:%d:%f" % (hrs, mins, secs))

What might be the bottleneck? Is there something terribly wrong with the code? (I began learning Python 4 months ago) Has anyone who has had similar issues with Sample/Extract Multi Values to Points/zonal Statistics as Table found a more efficient workaround for null values or missing records? Thank you in advance for any help you can offer.

share|improve this question
    
To clarify, only this section, with other unrelated processes commented out, resulted in the above times for each chunk of data. –  Barbarossa Jun 19 '13 at 21:17
3  
Do all of your points fall within the DEM in all cases? I have found that in order to efficiently and accurately use Extract Values to Points I needed to first clip my points by the DEM, otherwise I risked errors or long processing times. That being said, I have extracted values for thousands of points without issue. –  Nick Ochoski Jun 19 '13 at 21:22
    
@Radar, yes all points fall within the DEM. The points are actually derived from DEMs with the exact same extents. I had processed as many as 200,000 features at once with no problem until my newest project with nearly 10x that many. –  Barbarossa Jun 20 '13 at 13:52
1  
My first thought (and this may be wrong) is you could be encountering a memory issue with that many points. The tool might be designed to extract the values using memory before writing everything to disk. Other than that, make sure you clean up any artifacts from your loops (cursors) that may be eating up memory. –  Nick Ochoski Jun 20 '13 at 15:52
    
Are you open to using R? You will likely get much better performance using R's extract command in the raster package. –  Aaron Apr 27 at 21:28
add comment

1 Answer 1

I find that Sample crashes when run in memory, EMVTP crashes to desktop when multiple rasters are specified, and that Sample saves incorrect null values if either raster contains any null cells.

The solution that works for me is to run EMVTP on one layer at a time. The python console makes this very easy to do:

layers =

band-select layers, drag and drop into console, hit enter,

for layer in layers:
  arcpy.sa.ExtractMultiValuesToPoints('RandomPoints', layer, 'NONE')

hit enter

Has been much easier and more foolproof than other methods. However, my problem is too many rasters, not too many points, so I don't know if it will help you.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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