requesting 3D spTransform()
On Wed, 2 Jan 2013, Tom Roche wrote:
Is there packaged code to convert geographical coordinates (e.g., longitude, latitude, elevation) to Cartesian coordinates in 3-space?
[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