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
Solved! Go to Solution.
Here is my solution to create a customized mask from raster, that covers the entire surface of interest:
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)
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
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?
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?
Here is my solution to create a customized mask from raster, that covers the entire surface of interest:
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)