4

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
8

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)
6
  • 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 Mar 23 '20 at 11:41
  • The CRS passed in should be the CRS of the input shapefile. Is the CRS of the shapefile missing? – snowman2 Mar 23 '20 at 13:15
  • Yes, the CRS of the shapefile is missing. – Ehsan Mar 23 '20 at 14:10
  • Ah, got it. I updated the geopandas read_file command so it should have the CRS set. – snowman2 Mar 23 '20 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 Jul 7 '20 at 5:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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