[R] code to convert 3D geographical coordinates to Cartesian?

[Tom Roche]

summary: I'm looking for packaged, tested code to convert geographical
coordinates (e.g., longitude, latitude, elevation) to Cartesian
coordinates (x,y,z) in 3-space. I know of R code for 2-space and for
spherical <-> Cartesian. This can be extended (I attach a quick kludge
extending pracma::sph2cart), but a proper extension is non-trivial.

details:

https://stat.ethz.ch/pipermail/r-help/2012-December/332658.html
>>>> Is there packaged code to convert geographical coordinates (e.g.,
>>>> longitude, latitude, elevation) to Cartesian coordinates in 
>>>> 3-space? [Since] there's certainly scope for error, [I'd] prefer
>>>> to use tested, well-used code if available.

https://stat.ethz.ch/pipermail/r-help/2012-December/332666.html
>>> Have you checked the spatial stats task view on CRAN?

https://stat.ethz.ch/pipermail/r-help/2012-December/332667.html
>> I have, and can't believe that functionality this fundamental is not
>> API. But I have looked at several packages and am not seeing it.

https://stat.ethz.ch/pipermail/r-help/2012-December/332671.html
(rearranged)
> How did you search?

> install.packages("sos")
> require(sos)
>  > findFn("latitude longitude cartesian")

Thanks! I didn't know about that search mechanism. I was googling via
rseek.org and reading <package/>.pdf for packages that looked useful.
Unfortunately

> sphereplot::sph2car "Transforms 3D spherical coordinates to cartesian
> coordinates"

Note this does spherical <-> Cartesian, not geographical <-> Cartesian.
That can be extended (see a quick hack following my .sig to end of post,
extending Borchers' implementation), but the extension is {non-trivial,
error-prone}. Hence my hope to find code that is much better
tested/known/used.

> GEOmap::Lll2xyz "List Lat-Lon to cartesian XYZ"

GEOmap::Lll2xyz (and all the functions I saw in GEOmap) are 2D-input:
they take arguments=(lat, lon), not (lat, lon, elev). Again, this could
be extended; but, again, the same caveats apply.

> You have a very different notion of what would be considered core  
> capabilities in a statistics program than I have.

Geospatial information comes in a wide variety of formats. Doing
*geospatial* statistics with *real-world data* requires the ability to
transform those various formats. I consider geospatiality to be a core
competency for R, and geographical <-> Cartesian transforms to be
fundamental to that: YMMV.

thanks again, Tom Roche <Tom_Roche at pobox.com>
------------------sample hack follows to end of post------------------

# Convert geographic coordinates (lon, lat, elevation) to
# Cartesian coordinates (x, y, z) relative to earth center.
geo2cart <- function(lon.lat.elev) {
  stopifnot(is.numeric(lon.lat.elev))
  library(pracma) # for sph2cart
  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
  # to 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)