Nested conditional statement in ArcGIS on multiple rasters

Question

I have a list of rasters in ArcGIS and I want to create a conditional (Con) nested statement based on all the rasters of the list. I want this to be reiterate through many lists (the number of rasters in the lists can change). So far, I found a non-elegant way which build a string and then pass the string to the exec function, but I read the use of exec is higly discouraged. Any clever way to do it?

listRst = arcpy.ListRasters("*", "TIF")
initString = "outCon = "
for rst in listRst:
 rstString = " Con(( Raster(\"" + rst + "\") > 0) & ( Raster(\"" + rst + "\") < 367), 1,"
 initString = initString + rstString
conString = initString + " 0" + ")"*len(listRst)
exec (conString)
outCon.Save("ConTest.tif")

• Menno: If I understand it correctly, you want to evaluate every raster for if there are values between 0 and 367. After that it's not completely clear for me. Do you want to generate one raster which has value = 1 for where any of the iterated rasters holds a valid value (a large-scale boolean OR), or one raster which has value = 1 for where every of the iterated rasters holds a value value(a large-scale boolean AND)?
• Gianluca Franceschini: I want to evaluate the first raster, if the value is > 0 and < 367, if not I evaluate the second raster, if not the third and so on. If in any of the raster the value is > 0 and < 367, the value of 0 is assigned. Of course I could do single Con statements on each raster and the merge somehow the result, but I am looking for other approaches

Answer 1

You are right that exec is frowned upon and I think you will be better suited to use arrays and numpy to solve your problem than using the raster calculator.

I assume that all of the rasters in your list have the same number of rows and columns because you talk about merging the results at the end. The following will result in:

• 1 for cells where none of the rasters have values between 0 and 367 in the cell
• 0 for cells where at least 1 of the rasters has a value between 0 and 367 in the cell

I am fairly certain that is what you want.

This first section of code gets all of the variables you will initially need:

#set the lower and upper values
lower_break_value = 0
upper_break_value = 367
#get raster list, properties of first raster to create out_raster_data to hold results
list_rasters = arcpy.ListRasters('*', "TIF")
temp_raster = arcpy.Raster(list_rasters[0])
raster_rows = temp_raster.width
raster_columns = temp_raster.height
out_raster_data = numpy.ones((raster_columns, raster_rows), numpy.int) #in your question you only refer to integer values but this can be changed

The following section uses a loop on the rasters in your list to do the following for each raster:

1. create a numpy array from the raster
2. use numpy.where to classify values that are greater than your lower_break_value (0)
3. use numpy.where to classify the values that are less than your upper_break_value (367). numpy.where can only evaluate one condition so we have to break it up into two separate statements.
4. add the two resulting arrays from the numpy.where operations. Anything that meets both conditions will add to 0.
5. create a temp_mask that sets all values to True that do not equal 0 and False that equal 0

6. add temp_mask to a list bands_list

   bands_list = []
   for i in list_rasters:
    temp_array = arcpy.RasterToNumPyArray(i).astype(numpy.int)
    greater_than = numpy.where(temp_array > lower_break_value, 0, temp_array)
    less_than = numpy.where(temp_array < upper_break_value, 0, temp_array)
    temp_result = greater_than + less_than
    temp_mask = temp_result != 0
    bands_list.append(temp_mask)
   

In the final portion of code we loop over the masked arrays in band_list and multiply them by out_raster_data which is just an array of ones. Any cell in any of the arrays\rasters that was between 0 and 367 is a 0, therefore when we multiply it by the other values it will result in a 0 in out_raster_data:

for i in range(0, len(bands_list)):
 out_raster_data *= bands_list[i]
#convert out_raster_data to a raster and save the results
result_raster = arcpy.NumPyArrayToRaster(out_raster_data)
result_raster.save(r'') #enter output raster path here