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

[Tom Roche]

https://stat.ethz.ch/pipermail/r-sig-geo/2012-December/017102.html
> Is there packaged code to convert geographical coordinates (e.g.,
> longitude, latitude, elevation) to [simple] Cartesian coordinates
> in 3-space?

Just to clarify: I'm looking for code like that which follows my .sig
(to end of post), but which

* is much better tested/known/used (I just hacked the following this
  evening)

* _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.

This is such low-level functionality that I can't believe it's not API
somewhere, but I have looked at several packages and am not seeing it.
Is it so low-level that it's just not exposed as API? or am I just not
seeing it in package docs? or did I fail a spatial-intelligence test?

TIA, Tom Roche <Tom_Roche at pobox.com>----sample code follows to EOP----

# Convert geographic coordinates (lon, lat, elevation) to
# Cartesian coordinates (x, y, z) relative to earth center.
geo2cart <- function(lon.lat.elev) {
  library(pracma) # for sph2cart
  stopifnot(is.numeric(lon.lat.elev))
  return(sph2cart(geo2sph(lat.lon.elev)))
} # end function geo2cart(lon.lat.elev)

# Convert geographic coordinates (lon, lat, elevation) to
# spherical coordinates (azimuth/Θ, polar angle/φ, radial distance/r),
# relative to earth center
geo2sph <- function(lon.lat.elev) {
  stopifnot(is.numeric(lon.lat.elev))

  # Use convention of package=pracma: see http://cran.r-project.org/web/packages/pracma/
  # same as described @ http://en.wikipedia.org/wiki/File:3D_Spherical_2.svg
  # convert
  # * longitude -> azimuth (Θ)
  # * latitude -> polar angle (φ)
  # * elevation -> radial distance from earth center (r)
  if (is.vector(lon.lat.elev) && length(lon.lat.elev) == 3) {
    theta <- lon2azi(lon.lat.elev[1])
    phi <- lat2pol(lon.lat.elev[2])
    r <- elev2rdist(lon.lat.elev[3])
    m <- 1
  } else if (is.matrix(lon.lat.elev) && ncol(lon.lat.elev) == 3) {
    theta <- lon2azi(lon.lat.elev[,1])
    phi <- lat2pol(lon.lat.elev[,2])
    r <- elev2rdist(lon.lat.elev[,3])
    m <- nrow(lon.lat.elev)
  } else {
    stop('geo2sph: ERROR: input must be a vector of length 3 or a matrix with 3 columns.')
  }

  if (m == 1) {
    tpr <- c(theta, phi, r)
  } else {
    tpr <- cbind(theta, phi, r)
  }
  return(tpr)
} # end function geo2sph(lon.lat.elev)

# Filter and convert longitudes to radians.
# Much-too-simple tests:
# > lon2azi(c(seq(0, 180), seq(-179, -1)))/pi
# > lon2azi(c(0:359))
lon2azi <- function(lon) {
  library(pracma) # for deg2rad

  stopifnot(is.numeric(lon))
  # only handle vectors
  if (!is.vector(lon)) {
    lon.vec <- c(lon) 
  } else {
    lon.vec <- lon
  }

  # only take -180 <= lon.vec <= 180
  # TODO: better error messages
  stopifnot(length(lon.vec[lon.vec < -180]) == 0)
  stopifnot(length(lon.vec[lon.vec > 180]) == 0)

  # Convert to azimuth=Θ in radians.
  # Convert to 0 <= lon.vec <= 360 (hmm ...) e.g., -179 -> 181.
  lon.vec[lon.vec < 0] <- 360.0 + lon.vec[lon.vec < 0]
  theta.vec <- deg2rad(lon.vec)

  if (!is.vector(lon)) {
    return(theta.vec[1])
  } else {
    return(theta.vec)
  }
} # end function lon2azi(lon)

# Filter and convert latitudes to radians.
# Much-too-simple tests:
# > lat2pol(c(seq(90, 0), seq(-1, -90)))/pi
# > lat2pol(c(-100:100))
lat2pol <- function(lat) {
  library(pracma) # for deg2rad

  stopifnot(is.numeric(lat))
  # only handle vectors
  if (!is.vector(lat)) {
    lat.vec <- c(lat) 
  } else {
    lat.vec <- lat
  }

  # only take -90 <= lat.vec <= 90
  # TODO: better error messages
  stopifnot(length(lat.vec[lat.vec < -90]) == 0)
  stopifnot(length(lat.vec[lat.vec > 90]) == 0)

  # Convert to polar angle=φ in radians.
  # northern hemisphere + equator:
#  lat.vec[lat.vec >= 0] <- 90.0 - lat.vec[lat.vec >= 0]
  # southern hemisphere:
#  lat.vec[lat.vec < 0] <- 90.0 - lat.vec[lat.vec < 0]
  phi.vec <- deg2rad(90.0 - lat.vec)

  if (!is.vector(lat)) {
    return(phi.vec[1])
  } else {
    return(phi.vec)
  }
} # end function lat2pol(lat)

# Convert elevation to radial distance
# Much-too-simple tests:
# > elev2rdist(c(seq(1000, by=1000, length=10)))
# > elev2rdist(c(seq(-1000, by=1000, length=10)))
elev2rdist <- function(elev) {
  r.earth.meters <- 6370000.0 # constant, kludge, see below

  stopifnot(is.numeric(elev))
  # only handle vectors
  if (!is.vector(elev)) {
    elev.vec <- c(elev) 
  } else {
    elev.vec <- elev
  }

  # Only take non-negative elevations.
  stopifnot(length(elev.vec[elev.vec < 0]) == 0)

  # Convert elevation to radial distance (from earth center).
  # TODO: use coordinate reference system.
  # Kludge: spherical earth with radius defined above.
  r.vec <- r.earth.meters + elev

  if (!is.vector(elev)) {
    return(r.vec[1])
  } else {
    return(r.vec)
  }
} # end function elev2rdist(elev)