Handle Curvilinear/rotated grid NetCDF file in R

Question

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 not 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[[1]])

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

Second, the resulting figure shows that the data were projected to the wrong location

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 rlon and rlat variables (instead of lon and 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 (CDO, NCO, 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.

Answer 1

After searching for a while, I found two possible methods to handle this problem.

Method 1

Using GDAL to reproject the data, using lat and lon variables, to a 2D regular grid. This method was mentioned before in here and the following command can be used:

gdalwarp -geoloc NETCDF:"infile.nc":required_variable out_proj.tif

The output file is actually a multiband raster that can be imported to any GIS software and will project properly. However, this solution removes the meta data of the NetCDF file. In other words, there will be no information related to the variable and its units, time steps, levels, etc.

---

Method 2

NetCDF files on a rotated grid include lat and lon as variables (not dimensions to map the file). The 2D lat and lon variables can be retrieved from the file. Then, they can be used to get the index of the different cells within the mask layer (using either R or GIS), to extract the values for those cells using R or any command-line tool. The following is a sample R script to handle this method

library('raster')
library('ncdf4')
library('ncdf4.helpers')
library('maptools')
nc file name
nc_file <- 'infile.nc'
get variables from the NetCDF file in R
nc_fo <- nc_open(nc_file)
names(nc_fo$dim) #display dimensions
names(nc_fo$var) #display variables
req_var <- names(nc_fo$var)[3] #precip, temp, pressure, etc.
lat <- ncvar_get(nc_fo, "lat") #2D lat values
lon <- ncvar_get(nc_fo, "lon") #2D lon values
ts <- as.POSIXct(nc.get.time.series(nc_fo))
Read the values of the variable at the first time step over all grid cells to be used for converting linear index to 2D (row and column) index.
Count of -1 mean all values along that dimension
dum_var <- ncvar_get(nc_fo, req_var,start=c(1,1,1),count = c(-1,-1,1)) 
#close the nc file
nc_close(nc_fo)

Now, the 2D lon and lat values can be exported as a csv file and imported to any GIS software or can be converted to a spatial point data frame in R.

coord <- data.frame(id=1:length(lon),lon=as.vector(lon), lat=as.vector(lat))
write.csv(coord, file= 'coordinates.csv')

Note, the coord$id is the linear index of the grid cells.

The cells and their index that lie in the mask shapefile can be identified from any GIS software. The index can be used to extract the data in either R or CDO.

An R example

nc_fo <- nc_open(nc_file)
data_out <- rep(0,length(ts)) #create zeros vector to output the average value
id <- c(16811,18526,157891) # index of cells within the mask shapefile
for (i in 1:length(id)) {
pt_id <- id[i]
  #convert linear ind to r and c index
  rc <- arrayInd(pt_id,dim(dum_var)) #get row and column index
#get the variable at the required location for all time steps.
pt_data <- ncvar_get(nc_fo, req_var,

                     start=c(rc[1],rc[2],1),count=c(1,1,-1)) #-1 read all time_steps
  data_out <- data_out+pt_data #sum data to be averaged after the loop
}
average value spatially
data_out <- data_out/length(id)
nc_close(nc_fo)

A CDO example using id=c(16811,18526,157891)

cdo -O selgridcell,16811,18526,157891 infile.nc extracted.nc #extract data for specific indexes
cdo fldmean,weights=FALSE extracted.nc avg_mask.nc #get the zonal average over the area

---

Other Methods

Other possible solutions to this problem are listed here