AnsweredAssumed Answered

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

Question asked by daniel.interrante.phl on Jun 20, 2019
Latest reply on Jun 21, 2019 by SWalbridge-esristaff

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)

Outcomes