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.
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)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 **).
## <Result '18'>## <Result '6'>## <Result '72'>## <Result '2'>## <Result '2985984'>## <Result '2'>## <Result '-12'># compound statements work too
raster3 = (raster1 / raster2) * (raster2 - 3.0)
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")## <Result '6'>Logical operators, including compound operations, work on rasters too.
## <Result '1'>## <Result '0'>## <Result '1'>## <Result '0'>## <Result '1'>## <Result '1'># compound statements work as expected too
raster3 = (raster1 > 10) & (raster2 > 10)
arcpy$management$GetRasterProperties(raster3, "MAXIMUM")## <Result '0'>## <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'>