Tom,
A RasterBrick is 3D. There is a 4D object RasterQuad designed to work with
the type of data you have, but I have not implemented any methods for it,
so it is currently of no use. I think all you can do with the raster
package is process the data level by level and then combine everything into
a single NCDF file by directly using functions from the ncdf(4) package.
Robert


On Thu, Apr 4, 2013 at 9:21 PM, Tom Roche <Tom_Roche@pobox.com> wrote:

>
> summary: When regridding a 4D Brick (e.g., the netCDF data variable
> Fraction_of_Emissions(TSTEP, LAY, ROW, COL) in
>
> > netcdf GFED-3.1_2008_N2O_3hourly_fractions {
> > dimensions:
> >   TSTEP = 12 ;
> >   LAY = 8 ;
> >   ROW = 360 ;
> >   COL = 720 ;
> > variables:
> >   float Fraction_of_Emissions(TSTEP, LAY, ROW, COL) ;
> >     Fraction_of_Emissions:units = "unitless" ;
> >     Fraction_of_Emissions:long_name = "3hourly fraction of monthly N2O
> emission" ;
> >     Fraction_of_Emissions:_FillValue = 9.96921e+36f ;
> >   int TSTEP(TSTEP) ;
> >     TSTEP:long_name = "index of month in 2008" ;
> >     TSTEP:units = "month" ;
> >   int LAY(LAY) ;
> >     LAY:units = "4-digit hour" ;
> >     LAY:long_name = "time of day starting 3-hour-long period" ;
> >   float ROW(ROW) ;
> >     ROW:units = "degrees_north" ;
> >     ROW:long_name = "latitude" ;
> >   float COL(COL) ;
> >     COL:units = "degrees_east" ;
> >     COL:long_name = "longitude" ;
>
> ) I want to keep all 4 dimensions in the output. But I am only able to
> write 3D output, since I lose dimension=LAY:
>
> > netcdf GFED-3.1_2008_N2O_3hourly_fractions_regrid {
> > dimensions:
> >   COL = 459 ;
> >   ROW = 299 ;
> >   TSTEP = UNLIMITED ; // (12 currently)
> > variables:
> >   double COL(COL) ;
> >     COL:units = "meter" ;
> >     COL:long_name = "COL" ;
> >   double ROW(ROW) ;
> >     ROW:units = "meter" ;
> >     ROW:long_name = "ROW" ;
> >   int TSTEP(TSTEP) ;
> >     TSTEP:units = "month" ;
> >     TSTEP:long_name = "TSTEP" ;
> >   float Fraction_of_Emissions(TSTEP, ROW, COL) ;
> >     Fraction_of_Emissions:units = "unitless" ;
> >     Fraction_of_Emissions:_FillValue = -3.4e+38 ;
> >     Fraction_of_Emissions:missing_value = -3.4e+38 ;
> >     Fraction_of_Emissions:long_name = "3hourly fraction of monthly N2O
> emission" ;
> >     Fraction_of_Emissions:projection = "+proj=lcc +lat_1=33 +lat_2=45
> +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000" ;
> >     Fraction_of_Emissions:projection_format = "PROJ.4" ;
>
> ) Can one
>
> + create a 4D RasterBrick or RasterStack
> + instruct projectRaster
>
> so as to get 4D regridded output? Or must one
>
> - create 3D RasterLayer's (i.e., one per LAY) from the 4D Brick
> - regrid each 3D Layer
> - reassemble the 4D Brick
>
> ? or Something Completely Different?
>
> details:
>
> As part of
>
> https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na
>
> I need to regrid all 4 dimensions of the data variable
> Fraction_of_Emissions(TSTEP, LAY, ROW, COL) in a netCDF file:
>
> > netcdf GFED-3.1_2008_N2O_3hourly_fractions {
> > dimensions:
> >   TSTEP = 12 ;
> >   LAY = 8 ;
> >   ROW = 360 ;
> >   COL = 720 ;
> > variables:
> >   float Fraction_of_Emissions(TSTEP, LAY, ROW, COL) ;
> >     Fraction_of_Emissions:units = "unitless" ;
> >     Fraction_of_Emissions:long_name = "3hourly fraction of monthly N2O
> emission" ;
> >     Fraction_of_Emissions:_FillValue = 9.96921e+36f ;
> >   int TSTEP(TSTEP) ;
> >     TSTEP:long_name = "index of month in 2008" ;
> >     TSTEP:units = "month" ;
> >   int LAY(LAY) ;
> >     LAY:units = "4-digit hour" ;
> >     LAY:long_name = "time of day starting 3-hour-long period" ;
> >   float ROW(ROW) ;
> >     ROW:units = "degrees_north" ;
> >     ROW:long_name = "latitude" ;
> >   float COL(COL) ;
> >     COL:units = "degrees_east" ;
> >     COL:long_name = "longitude" ;
>
> In
>
> https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/src/6a214b2076dee1ce4cae5ee44779642b33481d81/vis_regrid_vis.r?at=master
> I load the data like (substituting and reformatting for mail)
>
> hour3ly.in.raster <- raster::brick(
>   'GFED-3.1_2008_N2O_3hourly_fractions.nc',
>   varname='Fraction_of_Emissions')
> hour3ly.in.raster@crs <- global.crs
> hour3ly.in.raster
> # class       : RasterBrick
> # dimensions  : 360, 720, 259200, 12  (nrow, ncol, ncell, nlayers)
> # resolution  : 0.5, 0.5  (x, y)
> # extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> # coord. ref. : +proj=longlat +ellps=WGS84
> # data source : .../GFED-3.1_2008_N2O_3hourly_fractions.nc
> # names       : X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12
> # month       : 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12
> # varname     : Fraction_of_Emissions
> # level       : 1
>
> and then regrid with
>
> hour3ly.out.raster <-
>   raster::projectRaster(
>     from=hour3ly.in.raster, to=template.raster, crs=out.crs,
>     method='bilinear', overwrite=TRUE, format='CDF',
>     # args from writeRaster
>     varname=hour3ly_out_datavar_name,
>     varunit=hour3ly_out_datavar_units,
>     longname=hour3ly_out_datavar_long_name,
>     xname=hour3ly_out_x_var_name,
>     yname=hour3ly_out_y_var_name,
>     zname=hour3ly_out_z_var_name,
>     zunit=hour3ly_out_z_var_units,
>     filename=hour3ly_out_fp)
> hour3ly.out.raster
> # class       : RasterBrick
> # dimensions  : 299, 459, 137241, 12  (nrow, ncol, ncell, nlayers)
> # resolution  : 12000, 12000  (x, y)
> # extent      : -2556000, 2952000, -1728000, 1860000  (xmin, xmax, ymin,
> ymax)
> # coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97
> +a=6370000 +b=6370000
> # data source : /tmp/gfed-3.1_global_to_aqmeii-na/
> GFED-3.1_2008_N2O_3hourly_fractions_regrid.nc
> # names       : X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12
> # month       : 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12
> # varname     : Fraction_of_Emissions
>
> The output appears to horizontally-regrid correctly, but loses needed
> dimension=LAY:
>
> # netcdf GFED-3.1_2008_N2O_3hourly_fractions_regrid {
> # dimensions:
> #   COL = 459 ;
> #   ROW = 299 ;
> #   TSTEP = UNLIMITED ; // (12 currently)
> # variables:
> #   double COL(COL) ;
> #     COL:units = "meter" ;
> #     COL:long_name = "COL" ;
> #   double ROW(ROW) ;
> #     ROW:units = "meter" ;
> #     ROW:long_name = "ROW" ;
> #   int TSTEP(TSTEP) ;
> #     TSTEP:units = "month" ;
> #     TSTEP:long_name = "TSTEP" ;
> #   float Fraction_of_Emissions(TSTEP, ROW, COL) ;
>
> note the input (above) is
>     float Fraction_of_Emissions(TSTEP, LAY, ROW, COL) ;
>
> #     Fraction_of_Emissions:units = "unitless" ;
> #     Fraction_of_Emissions:_FillValue = -3.4e+38 ;
> #     Fraction_of_Emissions:missing_value = -3.4e+38 ;
> #     Fraction_of_Emissions:long_name = "3hourly fraction of monthly N2O
> emission" ;
> #     Fraction_of_Emissions:projection = "+proj=lcc +lat_1=33 +lat_2=45
> +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000" ;
> #     Fraction_of_Emissions:projection_format = "PROJ.4" ;
>
> Based on the latest version of
>
> http://cran.r-project.org/web/packages/raster/raster.pdf
>
> and on
>
> https://stat.ethz.ch/pipermail/r-sig-geo/2012-September/016175.html
>
> I thought use of the parameters `level` or `lvar` might allow me to
> usefully tell `raster` of my desire to get dimension=LAY in the input
> to, and the output from, the regridding process. However my attempts
> to date have failed :-( Given that in
>
> Fraction_of_Emissions(TSTEP, LAY, ROW, COL)
>
> the horizontal dimensions are in positions=[2,3] in the dimension
> array, and guessing at parameter values from the above documents
> (which IMHO could be clearer regarding the semantics of `level` and
> `lvar`, but perhaps I am missing something), I tried
>
> hour3ly.in.raster <- raster::brick(
>   hour3ly_in_fp, varname=hour3ly_in_datavar_name, level=0, lvar=1)
>
> drawing the warning
>
> > Warning message:
> > In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
> >   "level" set to 1 (there are 720 levels)
>
> so I instead tried
>
> hour3ly.in.raster <- raster::brick(
>   hour3ly_in_fp, varname=hour3ly_in_datavar_name, level=1, lvar=0)
> hour3ly.in.raster@crs <- global.crs
> hour3ly.out.raster <-
>   raster::projectRaster(
>     from=hour3ly.in.raster, to=template.raster, crs=out.crs,
>     method='bilinear', overwrite=TRUE, format='CDF',
>     # args from writeRaster
>     varname=hour3ly_out_datavar_name,
>     varunit=hour3ly_out_datavar_units,
>     longname=hour3ly_out_datavar_long_name,
>     xname=hour3ly_out_x_var_name,
>     yname=hour3ly_out_y_var_name,
>     zname=hour3ly_out_z_var_name,
>     zunit=hour3ly_out_z_var_units,
>     filename=hour3ly_out_fp)
>
> This proceeded without error or warning, but the output from
>
> system(sprintf('ncdump -h %s', hour3ly_out_fp))
>
> remains unchanged; i.e., still no dimension=LAY in the output file.
> Similarly I tried
>
> hour3ly.in.raster <- raster::brick(
>   hour3ly_in_fp, varname=hour3ly_in_datavar_name, level=2, lvar=3)
>
> hour3ly.in.raster <- raster::brick(
>   hour3ly_in_fp, varname=hour3ly_in_datavar_name, level=3, lvar=2)
>
> and got the same results. What am I doing wrong?
>
> TIA, Tom Roche <Tom_Roche@pobox.com> <Roche.Tom@epa.gov>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]

