- Home
- :
- All Communities
- :
- Developers
- :
- Python
- :
- Python Questions
- :
- Spatial Analyst appears to be inhibiting iterative...
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page
Spatial Analyst appears to be inhibiting iterative analysis when run in For Loop
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + 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:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])
print(pymsg)
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Dani
can you format your code
/blogs/dan_patterson/2016/08/14/script-formatting
and can you move your 'defs' outside the try except block
put them underneath your imports and outside the block
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Dan,
Thank you for sharing the formatting feature with me, is this more readable now? I appreciate you taking a look.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
