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

Tom Roche Tom_Roche at pobox.com
Fri Apr 5 06:21:32 CEST 2013


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 at 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 at 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 at pobox.com> <Roche.Tom at epa.gov>



More information about the R-sig-Geo mailing list