1

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
  • 2
    Have you timed each section/function to see which part takes up most time? What was the result? – Martin Apr 7 '16 at 12:11
  • Thanks for the hint. I identified the most computation intensive parts and added the Information to the question (see EDIT). – the_chimp Apr 7 '16 at 12:43
  • Your conversions to vector format are going to be horribly slow. Most calculations involving raster data can be performed entirely in raster format, with a huge speedup in processing time. (The calculations you describe ought to take tiny fractions of a second rather than minutes.) If you could explain what you are attempting to do at a sufficient level of abstraction--that is, without referring to a particular algorithm--then readers might be able to suggest much faster data structures and algorithms for achieving your ends. – whuber Apr 7 '16 at 21:58
  • Is it also possible to stay with rasters when intersecting a raster (my approach: raster that is converted to a poylgon) with a Point dataset? The calculateArea(raster) function should take a raster as a parameter and calculate its size, where I tried to count the cells and multiply it with the cell size. Any other approaches and ideas are welcome. – the_chimp Apr 8 '16 at 6:00
  • If I have understood what you mean by "intersect" a raster with points, then the answer is yes. All raster GISes support efficient mechanisms to query the values of rasters at point locations (without converting the entire raster back to vector format). It's unclear what you mean by the area of a raster: if it's just the count of cells in it, simply multiply the product of its dimensions by the squared cellsize. If it's the count of non-NoData cells, you can rapdily construct the non-NoData indicator grid: its attribute table will have a single row giving the count. – whuber Apr 8 '16 at 15:53
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)
  • I tried to use your code, but the code exits trying to calculate the numberOfCells. I get the error "'in_table' is not a table or a featureclass " - although the Int() calculation works fine. Any ideas? The cellSize can be queried, so I think the conversion to the integer raster works? – the_chimp Apr 8 '16 at 5:26
  • Conversion works. Calculate cells your way, it is fast anyhow – FelixIP Apr 8 '16 at 7:01
  • Either way I try to calculate the number of cells fails, as the created raster after the Int(raster) conversion cannot be read - neither with your suggestion nor with a simple cursor and row access. And if that does not work, I cannot convert it with the Int(raster) but have to stick with the Int_3d(raster) which takes multiple times longer. – the_chimp Apr 8 '16 at 7:14
  • Last try. Make raster layer out of integer raster and see if it will make table accessible – FelixIP Apr 8 '16 at 10:37
1

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.

  • Thanks a lot. I currently try to use your function, but I run into Problems, as the cell size is determined to be 1 and cell Count results in 47488 instead of 24375 - any ideas? But except these issues, the function only takes a few seconds - so it would be a huge improvement regarding computation time. I would like to stick with your Approach. – the_chimp Apr 8 '16 at 5:01
  • Hey the_chimp - I think the tuple returned by floored_array.shape[x] might actually be 3 elements long in the case of a single-band raster, which might be causing some of this problem. floored_array.shape[x] might be indexing a tuple like (1, 550, 750) - we are interested in only the last two dimensions of the array. I think the first dimension relates to the number of bands in the raster. Additionally, @dslamb's raster.meanCellWidth * ...Height-based method is a great alternative to .GetRasterProperties_man... Upvoted @dslamb. – Jim Apr 8 '16 at 15:08
1

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

  1. 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])

  1. 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.

  1. 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.

  • If the extent of the raster is not squared, but e.g. like the outline of a city, do the raster.width and raster.height refer to the bounding box? If so, I think then I could not use it in this way, as I need to get the exact number of cells covering the study area extent. – the_chimp Apr 8 '16 at 5:31
  • Good point. I also realized I left off that I needed to multiply by the cell area not just the length of side of the cell. I've corrected that above, and added the code to subtract out the NoData cells so the area correctly matches an irregular shaped polygon. It still makes use of NumPy but not of the Int_3D conversion. I tested this on a floating point raster too. It looks like you already solved your problem in another comment, so I've just added this as another option. – dslamb Apr 8 '16 at 13:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.