Issue working with thredds server and creating raster layers

[RMORSEcode]

Hi Mike, I am trying to use angstroms on a THREDDS server (not on a specific .nc file) and functions where a raster is returned do not seem to have the correct dimensions on my system.

For example, reading the depth layer 'h' from a roms file should give a raster with dimensions similar to [x, y, 1] , in my case is [242, 106, 1], but instead I get a raster with dimensions of [4,2,1]. I can work around this by using the ncdf4 library and calling nc_open on the file location on a case by case basis, but I cannot seem to modify the base functions (...or more likely, I am doing it improperly). I can send you the full thredds link privately if you want to try it yourself.

using default angstroms code

roms_file='http://tds.marine.rutgers.edu/.../2007-2016/avg' h=raster(roms_file, varname='h') dim(h) [1] 4 2 1 h class : RasterLayer band : 1 (of 2 bands) dimensions : 4, 2, 8 (nrow, ncol, ncell) resolution : 1, 1 (x, y) extent : 0, 2, 0, 4 (xmin, xmax, ymin, ymax) coord. ref. : NA data source : in memory names : avg values : 0, 5.787037e-06 (min, max)

using modified work-around

roms_file='http://tds.marine.rutgers.edu/.../2007-2016/avg' nc=nc_open(roms_file) h=raster(ncvar_get(nc, varid='h')) dim(h) [1] 242 106 1

within angstroms workaround

roms_file='http://tds.marine.rutgers.edu/.../2007-2016/avg' h=raster(ncget(x, varname = 'h')) dim(h) [1] 242 106 1 h class : RasterLayer dimensions : 242, 106, 25652 (nrow, ncol, ncell) resolution : 0.009433962, 0.004132231 (x, y) extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax) coord. ref. : NA data source : in memory names : layer values : 5, 5015.582 (min, max)

However, plot(h) shows that this is rotated at least 90 degrees and maybe inverted as well (???).

Partial Solution: Replacing raster(x, varname = depth) with raster(ncget(x, varname = depth)) in functions might fix the problem. Also applies to functions calling brick or stack. Thanks for your work on this! Cheers, Ryan

[mdsumner]

Please do send the thredds link. I'm surprised raster gets the wrong dimensions here, but I'll explore it.

The transpose/flip aspect is a perennial problem, NetCDF tends to be 0,0 based with Y pointing "up" - unlike most geo-spatial grid applications (it's because it's actually n-dimensional so why would Y be a special case ...) - but it's not always that simple. I think angstroms::romsdata uses the wrong default, but I'll see if I can figure out why you're getting a 4x2 grid from 'h'.

I tried what seems a similar source, and it's a problem with raster getting the wrong pair of dimensions, as if it is ignoring varname (it doesn't return the list of available varnames if you don't specifiy it either, which is unexpected).

ur <- "http://tds.marine.rutgers.edu/thredds/dodsC/roms/gom/g9/daily_avg"
library(raster)
brick(ur, varname = "h")
# class       : RasterBrick 
# dimensions  : 4, 2, 8, 2  (nrow, ncol, ncell, nlayers)
# resolution  : 1, 1  (x, y)
# extent      : 0, 2, 0, 4  (xmin, xmax, ymin, ymax)
# coord. ref. : NA 
# data source : http://tds.marine.rutgers.edu/thredds/dodsC/roms/gom/g9/daily_avg 
# names       : daily_avg.1, daily_avg.2 

library(ncdf4)
nc <- nc_open(ur)
h <- ncvar_get(nc, "h")
dim(h)
[1] 160 100

So, I think this is a raster/dods-variant problem - you are right about the workarounds but the fix (for now) belongs in the raster package. I'll try to investigate more. Thanks!

[mdsumner]

(But, why can't you share the full link here?)

[mdsumner]

Note also that raster(ncget(x, varname = 'h')) will give a raster in extent(0, 1, 0, 1) which is not very helpful, romsdata() goes an extra step and makes sure the index space is the extent.

[mdsumner]

Here's what I think is happening, I think the raster package logic for a NetCDF/Thredds source is being ignored - and the logic is passed to rgdal::readGDAL, and that's what's being returned. Raster triggers on ".nc" and not much else I think, so this is a rare case where your rgdal (and mine) is built for Thredds-NetCDF and the raster support is bypassed. (It's otherwise very hard to over-ride the use of ncdf4 and get rgdal instead without removing its installation - this will now occur on Windows too as the CRAN binaries can support NetCDF).

If you can get a source URL that has ".nc" in the path then it might work the usual way. I'll investigate more and let the raster author know (he's back, having indicated that ongoing support will occur).

[RMORSEcode]

Thanks Mike, I sent you an email with the link, but it sounds like this is a raster issue...

[mdsumner]

(I might move from using raster in favour of RNetCDF since that now has Thredds for Windows, so it's helpful)

[mdsumner]

Actually .... raster has an argument ncdf = TRUE and angstroms needed that, fixed in fb340cea8e219c6af1150dda6952e425d6146bd1

I also made "varname" required, as it was trying to be too clever. Your example should now work like this:

library(raster)
romsdata(ur, varname = "h")

Thanks!

[mdsumner]

I think we're good now.