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,:,:])
My problem is here: II 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,:,:])
