Raster math with arcpy

The arcpy package, via reticulate, supports Raster math operations. This vignette demonstrates basic usage of arcpy Raster Algebra from R. First, we'll set up a workspace and create some temporary rasters.

library(arcpy)
arcpy$env$workspace = tempdir()
arcpy$env$scratchWorkspace = tempdir()

cellSize = 2
outExtent = arcpy$Extent(0, 0, 250, 250)

# Check out the ArcGIS Spatial Analyst extension license
arcpy$CheckOutExtension("Spatial")
## [1] "CheckedOut"
# Execute CreateConstantRaster
raster1 = arcpy$sa$CreateConstantRaster(12, "FLOAT",
  cellSize, outExtent)

raster2 = arcpy$sa$CreateConstantRaster(6, "FLOAT",
  cellSize, outExtent)

Math operations

All of the basic mathematical operations are supported. Usage follows the R operator symbology, rather than the Python operator symbology or arcpy operator symbology (e.g., Modulus is achieved using %%, not %, and exponentiation is achieved using ^ rather than pow() or **).

raster3 = raster1 + raster2
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '18'>
raster3 = raster1 - raster2
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '6'>
raster3 = raster1 * raster2
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '72'>
raster3 = raster1 / raster2
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '2'>
raster3 = raster1 ^ raster2
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '2985984'>
raster3 = raster1 %/% raster2
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '2'>
raster3 = -raster1
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '-12'>
# compound statements work too
raster3 = (raster1 / raster2) * (raster2 - 3.0)
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '6'>

Logical operators

Logical operators, including compound operations, work on rasters too.

raster3 = raster1 > raster2
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '1'>
raster3 = raster1 < raster2
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '0'>
raster3 = raster1 != raster2
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '1'>
raster3 = raster1 == raster2
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '0'>
raster3 = raster1 >= 12
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '1'>
raster3 = raster2 <= 6
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '1'>
# compound statements work as expected too
raster3 = (raster1 > 10) & (raster2 > 10)
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '0'>
raster3 = (raster1 > 10) | (raster2 > 10)
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '1'>
# The not operator works in R too
raster3 = !(raster2 > 10)
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '1'>
# or you can do
raster3 = arcpy$sa$BooleanNot(raster2 > 10)
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")
## <Result '1'>