Mask xarray dataset using a shapefile [closed]

Question

I have used xarray on some satellite data that I'd like to mask it using a shapefile. I have previously been doing this using a combination of PIL, numpy and gdal, but it would be cleaner to do this directly on the xarray dataset with Rasterio and GeoPandas.

I have the following xarray Dataset:

Which looks as follows:

I've opened my shapefile in GeoPandas using:

sf = geopandas.read_file('GRI_jergetal.shp')

I then want to mask my xarray dataset using the shapefile in GeoPandas, and it seems like rasterio.mask.mask has this capability. However, the code that I have used gives the following error (adding in a nodata and a transform attribute to my xarray as it did not contain those parameters originally).

sf = sf.to_crs('EPSG:32643')
ndvi.attrs['nodata'] = np.nan
ndvi.attrs['transform'] = good_ds.affine
rasterio.mask.mask(ndvi.isel(time=0), sf.loc[0, 'geometry'])

[Output]:
`TypeError: 'Polygon' object is not iterable`.

I've tried a method posted in a previous post Python mask NetCDF data using shapefile (xarray & GeoPandas) & Mask Rasterio raster with GeoPandas shapefile, but the method doesn't seem right for my present case. Have I made a mistake in invoking the rasterio.mask.mask() function or should I try a different method for getting this to work?

Answer 1

Masking NetCDF time series data from shapefile using Python

You can use rioxarray. Here is an example: https://corteva.github.io/rioxarray/stable/examples/clip_geom.html

import rioxarray
import geopandas

geodf = geopandas.read_file(...)
xds = rioxarray.open_rasterio(...)
clipped = xds.rio.clip(geodf.geometry.values, geodf.crs)

Answer 2

You can to use the rasterio.features.geometry_mask function to do this:

ShapeMask = rasterio.features.geometry_mask(sf.iloc[0],
                                      out_shape=(len(ndvi.y), len(ndvi.x)),
                                      transform=ndvi.geobox.transform,
                                      invert=True)
ShapeMask = xr.DataArray(ShapeMask , dims=("y", "x"))

# Then apply the mask
NDVImasked = ndvi.where(ShapeMask == True)

I think from your code snippet above sf is a geopandas dataframe. If that's the case, then you need to grab just the polygon you want to apply as the mask (as per sf.iloc[0] above).

Answer 3

You can try using Digital Earth Australia's xr_rasterize function to convert your geopandas geodataframe into an xarray object, and then use xarray's .where() method to mask you're array.

e.g.

mask = xr_rasterize(gdf, da)
masked_da = da.where(mask)

If you would prefer to use rasterio.features.geometry_mask, then the following code should work. You may need to adjust how you grab the transform information if your dataset does not have a ds.geobox attribute (an open data cube query object will have this attribute)

import xarray as xr
import geopandas as gpd

#open your shapefile and xarray object
ds = xr.open_dataset('your_dataset')
gdf = gpd.read_file('your_shape.shp')

#convert the geometry of a single polygon in your gdf (adjust the index, [0], to match the row of the dataframe you care about)
geom = geometry.Geometry(
            gdf.geometry.values[0].__geo_interface__, geometry.CRS(
                'epsg:4326'))

# create polygon mask
mask = rasterio.features.geometry_mask(
            [geom.to_crs(ds.geobox.crs) for geoms in [geom]],
            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 == False)