rasterio - mask raster against shape, with all pixels completely inside shape

Question

I'm cutting a tiff file against a given shape, to filter the pixels inside.

I'm using rasterio to do it like this:

mask = rasterio.mask.mask(tiff, geoms, crop=True, all_touched=False)

This works but doesn't yield exactly what I want (which is: give me only the pixels that are completely inside the shape).

see bigger picture

So, what all_touched=False does is, according to the docs, that it restricts the criteria from all_touched=True (fig. 1) by only including pixels that have its center inside the shape (fig. 2).

How can I select only the pixels that are completely inside the shape (like in fig. 3)?

Answer 1

You can create an inverted polygon mask, mask with all_touched=True to find the pixels you dont want. Then overlay this with your raster using np.where to extrakt the pixels that are fully inside your shapes.

import numpy as np
import geopandas as gpd
import rasterio
from rasterio import mask, plot
import matplotlib.pyplot as plt
from shapely.geometry import box

#Read raster
raster_file = r"C:\Users\bera\Desktop\gistest\nmd_sample2.tif"
dataset = rasterio.open(raster_file)
arr = dataset.read(1)
transform = dataset.transform

#Read vector geometries
df = gpd.read_file(r"C:/Users/bera/Desktop/gistest/mask_shapefile.shp")
fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(20,7))

#Plot the original data
plot.show(source=arr, ax=ax[0], cmap="tab10", transform=transform)
df.plot(ax=ax[0], color="None", hatch="X", edgecolor="black")
ax[0].set_title("Original data")
ax[0].set_axis_off()

#Create an inverted polygon covering the raster, with holes cut out
#  at the original polygons
extent = box(*dataset.bounds)
inverted_polygon = gpd.GeoSeries(extent, crs=dataset.crs).difference(
    df.dissolve())

#Mask the dataset using the inverted polygon
arr2, trans2 = mask.mask(dataset=dataset, 
                         shapes=inverted_polygon, all_touched=True,
                         filled=False)
plot.show(source=arr2, ax=ax[1], cmap="tab10", transform=transform)
inverted_polygon.plot(ax=ax[1], color="None", hatch="X", edgecolor="black")
ax[1].set_title("Inverted polygon")
ax[1].set_axis_off()


#arr2 above is a numpy masked array. Extract the mask to get a boolean
#  array with the same shape as the original array with value True,
#  where there are pixels not touching the inverted polygon above.
pixels_to_keep = np.ma.getmask(arr2)

#Keep the original array value where the mask is True, else np.nan
pixels_fully_within = np.where(pixels_to_keep, arr, np.nan)
plot.show(source=pixels_fully_within, ax=ax[2], cmap="tab10", transform=transform)
df.plot(ax=ax[2], color="None", hatch="X", edgecolor="black")
ax[2].set_title("Pixels fully within mask polygons")
ax[2].set_axis_off()