4

I've recently started working with spatial data in R. I am currently working on a project that uses data in both raster and netcdf formats, and I'd be happy to convert the latter to the former, however I have certain doubts about the procedure.

First I tried the very simple approach of just opening the netcdf using stack(), so basically just:

library(raster)
filename <- "C:/.../X.nc"
pol_st <- stack(filename)

This seems to work fine, I can plot the rasters etc. However, the resulting rasters do not have the proper extent:

> pol_st
class       : RasterStack 
dimensions  : 576, 992, 571392, 4475  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 0.5, 992.5, 0.5, 576.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
names       : X0, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14, ... 

Now, I know the actual extent should be somewhere around 73, 135, 18, 54. I suppose it wouldn't be a problem to change the extent of the raster (right?), but:

  • Can I be sure that the projection is correct? (Do I understand correctly that there is no projection in case of netcdf, just an array of longitudes, latitudes and values?)
  • I still need to get the exact extent from the netcdf data.

I was looking for some answers and found threads such as this and that, which made me think that actually the solution might be more complicated.

I've tried to recreate the methods shown in those threads, but have difficulties in creating rasters of lat and lon, that is, for:

library(raster)
filename <- "C:/.../X.nc"
vname <- "lat"
pol_st <- stack(filename, varname = vname)

I get:

Error in .varName(nc, varname) : 
  varname: lat does not exist in the file. Select one from:
NO2

I've tried different names for the variable too, but the "Select one from: NO2" makes it look like the data is missing lon and lat?

Unfortunately I am unable to share the data, but if I run:

library(ncdf4)
polnc <- nc_open(filename)
polnc

I get:

File C:/.../X.nc (NC_FORMAT_NETCDF4_CLASSIC):
 3 variables (excluding dimension variables):
    float latitude[lat]   
    float longitude[lon]   
    float NO2[lon,lat,time]   

 3 dimensions:
    lat  Size:576
    lon  Size:992
    time  Size:4475
2

Your Netcdf has three variables and three dimensions. The variables are latitude, longitude, and NO2:

float latitude[lat]   
float longitude[lon]   
float NO2[lon,lat,time]   

each of those is indexed by dimensions. latitude is indexed by lat, longitude by lon, and NO2 by three dimensions of lon,lat, and time. The dimensions section tells you how big these dimensions are:

3 dimensions:
    lat  Size:576
    lon  Size:992
    time  Size:4475

To get the value of the latitude and longitude coordinates you need to get the corresponding variable from the NetCDF, something like:

ncvar_get(polnc, "longitude")

might get you a vector of the longitude coordinates, which you can then assign to the raster - and similarly for latitude.

You might also try getting the dimensions from the $dim element:

polnc$dim["lon"]

raster will probably only read two-dimensional variables from NetCDFs, which is probably why it only suggests NO2 as a readable variable. In any case lat isn't a variable (its a dimension) - latitude is a variable.

  • Thank you for the answer! ncvar_get worked like a charm. Interestingly, $dim$vals is just a vector 1:992. So I suppose I do not need to worry about the projection, right? – yassem Feb 18 '18 at 2:43
  • I suspect the dimensions are all 1:N and the values are encoded in the variable with that dimension, for some reason. If the coords look like lat-long then they are probably EPSG:4326 geographic coordinates, but they could be other lat-long schemes based on different earth spheroid models or datums. Check with the source if you need high precision! – Spacedman Feb 18 '18 at 12:29
3

To add a few notes:

The problem is with the file content, the coordinates are not provided in the standard way. longitude and latitude should normally be dimensions, not variables. lon and lat are then redundant. I would ask the source of the files to fix them.

The coordinate reference system (crs) is probably not provided in the file, so you get a default guess; perhaps questionable. Many NetCDF files do provide the coordinate reference system; in most cases it is longitude/latitude.

If you happend to know the actual extent and crs you can set these.

extent(pol_st) <- c( , , ,  )
crs(pol_st) <- " "

Also, instead of using stack, it is much more efficient to do

pol_st <- brick(filename)
  • Just to clarify - when I use ncvar_get I can see the actual coords, it's just $dim$vals that returns 1:992. So the NetCDF does provide actual coordinates. In that case could I use a method like one of the previous threads? Is there even any need for that? As for the stack, I just assumed that since it requires the rasters to be of the same extent it is somehow more efficient, but good to know, thanks! – yassem Feb 20 '18 at 16:25
  • You are right, and I edited my answer a bit to reflect that. You should be able to use these files by creating Raster objects "by hand" – Robert Hijmans Feb 20 '18 at 18:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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