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



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))

The way I would like images to be paired

1 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:])

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.

  • Thank you so much! I'm getting an error saying "NameError: name 'con' is not defined". Does 'con' even need to be defined since it's part of Raster calculator?
    – seak23
    Aug 23, 2019 at 19:45
  • 1
    Oh yeah, I didn't notice. You have to import it first. Write this line at the top of your script from arcpy.sa import * and change con for Con. Here is the documentation with some examples . Aug 23, 2019 at 19:53
  • I think that worked but I got another error,Runtime error Traceback (most recent call last): File "<string>", line 19, in <module> File "C:\Python27\ArcGIS10.6\Lib\ntpath.py", line 67, in join p_drive, p_path = splitdrive(p) File "C:\Python27\ArcGIS10.6\Lib\ntpath.py", line 115, in splitdrive if len(p) > 1: TypeError: object of type 'Raster' has no len()
    – seak23
    Aug 23, 2019 at 20:02
  • Check the edit. It's probably because I used a Raster object instead of a string as the name of the resulting raster. Now it will match the name of the corresponding cvr raster. You can change that, of course. Aug 23, 2019 at 20:04
  • Sorry, I'm very new to this. Is a Raster object the automatic name given to a raster when you perform some sort of function? Does this mean I should change the cvr_ras and lss_ras names, or the cvr_rasters and lss_rasters?
    – seak23
    Aug 23, 2019 at 20:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.