ExtractByMask - Script to include pixels whose center falls outside of mask

259
5
Jump to solution
2 weeks ago
BenjaminBigarre
New Contributor II

Hi, 

I would like to know if there is a way to extract all pixels of a raster touched by a mask, not just the pixels whose center falls within the mask, even if it requires some scripting. 

I have been trying to clip a raster with a mask (tried raster and feature), however some pixels are excluded from the output since their center fall outside of the mask (see picture). I cannot change the pixel size of the input raster. 

I know gdalwarp has an option to keep all pixels that touch the mask, is there a way to code something similar with arcpy?

I attached my code for reference. 

Benjamin

0 Kudos
1 Solution

Accepted Solutions
BenjaminBigarre
New Contributor II

Here is my solution to create a customized mask from raster, that covers the entire surface of interest: 

  1. Convert raster to points
  2. Identify points of cells overlapping my surface of interest
  3. Select points, export to a table
  4. Use a custom function to extract cells of raster with list of points, and create the raster mask. 

Below is my code, screenshot of final result on arcgis is attached

import arcpy
from arcpy import env as e,SpatialReference
import os
from arcpy import da as d
from arcpy.sa import *
import pandas as pd
import numpy as np
from arcgis.features import *

projGDB = r"C:\Users\bigar\Documents\ArcGIS\Projects\MooreaHydro2024\MooreaHydro2024.gdb"
e.workspace = r"C:\Users\bigar\Documents\ArcGIS\Projects\MooreaHydro2024\MooreaHydro2024.gdb"
e.overwriteOutput = True
pd.set_option('display.max_rows', None)

def createmyMask(inRaster,pointlist):
    ## Creates a raster mask by extracting raster cells by coordinates.
    desc=arcpy.Describe(inRaster)
    rastermask=ExtractByPoints(in_raster=inRaster, points=pointlist,extraction_area="INSIDE")
    rastermask.save("C:/Users/bigar/Documents/ArcGIS/Projects/MooreaHydro2024/MooreaHydro2024.gdb/"+desc.baseName+"Mask")
    desc2=arcpy.Describe(rastermask)
    print("Mask {} was created".format(desc2.baseName))

def extractXY(gdbTable):
    df = pd.DataFrame(GeoAccessor.from_table(gdbTable))
    # df = df.astype({'X': int, 'Y':int})
    coordlist=[]
    for index, row in df.loc[:, ['X', 'Y']].iterrows():
        # print(type(int(row['X'])), type(int(row['Y'])))
        coordlist.append(arcpy.Point(int(row['X']),int(row['Y'])))
        print(coordlist)
    return coordlist

if __name__ == "__main__":
    XYlist = extractXY(r"C:\Users\bigar\Documents\ArcGIS\Projects\MooreaHydro2024\MooreaHydro2024.gdb\Opunohu_MODIS_CELL_COORD")

    createmyMask("C:/Users/bigar/Documents/ArcGIS/Projects/MooreaHydro2024/MooreaHydro2024.gdb/MOD16A3GF_061_ET_500m_doy2004001_aid0001_reproj", XYlist)

 

View solution in original post

5 Replies
DanPatterson
MVP Esteemed Contributor

It could only be done if the mask was 'buffered' using the Expand tool, expanding by one cell, however, it doesn't do it selectively in preferred areas

Alternately,

Boundary Clean (Spatial Analyst)—ArcGIS Pro | Documentation

maybe you could clean up the mask or the input data a bit


... sort of retired...
MarkBoucher
Occasional Contributor III

If you convert the polygon mask to a raster (with origin and cell size the same as the raster you want to extract) can you use that raster as a mask?  Would it include the all the edge cells you want?

BenjaminBigarre
New Contributor II

I resampled it to the same cell size as the input raster (500m), and I have the same issue. Resampling only keeps pixels whose center is within the original surface. 

btw I have more than a hundred of rasters to clip with same mask. 

I think a way for me to solve my issue could be to extract the pixels I want from the original raster by XY coordinates, then use that as a raster mask to clip the other rasters. Does it make sense?

0 Kudos
MarkBoucher
Occasional Contributor III
That makes sense. Also, if you convert the points to rasters then you would think that they would line up with the raster cells and be on the raster cells you are interested in. If you then converted those back to points, you would think they would be centered on the raster cells (if you needed to use the points to extract info from the raster).
BenjaminBigarre
New Contributor II

Here is my solution to create a customized mask from raster, that covers the entire surface of interest: 

  1. Convert raster to points
  2. Identify points of cells overlapping my surface of interest
  3. Select points, export to a table
  4. Use a custom function to extract cells of raster with list of points, and create the raster mask. 

Below is my code, screenshot of final result on arcgis is attached

import arcpy
from arcpy import env as e,SpatialReference
import os
from arcpy import da as d
from arcpy.sa import *
import pandas as pd
import numpy as np
from arcgis.features import *

projGDB = r"C:\Users\bigar\Documents\ArcGIS\Projects\MooreaHydro2024\MooreaHydro2024.gdb"
e.workspace = r"C:\Users\bigar\Documents\ArcGIS\Projects\MooreaHydro2024\MooreaHydro2024.gdb"
e.overwriteOutput = True
pd.set_option('display.max_rows', None)

def createmyMask(inRaster,pointlist):
    ## Creates a raster mask by extracting raster cells by coordinates.
    desc=arcpy.Describe(inRaster)
    rastermask=ExtractByPoints(in_raster=inRaster, points=pointlist,extraction_area="INSIDE")
    rastermask.save("C:/Users/bigar/Documents/ArcGIS/Projects/MooreaHydro2024/MooreaHydro2024.gdb/"+desc.baseName+"Mask")
    desc2=arcpy.Describe(rastermask)
    print("Mask {} was created".format(desc2.baseName))

def extractXY(gdbTable):
    df = pd.DataFrame(GeoAccessor.from_table(gdbTable))
    # df = df.astype({'X': int, 'Y':int})
    coordlist=[]
    for index, row in df.loc[:, ['X', 'Y']].iterrows():
        # print(type(int(row['X'])), type(int(row['Y'])))
        coordlist.append(arcpy.Point(int(row['X']),int(row['Y'])))
        print(coordlist)
    return coordlist

if __name__ == "__main__":
    XYlist = extractXY(r"C:\Users\bigar\Documents\ArcGIS\Projects\MooreaHydro2024\MooreaHydro2024.gdb\Opunohu_MODIS_CELL_COORD")

    createmyMask("C:/Users/bigar/Documents/ArcGIS/Projects/MooreaHydro2024/MooreaHydro2024.gdb/MOD16A3GF_061_ET_500m_doy2004001_aid0001_reproj", XYlist)