Rasterio geometry mask example for Stack Exchange

geometry_mask_example.py

import xarray as xr
import geopandas as gpd
import rasterio

# Open your shapefile and xarray object
ds = raster_mask
gdf = vector_mask

# Select shapefile feature you want to analyse
# and reproject to same CRS as xarray
gdf = gdf.iloc[[0]].to_crs(ds.crs)

# create polygon mask
mask = rasterio.features.geometry_mask(
            gdf.geometry,
            out_shape=ds.geobox.shape,
            transform=ds.geobox.affine,
            all_touched=False,
            invert=False)

mask = xr.DataArray(mask, dims=("y", "x"))

# mask ds with rasterized gdf
ds = ds.where(~mask)

Thanks a lot for this. I get the following error when applying a gdf on a xarray dataset:

AttributeError: 'Dataset' object has no attribute 'geobox'

Is geobox something specific to DigitalEarthAU?

Cheers!

Hi @e5k, yep: a GeoBox is a method that gets added to xarray data by Open Data Cube which allows you to easily obtain info on the spatial grid of the data (e.g. resolution, CRS, transform etc). Luckily there's now a really nice and lightweight odc-geo Python package that recreates the same functionality in a more generic way: try this code and see if it works better for you:

pip install odc-geo

import xarray as xr
import geopandas as gpd
import rasterio
import odc.geo.xr

# Open your shapefile and xarray object
ds = raster_mask
gdf = vector_mask

# Select shapefile feature you want to analyse
# and reproject to same CRS as xarray
gdf = gdf.iloc[[0]].to_crs(ds.odc.geobox.crs)

# create polygon mask
mask = rasterio.features.geometry_mask(
            gdf.geometry,
            out_shape=ds.odc.geobox.shape,
            transform=ds.odc.geobox.affine,
            all_touched=False,
            invert=False)

mask = xr.DataArray(mask, dims=("y", "x"))

# mask ds with rasterized gdf
ds = ds.where(~mask)

Thanks for the answer @robbibt. Unfortunately, using odc.geo.xr does not add the geobox function onto my xarray dataset. Is it ok to simply open it like this or do I have to link it to old somehow?

EVI_post = xr.open_dataset('~/myfile.nc')