Masking NetCDF data using shapefile (xarray & geopandas) [duplicate]

Question

I am using the following packages:

import pandas as pd
import numpy as np
import xarray as xr
import geopandas as gpd

I have the following objects storing data:

print(precip_da)

Out[]:
    <xarray.DataArray 'precip' (time: 13665, latitude: 200, longitude: 220)>
    [601260000 values with dtype=float32]
    Coordinates:
      * longitude  (longitude) float32 35.024994 35.074997 35.125 35.175003 ...
      * latitude   (latitude) float32 5.0249977 5.074997 5.125 5.174999 ...
      * time       (time) datetime64[ns] 1981-01-01 1981-01-02 1981-01-03 ...
    Attributes:
        standard_name:       convective precipitation rate
        long_name:           Climate Hazards group InfraRed Precipitation with St...
        units:               mm/day
        time_step:           day
        geostatial_lat_min:  -50.0
        geostatial_lat_max:  50.0
        geostatial_lon_min:  -180.0
        geostatial_lon_max:  180.0

This looks as follows:

precip_da.mean(dim="time").plot()

I have my shapefile as a geopandas.GeoDataFrame which represents a polygon.

awash = gpd.read_file(shp_dir+"/Export_Output.shp")

awash
Out[]:
  OID_         Name      FolderPath  SymbolID  AltMode Base  Clamped Extruded  Snippet PopupInfo Shape_Leng  Shape_Area  geometry
0     0 Awash_Basin Awash_Basin.kml         0        0  0.0       -1        0     None      None  30.180944    9.411263  POLYGON Z ((41.78939511000004 11.5539922500000...

Which looks as follows:

awash.plot()

Plotted one on top of the other they look like this:

ax = awash.plot(alpha=0.2, color='black')
precip_da.mean(dim="time").plot(ax=ax,zorder=-1)

My question is, how do I mask the xarray.DataArray by checking if the lat-lon points lie INSIDE the shapefile stored as a geopandas.GeoDataFrame?

So I want ONLY the precipitation values (mm/day) which fall INSIDE that shapefile.

I want to do something like the following:

masked_precip = precip_da.within(awash)

OR

masked_precip = precip_da.loc[precip_da.isin(awash)]

EDIT

I have thought about using the rasterio.mask module but I don't know what format the input data needs to be. It sounds as if it does exactly the right thing:

"Creates a masked or filled array using input shapes. Pixels are masked or set to nodata outside the input shapes"

Answer 1

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.apply(mapping), geodf.crs)

Answer 2

It seems that regionmask does what you want.

regionmask is a Python module that:

• contains a number of defined regions, including: countries, a landmask and regions used in the scientific literature.
• can plot figures of these regions with matplotlib and cartopy.
• can be used to create masks of the regions for arbitrary longitude and latitude grids with numpy and xarray
• arbitrary regions can be defined easily

Answer 3

It seems like Overlay from Geopandas should work as well, through intersection http://geopandas.org/set_operations.html https://nbviewer.jupyter.org/github/geopandas/geopandas/blob/master/examples/overlays.ipynb But you need first: 1. Convert your netcdf into a dataframe, 2. convert latitude and longitud into a polygon like in this example, all the way to the end https://medium.com/@Arbolmarket/working-with-geospatial-data-in-python-a5ad984c1161 Or even better, like this answer https://stackoverflow.com/questions/46332479/store-netcdf-data-in-geodataframe