5

I would like to find a method for plotting an exaggerated geoid (e.g. topography):

enter image description here

I primarily use R for my work but would be interested in other open source software options (e.g. QGIS w/ GRASS). I found this wonderful tool for MATLAB (link) but haven't discovered other options.

  • R's rgl could do this, can you add a link to the data and I will have a go. – mdsumner Mar 24 '14 at 19:09
  • @mdsumner - I had actually thought of that and posed a similar question on SO (stackoverflow.com/questions/22607406/…). There is some sample R data over there - maybe you could just add some noise to the z-dimension to make it "lumpy". Would be great if you could figure out a way! – Marc in the box Mar 24 '14 at 19:13
6

EDIT: I've updated this to do an actual surface.

It's interactive with rgl, and you can zoom in to see the closed surface but you'll need more work to respect the actual WGS84 datum and get your vertical exaggeration just right.

Download the files with R:

baseurl <- "http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/GIS/world_geoid"

fs <- file.path(baseurl, c("n45w180.zip", "n45w135.zip", "n45w90.zip", "n45w45.zip",    "n45e00.zip", "n45e45.zip", 
"n45e90.zip", "n45e135.zip", "n00w180.zip", "n00w135.zip", "n00w90.zip", 
"n00w45.zip", "n00e00.zip", "n00e45.zip", "n00e90.zip", "n00e135.zip", 
"s45w180.zip", "s45w135.zip", "s45w90.zip", "s45w45.zip", "s45e00.zip", 
"s45e45.zip", "s45e90.zip", "s45e135.zip", "s90w180.zip", "s90w135.zip", 
"s90w90.zip", "s90w45.zip", "s90e00.zip", "s90e45.zip", "s90e90.zip", 
 "s90e135.zip"))

for (i in seq_along(fs)) download.file(fs[i], basename(fs)[i], mode = "wb")

 ## unzip
for (i in seq_along(fs)) utils::unzip(basename(fs[i]))

## build the next string with R too
txt <- gsub(".zip", "", basename(fs))
 cat(paste(file.path(txt, txt, "w001001.adf"), collapse = " "), "\n")

Mosaic them all with GDAL:

gdalbuildvrt geoid.vrt n45w180/n45w180/w001001.adf n45w135/n45w135/w001001.adf      n45w90/n45w90/w001001.adf n45w45/n45w45/w001001.adf n45e00/n45e00/w001001.adf n45e45/n45e45/w001001.adf n45e90/n45e90/w001001.adf n45e135/n45e135/w001001.adf n00w180/n00w180/w001001.adf n00w135/n00w135/w001001.adf n00w90/n00w90/w001001.adf n00w45/n00w45/w001001.adf n00e00/n00e00/w001001.adf n00e45/n00e45/w001001.adf n00e90/n00e90/w001001.adf n00e135/n00e135/w001001.adf s45w180/s45w180/w001001.adf s45w135/s45w135/w001001.adf s45w90/s45w90/w001001.adf s45w45/s45w45/w001001.adf s45e00/s45e00/w001001.adf s45e45/s45e45/w001001.adf s45e90/s45e90/w001001.adf s45e135/s45e135/w001001.adf s90w180/s90w180/w001001.adf s90w135/s90w135/w001001.adf s90w90/s90w90/w001001.adf s90w45/s90w45/w001001.adf s90e00/s90e00/w001001.adf s90e45/s90e45/w001001.adf s90e90/s90e90/w001001.adf s90e135/s90e135/w001001.adf 

And do a quick decimation, choose your own percentage. (This is longlat data so we really need more sophisticated resampling to get data evenly distributed on the globe, doable in R but for another time.)

gdal_translate geoid.vrt geoid04.tif -outsize 4% 4% -co COMPRESS=LZW -co TILED=YES

Now back to R, read that simplified raster and convert to points in XYZ.

library(raster);library(rgdal)
r <- raster("geoid04.tif")

llh <- as.data.frame(r, xy = TRUE)

## just spherical
llh2xyz <- function(lonlatheight, rad = 500) {
cosLat = cos(lonlatheight[,2] * pi / 180.0)
sinLat = sin(lonlatheight[,2] * pi / 180.0)
cosLon = cos(lonlatheight[,1] * pi / 180.0)
sinLon = sin(lonlatheight[,1] * pi / 180.0)

x = rad * cosLat * cosLon
y = rad * cosLat * sinLon
z = lonlatheight[,3] + rad * sinLat

cbind(x, y, z)
}

## points on the geoid
xyz <- llh2xyz(llh, rad = 800)

## colour mapping scale
scl <- function(x, nn = 50) {
1 + (nn-1) * ((x[!is.na(x)] - min(x,na.rm = TRUE))/diff(range(x, na.rm = TRUE)))
}

Build geometry

library(geometry)
tbr = t(surf.tri(xyz, delaunayn(xyz)))

And plot

library(rgl)

library(maptools)
data(wrld_simpl)
m <- coordinates(as(as(wrld_simpl, "SpatialLines"), "SpatialPoints"))
m <- cbind(m, extract(r, m))

mxyz <- llh2xyz(m, rad = 820)
n <- 150
rgl.triangles(xyz[tbr,1], xyz[tbr,2], xyz[tbr,3],col = terrain.colors(n)[scl(llh[tbr, 3], nn = n)])

points3d(mxyz[,1], mxyz[,2], mxyz[,3], col = "black")

enter image description here

  • This is really excellent stuff! I'm learning a lot here. As far as I can tell, your the first to do this in R. Cheers – Marc in the box Mar 25 '14 at 11:46
4

Here's a few improvements on previous answer.

  1. Use randomCoordinates to more evenly distribute the points, and save creating any intermediate rasters.
  2. Just a 2D triangulation, simpler and more directly what we want, though there's more required, you need work to "seal" the surface properly by triangulating in a smarter way.
  3. To really get the exaggeration in that original plot you'd need to exaggerate the relief about the mean surface, not just exaggerate it straight as I have done.

(This has been great, thanks for the question, I learnt a lot. )

library(raster)
library(rgdal)
r <- raster("geoid.vrt")


library(geosphere)
llh <- data.frame(randomCoordinates(50000))
llh$h <- extract(r, llh[,1:2])
## just spherical
 llh2xyz <- function(lonlatheight, rad = 500, exag = 1) {
cosLat = cos(lonlatheight[,2] * pi / 180.0)
sinLat = sin(lonlatheight[,2] * pi / 180.0)
cosLon = cos(lonlatheight[,1] * pi / 180.0)
sinLon = sin(lonlatheight[,1] * pi / 180.0)

x = rad * cosLat * cosLon
y = rad * cosLat * sinLon
z = (lonlatheight[,3] * exag) + rad * sinLat

cbind(x, y, z)
}


## triangulate first in lonlat
library(geometry)

tbr <- delaunayn(llh[,1:2])
tri.ind <- as.vector(t(tbr))

# points on the geoid
myrad <- 6378137
myexag <- 8000
xyz <- llh2xyz(llh, rad = myrad, exag = myexag)

# colour mapping scale
scl <- function(x, nn = 50) {
1 + (nn-1) * ((x[!is.na(x)] - min(x,na.rm = TRUE))/diff(range(x, na.rm = TRUE)))
}

n <- 150
## those colours, closer to original
jet.colors <-
   colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                      "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

 rgl.triangles(xyz[tri.ind,1], xyz[tri.ind,2], xyz[tri.ind,3],col = jet.colors(n)[scl(llh[tri.ind, 3], nn = n)])

library(maptools)
data(wrld_simpl)
m <- coordinates(as(as(wrld_simpl, "SpatialLines"), "SpatialPoints"))
m <- cbind(m, extract(r, m))

mxyz <- llh2xyz(m,rad = myrad + 2000, exag = myexag)
n <- 150

points3d(mxyz[,1], mxyz[,2], mxyz[,3], col = "black")

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.