Geo Computation

José Martín-García
Algorithms R&D
Wolfram Research

Join the Conversation #WolframTechConf

Geo Computation

• There are many computational aspects in the WL Geo project. Some fundamental aspects are:

- Geo measures and relations (distances, angles, areas, membership, ...)

- Cartographic geo projections and datums (mapping the Earth onto a flat map)

- Internal computations in GeoGraphics (estimates of ranges, zoom, projection, image size, ...)

- Geophysical modeling.

• Functionality has improved steadily since 10.0, both in coverage and speed. Main recent driver: GeoImage.

• Plan of the talk:

- Summary of existing functionality.

- Recent additions.

- Answers to frequent questions.

Measuring the Earth

Geo locations

• GeoPosition. Flexible specification of geo locations:

Here

GeoPosition[{lat, lon}]				Angles always in degrees
	GeoPosition[{lat, lon, height}]		Elevation in meters, default is zero
	GeoPosition[{lat, lon, height, date}]	Default date is now
	GeoPosition[{pos1, pos2, ...}]		Arrays of any structure, but consistent dimension
	GeoPosition[pos, datum]			Different datums for the Earth, or different celestial globes

• GeoPosition can also parse / extract / lookup information (convenient for pre-processing input):

GeoPosition[{{40.2,-91.51},{"30d20'51.4'' N","102deg, 31min, 45.02sec W"},{(Quantity[1, "rad", radians, "Radians"]),(QuantityPostfix[150, "°", degrees, "AngularDegrees"])},(Entity["London", Entity["City",{"London","GreaterLondon","UnitedKingdom"}], "Entity[\"City\", {\"London\", \"GreaterLondon\", \"UnitedKingdom\"}]", "city"])}]//InputForm

GeoPosition[]

GeoPosition[(Entity["International Space Station", Entity["Satellite","25544"], "Entity[\"Satellite\", \"25544\"]", "satellite"])]

GeoPosition[(Entity["Olympus Mons", Entity["SolarSystemFeature","OlympusMonsMars"], "Entity[\"SolarSystemFeature\", \"OlympusMonsMars\"]", "solar system feature"])]

• Other geo location specifications:

	GeoPositionXYZ[{x, y, z}]		Cartesian coordinates wrt the Earth-centered Greenwich/polar frame
	GeoPositionENU[{e, n, u}, loc]	Cartesian coordinates wrt a frame centered at loc
	GeoGridPosition[coords, proj]	Coordinates in any geo projection

Geodetic computation

• Geodesics: GeoDistance, GeoDistanceList, GeoDirection, GeoDisplacement & GeoDestination. (Local computations.)

GeoDistance[Here,GeoPosition[{0,0}]]

GeoDirection[Here,GeoPosition[{0,0}]]

GeoGraphics[{Arrow@GeoPath[{Here,%%,%}],GeoMarker[{0,0}]},GeoGridLines->Automatic,GeoProjection->"Mercator"]

• Other measures: GeoLength, GeoArea, GeoBounds. (Local computations.)

• Travel paths: TravelDirections, TravelDistance, TraveDistanceList, TravelTime. (Computations performed in the geo server.)

• Distance relations: GeoNearest, GeoEntities, GeoIdentify, GeoWithinQ. (Computations performed in the geo server.)

• Recent extension of Nearest. (Local computations.)

capitals=GeoPosition[(EntityClass["all US states with District of Columbia", EntityClass["AdministrativeDivision","AllUSStatesPlusDC"], "EntityClass[\"AdministrativeDivision\", \"AllUSStatesPlusDC\"]", "administrative divisions"])[(EntityProperty["capital city", EntityProperty["AdministrativeDivision","CapitalCity"], "EntityProperty[\"AdministrativeDivision\", \"CapitalCity\"]"])]]

Nearest[capitals,Here]

GeoIdentify["City",%]

Nearest is faster than GeoNearest, but has more reduced geo applicability:

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

Nearest[points,Here]//AbsoluteTiming

GeoNearest[Thread@points,Here]//AbsoluteTiming

GeoGraphics/GeoImage & geo projections

• GeoImage supports much of the power of GeoGraphics. Five important differences:
	- The output type is Image, and not (Geo)Graphics. Output size is now controlled by the image dimensions.
	- The default geo styling is “Satellite”, and not “StreetMap”.
	- The default geo styling is “Mercator” at all scales.
	- Automatic switch to the “DigitalGlobe” provider for zoom levels beyond 11.
	- The GeoBackground option is converted into a second argument.

GeoImage["World","SatelliteWithLabels",GeoServer->"DigitalGlobe",GeoZoomLevel->2]

Head[%]

• GeoImage has motivated a restructuring of part of the GeoGraphics code to avoid the dependence on Texture. It can handle the same geo input (ignoring the primitives). It can handle projections as well:

GeoImage["World",GeoProjection->"Mollweide"]

Recent: GeoAntipode

• Compute the antipodal position or antipodal geo primitive.

• GeoAntipode on explicit locations:

Here

GeoAntipode[Here]

• Like GeoPosition does, geo entities will be collapsed to a GeoPosition object:

GeoAntipode[(Entity["Madrid", Entity["City",{"Madrid","Madrid","Spain"}], "Entity[\"City\", {\"Madrid\", \"Madrid\", \"Spain\"}]", "city"])]

That location is in New Zealand:

GeoIdentify["Country",%]

• On geo primitives:

GeoAntipode[DayHemisphere[Now]]

GeoAntipode[GeoDisk[Here]]

LinguisticAssistant["Polygon"]

GeoAntipode[%]

GeoGraphics[%]

• Draw the antipode of Antarctica around the North Pole:

GeoGraphics[GeoAntipode[Polygon[(Entity["Antarctica", Entity["GeographicRegion","Antarctica"], "Entity[\"GeographicRegion\", \"Antarctica\"]", "geographic region"])]],GeoGridLines->Automatic]

• We can do the same thing for the whole world:

pol=(Entity["world", Entity["GeographicRegion","World"], "Entity[\"GeographicRegion\", \"World\"]", "geographic region"])["Polygon"];

GeoGraphics[{pol,Red,GeoAntipode[pol]},GeoRange->"World",GeoBackground->None]

Soon: GeoRotate

• Rotate a geo primitive in multiple ways.

uspol=(Entity["United States", Entity["Country","UnitedStates"], "Entity[\"Country\", \"UnitedStates\"]", "country"])["Polygon"]

GIS`GeoRotate[uspol,Here->(Entity["Paris", Entity["City",{"Paris","IleDeFrance","France"}], "Entity[\"City\", {\"Paris\", \"IleDeFrance\", \"France\"}]", "city"])]

GeoGraphics[{%,Arrow@GeoPath[{Here,(Entity["Paris", Entity["City",{"Paris","IleDeFrance","France"}], "Entity[\"City\", {\"Paris\", \"IleDeFrance\", \"France\"}]", "city"])}]},GeoRange->"World"]

• Interact with the rotated object, in any projection:

Manipulate[DynamicModule[{pt=First[GeoGridPosition[Here,proj]]},LocatorPane[Dynamic[pt],Dynamic@GeoGraphics[Point@@GIS`GeoRotate[uspol,Here->GeoGridPosition[pt,proj]],GeoRange->"World",GeoProjection->proj]]],{{proj,"Equirectangular"},GeoProjectionData["Spherical"]},SaveDefinitions->True]

• A 3D rotation has three degrees of freedom. One of them is being computed automatically.

graticule[pos_]:={
Table[GeoCircle[pos,r],{r,(Quantity[0, "Mm", megameters, "Megameters"]),(Quantity[20, "Mm", megameters, "Megameters"]),(Quantity[2, "Mm", megameters, "Megameters"])}],Table[GeoPath[{pos,(Quantity[20, "Mm", megameters, "Megameters"]),α}],{α,(QuantityPostfix[0, "°", degrees, "AngularDegrees"]),(QuantityPostfix[360, "°", degrees, "AngularDegrees"]),(QuantityPostfix[30, "°", degrees, "AngularDegrees"])}]
}

Manipulate[GeoGraphics[{Gray,graticule[{lat,lon}],GeoMarker[{lat,lon}],Black,Point@@GIS`GeoRotate[uspol,{{lat,lon},-α}]},
GeoRange->"World",
GeoProjection->"Mercator"
],
{{lat,90},-90,90},
{{lon,0},-180,180},
{{α,0},-90,270},
SaveDefinitions->True
]

Q1: How to get a map with ticks in miles?

• The key parameter is “ReferenceModel”, by default typically 1, or Pi/180 for cylindrical maps with longitude ticks in degrees.

GeoProjectionData["Orthographic"]

GeoProjectionData["Mercator"]

Tally[GeoProjectionData[#,"ReferenceModel"]&/@GeoProjectionData["Spherical"]]

• To get a map in miles, use “ReferenceModel” -> 3958.761

QuantityMagnitude[(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"])["Radius"],"Miles"]

GeoGraphics[GeoRange->"World",GeoProjection->{"Orthographic","ReferenceModel"->3958.761},Frame->True,GeoGridLines->Automatic,GeoBackground->"Satellite"]

• No projection represents all distances faithfully. The best we can get is having regions (frequently around some preferred line) with good approximation. Some projections have standard parallels as parameters, which are lines along which distances are properly represented.

GeoProjectionData["Albers"]

Manipulate[
GeoGraphics[{GeoPath[{"Parallel",lat1}],GeoPath[{"Parallel",lat2}]},GeoRange->"World",GeoProjection->{"Albers","ReferenceModel"->3958.761,"StandardParallels"->{lat1,lat2}},Frame->True,GeoGridLines->Automatic],
{{lat1,0},-90,90},
{{lat2,40},-90,90}
]

Q2: How to get a square map?

• Take a 6 x 6 degree geo range. The map will be only truly square in the “Equirectangular” projection:

GeoGraphics[GeoRange->{{40,46},{10,16}},GeoProjection->"Equirectangular",GeoGridLines->(QuantityPostfix[1, "°", degrees, "AngularDegrees"]),Frame->True]

• The map is not square in most projections, including the Mercator projection:

GeoGraphics[GeoRange->{{40,46},{10,16}},GeoProjection->"Mercator",GeoGridLines->(QuantityPostfix[1, "°", degrees, "AngularDegrees"]),Frame->True]

• Extract the projection, and use it together with a specification of a square region in the projected map:

projOption=Options[%,GeoProjection]

GeoGraphics[PlotRange->{{10,16},{44,50}},projOption]

• For any other desired aspect ratio, use a plot range that corresponds to that aspect ratio:

GeoGraphics[PlotRange->{{7,19},{44,50}},projOption]

• Same thing can be done with any other projection:

GeoGraphics[GeoRange->{{40,46},{10,16}},GeoProjection->"WinkelTripel",GeoGridLines->(QuantityPostfix[1, "°", degrees, "AngularDegrees"]),Frame->True]

projOption=Options[%,GeoProjection]

GeoGraphics[PlotRange->{{-0.04,0.04},{0.70,0.78}},projOption]

Q3: Difference between Line and GeoPath?

• Take three different “lines” between Oslo and Chicago:

prims={
Red,GeoPath[{(Entity["Oslo", Entity["City",{"Oslo","Oslo","Norway"}], "Entity[\"City\", {\"Oslo\", \"Oslo\", \"Norway\"}]", "city"]),(Entity["Chicago", Entity["City",{"Chicago","Illinois","UnitedStates"}], "Entity[\"City\", {\"Chicago\", \"Illinois\", \"UnitedStates\"}]", "city"])},"Geodesic"],
Green,GeoPath[{(Entity["Oslo", Entity["City",{"Oslo","Oslo","Norway"}], "Entity[\"City\", {\"Oslo\", \"Oslo\", \"Norway\"}]", "city"]),(Entity["Chicago", Entity["City",{"Chicago","Illinois","UnitedStates"}], "Entity[\"City\", {\"Chicago\", \"Illinois\", \"UnitedStates\"}]", "city"])},"RhumbLine"],
Blue,Line[{(Entity["Oslo", Entity["City",{"Oslo","Oslo","Norway"}], "Entity[\"City\", {\"Oslo\", \"Oslo\", \"Norway\"}]", "city"]),(Entity["Chicago", Entity["City",{"Chicago","Illinois","UnitedStates"}], "Entity[\"City\", {\"Chicago\", \"Illinois\", \"UnitedStates\"}]", "city"])}]
};

• Azimuthal projections make central geodesics appear almost straight. Note both the Line and the geodesic go through Iceland:

GeoGraphics[{prims},GeoGridLines->Automatic]

• Line always produces a straight result, but different for each projection. Note that Line does not go through Iceland now. GeoPath produces consistent results on the surface of the Earth.

GeoGraphics[{prims},GeoProjection->#,GeoGridLines->Automatic]&/@{"Mercator","CraigRetroazimuthal","Wiechel"}

• Same difference with the interior regions. Polygons do not even intersect in the equirectangular projection:

points=GeoPosition[{{85,-90},{60,-35},{0,0},{85,-90}}];

GeoGraphics[{
PointSize[Large],Point[points],
Red,FilledCurve[GeoPath[points]],
Blue,FilledCurve[Line[points]]
}]

However, they fully intersect in the orthographic projection, with appropriate centering:

GeoGraphics[{
PointSize[Large],Point[points],
Red,FilledCurve[GeoPath[points]],
Blue,FilledCurve[Line[points]]
},
GeoProjection->{"Orthographic","Centering"->{50,-10}},
GeoGridLines->Automatic]

• However: GeoPath is internally converted into Line with many more additional points. If points of the line are close enough, it is better to use Line (and similarly Polygon). For example, all our country polygons are given as Polygon[...] or FilledCurve[{..., Line[...], ...}] objects.

polygon=LinguisticAssistant["Polygon"]

GeoGraphics[{Red,Point@@polygon}]

Q4: GeoCenter, “Centering” and “GridOrigin”

• “GridOrigin”: constant shift of resulting map.
  GeoCenter: center of the included geo range (region of Earth’s surface).
  “Centering”: point on Earth that will be mapped to the grid origin.

Manipulate[
GeoGraphics[{GeoMarker[{0,gclon}],GeoMarker[{0,lon0},"Color"->Blue]},
GeoRange->"World",GeoCenter->{0,gclon},PlotRange->{{-4,4},{-1.5,1.5}},GeoGridLines->Automatic,GeoProjection->{"Mollweide","Centering"->{0,lon0},"GridOrigin"->{origx,0}},Frame->True],
{{gclon,0},-180,180},{{lon0,0},-180,180},
{{origx,0},-3,3}]

In most applications the “Centering” value is automatically chosen to agree with the GeoCenter value

Options[GeoGraphics[(Entity["United States", Entity["Country","UnitedStates"], "Entity[\"Country\", \"UnitedStates\"]", "country"]),GeoProjection->"Mollweide"],{GeoCenter,GeoProjection}]

• Example: How far are we from North Korea?

Manipulate[GeoGraphics[{GeoMarker[Here],Table[GeoCircle[(Entity["North Korea", Entity["Country","NorthKorea"], "Entity[\"Country\", \"NorthKorea\"]", "country"]),d],{d,(Quantity[0, "km", kilometers, "Kilometers"]),(Quantity[20000, "km", kilometers, "Kilometers"]),(Quantity[2000, "km", kilometers, "Kilometers"])}]},GeoProjection->{proj,"Centering"->(Entity["North Korea", Entity["Country","NorthKorea"], "Entity[\"Country\", \"NorthKorea\"]", "country"])},GeoRange->"World",GeoCenter->-150],{proj,{"Equirectangular","LambertAzimuthal","AzimuthalEquidistant"}}]