ArcPy Script Tool takes 30min - Ideas to accelerate the (raster) process

Question

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

Answer 1

I'd say importing SA functions, setting workspace to fastest media possible using TableToNumpyArray might help as well:

1. from arcpy.sa import *
2. from arcpy import env
3. env.workspace='in_memory'

If in_memory doesn't work set it to folder (not FGDB) on fastest disk. This is where ArcGIS stores temp rasters

    RasterInt=Int(raster)

    numberOfCells =sum([row[0] for row in arcpy.da.TableToNumPyArray(RasterInt,"COUNT")])

#

To get cell size simply use:

cellSize = RasterInt.meanCellHeight

NOTE: when you imported Spatial Analyst functions once, it is much easier and faster to use them, e.g.:

extracted_cells = Con(my_raster >= mean+std, 1)

Answer 2

The following might be more of a set of a partial suggestions for your situation rather than a direct and complete answer, but have you considered performing the arcpy.Int_3d() , numberOfCells and area work using numpy arrays? I've tested the function below on some small (6k x 2.5k pixels) single-band rasters and I hope it will reduce your running time. I could of course be completely wrong here; maybe the arcpy.Int_3d function uses similar numpy objects already?

def Int_3d_numpy_version(input_raster):
    'an attempt at arcpy.Int_3d() using numpy'
    input_array = arcpy.RasterToNumPyArray(input_raster)
    floored_array = input_array // 1 # Floor division cell-by-cell (element-wise)
    output_raster = arcpy.NumPyArrayToRaster(floored_array)
    cell_count = floored_array.shape[0] * floored_array.shape[1] # Assumes a 1-band raster, floored_array.size is another option
    cellSize = float(arcpy.GetRasterProperties_management(output_raster, 'CELLSIZEX').getOutput(0))
    area = cellSize * cellSize * cell_count
    return output_raster, cell_count, area

Best Luck.

Answer 3

I second that NumPy is fast for calculations on rasters (although converting back to a raster file can seemingly take a long time for the few times I've used it).

1. You could also calculate your area using the Raster properties. Something like this without the need of Int_3d: Edited to correct calculation

area = (my_raster.width * my_raster.height) * my_raster.meanCellWidth*my_raster.meanCellHeight

2. Also, you can shorten these to a single line of code (althought it probably won't affect speed):

count_intersected_features = arcpy.GetCount_management(intersectedFeatures) value4 = float(count_intersected_features.getOutput(0))

Instead do:

value4 = float(arcpy.GetCount_management(intersectedFeatures)[0])

3. Finally, using "in_memory" paths or using feature/raster layers can help speed things up too. This may cause memory issues if you don't delete them using Delete_management or overwrite them.

Instead of doing an intersect and writing to disk (which takes time) you can do something like this:

arcpy.MakeFeatureLayer_management(input_points, "points_layer")

arcpy.SelectLayerByLocation_management('points_layer', 'intersect', my_raster_poly)

matchcount = int(arcpy.GetCount_management('points_layer')[0])

(example above modified from: http://pro.arcgis.com/en/pro-app/tool-reference/data-management/select-layer-by-location.htm)

You can just create the points layer once, and then use it different times outside of the loop.

4. I've also found that running a toolbox in ArcCatalog instead of ArcMap can help speed up things a little. If you can, even just run the script in IDLE/some IDE/Commandline completely outside of ArcGIS and there will be some improvements.

Hope those help,

David

Edited

For irregular mask on a raster you can subtract out the no data cells.

Get a raster where NoData is 1 and Values are 0:

my_raster_null =arcpy.sa.isNull(my_raster)

Numpy gets you the total null cells:

nodata_cellcount = np.sum(arcpy.RasterToNumPyArray(my_raster_null))

Calculate area:

area = ((my_raster.height*my_raster.width)-nodata_cellcount) * (my_raster.meanCellWidth*my_raster.meanCellHeight)

This may not match the area of your polygon because of cell size.