Intersect issue in an iterative python script tool

[LisaNelson]

Hello,

I have a script tool (10.0 SP5, Python 2.6) that intersects a fishnet polygon and a line feature class.  The script uses an iterative loop to process individual single-part area of interest (i.e. extent) polys, running the intersect operation for each area of interest.  The code works great if no iteration occurs, but blows up when iterating the second feature (the first runs correctly).  The intersect operation is the point of failure - no output features are produced for the second output.  However, the same code executes correctly as snippets when run from the Geoprocessing --> Python window using the same inputs as the script tool.  Also, the toolbox Intersect tool works correctly as well.  It's just the script tool version that blows up.  Here's the snippet: 

arcpy.CopyFeatures_management(AOA_Layer,"in_memory\\temp_aoa"); showGpMessage() rows = arcpy.SearchCursor("in_memory\\temp_aoa"); showGpMessage()  for row in rows:            ...# variable set-up, etc.            ...# where Input_Roads_Features and Fishnet_Output vars are set by single-part area of interest extent            ...# repair geometry is run on Input_Roads_Features and Fishnet_Output prior to intersect attempt             # Intersect             inFeatures = []             inFeatures.append(Input_Roads_Features)             inFeatures.append(Fishnet_Output)             arcpy.Intersect_analysis(inFeatures, Intersect_Fishnet, "ALL", "", "INPUT"); showGpMessage()

For the second iteration, Intersect_Fishnet is produced but does not contain features.

Since the Python window and Intersect toolbox versions do work, my suspicion was that the 'inFeatures' list variable wasn't getting cleared/re-populated, but I think that does happen correctly since I re-initialize the list. 

Any insight/tips on iterative processing appreciated!

Lisa

[MathewCoyle]

You are going to need to post your entire code because what you posted here makes little sense without the context of where your variables are coming from.

Edit: All the examples seem to use lists fine.

[LisaNelson]

Good point!  Didn't want to spam the thread with hundreds of lines of code...  All inputs and outputs are file geodatabase objects.  So, here's 90 lines:

 # Generate distance rasters, iterating through AOA features if necessary:
        arcpy.CopyFeatures_management(AOA_Layer,tempFC); showGpMessage()
        rows = arcpy.SearchCursor(tempFC); showGpMessage()

        for row in rows:
            aoaFeature = row.shape
            aoaExt = aoaFeature.extent
            aoaFeatureArea = aoaFeature.area
            arcpy.env.extent = aoaExt
            fields = arcpy.ListFields(tempFC,aoaFieldName,"String"); showGpMessage()
            for field in fields: # only one iteration
                aoaName = row.getValue(field.name); showGpMessage()
                arcpy.AddMessage("\naoaName: "+aoaName)

            # Local variables
            arcpy.env.workspace = Output_Folder + "\\" + Output_GDBName
            
            Output_RoadDensity = Output_FileGDB + "\\src_" + aoaName.replace(" ", "_")  # Example: src_THRO_Central
            Intersect_Fishnet = Output_FileGDB + "\\src_" + aoaName.replace(" ", "_") + "_if"
            Constant_Raster = Output_FileGDB + "\\src_" + aoaName.replace(" ", "_") + "_int"
            Unique_IntegerRaster = Output_FileGDB + "\\src_" + aoaName.replace(" ", "_") + "_integer"
            Unique_Points = Output_FileGDB + "\\src_" + aoaName.replace(" ", "_") + "_pts"
            Fishnet_Output = Output_FileGDB + "\\src_" + aoaName.replace(" ", "_") + "_fishnet"  #Example: src_THRO_Central_fishnet
            Fishnet_Layer = "Fishnet_Layer"
            Sum_Table =  Output_FileGDB + "\\src_" + aoaName.replace(" ", "_") + "_sum"
   # Set up sum field name to use when calculating cell sums
            Calculated_Sum_Field = "" 
            Calculated_Sum_Field = "[src_" + aoaName.replace(" ", "_") + "_sum.SUM_CALC_LENGTH_KM]" # Example: src_THRO_Central_sum.SUM_CALC_LENGTH_KM
            Raster_Value_Source_Field = "src_" + aoaName.replace(" ", "_") + "_if.CELL_LENGTH_KM"
            
            #Delete existing files
            Files2Delete = [Fishnet_Output, Intersect_Fishnet, Fishnet_Layer,
                            Sum_Table, Constant_Raster, Unique_IntegerRaster, 
                            Unique_Points]
            for dfile in Files2Delete:
                try:
                    if arcpy.Exists(dfile): arcpy.Delete_management(dfile); showGpMessage()
                except: pass                
       
            arcpy.Buffer_analysis(aoaFeature,"in_memory\\temp_buffer", 1000 * .71); showGpMessage()
                 
            constRaster = arcpy.sa.CreateConstantRaster(1, "INTEGER", 1000, "in_memory\\temp_buffer"); showGpMessage()
            constRaster.save(Constant_Raster); showGpMessage()
            arcpy.DefineProjection_management(Constant_Raster, Output_Projection); showGpMessage()

            # Generate unique value raster to use as fishnet polygon source
            # Convert constant raster to points, use point ID to generate unique value integer raster
            constDesc = arcpy.Describe(Constant_Raster); showGpMessage()
            arcpy.env.extent = Constant_Raster
            arcpy.env.snapRaster = Constant_Raster
            arcpy.RasterToPoint_conversion(Constant_Raster, Unique_Points); showGpMessage()
            arcpy.Delete_management(Constant_Raster); showGpMessage()
            arcpy.PointToRaster_conversion(Unique_Points, "pointid", Unique_IntegerRaster, "MOST_FREQUENT", "", 1000); showGpMessage()
            arcpy.DefineProjection_management(Unique_IntegerRaster, Output_Projection); showGpMessage()
            arcpy.BuildRasterAttributeTable_management(Unique_IntegerRaster, "NONE"); showGpMessage()        

            # Raster to Polygon
            arcpy.RasterToPolygon_conversion(Unique_IntegerRaster, Fishnet_Output, "NO_SIMPLIFY", "VALUE"); showGpMessage()

            # Define Projection
            arcpy.DefineProjection_management(Fishnet_Output, Output_Projection); showGpMessage()
            arcpy.RepairGeometry_management(Fishnet_Output, "DELETE_NULL"); showGpMessage()
            
            # Intersect, attempting to re-initialize inFeatures list
            inFeatures = []
            inFeatures.append(Input_Roads_Features)
            inFeatures.append(Fishnet_Output)

            ##### Critical intersect operation: ########
            arcpy.Intersect_analysis(inFeatures, Intersect_Fishnet, "ALL", "", "INPUT"); showGpMessage()

            # Repair Geometry
            arcpy.RepairGeometry_management(Intersect_Fishnet, "DELETE_NULL"); showGpMessage()

            # Add CELL_LENGTH_KM Field (= sum of segment lengths within 'grid cell'/intersected fishnet poly)
            arcpy.AddField_management(Intersect_Fishnet, "CELL_LENGTH_KM", "FLOAT", "", "", "", "", "NULLABLE", "NON_REQUIRED", ""); showGpMessage()

            # Make Feature Layer
            arcpy.MakeFeatureLayer_management(Intersect_Fishnet, Fishnet_Layer); showGpMessage()
            
            # Re-Calculate CALC_LENGTH_KM Field: accounts for the intersection operation
            # Calculates segment lengths within 'grid cells' (i.e., within intersected fishnet poly)
            arcpy.CalculateField_management(Fishnet_Layer, "CALC_LENGTH_KM", "!shape.length@kilometers!", "PYTHON", ""); showGpMessage()

            # Summary Statistics: find the sum of segments within each fishnet polygon; unique case field = grid_code which came from Unique_IntegerRaster via Fishnet_Output
            arcpy.Statistics_analysis(Fishnet_Layer, Sum_Table, "CALC_LENGTH_KM SUM", "grid_code"); showGpMessage()
            # Add Join (inner join logic used)
            arcpy.AddJoin_management(Fishnet_Layer, "grid_code", Sum_Table, "grid_code", "KEEP_ALL"); showGpMessage()

            # Calculate CELL_LENGTH_KM Field = re-calculated, summed CALC_LENGTH_KM value from joined field
            
            #  Example Calculated_Sum_Field variable:     src_THRO_Central_sum.SUM_CALC_LENGTH_KM
            ##### Blows up here since SUM_CALC_LENGTH_KM is null ###########
            arcpy.CalculateField_management(Fishnet_Layer, "CELL_LENGTH_KM", Calculated_Sum_Field, "VB"); showGpMessage()  

It crashes ArcMap completely, so I don't get a record of output messages, which might help in further diagnosing the issue.  However, when the script bails, the Fishnet_Layer (i.e. the intersected fishnet polys and input road features) feature class is empty.  Obviously, this is responsible for the null/missing attributes that cause arcpy.CalculateField_management to fail.

Still can't figure out why the intersect op fails for the second feature, but not the first...

Thanks for any assistance!

Lisa

[MathewCoyle · Accepted solution]

I don't really like how this is used in a cursor loop. I imagine it would be cleaner to assemble a list of unique values and loop through that creating feature layers with definition queries to process on.

Is your aoaName variable a unique value? Also, you don't seem to be calculating the same field you create.

arcpy.AddField_management(Intersect_Fishnet, "CELL_LENGTH_KM", "FLOAT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "") arcpy.CalculateField_management(Fishnet_Layer, "CALC_LENGTH_KM", "!shape.length@kilometers!", "PYTHON", "")

Probably not related, but also post how you create your Output_FileGDB variable.

[LisaNelson]

Hello - yes, you were correct about the 'layer' vs. 'list' approach. I realized I was trying to re-iterate using the list when all I needed to do was an explicit select followed by the fishnet generation. 

Thanks for taking a look at this and for responding!

Lisa