2

Really simple problem with the raster package, also using ncdf4 to load in an ECMWF Era-Interim Netcdf file.

Simply doing this:

a <- nc_open("SSTs.nc")
B <- brick(a, varname="sst")

Returns this:

    Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘brick’ for signature ‘"ncdf4"’

The file is just SST data over the whole globe, for 1 month (Jan2016). When I convert it into an array (i.e. extract dimensions/variable, and convert time to UTC, shove it into an array) I don't get the same error, but the raster package says it supports .nc files straight in (so long as they're cf-1 compatible, which Era-Interim .nc's are)

Any help much appreciated, have tried this with many Netcdf files (non-Era Interim too).

Ndharwood
  • 123
  • 3
  • 11
  • Could you give the sample dataset, also you can try using `stack` or using previous `ncdf package` – Eko Susilo Mar 24 '17 at 09:27
  • the same question [Read Netcdf sub categories and convert to grid] http://stackoverflow.com/questions/33784940/read-netcdf-sub-categories-and-convert-to-grid – Eko Susilo Mar 24 '17 at 13:35
  • Thanks for the reply. Using stack returned this: Error in data.frame(values = unlist(unname(x)), ind = factor(rep.int(names(x), : arguments imply differing number of rows: 1654, 16 Sample dataset is this: https://drive.google.com/open?id=0Bz0W7Ut_SNfjcjg1ODVrc2FhN2s – Ndharwood Mar 24 '17 at 13:40
  • @Ndharwood There shouldn't be any problem if you use `s <- stack("SSTs.nc")` – Geo-sp Apr 04 '17 at 19:36

1 Answers1

1

thank Renaud Lancelot, who give clearly source code. I have modified his code to fit with your data

 # load package
 library(sp)
 library(raster)
 library(ncdf4)

 # read ncdf file
 nc<-nc_open('D:/SSTs.nc')

 # extract variable name, size and dimension
 v <- nc$var[[1]]
 size <- v$varsize
 dims <- v$ndims
 nt <- size[dims]              # length of time dimension
 lat <- nc$dim$latitude$vals   # latitude position
 lon <- nc$dim$longitude$vals  # longitude position

 # read sst variable
 r<-list()
 for (i in 1:nt) {
   start <- rep(1,dims)     # begin with start=(1,1,...,1)
   start[dims] <- i             # change to start=(1,1,...,i) to read    timestep i
   count <- size                # begin with count=(nx,ny,...,nt), reads entire var
   count[dims] <- 1             # change to count=(nx,ny,...,1) to read 1 tstep

   dt<-ncvar_get(nc, varid = 'sst', start = start, count = count)

   # convert to raster
   r[i]<-raster(dt)
 }

 # create layer stack with time dimension
 r<-stack(r)

 # transpose the raster to have correct orientation
 rt<-t(r)
 extent(rt)<-extent(c(range(lon), range(lat)))

 # plot the result
 spplot(rt)
Eko Susilo
  • 250
  • 2
  • 15
  • Thanks Eko - this is brilliant. Especially thankful for the orientation correction. Could I possibly ask you to clarify what the steps of the For loop do? This would be helpful so I can start messing around with it as the comments suggest. – Ndharwood Apr 11 '17 at 14:39
  • Regarding your ncdf have 3 dimensions (longitude, latitude, and time), For loop will read time series data (depend on time dimension) and assign the value into r[i], the result will be re-stack to get all data in a raster fil – Eko Susilo Apr 16 '17 at 05:12
  • Hi Eko, thanks again - can you tell me what lines 2 & 4 inside the For loop do? The 1st line replicates timestep 1 BY the number of dimensions in the variable, but why does the 2nd line subset dims - i.e. start[dims]? Again, the 3rd is to tell the ncvar_get function to count for the size of the variable, but what does the 4th line ( count[dims] ) do? Many thanks – Ndharwood Apr 20 '17 at 16:22
  • We need to define the `start` and `count` argument first: `start` A vector of indices indicating where to start reading the passed values (beginning at 1), the length must equal the number of dimensions (order is X-Y-Z-T (i.e., time is last) `count` A vector of integers indicating the count of values to read along each dimension (order is X-Y-Z-T), the length must equal the number of dimensions the variable – Eko Susilo Apr 21 '17 at 03:14
  • Eko - In this code there's no extraction of the time variable, which means the layers are named "layer 1" etc instead of the date. Is there anyway of setting the date as the layer name? I've tried setNames() using a 'date' object which I extracted in POSIXlt format but it didn't work. – Ndharwood Apr 27 '17 at 13:46
  • here additional command to names your raster (1) retreived time information `dt<-as.Date(nc$dim$time$vals/24, origin='1900-01-01 00:00:00')` and (2) rename layer `names(rt)<-dt` – Eko Susilo Apr 28 '17 at 03:28
  • Thanks Eko, you've been immensely helpful! The layer names have 'X' prepended onto them automatically (e.g. 'X2016-01-01') but I'm sure I'll figure out how to remove this. These replies did everything I need! – Ndharwood May 16 '17 at 15:16