Progress in Geo Computation

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

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

OUTLINE

• GeoGraphics Updates

New Geo Background Default

New Geo Range Default

Improved Handling of Labels

Support for Epilog

Multiposition Visible Region

• GeoOrientationData

Orientation of Earth in Space

Rotations

GeoGraphics Updates

1. New GeoBackground Default

Raster Backgrounds with Vector Labels

• We now use a mixture of raster street map images and vector labels:

GeoGraphics[(Entity["Italy", Entity["Country","Italy"], "Entity[\"Country\", \"Italy\"]", "country"]),GeoGridLines->Automatic]

• This is used with all geo projections, without distortion:

GeoGraphics[(Entity["Italy", Entity["Country","Italy"], "Entity[\"Country\", \"Italy\"]", "country"]),GeoProjection->"TransverseMercator",GeoGridLines->Automatic]

• In case the former behavior is needed, use:

GeoGraphics`$GeoBackgroundType="Legacy";

GeoGraphics[(Entity["Italy", Entity["Country","Italy"], "Entity[\"Country\", \"Italy\"]", "country"]),GeoProjection->"TransverseMercator",GeoGridLines->Automatic]

• Choose full vector maps:

GeoGraphics`$GeoBackgroundType="Vector";

GeoGraphics[(Entity["Italy", Entity["Country","Italy"], "Entity[\"Country\", \"Italy\"]", "country"]),GeoProjection->"TransverseMercator",GeoGridLines->Automatic]

• Restore the original value:

GeoGraphics`$GeoBackgroundType=.

Lists of Geo Backgrounds

• Combine geo backgrounds of different types:

GeoGraphics[(Entity["Austria", Entity["Country","Austria"], "Entity[\"Country\", \"Austria\"]", "country"]),ImageSize->Large,GeoBackground->{"StreetMapNoLabels","VectorLabels"}]

• Control parameters independently:

GeoGraphics[(Entity["Austria", Entity["Country","Austria"], "Entity[\"Country\", \"Austria\"]", "country"]),ImageSize->Large,GeoBackground->{
GeoStyling["StreetMapNoLabels",GeoZoomLevel->5],
GeoStyling["VectorLabels",GeoZoomLevel->8]
}]

• Opacity may be needed:

GeoGraphics[(Entity["Austria", Entity["Country","Austria"], "Entity[\"Country\", \"Austria\"]", "country"]),GeoBackground->{
"ReliefMap",
GeoStyling["StreetMapNoLabels",Opacity[0.6]],
"VectorLabels"
},ImageSize->Large]

2. New GeoRange Default

World Maps

• Since its beginning in WL10 (2014), GeoGraphics[ ] showed a local map.

• Now it shows a world map, still taking into account your geo location:

GeoGraphics[]

• For other projections too:

GeoGraphics[GeoProjection->"LambertAzimuthal"]

• If pointlike information is provided, the map is still local:

Here

GeoGraphics[%]

• Extended data, like polygons or entities with polygons, determine the geo range:

GeoGraphics[(Entity["Champaign", Entity["City",{"Champaign","Illinois","UnitedStates"}], "Entity[\"City\", {\"Champaign\", \"Illinois\", \"UnitedStates\"}]", "city"])]

Consequence for Invalid Maps

• Invalid maps with no location information will default to world maps:

GeoGraphics[GeoPath[foo]]

• GeoGraphics is now more lenient with incorrect parts of the input:

GeoGraphics[{Red,GeoPath[{"Parallel",Here}],Blue,GeoPath[foo],Green,GeoPath[{"Meridian",Here}]}]

3. Improved Handling of Labels

Labels Follow Lines

• Compute label collisions at the label of the individual letters:

GeoGraphics[(Entity["Monaco", Entity["Country","Monaco"], "Entity[\"Country\", \"Monaco\"]", "country"]),GeoRange->Quantity[3,"Miles"],GeoBackground->"VectorMarketing",ImageSize->1000]

Label Scaling

• Rescale labels together with the map:

GeoGraphics[Here,GeoBackground->"VectorDark"]

4. Support for Epilog

• Add geo primitives to a map, but keeping estimation of ranges and projection:

GeoGraphics[{Red,Polygon[(Entity["France", Entity["Country","France"], "Entity[\"Country\", \"France\"]", "country"])]},Epilog->{Blue,Polygon[(Entity["Germany", Entity["Country","Germany"], "Entity[\"Country\", \"Germany\"]", "country"])]}]

• Compare with using both sets of geo primitives together:

GeoGraphics[{{Red,Polygon[(Entity["France", Entity["Country","France"], "Entity[\"Country\", \"France\"]", "country"])]},{Blue,Polygon[(Entity["Germany", Entity["Country","Germany"], "Entity[\"Country\", \"Germany\"]", "country"])]}}]

5. Multiposition Visible Region

• We were able to show the visible region from a point:

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

• Compute points along the path of a flight:

locs=FlightData[(Entity[(RowDefault["JetBlue Airways", " flight ", "224", " (departing ", (RowWithSeparators[,, ",", "11 Oct 2021", (RowDefault["7:07 am", " ", "PDT"])]), ")"]), Entity["Flight","202110110000099785"], "Entity[\"Flight\", \"202110110000099785\"]", "flight"]),"FlightPath"]

• Show the visibility regions for the respective points of the path:

GeoGraphics[{Blue,Line[locs],Opacity[0.15,Pink],GeoVisibleRegionBoundary[locs]}]

GeoGraphics[{Blue,Line[locs],GeoStyling[Opacity[0.01,Pink]],GeoVisibleRegion[locs]}]

GeoOrientationData

Orientation of Earth in Space

Goal

• Transition from precise positioning on Earth to precise positioning in Astronomy, and viceversa.

• What is the rotation between an Earth-fixed frame and a celestial frame at a given time?

• Frames on Earth: Datums. Choose the ITRS (ITRFnn realizations).

GeodesyData["Datum"]

• Celestial frames: Choose the ICRS, determined by distant quasars.

• Information provided daily by the IERS [web].

Precision Level

• Geo computations on Earth are precise to less than a millimeter.

• We want astro computations to be precise to less than a milliarcsecond.

(Entity["Moon", Entity["PlanetaryMoon","Moon"], "Entity[\"PlanetaryMoon\", \"Moon\"]", "planetary moon"])[(EntityProperty["distance from Earth", EntityProperty["PlanetaryMoon","DistanceFromEarth"], "EntityProperty[\"PlanetaryMoon\", \"DistanceFromEarth\"]"])]Quantity[1,"Milliarcseconds"/"Radians"]//UnitConvert

(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"])[(EntityProperty["average radius", EntityProperty["Planet","Radius"], "EntityProperty[\"Planet\", \"Radius\"]"])]Quantity[1,"Milliarcseconds"/"Radians"]//UnitSimplify

It takes 10 minutes for the Earth to precess 1 mas:

Quantity[25800.,"Years"]Quantity[1,"Milliarcseconds"]/Quantity[2Pi,"Radians"]//UnitSimplify

It takes about 0.1 milliseconds for the Earth to rotate 1 mas:

Quantity[1.,"Days"]Quantity[1,"Milliarcseconds"]/Quantity[2Pi,"Radians"]//UnitSimplify

Decomposition in Three Rotations

1) Motion of Earth’s rotation axis in space

• Obliquity angle ϵ :

GeoOrientationData[Now,"Obliquity"]

dates=DateRange[DateObject[{-10000}],DateObject[{30000}],Quantity[10,"Years"]];

DateListPlot[Transpose[{dates,GeoOrientationData[#,"Obliquity"]&/@dates}]]

• Bias / Precesion / Nutation:

rotationAngleInAS[rot_]:=3600Max[Im[Log[Eigenvalues[rot]]]]/Degree

GeoOrientationData[Now,"BiasMatrix"]

rotationAngleInAS[%]

GeoOrientationData[Now,"PrecessionMatrix"]

rotationAngleInAS[%]

GeoOrientationData[Now,"NutationMatrix"]

rotationAngleInAS[%]

GeoOrientationData[Now,"BiasPrecessionNutationMatrix"]

rotationAngleInAS[%]

2) Rotation of Earth around its rotation axis

• Earth’s rotation is slowing down due to “tidal break”.

GeoOrientationData[DateObject[{2020,10,13}],"DayDuration"]

%//InputForm

GeoOrientationData[DateObject[{2020,10,13}],"DayDurationExcess","Around"]

• During a period of time:

DateInterval[{DateObject[{2018,1,1}],Today}]

GeoOrientationData[%,"DayDurationExcess"]

DateListPlot[%]

• Before 1960 we have data fitting observations, not direct data:

GeoOrientationData[DateInterval[{DateObject[{1900,1,1}],Today}],"DayDurationExcess"]//DateListPlot

GeoOrientationData[DateInterval[{{1000},{1900}}],"DayDurationExcess"]//DateListPlot

3) Motion of Earth’s rotation axis on Earth’s surface

• Position, from geodetic North, of the Earth’s rotation axis:

GeoOrientationData[Now,"PolarMotion"]

• Rotation axis between 2010 and 2015:

dates=DateRange[DateObject[{2010,1,1}],DateObject[{2015,1,1}],"Day"];

pm=GeoOrientationData[dates,"PolarMotion"]

The standard convention makes the x axis point toward Greenwich on the left and the y axis toward 90 degrees west on the top:

{ListPlot[pm,AspectRatio->Automatic,ScalingFunctions->{"Reverse",None}],}

Convert to a representation in meters:

pmM=pm/Quantity[1,"Radians"]GeodesyData["ITRF00","SemiminorAxis"]//UnitConvert

ListPlot[pmM,AspectRatio->Automatic,ScalingFunctions->{"Reverse",None}]

• This shows how the average position has moved in the last century (code in the refpage of GeoOrientationData):

GeoGraphics[GeoRange->{{60,90},All},GeoCenter->{90,0},GeoProjection->"Orthographic",GeoGridLines->Automatic,GeoZoomLevel->3]

Conclusions

• We have shown several improvements and updates in GeoGraphics, in our progress towards full vector maps.

• New function GeoOrientationData, establishing a connection between the Geo and Astro projects.