Obviously it's difficult to reproduce a netCDF here on Stackexchange, especially a curvilinear/rotated grid type.
So let me start with what I already know. I know how to mask a regular netCDF with rectangular/square grid with a shapefile.
For this I use NCL programming and it's routine "shapefile_utils.ncl". This is my code:
;-- copy shapefile_utils.ncl from
;-- https://www.ncl.ucar.edu/Applications/Scripts/shapefile_utils.ncl
load "shapefile_utils.ncl"
shpname = "gadm36_GER_0.shp"
maskname = "germany_mask.nc"
print_shapefile_info(shpname)
;-- open data file to get the grid to be used
f = addfile("Europe.nc","r")
;-- read variable
var = f->tp(0,:,:)
var@lat2d = f->latitude
var@lon2d = f->longitude
;-- 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 = "y"
mask_array!1 = "x"
mask_array@coordinates = "latitude longitude"
;-- create new netCDF file and write mask array
system("rm -f " + maskname)
fout = addfile(maskname,"c")
fout->mask_array = mask_array
fout->latitude = f->latitude
fout->longitude = f->longitude
So in this specific example of NCL code, gadm36_GER_0.shp
is the shapefile of Germany, then the mask nc file I want to create is germany_mask.nc
and the nc file I am going to mask is Europe.nc
.
Then finally in I use CDO in Linux to mask the netCDF.
cdo div inputfile maskfile outputfile
Or in this example it the code would be:
cdo div Europe.nc germany_mask.nc Germany_masked.nc
My question would be: how to do this with a netCDF with a curvilinear/rotated grid.
All I can really do is show you some output of curvilinear/rotated grid nc file in R.
5 variables (excluding dimension variables):
double time_bnds[bnds,time] (Chunking: [2,1])
double lon[rlon,rlat] (Contiguous storage)
standard_name: longitude
long_name: longitude
units: degrees_east
_CoordinateAxisType: Lon
double lat[rlon,rlat] (Contiguous storage)
standard_name: latitude
long_name: latitude
units: degrees_north
_CoordinateAxisType: Lat
char rotated_pole[] (Contiguous storage)
grid_mapping_name: rotated_latitude_longitude
grid_north_pole_latitude: 39.25
grid_north_pole_longitude: -162
long_name: coordinates of the rotated North Pole
float pr[rlon,rlat,time] (Chunking: [415,423,1])
standard_name: precipitation_flux
long_name: Precipitation
units: kg m-2 s-1
grid_mapping: rotated_pole
coordinates: lat lon
_FillValue: 1.00000002004088e+20
missing_value: 1.00000002004088e+20
cell_methods: time: mean
4 dimensions:
time Size:1 *** is unlimited ***
standard_name: time
long_name: time
bounds: time_bnds
units: days since 1969-12-01T00:00:00Z
calendar: proleptic_gregorian
axis: T
bnds Size:2 (no dimvar)
rlon Size:415
standard_name: projection_x_coordinate
long_name: longitude in rotated pole grid
units: degrees
axis: X
rlat Size:423
standard_name: projection_y_coordinate
long_name: latitude in rotated pole grid
units: degrees
axis: Y
When I directly apply the NCL code shown above on the this specific curvilinear/rotated grid netCDF, it doesn't work. It doesn't throw an error, but the entire masked nc file has NA values. It basically masks the entire matrix and not the region I want to mask.
My question: Does anybody have any ideas on how to mask a curvilinear/rotated grid with a shapefile?