[raster] problem with large netCDF file

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130210/07888461/attachment.pl>

Tom Roche Sun, 10 Feb 2013 14:40:05 -0500

$ ls -alh ./METCRO3D_080101
-r--r--r-- 1 me mygroup 5.3G Jan  8 11:32 ./METCRO3D_080101
# extension used below
$ ln -s ./METCRO3D_080101 ./METCRO3D_080101.nc

...

raster('./METCRO3D_080101.nc', varname='ZF')

...

Error in .local(.Object, ...) : 
  `[METCRO3D_080101]' not recognised as a supported file format.

Robert J. Hijmans Sun, 10 Feb 2013 13:53:56 -0800

ncdf files are indeed recognized by their extension (".NC",".CDF",

I was able rename the above file and subsequently load it with
`raster`, but only after copying the above file to my homespace :-(
which is generally unfeasible due to size.

I might add further checks that do not depend on the extension.

Or perhaps merely add a flag allowing the user to {assert, vouch for}
the file's netCDF-ness? Why I ask:

Files such as the above are generated by a production model of long
standing (~20 years), which moreover generates many other netCDF files
with the same naming convention ("<8 char description/>_<timestamp/>"
with no extension). Hence a request by me for a more appropriate
naming convention seems unlikely to be adopted, and would not help
work with many years of legacy output. Moreover these outputs are
usually quite large (i.e., tens of GB or more), and usually generated
to shared directories with restricted permissions. Hence it would be
quite desirable to enable current `raster` code to either

* respect the extension of a symlink

* take a runtime flag, e.g. 'format=CDF'

If there is somewhere to make a feature request, please let me know.

TIA, Tom Roche <Tom_Roche at pobox.com>

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130211/eccbfc17/attachment.pl>

summary: what are the semantics of 'band', 'level', and 'lvar' in
`raster:::.rasterObjectFromCDF` (et al)?

details:

https://stat.ethz.ch/pipermail/r-sig-geo/2013-February/017441.html

system('ls -alh ./METCRO3D_080101')

-r--r--r-- 1 me mygroup 5.3G Jan  8 11:32 ./METCRO3D_080101

system('ls -alh ./METCRO3D_080101.nc')

lrwxrwxrwx 1 me mygroup 17 Feb 10 13:46 ./METCRO3D_080101.nc -> ./METCRO3D_080101

raster('./METCRO3D_080101', varname='ZF')

...

Error in .local(.Object, ...) : 
  `[...]/METCRO3D_080101' not recognised as a supported file format.

https://stat.ethz.ch/pipermail/r-sig-geo/2013-February/017454.html

For now, the workaround you can use is this function

raster:::.rasterObjectFromCDF(filename, varname='', band=NA,
type='RasterLayer', lvar=3, level=0, warn=TRUE, ...)

I'm apparently not understanding the semantics of 'band', and possibly
'level', and 'lvar' (all of which I'm guessing `raster` is designed to
protect me from :-) As noted above, my datavar='ZF' of interest is in
netCDF='./METCRO3D_080101', which is

system(sprintf('ncdump -h %s', template.raw.fp))

netcdf METCRO3D_080101 {
dimensions:
        TSTEP = UNLIMITED ; // (25 currently)

TSTEP is the time dimension

        DATE-TIME = 2 ;
        LAY = 24 ;

LAY is the vertical dimension

        VAR = 17 ;
        ROW = 299 ;
        COL = 459 ;

ROW and COL are the horizontal dimensions

variables:

... followed by 10 datavars (none of which are coordinate variables),
then ...

        float ZF(TSTEP, LAY, ROW, COL) ;
                ZF:long_name = "ZF              " ;
                ZF:units = "M               " ;
                ZF:var_desc = "full-layer height above ground" ;

... followed by 7 more datavars (none of which are coordinate
variables), then many global attributes (of which this convention,
IOAPI, is inordinately fond :-) So I tried

raster:::.rasterObjectFromCDF(

+   './METCRO3D_080101', varname='ZF',
+   band=11,  # datavar=ZF is 11th reported by `ncdump -h`
+   level=1,  # "4th dimension variable" in (TSTEP, LAY, ROW, COL)
+   lvar=2,   # "3rd dimension variable" in (TSTEP, LAY, ROW, COL)
+   type='RasterLayer', warn=TRUE)
class       : RasterLayer 
band        : 11 
dimensions  : 299, 459, 137241  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0.5, 459.5, 0.5, 299.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : ./METCRO3D_080101 
names       : ZF 
z-value     : 11 
zvar        : ZF 
level       : 1 

But if I do

system('mv ./METCRO3D_080101 ./METCRO3D_080101.nc')
raster('./METCRO3D_080101.nc', varname='ZF')

I get

class       : RasterLayer 
band        : 1 
dimensions  : 299, 459, 137241  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0.5, 459.5, 0.5, 299.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : [...]/METCRO3D_080101.nc 
names       : ZF 
z-value     : 1 
zvar        : ZF 
level       : 1 

which is almost the same, but not quite. It seems my 'band' value is
wrong, so I'd like to know how to determine the appropriate value; I'd
also like to know if my comments above correctly reflect the semantics
of 'level' and 'lvar'. Unfortunately I'm not quite getting this from
raster.pdf, and the source for raster:::.rasterObjectFromCDF seems quite
inscrutable (YMMV :-)

TIA, Tom Roche <Tom_Roche at pobox.com>