[R-sig-Geo] function to convert 3D geographical coordinates to Cartesian?

[Tom Roche]

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?
>>> The spherical-to-Cartesian math is straightforward enough, but
>>> there's certainly scope for error, so I'd prefer to use tested,
>>> previously-used code if available.

https://stat.ethz.ch/pipermail/r-sig-geo/2013-January/017110.html
(rearranged)
> your apparent elevation at the geographical coordinate point is
> crucially dependent on the datum.

Fortunately, the atmospheric model with which I'm working assumes a
spherical earth. So (If I Understand Correctly--please correct if not),
I have no such dependency; hence proposed code @ end of

https://stat.ethz.ch/pipermail/r-sig-geo/2012-December/017105.html

In the general case, the user should provide a CRS (more below).

> If you look at 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.

So spTransform() is, by design, 2D-only? <sigh/>

> The most robust route forward might be to permit 3D SpatialPoints to
> be transformed directly, and choose a geo-centric target projection.

That's also my impression:

https://stat.ethz.ch/pipermail/r-sig-geo/2012-December/017105.html
>> I'm looking for code [which]
...
>> * _properly_ converts elevation to radial distance (from earth
>>   center), i.e., using a user-specified CRS. My function
>>   elev2rdist below is a blatant kludge, though it should work for
>>   the model for which I'm assimilating data.

> If you have a data set case with both the input and matching output
> data, I could take a look.

I don't have output data yet, nor even data in the form my output must
match (though I should have the latter soon, and I believe I know its
form--more below). I do have input data, with which I'm working @

https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na

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

(The input is, and the output must be, netCDF.) After considerable
research, ISTM the current best tool for this job is NCL

http://www.ncl.ucar.edu/

since it has API for

* converting HσP to elevation

* calling ESMF routines that (apparently) most conservatively 3D-regrid

But NCL and ESMF consider LCC (at least in the form in which my output
uses it) to be a Cartesian grid, so before I can use that, I must (IIUC)

1. horizontally crop the input to the output's bounds. Fortunately R
   package=raster makes this easy.

2. convert the cropped input grid's vertical coordinates from HσP to
   elevation. I can do this with NCL.

3. convert the cropped input grid's from geographical coordinates to
   Cartesian coordinates.

Step 3 motivates my request. Fortunately I'm working with spheres, and
there are several available 3D spherical-to-Cartesian routines (and even
I can handle that math :-) so for now I can do something like the code @

https://stat.ethz.ch/pipermail/r-sig-geo/2012-December/017105.html

But, longer-term, It Would Be Nice if someone provided a proper (i.e.,
taking a user-provided CRS argument) 3D geographical-to-Cartesian
conversion.

Your assistance is appreciated (even if only useful in future),
Tom Roche <Tom_Roche at pobox.com>