1

I'm attempting to read a netCDF file from ERA5 using the raster package in R. An example file is downloadable here.

As ERA5 detail here,

global ERA5 data was interpolated to a regular lat/lon grid with resolution r = 0.25 deg. The 'top left' ERA5 data point is always at Longitude=0 ; Latitude=90, with further grid points spaced by r, and the 'bottom right' grid point at Longitude=360-r ; Latitude=-90. Many software applications by default visualise regularly spaced data as a continuous tiled surface, as in (b). If you use this visualisation, think of coordinates as referencing the centroids of the tiles.

So essentially, ERA5 data is created like this: enter image description here

But the raster package in R reads a 0.25° data file as this:

library(raster)
stack(sample_ERA5.nc)

enter image description here

And therefore gives me this output:

class      : RasterStack 
dimensions : 721, 1440, 1038240, 249  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : -0.125, 359.875, -90.125, 90.125  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 

Note the wrong extent and the 721 rows, where there should only be 720. Of course I could resample() it to the right extent and resolution, but this distorts data, which is not what I want to do!

Any ideas on how to create an R raster of this data? Perhaps via loading them as points and then rasterizing it? I definitely want to keep the data in it's "original form".

2
  • 721 seems correct, its the length of the sequence from -90 to +90 by 0.25.
    – Spacedman
    Apr 4, 2022 at 20:36
  • not resample, just do extent(r) <- extent(0, 360, -90, 90) - might look at exactly why this nc gets interpreted this way, but netcdf is a huge ball of fun with infinite possibilities so it's not unsurprising
    – mdsumner
    Apr 5, 2022 at 7:03

2 Answers 2

1

I don't see any problems here. The extent as a raster is according to the interpretation in the ERA5 docs, such that the extent overlaps on the left and at the top and bottom, but doesn't reach +360, with the overlaps being all 1/2 the grid size:

> era = raster("./sample_ERA5.nc")
> extent(era)
class      : Extent 
xmin       : -0.125 
xmax       : 359.875 
ymin       : -90.125 
ymax       : 90.125 

If you want this as points, then rasterToPoints returns the point-based interpretation, with the coordinates being the correct rounded coordinates given in the ERA5 documentation:

> head(rasterToPoints(era))
        x  y       tp
[1,] 0.00 90 9.336591
[2,] 0.25 90 9.336591
[3,] 0.50 90 9.336591
[4,] 0.75 90 9.336591
[5,] 1.00 90 9.336591
[6,] 1.25 90 9.336591
2
  • Thanks Spacedman! So if I understand this correctly, then to get a raster that resembles raster(res = 0.25) , then I'm left only with resampling? Apr 4, 2022 at 21:35
  • If you want a raster with extent (0,360) and (-90,90) with 0.25 cells, then (because the definition of "extent" includes the cell sizes) you will have to interpolate 721 to 720 rows and interpolate across half a column width as well to shift the grid over. The docs seem to read that you should consider the data to have "point support" at the given coordinates, rather than as values over any grid cells.
    – Spacedman
    Apr 4, 2022 at 22:06
0
cdo -f nc4 -z zip_1 sellonlatbox,-180,180,-90,90 infile.nc outfile.nc

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.