Source

MOZART-global to AQMEII-NA / cartesianize.elevations.r

Full commit
# R code, driven by cartesianize_elevations.sh, to
# * read a matrix of points with schema=(longitude, latitude, elevation above mean sea level (AMSL))
# * "Cartesianize" those points to (x,y,z) relative to center of CMAQ spherical earth
# * write those points to CSV for pass to NCL
# TODO: pass parameters via commandline, not environment variables

# constants-----------------------------------------------------------

# all the following env vars must be set and exported in driver script (TODO: fail if not)
this.fn <- Sys.getenv('CALL_R_FN')  # FIXME: gets overidden by `source` below
# cat(sprintf('this.fn=%s\n', this.fn)) # debugging

# local constants, from environment-----------------------------------
# taking from Rscript call failed, annoyingly and mysteriously :-(

# cat(sprintf('%s: Setting constants from environment ...\n', this.fn)) # debugging

# source grids (for this script)
sg.ascii.fp <- Sys.getenv('SG_ASCII_FP') # path to input CSV file
col.names.str <- Sys.getenv('SG_ASCII_COL_NAMES_CSV_ORDERED')
col.names.vec <- c(unlist(strsplit(col.names.str, ',')))

# target grids for this script, intermediate grids from driver perspective
ig.ascii.fp <- Sys.getenv('IG_ASCII_FP') # path to output CSV file

# note: failing to quote following "horizontal line" produced
# intractable errors="invalid argument to unary operator" in following statements
# ----------------------------------------------------------------------

# # start debugging
# # TODO: check paths
# cat(sprintf('sg.ascii.fp==%s\n', sg.ascii.fp))
# cat(sprintf('col.names.str==%s\n', col.names.str))
# cat(sprintf('ig.ascii.fp==%s\n', ig.ascii.fp))
# # cat(sprintf('%%%==%s\n', %%%))
# #   end debugging

# functions-----------------------------------------------------------

# 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(pracma::sph2cart(geo2sph(lon.lat.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)

# code----------------------------------------------------------------

# cat(sprintf('%s: Loading package=pracma ...\n', this.fn)) # debugging
library(pracma)

# 1. Load the source grid from CSV

lon.lat.z <- as.matrix(
  read.csv(sg.ascii.fp, header=FALSE,col.names=col.names.vec))
# start debugging
cat('\n') # separator
cat(sprintf('%s: lon.lat.z: head, tail:\n', this.fn))
head(lon.lat.z)
tail(lon.lat.z)
#   end debugging
# > head(lon.lat.z)
#       lon      lat        z
# [1,]   -130 59.68421 1157.416
# [2,]   -130 59.68421 1281.332
# [3,]   -130 59.68421 1399.451
# [4,]   -130 59.68421 1518.972
# [5,]   -130 59.68421 1640.001
# [6,]   -130 59.68421 1762.510
# > tail(lon.lat.z)
#            lon      lat        z
# [35905,]  -57.5 21.78947 36905.33
# [35906,]  -57.5 21.78947 38412.34
# [35907,]  -57.5 21.78947 39982.21
# [35908,]  -57.5 21.78947 41617.34
# [35909,]  -57.5 21.78947 43320.32
# [35910,]  -57.5 21.78947 45093.82

# 2. Cartesianize using functions above

x.y.z <- geo2cart(lon.lat.z)
# start debugging
cat('\n') # separator
cat(sprintf('%s: x.y.z: head, tail:\n', this.fn))
head(x.y.z)
tail(x.y.z)
#   end debugging
# > head(x.y.z)
#             x        y       z
# [1,] -3535295 -4213201 3215941
# [2,] -3535364 -4213283 3216003
# [3,] -3535429 -4213361 3216063
# [4,] -3535496 -4213440 3216123
# [5,] -3535563 -4213520 3216184
# [6,] -3535631 -4213601 3216246
# > tail(x.y.z)
#                x        y       z
# [35905,] 1277820 -2005775 5949158
# [35906,] 1278120 -2006247 5950557
# [35907,] 1278433 -2006738 5952015
# [35908,] 1278760 -2007250 5953533
# [35909,] 1279099 -2007783 5955114
# [35910,] 1279453 -2008339 5956761

# 3. Write Cartesianized matrix back to CSV for consumption by NCL
write.table(x=x.y.z, file=ig.ascii.fp,
  sep=',', append=FALSE, quote=FALSE, row.names=FALSE, col.names=FALSE)
# start debugging
ls.cmd <- sprintf('ls -alh %s', ig.ascii.fp)
cat('\n') # separator
cat(sprintf('%s: %s==\n', this.fn, ls.cmd))
system(ls.cmd)
cat('\n') # separator
cat(sprintf('%s: `head %s`==\n', this.fn, ig.ascii.fp))
system(sprintf('head %s', ig.ascii.fp))
# -3535295.18688108,-4213200.73984654,3215940.65364899
# -3535363.9466358,-4213282.68453125,3216003.20211492
# -3535429.48965546,-4213360.79566041,3216062.82442377
# -3535495.81082981,-4213439.8341582,3216123.1545941
# -3535562.9684684,-4213519.86951526,3216184.24566817
# -3535630.94797151,-4213600.88433232,3216246.08436512
# -3535699.77100706,-4213682.90443219,3216308.69039551
# -3535769.50371912,-4213766.0086423,3216372.12392843
# -3535840.16719796,-4213850.22209708,3216436.40414898
# -3535911.74903902,-4213935.53001332,3216501.51977315

cat('\n') # separator
cat(sprintf('%s: `tail %s`==\n', this.fn, ig.ascii.fp))
system(sprintf('tail %s', ig.ascii.fp))
# 1276731.50100593,-2004067.02298057,5944091.3801341
# 1276988.5531689,-2004470.51405335,5945288.14041217
# 1277254.36222637,-2004887.75069712,5946525.67024989
# 1277531.22410427,-2005322.33679373,5947814.65881221
# 1277819.71946488,-2005775.18380034,5949157.80949412
# 1278120.2830282,-2006246.97409069,5950557.14622586
# 1278433.38437031,-2006738.44475157,5952014.85521734
# 1278759.5020598,-2007250.34698521,5953533.16454511
# 1279099.1498614,-2007783.48724068,5955114.46614859
# 1279452.86390567,-2008338.7070744,5956761.25609616
#   end debugging