[R-sig-Geo] extract netCDF layer then export

[Tim Howard]

Robert,
Yes, raster thinks my time steps are 'levels' rather than 'bands' (600 time steps in these examples):
> s <- stack("Extraction_Tavg.nc", bands=pos)
Error in .local(x, ...) : no valid bands supplied
In addition: Warning message:
In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
  "level" set to 1 (there are 600 levels)

> s <- stack("Extraction_Tavg.nc", levels=pos)
Error in .local(x, ...) : 
  Arguments should be Raster* objects or filenames
In addition: Warning messages:
1: In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
  "level" set to 1 (there are 600 levels)
2: In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
  "level" set to 1 (there are 600 levels)
 

It looks like a call to raster works, for a single time step:
 
> r <- raster("Extraction_Tavg.nc", level=535)
> r
class       : RasterLayer 
dimensions  : 49, 97, 4753  (nrow, ncol, ncell)
resolution  : 0.125, 0.125  (x, y)
extent      : -82.125, -70, 40, 46.125  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
values      : C:\mypath\Extraction_Tavg.nc 
zvar        : Tavg 
level       : 535 
 
 
The only thing that's confusing to me, is that ncdf *does* seem to know that the levels are time steps:

> ncTemp <- open.ncdf("Extraction_Tavg.nc")
> ncTemp
[1] "file Extraction_Tavg.nc has 4 dimensions:"
[1] "time   Size: 600"
[1] "latitude   Size: 49"
[1] "longitude   Size: 97"
[1] "projection   Size: 1"
[1] "------------------------"
[1] "file Extraction_Tavg.nc has 1 variables:"
[1] "float Tavg[longitude,latitude,time,projection]  Longname:Tavg Missval:1.00000002004088e+20"
> 
 
At any rate, I can work with this (via multiple calls to raster in a loop). Your code and examples are VERY helpful and I'm very happy with what I've got now. 
Thank you so much for the quick help. 
Best, 
Tim
 

>>> "Robert J. Hijmans" <r.hijmans at gmail.com> 4/9/2012 1:59 PM >>>
Tim, 

Currently you can have only one level (the fourth dimension) in a single RasterBrick, so for each level you need to make a new object (whether RasterLayer, Brick, or Stack). But you can make a RasterStack for only a few time steps (the third dimension, "bands"). 

s <- stack("Extraction_Tavg.nc", bands=c(1, 13, 15)) #grab only time steps 1, 13, 15

You can change what you consider a 'level' vs 'time' and choose the level with the 'level' and 'lvar' arguments (which do not appear to be documented). The defaults are lvar=3, level=0. 

But if you are comfortable with directly using the ncdf functions, that may of course work as well or better. 

Robert

On Mon, Apr 9, 2012 at 10:48 AM, Tim Howard <tghoward at gw.dec.state.ny.us> wrote:


I'm actually doing more than just grabbing a single layer/band in R. I've got a script set up to extract specific layers from the time dimension, get the average, then write it out. Ideally there would be a structure similar to this (is this possible?): 
b <- brick("Extraction_Tavg.nc", band = c(1, 13, 15)) #grab only time steps 1, 13, 15
This was why I resorted to using "get.var.ncdf" 
Perhaps there's a way I can initialize a rasterbrick or stack using your suggestions and then return to get.var.ncdf to get the layers I want. 
Tim

>>> Tim Howard 4/9/2012 1:29 PM >>>
Robert,
Thank you. That does pull the extent and resolution information over, but perhaps the netCDF files I am using are structured funny as it throws an error for the levels:
> b <- brick("Extraction_Tavg.nc")
Warning message:
In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
"level" set to 1 (there are 3 levels)
> projection(b) <- "+proj=longlat +datum=WGS84"
> r3 <- b[[3]]
Error in .local(x, ...) : not a valid subset
> r3 <- b[[2]]
Error in .local(x, ...) : not a valid subset
> b
class : RasterBrick 
dimensions : 2, 2, 4, 1 (nrow, ncol, ncell, nlayers)
resolution : 0.125, 0.125 (x, y)
extent : -74.25, -74, 42, 42.25 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
values : <<mypath>>Extraction_Tavg.nc 
zvar : Tavg 
level : 1 
Tim

>>> "Robert J. Hijmans" <r.hijmans at gmail.com> 4/9/2012 1:19 PM >>>
Tim,

You can do:

library(raster)
b <- brick("Extraction_Tavg.nc")
projection(b) <- '+proj=longlat +datum=WGS84'

# multilayer object
b <- writeRaster(b, 'boo.tif') 

# for single layer (time 3)
r3 <- b[[3]]
r3 <- writeRaster(r3, 'boo3.tif') 


# or
r2 <- raster("Extraction_Tavg.nc", band=2)
projection(r2) <- '+proj=longlat +datum=WGS84'
r2 <- writeRaster(r2, 'boo3.tif') 

Or something similar with writeRaster following you approach. 
Oh, and I thought that ArcGIS can now also open ncdf files

Robert



On Mon, Apr 9, 2012 at 10:02 AM, Tim Howard <tghoward at gw.dec.state.ny.us> wrote:


I've got most of the steps for manipulating netCDF files using packages ncdf and raster but I'm currently stuck getting the correct projection information (lat/long) over to the raster and then exporting it.

Here's a VERY simply netCDF raster (2 cells lat X 2 cells long X 3 time periods) and how I might grab one of those time periods (I tried using dput to actually provide the data, but couldn't get dget to re-load it, I'd be happy to supply data if someone can suggest how).

> library(ncdf)
> library(raster)
> ncTemp2 <- open.ncdf("Extraction_Tavg.nc")
> ncTemp2
[1] "file Extraction_Tavg.nc has 4 dimensions:"
[1] "time Size: 3"
[1] "latitude Size: 2"
[1] "longitude Size: 2"
[1] "projection Size: 1"
[1] "------------------------"
[1] "file Extraction_Tavg.nc has 1 variables:"
[1] "float Tavg[longitude,latitude,time,projection] Longname:Tavg Missval:1.00000002004088e+20"
> str(ncTemp2)

<<snip to save space>>>>

My Question: How do I now export the single raster, m, so that I can bring along all the relevant spatial information contained in the netCDF (which would be the lat long coordinates?). I'd like to export the raster into any format accessible by esri products (boo!).

Thanks in advance,
Tim




[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org 
https://stat.ethz.ch/mailman/listinfo/r-sig-geo