[R-sig-Geo] reboxing with package=gstat

[Tom Roche]

Apologies if the following reveals a profound lack of understanding of
spatiality or geostatistics (esp kriging), which I'm learning slowly and
"on-the-fly":

summary: for "reboxing" from a regular global unprojected grid to a
regular local projected grid (both 3D), should one use

* gstat::krige? If so, what model?

* gstat::idw? If so, what formula?

or something completely different?

details: 

https://stat.ethz.ch/pipermail/r-sig-geo/2012-November/016573.html
>>> [I seek to] incorporate N2O emissions inventories (EIs), initial
>>> conditions (ICs), and boundary conditions (BCs) from relatively
>>> coarse-scaled global inventories (with horizontal grids on the order
>>> of degrees lon-lat) into a [finer-scaled] regional atmospheric model
>>> (12-km LCC horizontal). The [EIs] are [surface/2D] and I now know
>>> how to "regrid," or interpolate from one 2D grid to another.

>>> However the model and [its] IC/BCs are 3D[, and] have [both]
>>> different vertical layering (i.e., [number and height of layers])
>>> [and] different horizontal grids. I'd like to know: is there R code
>>> available that will "rebox," i.e., interpolate vertically as well as
>>> horizontally?

To further specify: my global IC/BC estimates mass concentration for
each box/voxel in a regular, unprojected grid 1.875° x 2.5° x 56 levels.
I want to "rebox" from global IC/BCs to model-ready IC/BCs, where the
latter must estimate mass concentration over a regular grid that
horizontally projects to LCC over North America at 12-km resolution, and
which has 34 vertical levels. (Both the vertical and horizontal extents
of the model-ready IC/BC are subsets of the global extents.) I can
compute the sizes of the boxes, so I can compute the mass in each input
box, and recompute the output concentrations, presuming I can learn how
to rebox the masses.

https://stat.ethz.ch/pipermail/r-sig-geo/2012-November/016574.html
(rearranged)
>> [A] typical issue is [whether] variability, or spatial correlation,
>> in x/y differs from that in z for the same distance unit[.]

I believe the vertical and horizontal correlations of N2O concentrations
differ strongly:

* vertical distribution of N2O is governed by atmospheric processes
  (advection and diffusion) very different from the mostly biological
  pedospheric processes that govern the horizontal distribution of its
  emissions (though the two sets of processes are coupled via
  precipitation and temperature)

* N2O in the troposphere (which unfortunately is all my group can now
  model) has a long life (scale ~= century) relative to the temporal
  scale of the processes that produce/emit it (which vary diurnally and
  seasonally, with soil moisture and temperature)

That being said, N2O datasets are quite appallingly sparse, given its
current significance as a depleter of stratospheric O3 and as a GHG.

>> Package gstat provides 3D interpolation:

>> library(gstat)
>> demo(gstat3D)

>> in this case, if you'd leave out the variogram model, inverse distance
>> is used for interpolation.

I note

R session
> demo(gstat3D)
...
> # $Id: gstat3D.R,v 1.5 2007-02-23 13:34:07 edzer Exp $
> # simple demo of 3D interpolation of 50 points with random normal
> # values, randomly located in the unit cube
> n <- 50
> data3D <- data.frame(x=runif(n), y=runif(n), z=runif(n), v=rnorm(n))
> coordinates(data3D)=~x+y+z
> range1D <- seq(from=0, to=1, length=20)
> grid3D <- expand.grid(x=range1D, y=range1D, z=range1D)
> gridded(grid3D)=~x+y+z
> res3D <- krige(formula=v ~ 1, data3D, grid3D, model=vgm(1, "Exp", .2))

So for my usecase,

1 'coordinates(data3D)' would be the coordinates of the centers of the
  boxes defined by the global/unprojected grid

2 'data3D' would be the pairs of coordinates and mass concentrations
  for the global/unprojected grid

3 'grid3D' would be the coordinates of the centers of the
  local/projected gridboxes

4 'res3D' would be the pairs of coordinates and mass concentrations
  for the global/unprojected grid

correct? (Though I note that, in the demo

> names(res3D at data)
[1] "var1.pred" "var1.var" 

and I don't know why there are 2 data columns rather than 1.) If so, I
can create data3D and get its coordinates from my input grid, and
obtain grid3D by projecting coordinates(data3D) using R package=M3.
But to get res3D from grid3D, for my usecase, should I use

* gstat::krige? If so, what model? (or no model?)

* gstat::idw? If so, what formula?

* or something else?

Your assistance is appreciated! Tom Roche <Tom_Roche at pobox.com>