Cartography

Joshua R. Davis, jdavis@carleton.edu, revised 2020 October 6

Spherical coordinates

This parametrization of the unit sphere is based on spherical coordinates, with x playing the role of phi and y playing the role of theta. The following manipulator shows the image of the parametrization, with a single point highlighted. Moving the sliders changes the x and y used to make the highlighted point. (If the manipulator is too sluggish, try replacing “Quality” with “Speed”.)

Clear[x, y];
f = {Sin[x] Cos[y], Sin[x] Sin[y], Cos[x]};
xs = {0, Pi};
ys = {0, 2 Pi};
Manipulate[
Show[
ParametricPlot3D[f, {x, xs[[1]], xs[[2]]}, {y, ys[[1]], ys[[2]]}, 
PlotRange->{{-1.1, 1.1}, {-1.1, 1.1}, {-1.1, 1.1}}, PerformanceGoal->"Quality"], 
Graphics3D[Sphere[f /. {x -> x0, y -> y0}, 0.05]]],
{x0, xs[[1]], xs[[2]]}, {y0, ys[[1]], ys[[2]]}]

Jacobian, metric, area distortion factor.

dfdx = Simplify[Transpose[{D[f, x], D[f, y]}]];
MatrixForm[dfdx]
g = Simplify[Transpose[dfdx].dfdx];
MatrixForm[g]
Simplify[Sqrt[Det[g]], {xs[[1]] < x < xs[[2]], ys[[1]] < y < ys[[2]]}]

Stereographic projection

This parametrization is the inverse of stereographic projection. I’ve restricted the domain to (-5, 5) x (-5, 5), but feel free to change that.

Clear[x, y];
f = {2x / (1 + x^2 + y^2), 2 y / (1 + x^2 + y^2), (-1 + x^2 + y^2) / (1 + x^2 + y^2)};
xs = {-5, 5};
ys = {-5, 5};
Manipulate[
Show[
ParametricPlot3D[f, {x, xs[[1]], xs[[2]]}, {y, ys[[1]], ys[[2]]}, 
PlotRange->{{-1.1, 1.1}, {-1.1, 1.1}, {-1.1, 1.1}}, PerformanceGoal->"Quality"], 
Graphics3D[Sphere[f /. {x -> x0, y -> y0}, 0.05]]],
{x0, xs[[1]], xs[[2]]}, {y0, ys[[1]], ys[[2]]}]

Jacobian, metric, area distortion factor.

dfdx = Simplify[Transpose[{D[f, x], D[f, y]}]];
MatrixForm[dfdx]
g = Simplify[Transpose[dfdx].dfdx];
MatrixForm[g]
Simplify[Sqrt[Det[g]], {xs[[1]] < x < xs[[2]], ys[[1]] < y < ys[[2]]}]

Mercator projection

This parametrization is the inverse of the Mercator projection, which is one of the most popular map projections in history. The x variable is longitude and the y variable is related to latitude. The true domain of y is R, but I restrict to a finite interval for the sake of the computer.

Clear[x, y];
f = {Sin[2 ArcCot[E^y]] Cos[x], Sin[2 ArcCot[E^y]] Sin[x], Cos[2 ArcCot[E^y]]};
xs = {0, 2 Pi};
ys = {-2, 2};
Manipulate[
Show[
ParametricPlot3D[f, {x, xs[[1]], xs[[2]]}, {y, ys[[1]], ys[[2]]}, 
PlotRange->{{-1.1, 1.1}, {-1.1, 1.1}, {-1.1, 1.1}}, PerformanceGoal->"Quality"], 
Graphics3D[Sphere[f /. {x -> x0, y -> y0}, 0.05]]],
{x0, xs[[1]], xs[[2]]}, {y0, ys[[1]], ys[[2]]}]

Jacobian, metric, area distortion factor. The metric doesn’t look like a scalar matrix, but it is. Try proving that its diagonal entries are equal.

dfdx = Simplify[Transpose[{D[f, x], D[f, y]}]];
MatrixForm[dfdx]
g = Simplify[Transpose[dfdx].dfdx];
MatrixForm[g]
Simplify[Sqrt[Det[g]], {xs[[1]] < x < xs[[2]], ys[[1]] < y < ys[[2]]}]

Archimedes’s cylindrical projection

This parametrization is the inverse of a cylindrical projection developed by Archimedes about 2,200 years ago. If I understand the history, it was reinvented by Lambert about 200 years ago. The idea behind the projection is simple. Wrap the sphere in a cylinder, whose axis passes through the north and south poles. Given a point q on the sphere (other than one of the poles), project it away from the axis and onto the cylinder. Then unbend the cylinder into a plane. In the end, x plays the role of longitude, and y is related to latitude.

Clear[x, y];
f = {Sqrt[1 - y^2] Cos[x], Sqrt[1 - y^2] Sin[x], y};
xs = {0, 2 Pi};
ys = {-1, 1};
Manipulate[
Show[
ParametricPlot3D[f, {x, xs[[1]], xs[[2]]}, {y, ys[[1]], ys[[2]]}, 
PlotRange->{{-1.1, 1.1}, {-1.1, 1.1}, {-1.1, 1.1}}, PerformanceGoal->"Quality"], 
Graphics3D[Sphere[f /. {x -> x0, y -> y0}, 0.05]]],
{x0, xs[[1]], xs[[2]]}, {y0, ys[[1]], ys[[2]]}]

Jacobian, metric, area distortion factor.

dfdx = Simplify[Transpose[{D[f, x], D[f, y]}]];
MatrixForm[dfdx]
g = Simplify[Transpose[dfdx].dfdx];
MatrixForm[g]
Simplify[Sqrt[Det[g]], {xs[[1]] < x < xs[[2]], ys[[1]] < y < ys[[2]]}]

Let’s see some maps already

Mathematica has built-in features for accessing geographic databases online and then doing the math needed to make maps of the Earth. Here is a cylindrical projection that is similar to Archimedes’s (if not identical). Notice how shapes are greatly distorted in the polar regions.

GeoGraphics[
GeoCenter->GeoPosition[{0,0}], (* Edge through anti-prime meridian. *)
GeoRange->{{-89.9,89.9},{-180,180}}, (* Tuned to GeoCenter and equal-angle below. *)
GeoRangePadding->None, (* Unnecessary? *)
GeoGridLines->{Range[-85,85,5],Range[-175,175,5]}, (* Tuned to GeoCenter. *)
GeoGridLinesStyle->None, (* Lots of options. *)
GeoZoomLevel->1, (* 1 for exploration; 4 for high quality; 5 too big. *)
GeoProjection->"CylindricalEqualArea", (* Mercator or CylindricalEqualArea. *)
GeoBackground->GeoStyling["ReliefMap"] (* ReliefMap, CountryBorders, Satellite, or StreetMapNoLabels. *)
]

Here’s Mercator’s cylindrical projection. Notice how areas are greatly distorted in the polar regions.

GeoGraphics[
GeoCenter->GeoPosition[{0,0}], (* Edge through anti-prime meridian. *)
GeoRange->{{-89.9,89.9},{-180,180}}, (* Tuned to GeoCenter. *)
GeoRangePadding->None, (* Unnecessary? *)
GeoGridLines->{Range[-85,85,5],Range[-175,175,5]}, (* Tuned to GeoCenter. *)
GeoGridLinesStyle->None, (* Lots of options. *)
GeoZoomLevel->1, (* 1 for exploration; 4 for high quality; 5 too big. *)
GeoProjection->"Mercator", (* Mercator or CylindricalEqualArea. *)
GeoBackground->GeoStyling["ReliefMap"] (* ReliefMap, CountryBorders, Satellite, or StreetMapNoLabels. *)
]