How difficult could it be two project a two dimensional image onto a three dimensional sphere in R? Well, it's not really that difficult, but it took me quite a while to figure out how to achieve this. Now why would we want to do this in the first place? Well, you could use such a technique to create a three dimensional model of Earth on which you can plot all sorts of interesting global phenomena. Well, I couldn't think of many other useful applications besides modelling fancy Christmas ornaments. For now we'll stick to the model of Earth.
This blog was rewritten several times, before publication, as I tried several strategies with varying degrees of success. I started with the persp function from the graphics package. Which just didn't work. This function can plot a single hemisphere, with limited control over the colours used for its surface.
Then I ended up using the spheresurf3D function from the plot3D package. Although I finally got it to work, it required a horrific workaround. The plot3D package uses a matrix of colour values that will be used in combination with a colour ramp to visualise on the surface. My workaround was to get all unique values from the target image and create a matrix with indices for those colours. I will not present this strategy in detail here due to this monstrosity and the other graphical limitations of this package. Don't get me wrong, I think plot3D is a great package, with many useful applications, it's just not the ultimate tool for what I want here.
The easiest way to create a three dimensional model of the globe was with the rgl package. But first, some code to prepare ourselves. The code below downloads a topographic map from NASA's website, which we will use to project onto a sphere. We convert the jpeg image into a png file, as the rgl package will only accept that file type. Note that free third party software (ImageMagick) is used for this step. We also specify a function to convert polar coordinates (longitude, latitude and altitude) into Cartesian coordinates. We also specify a matrix with Cartesian coordinates (‘world_coords’) which we will use to georeference the topographic image.
| require(rgl) | |
| # Download a topographic image of Earth from NASA's website | |
| # The file downloaded here is actually the lowest available | |
| # resolution. Higher resolutions are available. | |
| download.file("http://eoimages.gsfc.nasa.gov/images/imagerecords/73000/73751/world.topo.bathy.200407.3x5400x2700.jpg", "world.topo.jpg", mode = "wb") | |
| # Use third party software (ImageMagick) to convert the jpeg | |
| # into a png file (rgl can only work with png). | |
| # I also reduce the size and number of colours to achieve an | |
| # acceptable file size | |
| shell("convert world.topo.jpg -resize 50% -colors 64 world.topo.png") | |
| # function to convert polar to Cartesian coordinates, | |
| # assuming that the world is a perfect sphere | |
| polar_to_cart <- function(long, lat, alt = 1) | |
| { | |
| x <- outer(long*pi/180, lat*pi/180, function(x, y) alt*cos(y)*cos(x)) | |
| y <- outer(long*pi/180, lat*pi/180, function(x, y) alt*cos(y)*sin(x)) | |
| z <- outer(long*pi/180, lat*pi/180, function(x, y) alt*sin(y)) | |
| return (list(x = x, y = y, z = z)) | |
| } | |
| # Specify reference coordinates to project our image onto | |
| lat <- seq(90,-90, len=50) | |
| long <- seq(-180, 180, len=50) | |
| # get Cartesian coordinates from the specified longitudes and latitudes | |
| world_coords <- polar_to_cart(long, lat) |
Next, we set up an rgl device and our model. We create an environment with a black background (space), and set up a light to represent the sun. With the persp3d function we create a model of Earth. We use the ‘texture’ argument to project the topographic image onto the sphere. We add a slightly larger transparent light blue sphere to represent the Earth's atmosphere. Finally, we add a small red sphere hovering above my home country. That's it! We have a simple model of our planet. The nice thing of rgl devices is that they are interactive. You can rotate your model and zoom in and out. With the writeWebGL you can even save your interactive model as a web page. For me, this web page worked offline, but due to restrictions of blogspot, I could get it to work on my blog.
| # set up an rgl device with a black background | |
| open3d(windowRect = c(100, 100, 250, 250)) | |
| bg3d("black") | |
| # Remove all lights present and add a light to represent the sun | |
| clear3d(type = "lights") | |
| rgl.light(theta = 30, phi = 0, viewpoint.rel = T, ambient = "#FFFF66") | |
| # Create a sphere on which the world topography is projected | |
| persp3d(world_coords[["x"]], world_coords[["y"]], world_coords[["z"]], | |
| col="white", | |
| texture="world.topo.png", | |
| specular="#303030", axes=FALSE, box=FALSE, xlab="", ylab="", zlab="", | |
| normal_x = world_coords[["x"]], normal_y = world_coords[["y"]], normal_z = world_coords[["z"]]) | |
| # add a transparent sphere to model the Earth's atmosphere | |
| persp3d(world_coords[["x"]]*1.05, world_coords[["y"]]*1.05, world_coords[["z"]]*1.05, | |
| col = "#ADD8E6", specular = "black", alpha = 0.1, add = T) | |
| mark_pos <- polar_to_cart(long = 5.5, lat = 52, alt = 1.06) | |
| spheres3d(mark_pos[["x"]], mark_pos[["y"]], mark_pos[["z"]], | |
| 0.03, col = "red", add = T) |
In order to produce something that I can show here, I created an animation. Therefore I first changed the orientation of the model, and then animated a spinning world with the movie3d function.
| # Rotate the model with the north pole up | |
| rgl.viewpoint(zoom = .60, theta = 180, phi = 90) | |
| # Tilt the rotational axis with 23.5 degrees | |
| # with a slighlty better view of the Northern hemisphere | |
| U <- par3d("userMatrix") | |
| par3d(userMatrix = rotate3d(U, pi*23.5/180, -0.7,-0.3,0)) | |
| # Spin the model and save movie frames | |
| movie3d(spin3d(axis = c(0, 0, 1), rpm = 15), duration = 4, | |
| convert = F, dir = getwd()) | |
| # Convert movie frames to animated GIF | |
| shell("convert movie*.png -delay 20 -dither FloydSteinberg -colors 128 globe.gif") | |
| # remove individual animation frames after creating GIF | |
| file.remove(list.files()[grepl("^movie.*png$", list.files())]) |
Once you know how to do it, it's not that difficult to create an animation of a rotating planet in R. Of course there is always stuff that I could improve upon. For example, the Earth isn't a perfect sphere, but I had assumed so as a simplification. If you have any refinements for the approach shown here, please let me know in the comments.
No comments:
Post a Comment