I developed an arcpy script tool that takes as input multiple rasters, a criterion to select specific cells from the raster which are then intersected with a point shapefile. At the end a value is calculated and stored in a created table.
The script is quite simple, but inputting 7 raster files where each has 24 300 cells makes the script to run for 30minutes. Now I try to find a way how to accelerate the whole process, as I have to use this script more often in the future.
Is there a way to make the computation process faster? Or is it the case with raster files that they increase the computing time in this amount? I already tried to simplify the tool and as was mentioned in several posts that the Cursors could make a Problem I commented them out and did not write the resulting value anywhere, just running through the script, but with no effect. So I think the main part, where it Comes to processing time, is the one where the size of the rasters are calculated and the condition Statement is used as well as intersecting the raster with the points.
Measured computation time for the whole script:
As suggested, I used the
time functionality to find out, which parts of the code need the most
computation time - it is the function
calculateArea(raster). This function is called two times per each interation - once to calculate the whole area of the raster (takes ~70 seconds) and once to calculate the area of those cells that match the criterion (takes ~130 seconds). Other parts that need quite some time are the intersection (takes ~30 seconds) and the condition Statement (takes ~15 seconds). All other calculations need less than 1-2 seconds and are in this case negligible.
Measured computation time for the function calculateArea():
I focused on the function
calculateArea() and had a look which parts take what computation time:
- Create integer raster: 40-70 seconds (once 120 seconds)
- Get cell size: 10 seconds
- Calculate area (cellSize * cellSize * numberCells): 8 seconds
- Cursor initialization and query the number of cells: negligible, 0.009 and 0.8 seconds
The main probem is how to query the number of cells in the raster for which the size should be calculated. Using the
arcpy.Int_3d(raster) function takes extremely Long compared to the
arcpy.Int(raster) function, but the latter does not allow access afterwards to get the number of cells. Neither the aproach with the search cursor nor the suggested ones (numpy) lets me access the "COUNT" field in the raster attribute table (although this exists, I tried it) - it always gives an error that the 'Input_table is not a table or feature class'
Here is my code:
import arcpy import os import sys #params input_raster = arcpy.GetParameterAsText(0) criterion = arcpy.GetParameterAsText(1) input_points = arcpy.GetParameterAsText(2) output = arcpy.GetParameterAsText(3) #create table table_path = os.path.splitext(os.path.dirname(output)) table_name = os.path.splitext(os.path.basename(output)) arcpy.CreateTable_management(table_path, table_name + '.dbf') arcpy.AddField_management(output + '.dbf', 'Name', 'TEXT') arcpy.AddField_management(output + '.dbf', 'Result', 'DOUBLE') #function to calculate the area of a given raster def calculateArea(raster): rasterInt = arcpy.Int_3d(raster) try: cur = arcpy.da.SearchCursor(rasterInt, 'Count') numberOfCells = 0 for row in cur: numberOfCells += row cellSize = float(arcpy.GetRasterProperties_management(rasterInt, 'CELLSIZEX').getOutput(0)) area = cellSize * cellSize * numberOfCells return area except Exception: return 0 count_points = arcpy.GetCount_management(input_points) value1 = float(count_points.getOutput(0)) insertCur = arcpy.InsertCursor(output + '.dbf') for input in input_raster.split(';'): #for each inputted raster value2 = calculateArea(input) my_raster = arcpy.Raster(input) mean = my_raster.mean std = my_raster.standardDeviation if criterion == 'values above mean': extracted_cells = arcpy.sa.Con(my_raster >= mean, 1) else: extracted_cells = arcpy.sa.Con(my_raster >= mean+std, 1) value3 = calculateArea(extracted_cells) if value3 > 0: #if cells which match the criterion exist, convert to polygon and intersect with points my_raster_poly = arcpy.RasterToPolygon_conversion(extracted_cells, arcpy.Polygon, 'NO_SIMPLIFY') intersectedFeatures = arcpy.Intersect_analysis([my_raster_poly, input_points], arcpy.Point) count_intersected_features = arcpy.GetCount_management(intersectedFeatures) value4 = float(count_intersected_features.getOutput(0)) else: value4 = 0 #calculate a result value result = (value3/value1) / (value4/value2) #insert row row = insertCur.newRow() row.Name = os.path.splitext(os.path.basename(input_raster)) row.Result = result insertCur.insertRow(row) del row, insertCur