CDO

Hi Bhandari,

the problem is the wrong order of the dimensions in variable precip(time,lon,lat) which should be precip(time,lat,lon). This can be corrected in the NCL script, too.
Additional, I've added the cdo command to change the NaNs to a missing value for the variable attributes missing_value and _FillValue.

load "/usr/lib/ncarg/nclscripts/csm/shapefile_utils.ncl" 

shpname  = "gcs30mbbuf.shp" 
maskname = "Mahi11.nc" 

print_shapefile_info(shpname)

;-- open data file to get the grid to be used
infile = "MahiBasin_Box_India_Projected_Rainfall_cccma-CanESM2_historical_2001_2005.nc" 

;-- set NaN values to missing and _FillValue
cdo_cmd = "cdo -setattribute,precip@missing_value=1e20 -setmissval,1e20 -setmissval,nan " + \
           infile + " tmp.nc" 
system(cdo_cmd)

;-- read variable
f = addfile("tmp.nc", "r")

var  = f->precip(0,:,:)
var := var(lat|:, lon|:)    ;-- reorder dimensions
var@lat2d = f->lat
var@lon2d = f->lon

;-- shapefile mask resources
opt = True
opt@return_mask = True      ;-- this forces the return of a 0s and 1s mask array

;-- create the mask based on the given shapefile
mask_array = shapefile_mask_data(var, shpname, opt)
mask_array!0 = "lat" 
mask_array!1 = "lon" 

;-- create new netCDF file and write mask array
system("rm -f " + maskname)
fout = addfile(maskname,"c")

fout->mask_array = mask_array
fout->lat        = f->lat
fout->lon        = f->lon

system("rm tmp.nc")

-Karin

Hi Karin

I followed your NCL script, however, the mask file I created has only 0 value. Can you please take a look at my script? Thanks a lot!

load "shapefile_utils.ncl" 

shpname  = "./CHN_adm0.shp" 
maskname = "./CHN_mask.nc" 

;-- open data file to get the grid to be used
f = addfile("mask.nc","r")

;-- read variable
var                    =  f->T_2M(0,:,:)
var@lat2d              =  f->lat
var@lon2d              =  f->lon

;-- shapefile mask resources
opt             =  True
opt@return_mask =  True    

;-- create the mask based on the given shapefile
mask_array             =  shapefile_mask_data(var, shpname, opt)
mask_array!0           = "rlat" 
mask_array!1           = "rlon" 
mask_array@coordinates = "lat lon" 

;-- create new netCDF file and write mask array
system("rm -f " + maskname)
fout = addfile(maskname,"c")

fout->mask_array =  mask_array
fout->lat      =  f->lat
fout->lon      =  f->lon

Hi Ruolan,

I'm able to reproduce the problem but have not been able to solve it yet. The problem is the rotated curvilinear grid. I'll see if I can write a workaround and get back to you as soon as possible.

-Karin