Masking NetCDF time series data from shapefile using Python

Question

I have a 3-D time-series precipitation data (187 x 1800 x 3600), stored in a NetCDF file. I need to obtain the precipitation data for a shapefile.

import matplotlib.pyplot as plt
from netCDF4 import Dataset, num2date, 
from matplotlib.pyplot import figure
from datetime import datetime, date, timedelta
import numpy as np
import xarray as xr
import pandas as pd
import geopandas as gpd     

MSWEP_monthly = 'D:\G3P\DATA\Models\MSWEP\MSWEP_monthly.nc4'

MSWEP_monthly = Dataset(MSWEP_monthly, 'r')
Pre_MSWEP = MSWEP_monthly.variables['precipitation'][:]


MSWEP_monthly2 = xr.open_dataarray('D:\G3P\DATA\Models\MSWEP\MSWEP_monthly.nc4')

Lon_MSWEP = MSWEP_monthly2.lon
Lat_MSWEP = MSWEP_monthly2.lat

Africa_Shape = gpd.read_file('D:\G3P\DATA\Shapefile\Africa_SHP\Africa.shp')



from osgeo import gdal,osr,ogr

def makeMask(lon,lat,res):
    source_ds = ogr.Open(shapefile)
    source_layer = source_ds.GetLayer()
 
    # Create high res raster in memory
    mem_ds = gdal.GetDriverByName('MEM').Create('', lon.size, lat.size, gdal.GDT_Byte)
    mem_ds.SetGeoTransform((lon.min(), res, 0, lat.max(), 0, -res))
    band = mem_ds.GetRasterBand(1)
 
    # Rasterize shapefile to grid
    gdal.RasterizeLayer(mem_ds, [1], source_layer, burn_values=[1])
 
    # Get rasterized shapefile as numpy array
    array = band.ReadAsArray()
 

    mem_ds = None
    band = None
    return array

shapefile = 'D:\G3P\DATA\Shapefile\Africa_SHP\Africa.shp'
source_ds = ogr.Open(shapefile)

# calculate the cellsize
cellsize = Lon_MSWEP[:][1] - Lon_MSWEP[:][0]
 
# create the mask
mask = makeMask(Lon_MSWEP,Lat_MSWEP,cellsize)

Now if I implement the following code, the precipitation data for the first day (from the time-series) can be obtained with shape of 1800x3600:

precip = np.ma.masked_where(mask==0,Pre_MSWEP[0,:,:])

I tried to use a for loop to mask the precipitation data for the entire time series (time, lon, lat) over the area of interest. however, the below code gives me a 2-D data, probably for the last day.

Why is that?

for i in range(len(Pre_MSWEP)):
    precip = np.ma.masked_where(mask==0,Pre_MSWEP[i,:,:])
    

Answer 1

Here is an example of using rioxarray to mask out data with a shapefile: https://corteva.github.io/rioxarray/stable/examples/clip_geom.html

import geopandas
import rioxarray
import xarray
from shapely.geometry import mapping


MSWEP_monthly2 = xarray.open_dataarray('D:\G3P\DATA\Models\MSWEP\MSWEP_monthly.nc4')
MSWEP_monthly2.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True)
MSWEP_monthly2.rio.write_crs("epsg:4326", inplace=True)
Africa_Shape = geopandas.read_file('D:\G3P\DATA\Shapefile\Africa_SHP\Africa.shp', crs="epsg:4326")

clipped = MSWEP_monthly2.rio.clip(Africa_Shape.geometry.apply(mapping), Africa_Shape.crs, drop=False)