Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Last active December 8, 2016 15:10
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mdsumner/dc80283de50bb23ff7681b14768b9367 to your computer and use it in GitHub Desktop.
Save mdsumner/dc80283de50bb23ff7681b14768b9367 to your computer and use it in GitHub Desktop.
#' Title
#'
#' @param lonlatheight matrix or data.frame of lon,lat,height values
#' @param rad radius of sphere
#' @param exag exaggeration to apply to height values (added to radius)
#'
#' @return matrix
#' @export
llh2xyz <- function(lonlatheight, rad = 6378137.0, 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)
rad <- (exag * lonlatheight[,3] + rad)
x = rad * cosLat * cosLon
y = rad * cosLat * sinLon
z = rad * sinLat
cbind(x, y, z)
}
## download Blue Marble
u <- "http://bl.ocks.org/rasmuse/raw/75fae4fee3354ec41a49d10fb37af551/0a91d6f75b95b90bda15a074ced45de454bb038f/world.topo.bathy.200411.3x540x270.png"
download.file(u, basename(u), mode = "wb")
library(raster)
library(rgl)
library(quadmesh)
## read RGB and set geo-spatial transform
r <- setExtent(brick(basename(u)), extent(-180, 180, -90, 90))
## extract hex strings for each pixel (quad)
cols <- rgb(r[[1]][], r[[2]][], r[[3]][], max = 255)
qm <- quadmesh(r[[1]]) ## the topology
qm$vb[3,] <- 0 ## drop the z
qm$vb[1:3, ] <- t(llh2xyz(t(qm$vb[1:3, ]))) ## geometry
## vis
shade3d(qm, col = rep(cols, each = 4), specular = "black")
light3d(cbind(0, 0, 0))
bg3d("black")
## this time with textures
library(raster)
library(rgl)
library(quadmesh)
r <- setExtent(brick(basename(u)), extent(-180, 180, -90, 90))
## somewhat different is to forget about the geometry of the image
## and simply texture drape it on whatever surface
library(geosphere)
library(geometry)
pts <- randomCoordinates(2^12)
xyz <- llh2xyz(cbind(pts, 0))
tri <- geometry::convhulln(xyz)
library(rgl)
## rgl boilerplate:
tm <- tetrahedron3d()
tm$vb <- t(cbind(xyz, 1))
tm$it <- t(tri)
tcoords <- xyFromCell(setExtent(r, extent(0, 1, 0, 1)), cellFromXY(r, pts))
## has to be "not-black", which is default . . .
shade3d(tm, texture = basename(u), texcoords = tcoords[tm$it, ], col = "white", specular = "black")
light3d(cbind(0, 0, 0))
bg3d("black")
@mdsumner
Copy link
Author

mdsumner commented Aug 7, 2016

rgl_bm

@rasmuse
Copy link

rasmuse commented Aug 7, 2016

Interesting approach to actually represent it in 3D. Works well for orthographic and for "satellite" projections of course, but not for other projections (cylindrical, conical, etc) if I understand correctly?

Anyway, this approach could also be implemented with 3D WebGL I guess. But I haven't worked with WebGL myself.

Also, for me it's usually important to have coordinate transformations available for overlaying vector data, text, etc. This is a nice thing about the little d3 library I posted the other day -- it works seamlessly with d3's utilities for reprojecting vector data.

@mdsumner
Copy link
Author

@rasmuse it can be done with planar projections too, I'm keen to learn WebGL approaches, I just need to lock away these data structures in R before I focus on that. Feeding spatial stuff to rgl requires most of the same work as feeding it to .js (I'm 99% certain . . .). And, sorry I missed your comment here, I just don't see any notification so I'm looking into that.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment