Introduction

One of the activities that I use in my 100-level earthquakes and society class is an exploration of Earth’s topography using Mathematica. Mathematica provides easy access to elevation data and has some relatively simple to use commands for making rudimentary maps. They are not competitive with a nice GMT (Generic Mapping Tools) product, but they are useful for exploring concepts. I keep things simple and have students simply execute the notebook commands to produce topographic profiles along the paths following a constant latitude (e.g. 30N and 30S are interesting). This keeps the commands short and limits their focus on the code. For this note I will use a profile along the latitude of Penn State, which is located near the geographic center of Pennsylvania (latitude is 40.795N).

The Code

We start by creating a reference map that shows the location of the latitudinal profile using a red line. Here’s the Mathematica code:

(* create a map and show where the latitude "slide" is *)
latitude = 40.795;
p0 = Framed[
  GeoGraphics[GeoRange -> "World", GeoProjection -> "CylindricalEquidistant",
  GeoZoomLevel -> 2, GeoBackground -> "ReliefMap", ImageSize -> Full,
  GeoGridLinesStyle -> Directive[Red, Thick], 
  GeoGridLines -> {{latitude}, None}]]

and here is the map:

Mathematica world map with a red reference line at a latitude equal to 40.795N, the latitude University Park, PA.

Now that we have defined our profile, we have to extract the elevation data from the Wolfram data servers. The process is not difficult, you simply make a list of latitude and longitudes where you want samples and then ask for the values. Here’s the code to acquire the elevation data and make a very rudimentary plot.

(* Get the elevation data along the red line of constant latitude *)
elevations = 
  GeoElevationData[
   GeoPosition[{latitude, #}] & /@ Range[-180, 180, 0.25], 
   GeoZoomLevel -> 3];

(* plot the elevations along the line and a dashed line for sea level *)
p1 = ListPlot[elevations, Joined -> True, 
   Frame -> {False, False, False, False}, Axes -> False, 
   Filling -> Bottom, ImageSize -> Full, AspectRatio -> 0.2, 
   Prolog -> {Dashed, Line[{{0, 0}, {Max[elevations][[1]], 0}}]}]

and here is the profile:

Profile across the Earth at a latitude of 40.795N. The dashed line represents sea level. The vertical exaggeration is extreme.

The command to get an elevation value is GetElevationData, the command below makes a list of latitude and longitude values (GeoPositions in Mathematica) that are handed to GetElevationData, and the elevation values are returned. The list of positions is constructed by this part of the first command,

GeoPosition[{latitude, #}] & /@ Range[-180, 180, 0.25],

which creates a list of elevation data at the latitude of interest from -180 to +180 longitude, with a sample interval of 0.25 degrees. Elevation values are returned as a Structured Array with elevation values in feet, but can easily be converted to meters. Students can right-click on the image in Mathematica to get a “live” cursor to show them the elevation in feet.

If you want elevations returned as a list of {latitude, longitude, elevation} in meters, use

(* Get the elevation data along the red line of constant latitude shown above *)

elevations = 
  GeoElevationData[
   GeoPosition[{latitude, #}] & /@ Range[-180, 180, 0.25], Automatic, 
   "GeoPosition", GeoZoomLevel -> 3];

But you will have to change the plot command to show that (not by much).

Summary

Mathematica makes it reasonably simple to extract and display topographic data.I use the simple charts to discuss the surface features on the planet, including ridges, islands, trenches, continents, mountains, etc. The examples above are simple, but in an introductory class, I simply want to show students one of the many things that geoscientists do to complete their work – use computers to explore the Earth.

Credits

I used the syntax highlighting tool at https://andrewsun.com/tools/syntax-highlighter/ to typeset the Mathematica code.