I need to multiply multiple rasters in a folder with their corresponding rasters in another folder.

2445
48
03-21-2019 12:12 AM
LochlanJones
New Contributor III

I have one folder (folder 1) with 200 or more rasters and another folder (folder 2) with the same number. The rasters in folder 1 are all named genus_species_01 (e.g. acan_rufo_01, lito_bico_01, and so on). The rasters in folder 3 are the same but end with _03 instead. Each raster needs to be multiplied with it's corresponding raster from the other folder and the output must retain the species name and be saved in another folder. I am completely inexperienced when it comes to Python but I believe it is the only way to achieve this task.

The link below describes a possible solution to my problem only I am not proficient enough at formatting code (or understanding it) to be able to adapt it to my problem.

python - How to use the math tool ("times") in ArcGIS modelbuilder to multiply two raster sets while... 

I understand what most steps do but in some cases I'm unsure of what I need to change. In particular some of it would need to be changed based on the fact my rasters are located in separate folders rather than all within the same geodatabase as it is in the linked solution. Another part I am unsure about is prefix = r.split("_")[0] + "_" . To my understanding, this will split the raster name at the '_' character in order to retrieve the identifying first part of the name. But as far as I know the [0] refers to the first file in the specified workspace. So if my folder has over 200 rasters would I need to include r.split("_")[0] through to [200]?

Another potential issue is that I need both the abbreviate genus and species so it can only split the name at the second '_' and not the first.

Is anyone able to offer some guidance?

 

48 Replies
DanPatterson_Retired
MVP Esteemed Contributor

logistically, I would suggest collecting the raster names from each of the two folders first, make the pairing and do the math

flder0 = ['a_x01', 'a_x02', 'a_x03', 'a_x04']
flder1 = ['b_c_y01', 'b_c_y02', 'b_yc_03', 'b_yc_04']
xs = [i.rsplit('_')[1] for i in flder0]
ys = [i.rsplit('_')[-1] for i in flder1]
xsys = list(zip(xs, ys))
names = ["{}{}".format(*i) for i in xsys]

for i, nme in enumerate(names):
eq = "{} = {} * {}".format(*[nme, flder0[i], flder1[i]])
print(eq)

x01y01 = a_x01 * b_c_y01
x02y02 = a_x02 * b_c_y02
x0303 = a_x03 * b_yc_03
x0404 = a_x04 * b_yc_04‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Lines 1 and 2, would be a arcpy.ListRasters on each folder

Lines 3 and 4 would simply be splitting the names into their bits to 

line 5 is just part of the process of formatting the output names

When assembled, lines 8 and on are a 'fake' run of the process.

0 Kudos
JohannesBierer
Regular Contributor

If I do get it right, the name difference between the two raster to chose is name_01 and name_03 ?

Maybe this script can help you?

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

arcpy.env.overwriteOutput = True

# Paths to your rasterdata ?
arcpy.env.workspace = r'C:\pathToRaster1'
inws = arcpy.env.workspace
inws1 = r'C:\pathToRaster2'
outws = r'C:\outRaster'

# List all rasters
rasters = arcpy.ListRasters()

# 2) Loop through your raster list
for r in rasters:

prefix = os.path.splitext(r)[0]
suffix = os.path.splitext(r)[1]

prefix1 = prefix.rstrip("1")

raster2 = os.path.join(inws1, prefix1 + "3" + suffix)

print raster2

outraster = os.path.join(outws, prefix1 + "4" + ".tif") # Change the number

print outraster

# perform the raster algebra
mycalc = arcpy.sa.Raster(r) * arcpy.sa.Raster(raster2)

# Save to disk (write to .tif)
mycalc.save(outraster)

arcpy.CheckInExtension ("Spatial")
AnnaKreij
New Contributor III

I have the same challenge as Lochlan Jones, but I am attempting to iterate through and multiply two sets of corresponding rasters, named 'date_ndvi' (e.g. '20050122_ndvi') and 'date_combined'  (e.g.'20050122_combined'). Johannes Bierer's script works for me up until the raster algebra step, where I get Runtime error message:

Traceback (most recent call last):
File "<string>", line 1, in <module>
RuntimeError: ERROR 000732: Input Raster: Dataset D:\1_USGS\6_2005\NDVI_Clips_scale\clouds_null\20050318_ndvi3.tif does not exist or is not supported

It is the correct path, and the raster exist. Would anyone know where I might have gone wrong?

Very grateful for any advice, as I have not found a solution for this anywhere else. 

0 Kudos
DanPatterson_Retired
MVP Esteemed Contributor

Anna...

The errors can vary depending on your version of python and how programs trap errors.... but your file paths aren't the greatest. 

You need to 'raw' encode them, use double backslashes or use forward slashes. 

For future reference, in the world of python 3 and Unicode, this will become a more important issue.

p = "D:\1_USGS\6_2005\NDVI_Clips_scale\clouds_null\20050318_ndvi3.tif"
File "<ipython-input-3-6b3231a93a98>", line 1
p = "D:\1_USGS\6_2005\NDVI_Clips_scale\clouds_null\20050318_ndvi3.tif"
^
SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 16-17: malformed \N character escape

# ---- use 'raw' encoding for your paths

p = r"D:\1_USGS\6_2005\NDVI_Clips_scale\clouds_null\20050318_ndvi3.tif"

print(p)
D:\1_USGS\6_2005\NDVI_Clips_scale\clouds_null\20050318_ndvi3.tif
AnnaKreij
New Contributor III

Hi Dan Patterson‌, thank you for this! I am a complete novice to python, so I really appreciate being pulled up on these things early on. I will raw encode and try again! 

0 Kudos
LochlanJones
New Contributor III

Hi Anna, could you provide your code that you used? I can’t remember right now and I don’t want to make promises but I think I encountered a similar issue and I can look at it when I get home later today (approx 9 hours from now).

Also is the error file the first one in the list?

Get Outlook for iOS<https://aka.ms/o0ukef>

0 Kudos
AnnaKreij
New Contributor III

Thank you Lochlan Jones‌! Dan (see above) highlighted that my file paths were not great and needed recoding, so I will do that first to see if it changes anything. It was the first error in the list, but I would be very interested to see if the script worked for you in the end. It has been very hard to get online support on this, so I really appreciate your efforts! I'll post my code with updated file paths and result in the next hour.

0 Kudos
AnnaKreij
New Contributor III

I have not had a chance to raw encode my file paths, so this may still be the issue, but if anyone can see other errors, please let me know. I had some indent errors, but when fixed the code ran until error 00732. The files ndvi_3 and ndvi_4 has also not been printed - i.e I cannot see the output in explorer/ArcCatalog... Again, I'm very new to any programming, so happy to take any feedback.

>>> import arcpy, os
>>> from arcpy.sa import*
>>> arcpy.CheckExtension("Spatial")
u'Available'
>>> arcpy.env.overwriteOutput = True
>>> arcpy.env.workspace = r 'D:\1_USGS\6_2005\NDVI_Clips_scale'
Parsing error SyntaxError: invalid syntax (line 1)
>>> arcpy.env.workspace = r'D:\1_USGS\6_2005\NDVI_Clips_scale\ndvi'
>>> inws = arcpy.env.workspace
>>> inws1 = r'D:\1_USGS\6_2005\NDVI_Clips_scale\clouds_null'
>>> outws = r'D:\1_USGS\6_2005\NDVI_outputs'
>>>
>>> rasters = arcpy.ListRasters()
>>> for r in rasters:
... prefix = os.path.splitext(r)[0]
... suffix = os.path.splitext(r)[1]
...
>>> prefix1 = prefix.rstrip("1")
...
Parsing error IndentationError: unexpected indent (line 1)
>>> prefix1 = prefix.rstrip("1")
>>> raster2 = os.path.join(inws1, prefix1 + "3" + suffix)
>>> print raster2
D:\1_USGS\6_2005\NDVI_Clips_scale\clouds_null\20050318_ndvi3.tif
>>> outraster = os.path.join(outws, prefix1 + "4" + ".tif")
>>> print outraster
D:\1_USGS\6_2005\NDVI_outputs\20050318_ndvi4.tif
>>> mycalc = arcpy.sa.Raster(r) * arcpy.sa.Raster(raster2)
Runtime error
Traceback (most recent call last):
File "<string>", line 1, in <module>
RuntimeError: ERROR 000732: Input Raster: Dataset D:\1_USGS\6_2005\NDVI_Clips_scale\clouds_null\20050318_ndvi3.tif does not exist or is not supported
0 Kudos
LochlanJones
New Contributor III

Hi Anna, please see my (hopefully helpful) explanation of the code I successfully used below with my data. Something of note is that how your files are named is extremely important to this code working which is why it is not working for you.
For my answer I'll assume your first raster set is named date_ndvi.tif and your second raster set is named date_combined.tif. Because I don't know what you want your output rasters to be called I will assume it will be date_multiply.tif.
Please read my hastily put together explanation in the attached image below before continuing.


Okay, now make sure you have the correct filepaths specified in your code and then replace the following lines:

  • Replace prefix1 = prefix.rstrip("1") with prefix1 = prefix.rstrip("ndvi")
  • Replace raster2 = os.path.join(inws1, prefix1 + "3" + suffix) with raster2 = os.path.join(inws1, prefix1 + "combined" + suffix)
  • Replace outraster = os.path.join(outws, prefix1 + "4") with outraster = os.path.join(outws, prefix1 + "multiply" + suffix) 

Now it should hopefully work. 
The code wasn't working originally because it couldn't find your second set of rasters. In the raster2 line you added a "3" to your file name and then it couldn't find any rasters with that name. This is why it says 20050318_ndvi3.tif doesn't exist because you accidentally added a 3 to it.

I hope this helps. Let me know how it goes. I'm a beginner as well and I only know what this code means because I worked backwards and looked at it piece by piece after the code worked after having it given to me by the good samaritan Johannes Bierer.