Message Boards Message Boards

On making a 3D globe

Posted 8 years ago

A 3D globe

I'd always wanted to make a 3D globe in Mathematica. Having only recently figured out how to go about it, and surprised by how little code it took to do, I've elected to share this nice little demo.

The heart of this implementation is the nifty function GeoElevationData[] (new in version 10), which allows one to generate an array of elevation values for any point on the Earth. From this, one can either choose to generate a relief map for textures, or as in the image above, an actual 3D globe.

To start off, here's the raw data to be used by the subsequent plots:

gdm = Reverse[QuantityMagnitude[GeoElevationData["World", "Geodetic", GeoZoomLevel -> 3, UnitSystem -> "Metric"]]];

(Increase the setting of GeoZoomLevel, or change the UnitSystem as seen fit.)

From this, one can use produce a quick visualization using ReliefPlot[], which could also be used as a nice texture:

ReliefPlot[gdm, BoxRatios -> {2, 1, 1/2}, ColorFunction -> "HypsometricTints", ColorFunctionScaling -> False, 
           DataRange -> {{-180, 180}, {-90, 90}}, Frame -> False, PlotRangePadding -> None]

relief map

(I will note at this juncture that this is effectively what's being done behind the scenes when you use GeoStyling["ReliefMap"] in GeoGraphics[].)

But why be contented with a flat map, when we have actual elevation data! To make plotting slightly easier, build an interpolating function from the elevation data:

et = ListInterpolation[gdm, {{-90, 90}, {-180, 180}}];

With that, it is now a simple matter to generate a nice globe:

With[{s = 2*^5},
     ParametricPlot3D[(1 + et[?, ?]/s) {Cos[? °] Cos[? °],Cos[? °] Sin[? °], Sin[? °]}, {?, -180, 180}, {?, -90, 90},
                      Axes -> None, Boxed -> False, Mesh -> False, ColorFunctionScaling -> False,
                      ColorFunction -> (With[{r = Norm[{#1, #2, #3}]}, ColorData["HypsometricTints", s r - s]] &),
                      MaxRecursion -> 1, PlotPoints -> {1000, 500}, ViewPoint -> {-1.3, 2.4, 2.}]]

some globe

where I had chosen a scaling factor s that is smaller than the Earth's radius, to make the depressions and elevations slightly more prominent.

(Might be a bit hard to 3D-print, tho. Some smoothing and seam closing would be necessary)


Not being contented with the default ColorData["HypsometricTints"], I decided to search for other nice color gradients. I found a few nice ones from cpt-city, but ultimately decided to synthesize my own using pieces of other gradients in that archive. The result was the following:

myGradient1 = Blend[{{-8000, RGBColor["#000000"]}, {-7000, RGBColor["#141E35"]}, {-6000, RGBColor["#263C6A"]},
                     {-5000, RGBColor["#2E5085"]}, {-4000, RGBColor["#3563A0"]}, {-3000, RGBColor["#4897D3"]},
                     {-2000, RGBColor["#5AB9E9"]}, {-1000, RGBColor["#8DD2EF"]}, {0, RGBColor["#F5FFFF"]},
                     {0, RGBColor["#699885"]}, {50, RGBColor["#76A992"]}, {200, RGBColor["#83B59B"]}, {600, RGBColor["#A5C0A7"]},
                     {1000, RGBColor["#D3C9B3"]}, {2000, RGBColor["#D4B8A4"]}, {3000, RGBColor["#DCDCDC"]},
                     {5000, RGBColor["#EEEEEE"]}, {6000, RGBColor["#F6F7F6"]}, {7000, RGBColor["#FAFAFA"]},
                     {8000, RGBColor["#FFFFFF"]}}, #] &;

Using this as the replacement for "HypsometricTints" as the ColorFunction yields the image at the beginning of my post.

As a bonus image, I'd also found the colors needed to (more or less) reproduce ETOPO1; I'll omit them here, but show the resulting image (which some of you from Stack Exchange may have seen as my current Gravatar there):

ETOPO1

As a bonus bonus, one could also probably look up DEM data for other (terrestrial) planets and satellites, and render them using similar methods. I was able to generate a model of the Moon in this manner:

now that's amore

(I was prompted to post this, since Marco sorta kinda beat me to a nice and similar visualization in this thread. :D)

POSTED BY: J. M.
10 Replies

enter image description here - you earned "Featured Contributor" badge, congratulations !

This is a great post and it has been selected for the curated Staff Picks group. Your profile is now distinguished by a "Featured Contributor" badge and displayed on the "Featured Contributor" board.

POSTED BY: EDITORIAL BOARD

Thanks for sharing!

POSTED BY: Sander Huisman

Thanks for the great post. Similar to a previous post on screen flickering when using AnatomyPlot3D (thread), this ParametricPlot3D code consistently crashes Mathematica on my 2016 MacBook (have reported the crash to Apple and will submit a bug report to Wolfram).

POSTED BY: Timothy Ewing
Posted 8 years ago

Can you try with GeoZoomLevel -> 1 in GeoElevationData[] and a reduced setting of PlotPoints in ParametricPlot3D[]?

POSTED BY: J. M.

Hi JM,

yes indeed nice post. I woner whether we should combine our two plots and show how much elevation predicts the gravitational field, i.e. plot a normalised ratio g-strength/elevation. Plroblem is that your data has a much higher resolution.

By the way, the MacBook issue does not seem to go away using a reduced GeoZoomLevel. I don't think ist purely a performance issue. My old MacBook (2015) has. o problems, my new one (2016) does have serious issues under MMA11.

Cheers,

M.

PS Another nice representarion is when you "deform" the earth under atmospheric pressure, i.e higher pressure decreases the radius). Goucan make a little animation of how the earth deforms.

POSTED BY: Marco Thiel
Posted 8 years ago

Regarding

I wonder whether we should combine our two plots and show how much elevation predicts the gravitational field, i.e. plot a normalised ratio g-strength/elevation. Problem is that your data has a much higher resolution.

and

Another nice representation is when you "deform" the earth under atmospheric pressure, i.e higher pressure decreases the radius). Goucan make a little animation of how the earth deforms.

it should be doable, and if need be we can fall back on the (high order!) spherical harmonic expansions. Let me think about it.

POSTED BY: J. M.

Using GeoZoomLevel -> 1 with PlotPoints -> {1000, 500} works. Also, using GeoZoomLevel -> 3 with PlotPoints -> {500, 250} works. Rotating the resulting globes causes the previously noted screen flickering issue (thread). Edit: clarification that the code works without crashing Mathematica.

Thanks for the suggestions J.M.

POSTED BY: Timothy Ewing

Did tech support get back to you about this issue? I was told that they could not reproduce the problem, and suggested that the difficulty was with the complexity of the plot. This may be true, but I don't think that a graceful way to handle a complex plot is to crash to the finder.

I have this problem on both a 2 year old 13 Inch Retina MacBook pro, and a 3 year old 15 inch Retina macBook Pro with dedicated graphics. I doubt that newer hardware is that much faster/more capable.

Just curious.

I did not hear from tech support about this issue. I did get a response from them about the screen flickering issue. They had me install and test a beta version that seems to have fixed the flickering issue--still very slow to rotate the 3D anatomy models on the 2016 MacBook, but there are now no noticeable visual artifacts.

So, just ran this 3D globe again on the beta version and had the same problem when I reached the fourth code block: CPU 99%, RAM ~4GB, spinning beach ball with about 2 minutes hanging, and then a complete crash which generated a report to Apple dialog box. Sent all of this info on to the same support person as is covering the screen flickering issue.

POSTED BY: Timothy Ewing

Wonderful!

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract