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_ro...@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_ro...@pobox.com> <roche....@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]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo