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
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?