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

DaniInterrante — Sat, 11 Dec 2021 06:07:17 GMT

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)
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

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

DanPatterson_Retired — Fri, 21 Jun 2019 00:09:45 GMT

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

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

DaniInterrante — Fri, 21 Jun 2019 13:54:03 GMT

Hi Dan,

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

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

DanPatterson_Retired — Fri, 21 Jun 2019 14:17:20 GMT

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.

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

ShaunWalbridge — Fri, 21 Jun 2019 16:26:27 GMT

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

topic Re: Spatial Analyst appears to be inhibiting iterative analysis when run in For Loop in Python Questions

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

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

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

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

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