Re: [R-sig-Geo] [raster] conserving dimensions when regridding a 4D Brick or Stack

2013-04-09 Thread Robert J. Hijmans 
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=637 +b=637 ;
  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
 #

Re: [R-sig-Geo] [raster] conserving dimensions when regridding a 4D Brick or Stack

2013-04-09 Thread Tom Roche 

https://stat.ethz.ch/pipermail/r-sig-geo/2013-April/017947.html
 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

https://stat.ethz.ch/pipermail/r-sig-geo/2013-April/017973.html
 A RasterBrick is 3D [only].

also RasterStack, no? It Would Be Nice if the docs made this explicit--
or am I missing something?

 all you can do with the raster package is process the data level by
 level

I guessed:

# regridding, including the 4D 3hourlies
https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/src/f5e1fd0d607da8bbfe67bb18021dcf28e1f2/vis_regrid_vis.r?at=master
# reassembling the 4D 3hourlies
https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/src/f5e1fd0d607da8bbfe67bb18021dcf28e1f2/fix_3hourlies.ncl?at=master

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

[R-sig-Geo] [raster] conserving dimensions when regridding a 4D Brick or Stack

2013-04-04 Thread Tom Roche 

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=637 +b=637 ;
 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, 186  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=637 
+b=637 
# 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