[GIS] How to use the math tool (“times”) in ArcGIS modelbuilder to multiply two raster sets while using an iterator

iteratormodelbuilderpythonraster-calculator

I have two sets of rasters (in the same gdb), representing "loss" and "forestmask".

I would like to multiply each forestmask raster (value 1 = forest) with its corresponding loss raster (value 1=loss) to create a raster that shows deforestation (in which value 1 = deforestation).

Each loss raster should be linked to one corresponding forestmask raster (they have the same extent):
e.g. 

GFC2014_forestmask_00N_020E * GFC2014_loss_00N_020E = GFC2014_deforestation_00N_020E (=output filename)

I guess this should not be too hard to do, and I tried to do this with the ArcGIS modelbuilder, iterating through the "forest mask" rasters using a wildcard, but I got stuck when I wanted to make use of the Parse Path tool and in-line variable tricks for defining the second input file (i.e. the loss raster). Right now, I tried to add a new variable (=file) in the model, called loss below, with input:

"GFC_2014_loss_"+ "%Name%"[21:]

But this is not working: the Times tool does not accept the loss variable as it is as the second input. (not sure if the syntax is correct like this, and if this can be used in a "file" variable).

model screenshot

I think I know how to do this in R with the raster package, but I heard Python is faster for these kind of calculations, so any help is welcome. Preferably I would like to do this in the model builder so that I can easily link it to other models I created.

Best Answer

Since you have the Python tag, the following example is how you would accomplish this task using raster algebra in Python. Here I am assuming you have a variety of rasters in one GDB with a unique ID for each pair. The output is written to a separate folder in .tif format.

import arcpy, os
from arcpy.sa import *
arcpy.CheckOutExtension ("Spatial")

arcpy.env.workspace = r'C:\temp\temp.gdb'
inws  = arcpy.env.workspace
outws = r'C:\temp\out'

# List all of the "forestmask" rasters
rasters = arcpy.ListRasters("*forestmask*")

# 1) Assuming there is a "loss" counterpart with a unique ID
# 2) Loop through your raster list
for r in rasters:
    # Extract prefix and suffix for naming purposes later
    prefix = r.split("_")[0] + "_"  # Get just the prefix
    suffix = "_" + r.split("_")[2] + "_" + r.split("_")[3] # Get just the suffix

    # Define the "forestmask" and "loss" paths
    forestmask = os.path.join(inws, r)
    loss       = os.path.join(inws, r.replace("forestmask", "loss")) # Replace "forestmask" with "loss"

    # Convert these to raster objects
    r1 = Raster(forestmask)
    r2 = Raster(loss)

    # perform the raster algebra
    mycalc = r1 * r2

    # Save to disk (write to .tif)
    mycalc.save(os.path.join(outws, prefix + "deforestation" + suffix + ".tif"))

arcpy.CheckInExtension ("Spatial")
Related Question