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,:,:])
  • 1
    I too am facing a similar problem, detailed on Stackoverflow - link Commented Mar 12, 2021 at 18:31
  • The short answer is, because when you do Pre_MSWEP[i,:,:], you extract 2-D data, and mask is 2-D. A guess (without having tried the code), is that you could remove the loop and instead effectively making your mask 3-D by creating a new dimension with length 1 and rely on broadcasting of the mask along that dimension, by doing precip = np.ma.masked_where(np.expand_dims(mask, axis=0) == 0, Pre_MSWEP). Commented Dec 15, 2022 at 13:25
  • If that doesn't work or the mask istn't broadcasted, you may also have to repeat the mask manually: precip = np.ma.masked_where(np.repeat(np.expand_dims(mask, axis=0) == 0, Pre_MSWEP.shape[0], axis=0), Pre_MSWEP). Commented Dec 15, 2022 at 13:34
  • 1
    In general, using a variable that is created inside a loop, after the loop has finished, may indicate that you have made a mistake (unless your intention is to use the value the variable is assigned in the very last iteration of that loop, and you know that your loop will iterate at least once). In many other programming languages (in particular statically typed), such variables would go out of scope in the end of each iteration and thus not be accessible after the loop. Commented Dec 15, 2022 at 13:37

1 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)
  • Dear snowman2, Thanks for your answer. That works like a charm, very straightforward and super-fast, but I think you need to update your answer as follow, then I can confirm it as the correct answer. clipped = MSWEP_monthly2.rio.clip(Africa_Shape.geometry.apply(mapping), MSWEP_monthly2.rio.crs, drop=False)
    – Ehsan
    Commented Mar 23, 2020 at 11:41
  • The CRS passed in should be the CRS of the input shapefile. Is the CRS of the shapefile missing?
    – snowman2
    Commented Mar 23, 2020 at 13:15
  • Yes, the CRS of the shapefile is missing.
    – Ehsan
    Commented Mar 23, 2020 at 14:10
  • Ah, got it. I updated the geopandas read_file command so it should have the CRS set.
    – snowman2
    Commented Mar 23, 2020 at 14:44
  • I used your script to mask nc files from CHIRPS: chirps-v2.0.1981.days_p05.nc Source: data.chc.ucsb.edu/products/CHIRPS-2.0/global_daily/netcdf/p05 using shapefile. But experience an error: ``` File "Clip_NetCDF_with_SHP.py", line 10 CHIRPS_daily = xarray.open_dataarray('Z:\Temp\CHIRPS\NC\chirps-v2.0.1981.days_p05.nc') ^ SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 14-15: malformed \N character escape```
    – user97103
    Commented Jul 7, 2020 at 5:11

Not the answer you're looking for? Browse other questions tagged or ask your own question.