Thiago,

On Sat, Sep 29, 2012 at 9:15 AM, Thiago Veloso <thi_veloso@yahoo.com.br>wrote:

> Robert,
>
> Thank you very much for the tip. However, I still cannot load all the "z"
> layers in only one brick. Please see the comments below:
>
>
That is because this file has levels. It has lat/long/depth/time. The
default is for brick() to return all time steps (z) for a  level (depth) of
choice. This file is odd because it has a single, meaningless, time step
that should have been omitted.

The trick is to tell brick that depth is the 4th dimension (z) and that
time is the 3rd dimension (level), using the lvar (level variable)
argument. You also need to set the NA value.

library(raster)
f <- "soita.clay.nc"
b <- brick(f, lvar=4)
NAvalue(b) <- 9e+20
plot(b)

# less efficiently, to show how it works, you could also have done

r1 <- raster(f, level=1)
r2 <- raster(f, level=2)
r3 <- raster(f, level=3)
r4 <- raster(f, level=4)
r5 <- raster(f, level=5)
r6 <- raster(f, level=6)
s <- stack(r1, r2, r3, r4, r5, r6)
NAvalue(s) <- 9e+20


Best, Robert


> library(ncdf)
> > library(raster)
>
> # Reading file with ncdf package to examine its strucuture.
> > n <- open.ncdf ('~/Dropbox/web/soita.clay.nc')
> > n
> [1] "file ~/Dropbox/web/soita.clay.nc has 5 dimensions:"
> [1] "longitude   Size: 720"
> [1] "latitude   Size: 360"
> [1] "layer   Size: 6"
> [1] "time   Size: 1"
> [1] "lengthd   Size: 10"
> [1] "------------------------"
> [1] "file ~/Dropbox/web/soita.clay.nc has 1 variables:"
> [1] "double claypct[longitude,latitude,layer,time]  Longname:claypct
> Missval:1e+30"
>
> # The dimension I need to make a brick of is "layer", which contains data
> in six soil levels.
> # Therefore, I expected to create a brick with nlayers=6, but that doesn't
> happen. See below:
> > b <- brick('~/Dropbox/web/soita.clay.nc',level=3)
> Warning message:
> In .doTime(r, nc, zvar, dim3, ncdf4) : assuming a standard calender:0
>
> > b
> class       : RasterBrick
> dimensions  : 360, 720, 259200, 1  (nrow, ncol, ncell, nlayers)
> resolution  : 0.5, 0.5  (x, y)
> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84
> data source : ~/Dropbox/web/soita.clay.nc
> names       : X1958.01.01
> Date        : 1958-01-01
> varname     : claypct
> level       : 3
>
> It looks like raster assumes that the third level in the file contains
> only one layer. Is this right?
>
>

Yes, and that is a correct assumption. Each level as a single time layer




> In case there's a need to inspect the input file, it is 12MB:
> https://www.dropbox.com/s/hyalyip9ai8yidc/soita.clay.nc
>
> Best,
> Thiago.
>
>   ------------------------------
> *From:* Robert J. Hijmans <r.hijmans@gmail.com>
> *To:* Thiago Veloso <thi_veloso@yahoo.com.br>
> *Cc:* R-SIG list <r-sig-geo@r-project.org>
> *Sent:* Saturday, September 29, 2012 4:10 AM
> *Subject:* Re: [R-sig-Geo] Create a raster brick from a netcdf using all
> levels
>
>
> Thiago,
>
> You can try using the (perhaps undocumented) 'level' argument
>
> clay3 <- brick('/data/input-inland/soita.clay.nc', level=3)
>
> With a RasterBrick you can only one level at a time (i.e. x, y and z, but
> not a fourth dimension), but you can create several RasterBricks, one for
> each level.
>
> Robert
>
>
> On Thu, Sep 27, 2012 at 6:41 AM, Thiago Veloso <thi_veloso@yahoo.com.br>wrote:
>
>   Dear list,
>
>   I have a netcdf file that contains 6 levels. Below is the structure of
> the file:
>
> require(ncdf)
> require(raster)
> > nctemp <- open.ncdf('/data/input-inland/soita.clay.nc')
>
> > nctemp
> [1] "file /data/input-inland/soita.clay.nc has 5 dimensions:"
> [1] "longitude   Size: 720"
> [1] "latitude   Size: 360"
> [1] "layer   Size: 6"
> [1] "time   Size: 1"
> [1] "lengthd   Size: 10"
> [1] "------------------------"
> [1] "file /data/input-inland/soita.clay.nc has 1 variables:"
> [1] "double claypct[longitude,latitude,layer,time]  Longname:claypct
> Missval:1e+30"
> >
>
>   What I need to do is to load this netcdf to a brick raster object using
> all levels. When creating the brick file, specifying "nl=6" results in
> loading only the last level:
>
> > clay <- brick('/data/input-inland/soita.clay.nc', nl=6)
> Warning messages:
> 1: In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
>   "level" set to 1 (there are 6 levels)
> 2: In .doTime(r, nc, zvar, dim3, ncdf4) : assuming a standard calender:0
> > clay
> class       : RasterBrick
> dimensions  : 360, 720, 259200, 1  (nrow, ncol, ncell, nlayers)
> resolution  : 0.5, 0.5  (x, y)
> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84
> data source : /data/input-inland/soita.clay.nc
> names       : X1958.01.01
> Date        : 1958-01-01
> varname     : claypct
> level       : 1
>
>   I have also tested another approach, for example specifying "nl=1:6",
> however with no sucess at all. Does anyone know a workaround to load all
> leves of a netcdf file in a single raster brick object?
>
>   Many thanks in advance,
>   Thiago.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
>
>

	[[alternative HTML version deleted]]

