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

1334
15
Jump to solution
06-19-2021 04:56 PM
Labels (3)
LiYi1
by
New Contributor

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.

0 Kudos
1 Solution

Accepted Solutions
DavidPike
MVP Frequent Contributor

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)

View solution in original post

15 Replies
DavidPike
MVP Frequent Contributor

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

0 Kudos
LiYi1
by
New Contributor

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!

0 Kudos
DavidPike
MVP Frequent Contributor

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
by
New Contributor

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

0 Kudos
DavidPike
MVP Frequent Contributor

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)
0 Kudos
LiYi1
by
New Contributor

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)
0 Kudos
DavidPike
MVP Frequent Contributor

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)

 

0 Kudos
LiYi1
by
New Contributor

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

 

0 Kudos
DavidPike
MVP Frequent Contributor

That's the badger.

0 Kudos