Looping through two raster folders to perform raster calculation?

Question

I have two geodatabases with rasters containing information about "treecover" and "loss". The file names are slightly different but can be matched based on the last 8 characters (the coordinates of each raster granule).See image

I can perform this simple conditional statement using either CON or RASTER CALCULATOR to retrieve all "Loss" data situated on pixels which have a value for "treecover" greater than or equal to 50, with 9999 being assigned as the output value if the statement is FALSE. The output raster is exactly what I want. I just need to do it for 137 pairs of rasters

RASTER CALC STATEMENT: con("treecover raster" >= 50, "Loss raster", 9999)

I'm very much a beginner at using Python, but this is what I've managed to do so far based on similar questions posted here. Fair warning, I don't really know what some parts of the script do.

import arcpy, os, math

arcpy.CheckOutExtension("Spatial")


CVR =
r"C:\Aquifer_recalc\Aquifer_recalc_NA\Treecover_NA\Treecover_Inputs.gdb"

LSS = r"C:\Aquifer_recalc\Aquifer_recalc_NA\Lossyear_NA\Lossyear_Inputs.gdb"


outws = r"C:\Aquifer_recalc\Aquifer_recalc_NA\ras_calc_con"


arcpy.env.workspace = CVR


rasterlist  = arcpy.ListRasters()


for r in rasterlist:
    r1 = arcpy.sa.Raster(r)
    r2 = arcpy.sa.Raster(os.path.join(LSS,r))

    result =  con(r1 >=50,r2,999)

    result.save(os.path.join(outws, r))

Answer 1

The arcpy.ListRasters() function will give you a list with the filenames of all the rasters in a given workspace. This way, you can specify the first GDB as your workspace and get a list of all the rasters in it. Then, you can repeat this process for the second GDB.

CVR = r"C:\Aquifer_recalc\Aquifer_recalc_NA\Treecover_NA\Treecover_Inputs.gdb"
LSS = r"C:\Aquifer_recalc\Aquifer_recalc_NA\Lossyear_NA\Lossyear_Inputs.gdb"

# get first set of rasters
arcpy.env.workspace = CVR
cvr_rasters = arcpy.ListRasters()

# get second set of rasters
arcpy.env.workspace = LSS
lss_rasters = arcpy.ListRasters()

Because the order of your rasters match, you don't need to worry to order the rasters in a way pairs should match. This means the first name of cvr_rasters corresponds to the first name of lss_rasters (as they are sorted alphabetically). You can inspect these variables to make sure.

Having these two lists, you can take advantage of the zip() function, which will create a list of tuples for you to iterate through. The first tuple will have the first raster from CVR and the first raster from LSS, and so on. This way you can access two rasters at the same time:

for cvr_ras, lss_ras in zip(cvr_rasters, lss_rasters):
    r1 = arcpy.sa.Raster(cvr_ras)
    r2 = arcpy.sa.Raster(lss_ras)
    result = Con(r1 >=50, r2, 999)
    result.save(os.path.join(outws, cvr_ras))

Note that Raster is a class in arcpy that lets you create a Raster object from a raster file so you can perform raster calculator operations. It does not have anything to do with the name of your result. If you inspect either of cvr_rasters or lss_rasters you will see the oputput will be something like this:

['Hansen_GFC2015_treecover2000_60N_160W', ...]

This is a list containing one string for each raster you have in that GDB. However, you need something to tell arcpy you are dealing with rasters (hencer the Raster class). After this, you execute the Con() function and arcpy knows what you rasters actually are.

Regarding the name of the output, the Con() function returns a new Raster object. It has a .save() method that takes as a parameter the output name or path of the file. If you pass only a string and not a full path (e.g. 'my_raster' instead of r'C:\users\user01\my_raster.tif') it will try to save the raster in your current workspace. Right now, the name of one of your output file will be

"C:\Aquifer_recalc\Aquifer_recalc_NA\ras_calc_conHansen_GFC2015_treecover2000__60N_160W" or whatever numbers the first raster have. This is because it is joining outws with cvr_ras (which changes every iteration) using os.path.join().

I recommend you to add the following line before the for loop:

arcpy.env.workspace = r"C:\Aquifer_recalc\Aquifer_recalc_NA"

And then add the following after result = Con(r1 >=50, r2, 999)

out_name = "loss_{}".format(cvr_ras[-8:])
result.save(out_name)

In the case of the example raster, out_name will be "loss_60N_160W" and will be stored in C:\Aquifer_recalc\Aquifer_recalc_NA. It will take the last 8 characters of the current value of cvr_ras every loop and use them to create the output file name. You can play around with the out_name variable to get different names if you wanted to.