# RemoteSensing.io

Remote Sensing tips and tricks

## Using ArcPy for Raster Analysis

In order to do bulk raster analysis, or anything slightely bespoke, I will generally opt for using GDAL with Python. However, its possible to achieve similar things within an ESRI ArcGIS world if you want to, using ArcPy.

If all you need to do is a simple point operation (in this case, multiplying each pixel by 2), then this can very quickly be achieved using ArcPy’s raster calculator. The Spatial Analysis extension is required for this.

``````import arcpy
import numpy as np
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")

infile = "inputfile.tif"

# do any calculation on the raster object
x = Raster(infile) * 2

# save the new raster
x.save("outfile.tif")``````

If this is implemented using the Python prompt in Arcmap, as soon as `x` has been created it will appear as a layer in the “Table of Contents” labelled `x`. At this point it’s only saved in memory. If you quit Arcmap, the raster will be lost.

To do anything more complex with rasters such as filtering and neighbourhood functions, this method won’t cut it. In this case it’s necessary to start using NumPy arrays. Just as I have previously written about using NumPy arrays for raster processing with GDAL, below is a similar example using ArcPy.

Just for the purpose of this post, I’m only going to multiply the raster by 2 again. But obviously once data is in an array, anything in possible.

``````import arcpy
import numpy as np

infile  = "inputfile.tif"
in_arr  = arcpy.RasterToNumPyArray(infile, nodata_to_value=9999)

# NB: I find it useful to use the optional argument `nodata_to_value`.
# That way, ArcPy will read the no data value from the dataset, and
# set it to this defined value in the NumPy array. You then know exactly
# what number represent no data in the array.

# do some processing
out_arr = in_arr * 2

# gather some information on the original file
spatialref = arcpy.Describe(infile).spatialReference
cellsize1  = arcpy.Describe(infile).meanCellHeight
cellsize2  = arcpy.Describe(infile).meanCellWidth
extent     = arcpy.Describe(infile).Extent
pnt        = arcpy.Point(extent.XMin,extent.YMin)

# save the raster
out_ras = arcpy.NumPyArrayToRaster(out_arr,pnt,cellsize1,cellsize2)
arcpy.CopyRaster_management(out_ras,"outfile.tif")
arcpy.DefineProjection_management("outfile.tif", spatialref)``````
Spot an error raise a ticket