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.
EDIT:
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
EDIT 2
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))[0]
table_name = os.path.splitext(os.path.basename(output))[0]
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[0]
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))[0]
row.Result = result
insertCur.insertRow(row)
del row, insertCur