Possible to iterate both input raster and mask in "extract by mask" tool in model builder or ArcPy?

[LiYi1]

Hi, I think what I am trying to do is very straigitforward but I was stuck in two days and have exhausted all previous solutions here and stakcoverflow. So, I am posting here with the hope of getting some help from someone who is good at model builder or arcpy.

Basically, I have a folder of raster grids converted from polygons (N=x)

I have another folder of NDVI rasters (N=y)

I am trying to build a loop extract by mask tool so that the tool will take x*y combinations of raster grids ad the mask and NDVI as raster inputs so that each NDVI will be clipped by the raster grids.

Model builder does not allow two iterators, and the sub-model within the model cannot be used to bypass this since it only takes the last value of the sub-model in the main model. I see someone mention arcpy method but I have zero skill in python so no idea how this can be achieved. I am hoping for any comments on whether this is indeed available.

Any insights are greatly appreciated. 

Thanks a lot.

[DavidPike]

Why did you convert the polygons to rasters?

It seems like it can be done very easily using arcpy/python.  If you give some examples of your folder structures and naming conventions (both sample names and extensions of what you have already, and how you would like the output files named).

[LiYi1]

Hi David, appreciate your reply. I am converting mask poly to raster since I need to have the output raster at the same grid size of the kernel density grid that I later will multiply with raster calculator.

Thanks so much for helping with arcpy, so I have three geodatabases setup called MASK.gdb, INRASTER.gdb, OUTRASTER.gdb within a folder called Expo

and then within MASK.gdb I have MASK1, MASK2 ... stored as raster format

INRASTER i have NDVI1, NDVI2, ... stored as raster form

output OUTRASTER should have name of both MASK1_NDVI1 ... MASK1_NDVI2 ... MASK2_NDVI1...

Looking forward to your reply!

[DavidPike · Accepted solution]

This should do it

import arcpy
import os

# specify paths of your GDBs
inras_ws = r'C:\.......\INRASTER.gdb'
mask_ws = r'C:\.......\MASK.gdb'
outras_ws = r'C:\.......\OUTRASTER.gdb'

# Set the current workspace
arcpy.env.workspace = inras_ws
# create a list of in_raster paths
in_rasters = [os.path.join(inras_ws, raster) for raster in arcpy.ListRasters()]

# Set the current workspace
arcpy.env.workspace = mask_ws
# create a list of mask_raster paths
mask_rasters = [os.path.join(mask_ws, raster) for raster in arcpy.ListRasters()]

#iterate and save outputs
for raster in in_rasters:
    raster_name = os.path.basename(raster)
    for mask_raster in mask_rasters:
        mask_raster_name = os.path.basename(mask_raster)
        out_raster_name = raster_name + "_" + mask_raster_name
        out_path = os.path.join(outras_ws, out_raster_name)
        outExtractByMask = arcpy.sa.ExtractByMask(raster, mask_raster)
        outExtractByMask.save(out_path)

[LiYi1]

David, thanks very much. I modified the folder path and tried use IDLE (ArcGIS Pro) and then F5 run module. The Shell does not give me any error, but neither it does seem running. Do I run this code incorrectly? How should it be run usually? Excuse my rudimentary knowledge with ArcPy. Thanks, Li

[DavidPike]

I can't really say.  I've added some print statements which should help.

import arcpy
import os

# specify paths of your GDBs
inras_ws = r'C:\.......\INRASTER.gdb'
mask_ws = r'C:\.......\MASK.gdb'
outras_ws = r'C:\.......\OUTRASTER.gdb'

# Set the current workspace
arcpy.env.workspace = inras_ws
# create a list of in_raster paths
in_rasters = [os.path.join(inras_ws, raster) for raster in arcpy.ListRasters()]
print(in_rasters)
# Set the current workspace
arcpy.env.workspace = mask_ws
# create a list of mask_raster paths
mask_rasters = [os.path.join(mask_ws, raster) for raster in arcpy.ListRasters()]
print(mask_rasters)
#iterate and save outputs
for raster in in_rasters:
    raster_name = os.path.basename(raster)
    for mask_raster in mask_rasters:
        mask_raster_name = os.path.basename(mask_raster)
        out_raster_name = raster_name + "_" + mask_raster_name
        out_path = os.path.join(outras_ws, out_raster_name)
        print("running extract by mask - " + out_raster_name)
        outExtractByMask = arcpy.sa.ExtractByMask(raster, mask_raster)
        outExtractByMask.save(out_path)
        print("completed extract by mask - " + out_raster_name)

[LiYi1]

Hi David, thanks. I was able to make it work, now I need to take one step forward to use rastercalculator to multiply in_mask with the outextractbymask object. I took a stab here but it does not seem to work. Would you be so kind to look at it and let me know what I did wrong here? Much appreciated.

import arcpy
import os

# specify paths of your GDBs
inras_ws = r'P:\MADRES_GPS\Processed\Exposure\DayGPSExposure\in_raster.gdb'
mask_ws = r'P:\MADRES_GPS\Processed\Exposure\DayGPSExposure\in_mask.gdb'
outras_ws = r'P:\MADRES_GPS\Processed\Exposure\DayGPSExposure\out_raster.gdb'

# Set the current workspace
arcpy.env.workspace = inras_ws
# create a list of in_raster paths
in_rasters = [os.path.join(inras_ws, raster) for raster in arcpy.ListRasters()]

# Set the current workspace
arcpy.env.workspace = mask_ws
# create a list of mask_raster paths
mask_rasters = [os.path.join(mask_ws, raster) for raster in arcpy.ListRasters()]

#iterate and save outputs
for raster in in_rasters:
    raster_name = os.path.basename(raster)
    for mask_raster in mask_rasters:
        mask_raster_name = os.path.basename(mask_raster)
        out_raster_name =  mask_raster_name + "_" + raster_name 
        out_path = os.path.join(outras_ws, out_raster_name)
        print("running extract by mask - " + out_raster_name)
        outExtractByMask = arcpy.sa.ExtractByMask(raster, mask_raster)
        print("completed extract by mask - " + out_raster_name)
        outRC=RasterCalculator([outExtractByMask, mask_raster], ["x", "y"],
                               "x*y")
        outRC.save(out_path)
        print("completed raster calculator - " + out_raster_name)

[DavidPike]

I'd probably just do algebra on the Raster objects.  outExtractByMask should already be a raster object, however mask raster will need to be cast as one with arcpy.raster(...)

outRC = outExtractByMask ** arcpy.Raster(mask_raster)

[LiYi1]

I see, so in this way, if I wish to scale it to 0.0001 I just do the below right?

outRC = outExtractByMask ** arcpy.Raster(mask_raster) ** 0.0001

[DavidPike]

That's the badger.