Geo Computation in the Wolfram Language

José Martín-García, Wolfram Research

Join the Conversation #WolframTechConf

0. Summary

• 1. Coordinates on the maps. The role of projections.

• 2. Coordinates on the Earth. The role of datums.

• 3. Ellipsoidal geometry. GeoDistance, GeoDirection, GeoArea, etc. Compare to TravelDistance.

• 4. Low-level computational improvements in GeoGraphics v11.

• 5. Geo physics.

• 6. Future directions.

1.1. Geo Projections

• Earth’s surface is curved. A flat map will always introduce distortions. The mapping from {lat, lon} on the Earth to {x, y} on the map is called the “projection”. Every flat map uses a geo projection.

• “Equirectangular” distorts everything, mostly angles:

GeoPosition[(Entity["Barrow", Entity["City",{"Barrow","Alaska","UnitedStates"}], "Entity[\"City\", {\"Barrow\", \"Alaska\", \"UnitedStates\"}]", "city"])]

{GeoGraphics[(Entity["Barrow", Entity["City",{"Barrow","Alaska","UnitedStates"}], "Entity[\"City\", {\"Barrow\", \"Alaska\", \"UnitedStates\"}]", "city"]),GeoRange->Quantity[1,"Kilometers"]],GeoGraphics[(Entity["Barrow", Entity["City",{"Barrow","Alaska","UnitedStates"}], "Entity[\"City\", {\"Barrow\", \"Alaska\", \"UnitedStates\"}]", "city"]),GeoRange->Quantity[1,"Kilometers"],GeoProjection->"Equirectangular"]}

• “Mercator” preserves angles, but distorts angles:

GeoArea[(Entity["Africa", Entity["GeographicRegion","Africa"], "Entity[\"GeographicRegion\", \"Africa\"]", "geographic region"])]/GeoArea[(Entity["Greenland", Entity["Country","Greenland"], "Entity[\"Country\", \"Greenland\"]", "country"])]

With[{prims={Red,Polygon[(Entity["Africa", Entity["GeographicRegion","Africa"], "Entity[\"GeographicRegion\", \"Africa\"]", "geographic region"])],Green,Polygon[(Entity["Greenland", Entity["Country","Greenland"], "Entity[\"Country\", \"Greenland\"]", "country"])]}},{GeoGraphics[prims,GeoRange->"World",GeoGridLines->Automatic,Frame->True],
GeoGraphics[prims,GeoRange->"World",GeoProjection->"Mercator",GeoGridLines->Automatic,Frame->True]}]

1.2. Large variety of projections in WL

• They can be classified in multiple forms:

GeoProjectionData["Classes"]

a) Properties preserved:

GeoProjectionData["EqualArea"]

GeoGraphics[{Red,Polygon[(Entity["Africa", Entity["GeographicRegion","Africa"], "Entity[\"GeographicRegion\", \"Africa\"]", "geographic region"])],Green,Polygon[(Entity["Greenland", Entity["Country","Greenland"], "Entity[\"Country\", \"Greenland\"]", "country"])]},GeoRange->"World",GeoProjection->"Mollweide",GeoGridLines->Automatic]

GeoProjectionData["Conformal"]

GeoProjectionData["Equidistant"]

b) Method of projection:

GeoProjectionData["Cylindrical"]

GeoProjectionData["Azimuthal"]

GeoProjectionData["Conic"]

c) Ellipsoidal or spherical reference model

GeoProjectionData["Ellipsoidal"]

GeoProjectionData["Spherical"]

• Projections have parameters. GeoProjectionData gives their defaults:

GeoProjectionData["Albers","Classes"]

GeoProjectionData["Albers"]

Graphics3D[{Opacity[0.5],Cone[{{0,0,-1},{0,0,3}}],Blue,Sphere[]},Boxed->False]

“Centering” and “GridOrigin”:

GeoGraphics[GeoRange->"World",GeoProjection->"Albers",GeoGridLines->Automatic,Frame->True]

GeoGraphics[GeoMarker[{50,100}],GeoRange->"World",GeoProjection->{"Albers","Centering"->{50,100},"GridOrigin"->{1,1}},GeoGridLines->Automatic,Frame->True]

Modify the “StandardParallels” parameter:

GeoGraphics[GeoRange->"World",GeoProjection->{"Albers","StandardParallels"->{-10,20}},GeoGridLines->Automatic,Frame->True]

1.3. A collage

graticule[proj_]:=Framed[Labeled[
GeoGraphics["World",GeoProjection->{proj,"Centering"->{0,0,1.5}},GeoGridLines->(QuantityPostfix[10, "°", degrees, "AngularDegrees"]),GeoGridLinesStyle->{Red,Blue},GeoBackground->None,PlotRangePadding->Scaled[0.1]],
proj]]

graticule/@{"Armadillo","AmericanPolyconic","Bonne"}

plots=graticule/@GeoProjectionData["Spherical"];
images=Rasterize[#,ImageSize->200]&/@plots;
ImageCollage[images]

1.4. Future

• Oblique projections.

• Interrupted projections.

• Geo distortion functions.

(Legended[, (BarLegend)])

2.1. Assigning Coordinates on the Earth

• Take a relief model of the Earth (this a lot more distorted than the Earth really is!). How do we assign a pair {lat, lon} to each point on the surface? There is no natural way of drawing lines of latitude and longitude. It is not even clear where the equator should be.

a=100000;
b=90000;

{pos,elev}=GeoElevationData["World",Automatic,"GeoPositionElevation",GeoZoomLevel->2];
elevation=Interpolation[Flatten[Transpose[Append[Transpose[pos[[1]],{2,3,1}],elev+a],{3,1,2}],1]];
map=Image[GeoGraphics["World",GeoBackground->"ReliefMap",GeoZoomLevel->2]];
earthGraphics3D=Quiet@SphericalPlot3D[elevation[180/Pi (Pi/2-lat),180/Pi lon](b+(a-b)Sin[lat])/a,{lat,0,Pi},{lon,-Pi,Pi},PlotPoints->100,Mesh->None,PlotStyle->Texture[map],TextureCoordinateScaling->False,TextureCoordinateFunction->({(#5+Pi)/(2Pi),-#4/Pi}&),Ticks->False]

• Choose an oblate ellipsoid of revolution that approximates the Earth. The poles are uniquely defined, and the equator plane, but not the prime meridian.

ellipsoid={
Green,Tube[{{{0,0,-1.1a},{0,0,1.1a}}}],
Red,Opacity[0.5],Ellipsoid[{0,0,0},{a,a,b}]
};

ellipsoidGraphics3D=Graphics3D[Translate[Rotate[ellipsoid,Pi/100,{1,2,3}],{10,-20,30}]]

Show[earthGraphics3D,ellipsoidGraphics3D,PlotRange->All]

• We need to choose a prime meridian, and an orientation North/South.

primemeridianGraphics3D=Graphics3D[{Green,Arrow@Tube[Table[1.05{a Cos[θ],0,b Sin[θ]},{θ,-Pi/2,Pi/2,Pi/100}]]}];

Show[earthGraphics3D,ellipsoidGraphics3D,primemeridianGraphics3D,PlotRange->All]

• If a point of the ellipsoid has coordinates {lat, lon}, then a point of the surface of the Earth will have “geodetic coordinates” {lat, lon, h}, if it is on the local vertical (perpendicular to the ellipsoid), at a distance h on that vertical from the ellipsoid. Note that latitude is “geodetic”, and not “geocentric”:

Clear[a,b];
Manipulate[With[{ellip={1,b}},With[{ellipparam=With[{f=1-ellip[[2]]/ellip[[1]]},Join[ellip,{f,2f-(f^2)}]]},With[{xz=getXZ[lat,ellipparam],tlat=Tan[lat Degree]},
Graphics[{
Circle[{0,0},ellip],
Blue,Line[{{0,0},xz}],Circle[{0,0},.2,{0,geocentricLatitude[lat,ellipparam]Degree}],Text["Geocentric latitude"->geocentricLatitude[lat,ellipparam]HoldForm@Degree,{0,0},{1,-1.5}],
Red,Line[{{xz[[1]],xz[[2]]}+h{1,tlat},{0,xz[[2]]-Abs[xz[[1]]] Tan[lat Degree]}}],
Line[{xz-{tlat,-1}/4,xz+{tlat,-1}/4}],
PointSize[Large],Point[{xz[[1]],xz[[2]]}+h{1,tlat}],
Circle[{xz[[1]]-xz[[2]]/tlat,0},.22,{0,lat Degree}],Text["Latitude"->lat HoldForm@Degree,{xz[[1]]-xz[[2]]/tlat,0},{-1,Sign[xz[[2]]]}]
},Axes->True,PlotRange->{{-1.5,1.5},{-1,1}},ImageSize->600]]]]
,
{{b,1,"b"},0.5,1},
{{lat,30,"latitude"},-90,90},
{{h,0,"height"},-1,1},
SaveDefinitions->True,Initialization->{geocentricLatitude[lat_,ellipparam_]:=N@ArcTan[(1-ellipparam[[4]])Tan[lat  Degree]/.ComplexInfinity:>Sign[lat]Infinity]/Degree;
reducedLatitude[lat_,ellipparam_]:=N@ArcTan[(1-ellipparam[[3]])Tan[lat Degree]/.ComplexInfinity:>Sign[lat]Infinity]/Degree;
getXZ[lat_,ellipparam_]:=With[{reduced=reducedLatitude[lat,ellipparam]Degree},ellipparam[[1;;2]]{Cos[reduced],Sin[reduced]}];}]

GeodesyData[{1,0.5},{"GeocentricLatitude",Quantity[47.2,"AngularDegrees"]}]

GeodesyData[{1,0.5},{"FromGeocentricLatitude",%}]

2.2. Datums & Datum Transformations

• A “datum” is this combination:
	- An ellipsoid (values of a, b, or equivalent).
	- Location for the center of the ellipsoid.
	- Orientation of the polar axis of the ellipsoid.
	- Choice of the prime meridian.
   That is an ellipsoid plus a 3D Cartesian XYZ reference frame.

• This is what GeoPosition stores: geodetic coordinates:

GeoPosition[{lat,lon}]
GeoPosition[{lat,lon,h}],
GeoPosition[{lat,lon,h,t}]
GeoPosition[coords,datum]

Once we have a datum, latitude, longitude and height are well defined. The default datum is “ITRF00”:

GeodesyData["ITRF00","ReferenceEllipsoid"]

GeodesyData["GRS80","EllipsoidParameters"]

GeodesyData["GRS80","Semiaxes"]

There are other datums:

GeodesyData["Datum"]

And other ellipsoids:

GeodesyData["ReferenceEllipsoid"]

• The orientation of the polar axis may be given with respect to the stars, but it is not easy, and the polar axis moves. In practical terms it is even more difficult to specify the location of the center of the Earth. Usually a datum is given as a transformation from another datum. Nowadays the fundamental datums are those computed using satellites, like WGS84.

p=GeoPosition[{40.11,-88.24},"WGS72"]

The transformation is clearer in terms of XYZ vectors:

xyz=GeoPositionXYZ[p]

{translation,rot,scale}={"Translation","Rotation" /1000./3600. Pi/180,"Scale"/ 10.^8}/.GeodesyData[{"WGS72","NAD27"}]

matrix={{1,rot[[3]],-rot[[2]]},{-rot[[3]],1,rot[[1]]},{rot[[2]],-rot[[1]],1}};

This is called a (7-parameter) Helmert transformation:

translation+(1+scale) matrix.First[xyz]

newxyz=GeoPositionXYZ[xyz,"NAD27"]

%%===First[newxyz]

Now we can go back to lat, lon coordinates:

GeoPosition[newxyz]

There are also 15-parameter Helmert transformations, which take into account time derivatives:

GeodesyData[{"ITRF90","ITRF00"}]

2.3. GeoPosition arrays & QuantityArray

• GeoPosition[{lat, lon}] for a single point.

GeoPosition[(Entity["Paris", Entity["City",{"Paris","IleDeFrance","France"}], "Entity[\"City\", {\"Paris\", \"IleDeFrance\", \"France\"}]", "city"])]

• For collections of points, use GeoPosition[array], and not a list of GeoPosition objects.

geopos=GeoPosition[RandomReal[{-90,90},{1000,2}]]

Short[GeoGridPosition[geopos,"Mercator"]]//AbsoluteTiming

• Use Thread to convert to a list of GeoPosition objects.

Short[listpos=Thread[geopos]]

Use GeoPosition again to group a list of GeoPosition objects:

GeoPosition[listpos]

• Computations are then slower:

Short[GeoGridPosition[#,"Mercator"]&/@listpos]//AbsoluteTiming

• The Geo functions will try to preserve the numeric-array character, using QuantityArray if necessary:

LatitudeLongitude[geopos]//AbsoluteTiming

LatitudeLongitude/@listpos;//AbsoluteTiming

Or with functions that need to access the geo server:

GeoElevationData[geopos]//AbsoluteTiming

2.4. An Example

• A recent example in MSE.

pol=Polygon[{{115050.,485415.},{115614.,485596.},{115601.,485634.},{115887.,485725.},{116049.,485777.},{116052.,485768.},{116124.,485793.},{116136.,485798.},{116154.,485809.},{116162.,485816.},{116171.,485825.},{116262.,485750.},{116263.,485687.},{116263.,485666.},{116236.,485665.},{116206.,485662.},{116148.,485649.},{116122.,485642.},{115834.,485549.},{115846.,485509.},{115684.,485456.},{115676.,485481.},{115578.,485449.},{115633.,485279.},{115564.,485257.},{115494.,485229.},{115211.,485138.},{115146.,485116.},{115061.,485378.},{115050.,485415.}}];

proj={"Stereographic",
"ReferenceModel"->"Bessel1841",
"GridOrigin"->{155000,463000},
"Centering"->{52.1562,5.38764},
"CentralScaleFactor"->0.999908
};

GeoGraphics[{
EdgeForm[Red],Red,Polygon@GeoPosition[GeoGridPosition[pol[[1]],proj]]
},GeoProjection->proj,Frame->True]

• We need to use the fact that the {lat, lon} coordinates are not in our default “ITRF00” datum, but in the “Amersfoort” datum:

datum={"Helmert",
"ITRF00",
"Translation"->{-565.4171,-50.3319,-565.5524},
"Rotation"->{-398.957,343.988,-1877.40},
"Scale"->407.25
};

GeoGraphics[{
EdgeForm[Red],Red,Polygon@GeoPosition[GeoGridPosition[pol[[1]],proj]],
EdgeForm[Blue],Blue,Polygon@GeoPosition[GeoGridPosition[pol[[1]],proj,datum],"ITRF00"]
},GeoProjection->proj,Frame->True]

3.1. Ellipsoidal Geometry

• The fundamental problem is how to find geodesics on the surface of the ellipsoid.

• We have implemented the algorithms suggested by Charles F. F. Karney, 2011. This is Fig 1 from his article:

• By axisymmetry, there is an integration constant. This is Clairaut theorem 1735 (“conservation of angular momentum” in physics):

					Sin[(Subscript[α, 0])]==Cos[(Subscript[β, 1])] Sin[(Subscript[α, 1])]==Cos[(Subscript[β, 2])] Sin[(Subscript[α, 2])]

   with β being “reduced latitude”: Tan[β]==(Sqrt[1-(e^2)])Tan[ϕ].

• It is possible to map locally the ellipsoid onto an auxiliary sphere. Geodesics of the ellipsoid are mapped to geodesics (great circles) of the sphere. Latitude on the sphere is β, but the relation between longitudes is more complicated and involves elliptic integrals. Use series expansions, using the fact that the eccentricity of the Earth is small for fast evaluation. Use Clenshaw summation.

• There are two basic problems:

	- Direct problem: Given a point A, an azimuth (Subscript[α, 1]) and a distance (Subscript[s, 12]), find the point B. This is a matter of evaluating elliptic integrals.
	
	- Inverse problem: Given points A and B, find the azimuths (Subscript[α, 1]) and (Subscript[α, 2]), as well as the distance (Subscript[s, 12]). This is a root-finding problem. More complicated. Some cases are particularly difficult, like near-antipodal A, B, or when the geodesic passes near a pole.

3.2. GeoDistance, GeoDirection, GeoArea, ...

• Solve the direct geodetic problem:

p=Here

disp=GeoDisplacement[{Quantity[4000,"Miles"],Quantity[35,"AngularDegrees"]}]

q=GeoDestination[p,disp]

GeoGraphics[{Arrow@GeoPath[{p,disp}],GeoMarker[{p,q}]},GeoGridLines->Automatic]

• Solve the inverse geodetic problem:

GeoDisplacement[p,q]

{GeoDistance[p,q],GeoDirection[p,q]}//InputForm

• GeoLength can compute directly the length of a GeoPath object:

GeoLength[GeoPath[{p,q}]]

• The area between a geodesic and the equator can also be computed by evaluation of elliptical integrals.

outer=GeoPath[{{-50,-50},{50,-50},{50,50},{-50,50},{-50,-50}}];
inner=GeoCircle[{0,0},Quantity[2000,"Miles"]];

fc=FilledCurve[{{outer},{inner}}];

GeoGraphics[fc]

GeoArea[fc]

{GeoArea[FilledCurve[outer]],GeoArea[FilledCurve[inner]]}

Subtract@@%

3.3. Compare with Travel Directions

• Compute travel directions between two distant cities:

td=TravelDirections[{Entity["City",{"Lisbon","Lisboa","Portugal"}],Entity["City",{"Beijing","Beijing","China"}]}]

The geodesic (shortest path) goes much more North than the path given by TravelDirections:

GeoGraphics[{Thick,Red,GeoPath[td],Blue,GeoPath[{(Entity["Lisbon", Entity["City",{"Lisbon","Lisboa","Portugal"}], "Entity[\"City\", {\"Lisbon\", \"Lisboa\", \"Portugal\"}]", "city"]),(Entity["Beijing", Entity["City",{"Beijing","Beijing","China"}], "Entity[\"City\", {\"Beijing\", \"Beijing\", \"China\"}]", "city"])}]},GeoProjection->"Mercator",GeoGridLines->Automatic]

Again, the Mercator projection was misleading. An azimuthal projection shows much better how the paths actually compare:

Show[%,GeoProjection->"LambertAzimuthal",GeoZoomLevel->4]

• Another example of TravelDirections computation:

museums={(Entity["Museum of Modern Art", Entity["Museum","MuseumOfModernArt::nfkn8"], "Entity[\"Museum\", \"MuseumOfModernArt::nfkn8\"]", "museum"]), (Entity["American Museum of Natural History", Entity["Museum","AmericanMuseumOfNaturalHistory::gf4kv"], "Entity[\"Museum\", \"AmericanMuseumOfNaturalHistory::gf4kv\"]", "museum"])};

drive=TravelDirections[museums,TravelMethod->"Driving"]

walk=TravelDirections[museums,TravelMethod->"Walking"]

Compare both paths:

GeoGraphics[{GeoMarker[museums],Thick,Red,GeoPath[drive],Blue,GeoPath[walk]},GeoRange->Quantity[1,"Mile"]]

{walk["TravelDistance"],walk["TravelTime"]}

{drive["TravelDistance"],drive["TravelTime"]}

4.1. Version 11 Improvement: Clipping

• Longitude periodicity is difficult to handle with non-cylindrical projections

Old result:

Now:

GeoGraphics[{AbsoluteThickness[3],
{Red,GeoPath["Equator"]},
{Orange,GeoPath["NorthernTropic"],GeoPath["SouthernTropic"]},
{Blue,GeoPath["ArcticCircle"],GeoPath["AntarcticCircle"]},
{Green,GeoPath["GreenwichMeridian"],GeoPath["DateLineMeridian"]},
{Purple,GeoPath["DateLine"]}
},GeoProjection->"WinkelTripel"]

4.2. Version 11 Improvement: Geo Grid Lines

Old result:

Now:

GeoGraphics[GeoRange->"World",GeoProjection->{"Orthographic","Centering"->{40,140}},GeoGridLines->30,GeoGridLinesStyle->{Black,None},GeoBackground->None,Frame->True]

Random example:

proj=RandomChoice[GeoProjectionData["Spherical"]]
centering={RandomReal[{-90,90}],RandomReal[{-180,180}],RandomReal[{1,2}]}
n=RandomInteger[{1,100}]
GeoGraphics[{Thick,Red,GeoPath["Equator"],Blue,GeoPath["PrimeMeridian"]},GeoRange->"World",GeoProjection->{proj,"Centering"->centering},GeoGridLines->n,GeoGridLinesStyle->{Red,Blue},GeoBackground->None,Frame->True,PlotRangePadding->0.1]

5.1 GeogravityModelData

• The gravitational potential W(x, y, z), which includes the centrifugal term, is described as a sum

	W(x, y, z) == U(x, y, z) + T(x, y, z)

where U is the “normal potential” and T is the “disturbing potential”. U is an exact formula, even in the ellipsoidal case. It is possible to make the reference ellipsoid an equipotential surface of the normal potential with value (Subscript[W, 0]). The irregular surface that is equipotential for the full potential W is called the “geoid”. Usually W is expressed as a sum of many spherical harmonics, and both U and T have all terms. We use the EGM96 model of the Earth, that provides those harmonic coefficients. The gravitational vector field is the gradient of W.

• Together with the semiaxes of an ellipsoid, we need two physical parameters to know exactly the form of the normal potential. For example the “GRS80” ellipsoid is complemented with this product:

GM=(Quantity[3.986005*^14, (("m")^3)"/"(("s")^2), meters cubed per second squared, ((("Meters")^3)/(("Seconds")^2))]);

and this angular velocity, which is just Quantity[2 Pi, 1 / “SiderealDay”]:

ω=(Quantity[0.00007292115, "per ""second", reciprocal seconds, (1/("Seconds"))]);

• Example at a point:

GeogravityModelData[(Entity["New York City", Entity["City",{"NewYork","NewYork","UnitedStates"}], "Entity[\"City\", {\"NewYork\", \"NewYork\", \"UnitedStates\"}]", "city"])]

• Values for a region are computed as a QuantityArray object.

GeogravityModelData[{{40, -100},{45, -90}},"Potential",GeoZoomLevel->3]

5.2 GeomagneticModelData

• The magnetic field is calculated for the current time using the World Magnetic Model or "WMM". Other models are available.

GeomagneticModelData[GeoPosition[{45, -120, 100}]]

The model can give the positions of the magnetic poles:

GeoPosition["NorthGeomagneticPole"]

GeoPosition["NorthModelDipPole"]

dates=DateRange[(DateObject["1656", DateObject[{1656}]]),(DateObject["2016", DateObject[{2016}]]),(Quantity[10, "yr", years, "Years"])];

geomagneticNorth=GeomagneticModelData["NorthGeomagneticPole",dates];

dipNorth=GeomagneticModelData["NorthModelDipPole",dates];

GeoGraphics[{
Red,Point[geomagneticNorth],Blue,Point[dipNorth],
PointSize[Large],
Red,Point[Last@geomagneticNorth],Blue,Point[Last@dipNorth]
},GeoRange->{{65,90},All},GeoGridLines->Automatic]

5.3 GeoElevationData

• Three-argument elevation data:

GeoElevationData[loc,elevationtype,outputtype]

• The geoid surface has a geodetic elevation (i.e. wrt the reference ellipsoid) called undulation. We have EGM96 data for undulation. Elevation with respect to the reference ellipsoid is called “geodetic”, and with respect to the geoid is called “orthometric”.

• Possible elevation types are “Geodetic”, “Orthometric”, “Undulation”. Default is “Orthometric”.

• Possible output types are “QuantityArray”, “GeoPosition”, “GeoPositionXYZ”, and in v11 also “Region”. Default is Automatic.

GeoElevationData[Here,UnitSystem->"Metric"]

GeoElevationData[Here,"Geodetic",UnitSystem->"Metric"]

GeoElevationData[Here,Automatic,"GeoPosition"]

GeoElevationData[Here,"Undulation",UnitSystem->"Metric"]

• Francisco will talk about new elevation data coming soon...

6. Future

• New projection types.

• EPSG data.

• Geo distortion functions.

• Add algorithms to avoid having GeoWithinQ and other functions performing geo server calls if possible.

• Incorporate elevation and vertical datums into our computations.

• Continue improving the low-level functionality used by GeoGraphics and the rest of the Geo functionality.