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'
tileOutlineAllFeature = 'F:\\Lidar\\lidar_2018\\RearRoofTesting\\lidartesting\\BadHouse.gdb\\Tile_26875E244534N'
tileOutlineAll = 'tileOutlineAllLayer'
tileOutlineCurrentFeature = 'tileOutlineCurrent'
tileOutlineCurrent = 'tileOutlineCurrentLayer'
buildingFootprintsAllFeature = 'F:\\Lidar\\lidar_2018\\RearRoofTesting\\lidartesting\\BuildingFootprints.gdb\\TestBuildingsWithVios'
buildingFootprintsAll = 'buildingFootprintsAllLayer'
buildingFootprintsFeature = 'buildingFootprintsCurrent'
buildingPointsFeature = 'CurrentInputBuildingPoints'
buildingFootprints = 'buildingFootprintsLayer'
outputPoints = outputGbd + '\\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):
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'
outputGbd = 'F:\\Lidar\\lidar_2018\\RearRoofTesting\\lidartesting\\TestFinal.gdb'
outputPoints = outputGbd + '\\HoleFlags'
idw = Idw(binPoints, 'POINT_Z', '0.5', in_barrier_polyline_features="BIN_" + bin + "_Outline")
arcpy.Delete_management("BIN_" + bin + "_Outline")
print('Slope Analysis')
slope = Slope(idw, 'PERCENT_RISE', '', 'GEODESIC')
print('Con Analysis')
rawGroups = Con((slope < 90) & ~IsNull(minBRaster), 1)
print('Boundary Clean')
rawGroupsClean = BoundaryClean(rawGroups, '', 'ONE_WAY')
print('Region Group')
groups = RegionGroup(rawGroupsClean, 'FOUR', 'WITHIN')
print('Zonal Statistics')
ZonalStatisticsAsTable(groups, 'Value', idw, zonalStats)
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
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)
print('Pulling stats')
statCursor = arcpy.da.SearchCursor(zonalStats, ['Value', 'MIN', 'MAX', 'MEAN'])
statDict = {row[0]: list(row[1:]) for row in statCursor}
del statCursor
print('Feature Vertices to Points')
arcpy.FeatureVerticesToPoints_management(roofZones, roofZoneVertFeature)
print('Making Feature Layer')
arcpy.MakeFeatureLayer_management(roofZoneVertFeature, roofZoneVerts)
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')
arcpy.SelectLayerByAttribute_management(roofZoneVerts, 'REMOVE_FROM_SELECTION',
"gridcode = " + str(
row[0]) + " or gridcode in (" + ', '.join(
evalList) + ")")
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')
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
angleRange = angleMAX - angleMIN
if angleRange >= 180:
print('Found Possible hole')
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))
arcpy.Delete_management(zonalStats)
arcpy.Delete_management(roofZones)
arcpy.Delete_management(roofZoneVerts)
arcpy.Delete_management(roofZoneVertFeature)
arcpy.Delete_management(roofZoneVerts)
arcpy.Delete_management(roofZoneCurrent)
break
else:
arcpy.Delete_management(roofZoneCurrent)
arcpy.SelectLayerByAttribute_management(roofZones, 'CLEAR_SELECTION')
arcpy.SelectLayerByAttribute_management(roofZoneVerts, 'CLEAR_SELECTION')
del hCursor
print('Exited Loop')
arcpy.Delete_management(zonalStats)
arcpy.Delete_management(roofZones)
arcpy.Delete_management(roofZoneVerts)
arcpy.Delete_management(roofZoneVertFeature)
arcpy.Delete_management(roofZoneVerts)
arcpy.Delete_management(roofZoneCurrent)
del binPoints
del minB
del minBRaster
del zonalStats
del roofZoneFeature
del roofZones
del roofZoneVertFeature
del roofZoneVerts
del roofZoneCurrent
try:
print('Begining Tile Cursor')
for tile in tileCursor:
print(tile)
currentTileFeature = 'F:\\Lidar\\lidar_2018\\RearRoofTesting\\ExportForTesting2.gdb\\Tile_' + tile[0].rstrip('.las')
arcpy.env.outputCoordinateSystem = currentTileFeature
tileBINCursor = arcpy.da.SearchCursor(buildingFootprintsAll, 'BIN', where_clause="FileName = '" + tile[0] + "'")
tileBINList = ['1019616', '1064332']
del tileBINCursor
binCount = 0
for bin in tileBINList:
binCount += 1
if bin not in []:
print('Starting BIN ' + str(bin))
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
binPointsName = 'BIN_' + bin
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'
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)
"""
print('Establishing Minimum Bounds')
arcpy.MinimumBoundingGeometry_management(binPoints, minB, 'CONVEX_HULL')
arcpy.PolygonToRaster_conversion(minB, 'OBJECTID', minBRaster, 'MAXIMUM_AREA', cellsize='0.5')
arcpy.env.cellSize = minBRaster
arcpy.env.snapRaster = minBRaster
arcpy.env.extent = minBRaster
"""
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
"""
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)