[R-sig-Geo] requesting 3D spTransform()

[Roger Bivand]

On Wed, 2 Jan 2013, Tom Roche wrote:

>
> https://stat.ethz.ch/pipermail/r-sig-geo/2012-December/017099.html
>>> Is there packaged code to convert geographical coordinates (e.g.,
>>> longitude, latitude, elevation) to Cartesian coordinates in 3-space?
>
> https://stat.ethz.ch/pipermail/r-sig-geo/2013-January/017110.html
>> [In] the C code underlying spTransform(), you'll see that the call to
>> pj_transform() sets z to zero, that is to the datum surface, and the
>> returned values are discarded.
>
> ISTM "The Better Way" would be to add a boolean argument: e.g.,
> spTransform(..., 3d=false, ...)
>
> * make the current behavior the default
> * for 3d=true, don't discard the elevations
>
> Am I missing something?
>
>> The most robust route forward might be to permit 3D SpatialPoints to
>> be transformed directly, and choose a geo-centric target projection.
>
> Is there a way to request this? I note package=sp allows FRs
>
> https://r-forge.r-project.org/tracker/?atid=4032&group_id=1014&func=browse
>
> but I'm not seeing equivalent UI for package=rgdal.

It is much preferred for users to submit patches, or at least workable 
examples. Wishes tend only to vex developers, and the sp list is usually 
ignored - all the useful communication happens on R-sig-geo. I'm assuming 
that you can install rgdal from source, that is that you can check out the 
current R-forge rgdal SVN (revision 421), and install it. There you will 
find that a SpatialPoints object with a third dimension will now be 
transformed as you seem to think you need, since pj_transform() does take 
a z vector. SpatialLines and SpatialPolygons objects are limited to 2D and 
will not be accommodated.

I'm concerned that two of your steps make your workflow debatable, one 
that you are not actually dealing with height (which on the sphere is of 
no consequence); it seems that pj_transform is useful where z changes when 
datum transformation occurs. If both your cases are spherical, there is no 
datum transformation, so admitting z in spTransform() is a red herring. 
I'm also concerned that you appear to expect the input and output to be 
gridded, but after projection they will not be gridded without warping. In 
addition, you need to interpolate to a finer grid.

I'm willing to hazard the view that you would be better placed using 2D 
interpolation methods warping your model values to the finer projected 
grid after having projected the model input to irregular points. That is, 
3D interpolation is still a very different animal from what you need to 
do. Unfortunately both atmospheric and oceanographic models use their own 
dimensions in their own grids, none of which map easily into the real 
world. Please do try several approaches and let us know which actually 
works - 3D datum transformation is now implemented in rgdal, but I don't 
think it is what you really need. The real problem is change of support in 
3D, in addition to interpolation to a finer 2/3D grid.

Worry about the output format only once you reach a sensible outcome (in 
reply to your further NetCDF question). My guess would be that some form 
of hierarchical geostatistical model for interpolation is needed to carry 
out the change of support. Please wait before overwhelming people with 
your problems, I'd be grateful now for feedback on spTransform() only. If 
you cannot help with that, change your OS and get back when you are able 
to contribute, nobody here gets paid for answering!

Please also always quote verbatim your full thread - you've been dropping 
things, here starting a new thread (sigh!!), and wasting time because I 
have to scan several days' mail to find your workflow:

"My usecase is like this:

I need to conservatively 3D-regrid atmospheric concentration data from a
global inventory (aka "the input") which

- has lon-lat horizontal coordinates

- has hybrid sigma-pressure

http://en.wikipedia.org/wiki/Sigma_coordinate_system#hybrid_sigma-pressure

   (aka "HσP") vertical coordinates

- assumes a spherical earth

to a regional space (aka "the output") which

+ is Lambert conformal conic horizontal (at much finer resolution than
   the input)

+ is (I believe) HσP vertical (with different vertical levels than the
   input)

+ assumes a spherical earth"

OK?

Roger

>
> TIA, Tom Roche <Tom_Roche at pobox.com>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no