1

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.

1

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")
  • Thank you Aaron for your help! Running the 'Times' tool, I do get an error however RuntimeError: ERROR 000824: The tool is not licensed. Any idea how to fix this? Following this thread I tried importing arcinfo, but that does not solve the problem. Secondly, if outws is just the same gdb as the inws, can I save the results to the gdb with the same .save command, but then just excluding the ".tif" extension? – AstridB Apr 21 '15 at 11:41
  • edit: if I run the code in Python IDLE, I get the license error, but when I try to run the code in the python window within ArcGIS, it keeps on crashing... – AstridB Apr 21 '15 at 12:19
  • You need to have the Spatial Analyst extension for ArcGIS. Try adding arcpy.sa. in front of the Raster commands (e.g. arcpy.sa.Raster(forestmask) and arcpy.sa.Raster(loss). Regarding the second point, outws is not the same workspace as inws. However, you can change this to so that the output directory is the input directory. Yes, removing the .tif extension will be required to write to FGDB raster format. – Aaron Apr 21 '15 at 12:21
  • 1
    @yanes Sorry, my error arcpy.CheckInExtension allows both cases. Post edited to include your comments. Thanks. – Aaron Nov 3 '16 at 2:57
  • 1
    Actually, thanks. it is a good reminder for me. As- sometimes I just get stuck for not checking for case-sensitive commands, I always forget how case sensitive python could be. – yanes Nov 3 '16 at 2: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.