Create a raster brick from a netcdf using all levels

Thiago,?


On Sat, Sep 29, 2012 at 9:15 AM, Thiago Veloso <thi_veloso at 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 at gmail.com>
>To: Thiago Veloso <thi_veloso at yahoo.com.br> 
>Cc: R-SIG list <r-sig-geo at 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 at 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 at r-project.org
>>https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
>??