Spatial Analyst appears to be inhibiting iterative analysis when run in For Loop

Question

I am attempting to iterate a process using spatial analyst tools which seems to be causing some type of data lock type error that I am struggling to work around. I am attempting to perform an analysis on building roofs (identified by the ‘BIN’ (Building ID Number) value referenced in the code) using point elevation readings captured from lidar. The data throughout the project is stored in a series file geodatabases. Each iteration in the for loop is supposed to copy the source lidar points (which have been previously converted to vector) and polyline outline of the building footprint into a working geodatabase, where several geoprocessing tools are run. Once the iteration is complete, the database is cleared out (currently deletes entire GDB, but that was implemented as a failed work around) when finished. The analysis always seems to work as intended on the first iteration, but on the second iteration all the data copied into the workspace GDB ends up being empty geometry. The empty geometry then causes any following geoprocessing to fail. This has occurred in both python 2 and 3 as well as in ArcMap, ArcPro and stand-alone python (where it’s intended to be run). I’ve tried the initial copy into the workspace GDB via feature class to feature class, making a layer selecting and ‘copy’ and ‘copy feature’, copying xy coordinates and building the elevation points from scratch (which took too long to use), as well as copying each BIN to it’s own feature class and just being pointed via a python list. I’ve tested the script by commenting out the spatial analysis and the workspace copies work as intended, leading me to think it’s something about the spatial analyst tools causing it. I’ve also built python processes of this scale in the past using strictly vector, and never encountered this before. I’ve checked the source GDB to make sure there are no open cursors, and have found none. Has anyone encountered something like this before? For context, this is the project I'm working on Here’s my code (sorry it’s kind of a mess): import arcpy import sys import traceback from arcpy.sa import * arcpy.CheckOutExtension("Spatial") outputGbd = 'F:\Lidar\lidar_2018\RearRoofTesting\lidartesting\TestFinal.gdb' #TODO tileOutlineAllFeature = 'F:\Lidar\lidar_2018\RearRoofTesting\lidartesting\BadHouse.gdb\Tile_26875E244534N' tileOutlineAll = 'tileOutlineAllLayer' tileOutlineCurrentFeature = 'tileOutlineCurrent' tileOutlineCurrent = 'tileOutlineCurrentLayer' #currentTileFeature = 'F:\Lidar\lidar_2018\RearRoofTesting\ExportForTesting2.gdb\Tile_26875E244534N' #TODO create this-->This is the lidar points for current tile #currentTile = 'currentTileLayer' #buildingFootprintsAllFeature = 'F:\Lidar\LiDARAnalysis\Footprints.gdb\LI_BUILDING_FOOTPRINTS' #TODO update this to footprints layer used for lidar buildingFootprintsAllFeature = 'F:\Lidar\lidar_2018\RearRoofTesting\lidartesting\BuildingFootprints.gdb\TestBuildingsWithVios' #TODO Delete after testing buildingFootprintsAll = 'buildingFootprintsAllLayer' buildingFootprintsFeature = 'buildingFootprintsCurrent' buildingPointsFeature = 'CurrentInputBuildingPoints' #buildingFootprintsFeature = 'BIN_1009464_Footprint' buildingFootprints = 'buildingFootprintsLayer' outputPoints = outputGbd + '\HoleFlags' #outputPoints = 'HoleFlags' arcpy.MakeFeatureLayer_management(tileOutlineAllFeature, tileOutlineAll) arcpy.MakeFeatureLayer_management(buildingFootprintsAllFeature, buildingFootprintsAll) tileResult = arcpy.GetCount_management(tileOutlineAllFeature) tileCount = int(tileResult.getOutput(0)) tileCursor = arcpy.da.SearchCursor(tileOutlineAll, 'FileName') def rasterizeRoof(currentBuilding, bin, binPoints, workspace): # Working Variables binPointsName = 'BIN_' + bin minB = binPointsName + '_MinB' minBRaster = minB + '_Rast' zonalStats = binPointsName + '_ZonalStats' roofZoneFeature = binPointsName + '_RoofZones' roofZones = 'RoofZonesLayer' roofZoneVertFeature = binPointsName + '_RoofZoneVerts' roofZoneVerts = 'roofZoneVertsLayer' roofZoneCurrent = binPointsName + '_RoofZoneCurrent' # Output Variables outputGbd = 'F:\Lidar\lidar_2018\RearRoofTesting\lidartesting\TestFinal.gdb' outputPoints = outputGbd + '\HoleFlags' # TODO Temp idw = Idw(binPoints, 'POINT_Z', '0.5', in_barrier_polyline_features="BIN_" + bin + "_Outline") arcpy.Delete_management("BIN_" + bin + "_Outline") # Create Slope Raster Of Roof print('Slope Analysis') slope = Slope(idw, 'PERCENT_RISE', '', 'GEODESIC') # Use Con to create a raster representing areas of the idw that do not have a slope of 90 and are within the minimum bounds. # This Raster will only have a value of 1 as a constant for the Region Group Prcoess print('Con Analysis') rawGroups = Con((slope < 90) & ~IsNull(minBRaster), 1) # Clean the Boundaries to eliminate any small cells print('Boundary Clean') rawGroupsClean = BoundaryClean(rawGroups, '', 'ONE_WAY') # Region group all isolated sections of roof elevation print('Region Group') groups = RegionGroup(rawGroupsClean, 'FOUR', 'WITHIN') # Zonal Statistics to find elevation stats for each group print('Zonal Statistics') ZonalStatisticsAsTable(groups, 'Value', idw, zonalStats) # Convert to vector polygons and delineate the "Large" sections (over 25 ft) print('Converting back to vector') arcpy.RasterToPolygon_conversion(groups, roofZoneFeature, 'SIMPLIFY', 'Value') arcpy.MakeFeatureLayer_management(roofZoneFeature, roofZones) arcpy.AddField_management(roofZones, 'Large_Section', 'SHORT') idLargeCursor = arcpy.da.UpdateCursor(roofZones, ['Shape_Area', 'Large_Section']) for row in idLargeCursor: if row[0] >= 25: row[1] = 1 idLargeCursor.updateRow(row) del idLargeCursor # Delete all raster related datasets arcpy.Delete_management(minB) arcpy.Delete_management(minBRaster) arcpy.Delete_management(idw) arcpy.Delete_management(slope) arcpy.Delete_management(rawGroups) arcpy.Delete_management(rawGroupsClean) arcpy.Delete_management(groups) # Pull height stats for roof zones into a dictionary print('Pulling stats') statCursor = arcpy.da.SearchCursor(zonalStats, ['Value', 'MIN', 'MAX', 'MEAN']) statDict = {row[0]: list(row[1:]) for row in statCursor} del statCursor # Create feature class and layer of roof zone vertices print('Feature Vertices to Points') arcpy.FeatureVerticesToPoints_management(roofZones, roofZoneVertFeature) print('Making Feature Layer') arcpy.MakeFeatureLayer_management(roofZoneVertFeature, roofZoneVerts) # Iterate through all "small" roof zones (under 25 Ft), these "small" zones are believed to represent roof holes hCursor = arcpy.da.SearchCursor(roofZones, ['gridcode', 'Large_Section']) for row in hCursor: print('Roof Section ' + str(row[0])) if row[1] is None: print('Analysing small section') currentMIN = statDict[row[0]][0] currentMAX = statDict[row[0]][1] currentMEAN = statDict[row[0]][2] evalList = [] for k, v in statDict.items(): if k != row[0]: if currentMAX < v[0]: evalList.append(str(k)) if not evalList: print('Nothing to evaluate') continue arcpy.SelectLayerByAttribute_management(roofZones, 'NEW_SELECTION', "gridcode = " + str(row[0])) arcpy.SelectLayerByLocation_management(roofZoneVerts, 'WITHIN_A_DISTANCE', roofZones, '5 Feet', 'NEW_SELECTION') # Remove vertices from self and those that are lower in elevation arcpy.SelectLayerByAttribute_management(roofZoneVerts, 'REMOVE_FROM_SELECTION', "gridcode = " + str( row[0]) + " or gridcode in (" + ', '.join( evalList) + ")") # Near Analysis to find angles from current section's center to selected vertices angleCount = int(arcpy.GetCount_management(roofZoneVerts)[0]) if angleCount <= 1: print('Not enough angles') continue arcpy.FeatureToPoint_management(roofZones, roofZoneCurrent) arcpy.Near_analysis(roofZoneVerts, roofZoneCurrent, location='LOCATION', angle='ANGLE') # Calculate range between min and max angles ( > 180 representing a hole) rangeCursor = arcpy.da.SearchCursor(roofZoneVerts, 'NEAR_ANGLE', sql_clause=(None, "ORDER BY NEAR_ANGLE")) firstRow = True angleMIN = None angleMAX = None rangeCount = 0 for angle in rangeCursor: rangeCount += 1 if firstRow is True: angleMIN = angle[0] firstRow = False if rangeCount == angleCount: angleMAX = angle[0] del rangeCursor # Calculate range in angles angleRange = angleMAX - angleMIN if angleRange >= 180: print('Found Possible hole') # Add building point to output, delete all related features for analysis and continue to next building arcpy.FeatureToPoint_management(currentBuilding, 'in_memory\BIN_' + str(bin)) arcpy.Append_management('in_memory\BIN_' + str(bin), outputPoints, 'NO_TEST') arcpy.Delete_management('in_memory\CurrentBuilding') arcpy.Delete_management('in_memory\BIN_' + str(bin)) # Delete all analysis related feature classes arcpy.Delete_management(zonalStats) arcpy.Delete_management(roofZones) arcpy.Delete_management(roofZoneVerts) arcpy.Delete_management(roofZoneVertFeature) arcpy.Delete_management(roofZoneVerts) arcpy.Delete_management(roofZoneCurrent) # arcpy.Delete_management(binPoints) # arcpy.Delete_management(currentTile) # arcpy.Delete_management(currentBuilding) break else: # Clear selections and current zone centroid and continue on current building arcpy.Delete_management(roofZoneCurrent) arcpy.SelectLayerByAttribute_management(roofZones, 'CLEAR_SELECTION') arcpy.SelectLayerByAttribute_management(roofZoneVerts, 'CLEAR_SELECTION') del hCursor print('Exited Loop') # Delete all analysis related feature classes arcpy.Delete_management(zonalStats) arcpy.Delete_management(roofZones) arcpy.Delete_management(roofZoneVerts) arcpy.Delete_management(roofZoneVertFeature) arcpy.Delete_management(roofZoneVerts) arcpy.Delete_management(roofZoneCurrent) # arcpy.Delete_management(binPoints) del binPoints # arcpy.Delete_management(currentTile) # arcpy.Delete_management(currentBuilding) del minB del minBRaster del zonalStats del roofZoneFeature del roofZones del roofZoneVertFeature del roofZoneVerts del roofZoneCurrent try: #Scroll through tile outlines to perform analysis on tile by tile basis print('Begining Tile Cursor') for tile in tileCursor: print(tile) currentTileFeature = 'F:\Lidar\lidar_2018\RearRoofTesting\ExportForTesting2.gdb\Tile_' + tile[0].rstrip('.las') #TODO Update to production path arcpy.env.outputCoordinateSystem = currentTileFeature #TODO append Surrounding Tiles #Select out all BINS in tile for analysis tileBINCursor = arcpy.da.SearchCursor(buildingFootprintsAll, 'BIN', where_clause="FileName = '" + tile[0] + "'") #tileBINList = list(set([r[0] for r in tileBINCursor]))#TODO This will be pulled from the lidar source in production tileBINList = ['1019616', '1064332'] #TODO temp BIN input for testing del tileBINCursor binCount = 0 for bin in tileBINList: binCount += 1 if bin not in []: #TODO Use for already processed BINs in testing, remove line for prod print('Starting BIN ' + str(bin)) #Create Local Workspace for each iteration workspace = 'F:\Lidar\lidar_2018\RearRoofTesting\lidartesting\WorkspaceGDB.gdb' if arcpy.Exists(workspace): arcpy.Delete_management(workspace) arcpy.CreateFileGDB_management('F:\Lidar\lidar_2018\RearRoofTesting\lidartesting', 'WorkspaceGDB') arcpy.env.workspace = workspace arcpy.env.overwriteOutput = True #Local Variables binPointsName = 'BIN_' + bin #binPoints = 'F:\Lidar\lidar_2018\RearRoofTesting\lidartesting\BinLidar.gdb\' + binPointsName binPoints = binPointsName currentBuilding = 'F:\Lidar\lidar_2018\RearRoofTesting\lidartesting\BuildingFootprints.gdb\BIN_'+bin minB = binPointsName + '_MinB' minBRaster = minB + '_Rast' zonalStats = binPointsName + '_ZonalStats' roofZoneFeature = binPointsName + '_RoofZones' roofZones = 'RoofZonesLayer' roofZoneVertFeature = binPointsName + '_RoofZoneVerts' roofZoneVerts = 'roofZoneVertsLayer' roofZoneCurrent = binPointsName + '_RoofZoneCurrent' # Copy Building Outline and elevation points to workspace currentTile = 'currentTileLayer' arcpy.MakeFeatureLayer_management(currentTileFeature, currentTile) arcpy.FeatureToLine_management(currentBuilding, "BIN_" + bin + "_Outline") arcpy.SelectLayerByAttribute_management(currentTile, 'NEW_SELECTION', "BIN = '" + bin + "' AND BUILDING_PERIMETER IS NULL") arcpy.FeatureClassToFeatureClass_conversion(currentTile, workspace, binPoints) """ print('Starting Point Extraction') pointCursor = arcpy.da.SearchCursor(currentTile, ['POINT_X', 'POINT_Y', 'POINT_Z']) pointList = [row[:] for row in pointCursor] del pointCursor print('Inserting Points') pointCount = int(arcpy.GetCount_management(currentTile)[0]) if arcpy.exists(binPoints) and int(arcpy.GetCount_management(binPoints)[0]) == pointCount: print('Point Layer Exists') continue else: arcpy.CreateFeatureclass_management(workspace, binPoints, 'POINT', currentTile, has_z='SAME_AS_TEMPLATE') count = 0 breaks = [int(float(pointCount) * float(b) / 100.0) for b in range(10, 100, 10)] for point in pointList: count += 1 if count in breaks: print('Appending Points ' + str(count) + '% Complete') arcpy.Append_management(arcpy.PointGeometry(arcpy.Point(point[0], point[1], point[2])), binPoints, 'NO_TEST') del pointList print(arcpy.GetCount_management(binPointsName)[0]) desc = arcpy.Describe(currentTile) print('Bin Points Selected') print(desc.FIDset) #arcpy.MakeFeatureLayer_management(currentTile, binPoints)#TODO Testing ->Hoping to keep """ # Establish Minimum Bounds For Processing On BIN print('Establishing Minimum Bounds') arcpy.MinimumBoundingGeometry_management(binPoints, minB, 'CONVEX_HULL') arcpy.PolygonToRaster_conversion(minB, 'OBJECTID', minBRaster, 'MAXIMUM_AREA', cellsize='0.5') #Set Raster Environments based off minimum bounds raster arcpy.env.cellSize = minBRaster arcpy.env.snapRaster = minBRaster arcpy.env.extent = minBRaster ###Create Roof Zones to compare heights#### #Create IDW Raster of Roof Elevations #arcpy.SelectLayerByAttribute_management(buildingFootprintsAll, 'NEW_SELECTION', "BIN = '" + str(bin) + "'") #print('Bin Outline Selected') #desc = arcpy.Describe(buildingFootprintsAll) #print(desc.FIDset) #Extract XY Vertices from BIN Feature """ features = [] outlineCursor = arcpy.da.SearchCursor(buildingFootprintsAll, 'SHAPE@') for b in outlineCursor: for part in b: features.append(arcpy.Polyline(arcpy.Array([[point.X, point.Y] for point in part]))) del outlineCursor """ #arcpy.CopyFeatures_management(buildingFootprintsAll, currentBuilding) #arcpy.MakeFeatureLayer_management(buildingFootprintsAll, currentBuilding) #print(arcpy.GetCount_management(currentBuilding)[0]) print('Beginning Raster Analysis') rasterizeRoof(currentBuilding, bin, binPoints, workspace) arcpy.Delete_management(workspace) print('Script Complete') except arcpy.ExecuteError: tb = sys.exc_info()[2] tbinfo = traceback.format_tb(tb)[0] pymsg = "PYTHON ERRORS: Traceback info: " + tbinfo + " Error Info: " + str(sys.exc_info()[1]) msgs = arcpy.GetMessages(2) print(msgs) print(pymsg) except: tb = sys.exc_info()[2] tbinfo = traceback.format_tb(tb)[0] pymsg = "PYTHON ERRORS: Traceback info: " + tbinfo + " Error Info: " + str(sys.exc_info()[1]) print(pymsg) ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

ShaunWalbridge · Answer

Hello Dani,

Patrick Hammons‌ clued me into this issue, thanks for posting on GeoNet. As Dan mentions, your script has a lot going on which can make debugging complicated. The raster intermediary layers will be typically saved out to their own temporary files, though you're just referencing them as variables in Python. One approach would be to create explicitly new GDB files for each bin step. Check that these results are as expected, and copy the final output you're interested in into another GDB. If that approach is working, then you can just delete this GDB, which will clear out all of the files used by it, assuming no locks are preventing its deletion. For the SA outputs, I would recommend explicitly binding these to rasters you want by using `raster.save` on them, and storing those in the output GDB.

The other more in-depth approach is to do what Dan suggested -- start adding debugging information, and dig into why it is failing in the second iteration as-is. You can start by just sprinkling print statements after each step, and making sure that it is executing the way you're expecting. There are more complicated debugging solutions, like pdb, but I'd start simple and see if that's enough to figure out the issue.

Cheers,
Shaun

DanPatterson_Retired · Answer

Dani, throw in some print statements in and around line 232- to report the outcome of whether the delete actually did anything or not and whether the new workspace is correct.

Also, the overwriteoutputs=True line should be out of the loop as a main property at the top of the script.  turning it on and off …. hmmm not sure if that works, just turn it on for the whole script so that you can reuse files and if you want to backup a result, use 'Copy' to make it a new name

There is a lot going on there. It would be nice to terminate after one iteration and see what remains at that point.

DaniInterrante · Answer

Hi Dan,

Thank you for sharing the formatting feature with me, is this more readable now?  I appreciate you taking a look.