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)

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
2: In .getCRSfromGridMap4(atts) : cannot create a valid CRS

Second, the resulting figure shows that the data were projected to the wrong location enter image description here

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

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.

  • How did you define this index id <- c(16811,18526,157891)? I address this question because I am trying to follow the code in R in oder to convert a rotated netCDF file to regular grid/datafram. But I did not understand how you set the id index which you consider as mask Commented May 15, 2022 at 8:27
  • I have identified that from the coord data.frame. Once you define that, you can use your shapefile and identify/crop the points within it. You can use the cropped points id as your index id to extract the values for. Commented May 16, 2022 at 17:33

1 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


# 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

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)

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.