I have a NetCDF file with a curvilinear (rotated) grid that contains meteorological data.
I am trying to crop the NetCDF based on a mask shapefile in R, but the rotated grid does project the data in the correct coordinates.
I have tried the following in R:
library(ncdf4) # package for netcdf manipulation library(raster) # package for raster manipulation nc_file <- '2000010112.nc' nc_ras <- brick(nc_file) plot(nc_ras[])
First, when reading the NetCDF file as brick, I encountered the following error
Warning messages: 1: In .getCRSfromGridMap4(atts) : cannot process these parts of the CRS: long_name=coordinates of the rotated North Pole grid_north_pole_latitude=31.7583160400391 grid_north_pole_longitude=87.5970306396484 north_pole_grid_longitude=0 2: In .getCRSfromGridMap4(atts) : cannot create a valid CRS
This is not the correct location of the data as the extents should be
lon : -127.0499 to -87.86075 degrees_east and
lat : 42.03678 to 62.11168 degrees_north.
R uses the
rlat variables (instead of
lat) in the NetCDF file as the coordinates. These extent information were obtained from
CDO using the
sinfon command and the following was the result:
Grid coordinates : 1 : curvilinear : points=42147 (223x189) lon : -127.0499 to -87.86075 degrees_east lat : 42.03678 to 62.11168 degrees_north mapping : rotated_latitude_longitude rlon : 342.1528 to 362.1328 by 0.09000005 degrees rlat : -13.05 to 3.87 by 0.09 degrees
Apparently, I cannot use any shapefile as a mask to crop the data from the file because the shapefile will not be within the same extents as the NetCDF file.
I tried to find the
crs(nc_ras) of the NetCDF file but it returned the following:
CRS arguments: NA
I have also tried the
stars package that handles the curvilinear grid
library(stars) prec = read_ncdf(nc_file, curvilinear = c("lon", "lat"))
but I encountered the following error
Error in UseMethod("GPFN") : no applicable method for 'GPFN' applied to an object of class "rotated_latitude_longitude" Error in .check_curvilinear(c_v, var, meta$variable, curvilinear) : Specified curvilinear coordinates are not 2-dimensional.
I have tried several solutions using command-line tools (
NCL, which are reported on here, here, and other places) to convert the NetCDF file from the rotated grid to regular lonlat grid, but with no luck. I need to preserve the same values on the new grid (i.e., no interpolation or remapping)
So, my question, is there a way to handle this in R? In other words, handle the rotated curvilinear grid and crop data out of it based on another shapefile.