Geo Computation

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

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

Abstract

This talk reports on the status of several Geo projects, as well as discussing some points of the relation between Geo and the new Astro project in version 13.2 of Wolfram Language.

Introduction

Different fronts of progress in Geo:

Continue improving vector maps.

Add new Geo formats, by linking to the GDAL library.

Link to external libraries for geometric and projection computations.

Update of provider of satellite imagery.

Additions required by the Astro project, mainly “obliqueness”.

Connection between the Geo and Astro projects. Much in common.

(Advance part of the Astro talk to this one ☺.)

Outline

1) Vector maps

2) Geo formats

3) Geo computation

4) New satellite imagery

5) Obliqueness

6) Reference frames and time systems

Conclusions

1. Vector Maps

Use GeoBackground to access vector maps:

GeoGraphics[Here,GeoBackground->"VectorClassic"]

Compare with our current default with raster images:

GeoGraphics[Here]

2. Geo Formats

More formats via newer version of GDAL.

Example of vector format: OSM

File taken from openstreemap.org:

Import["~/Downloads/map.osm"]

Example of raster format: GXF

Bouguer gravitational anomaly in the area of Lake Superior [USGS]:

image=Import["~/Downloads/LakeSupRegBouguer.gxf"]

The projection is “UTMZone16”:

proj=Import["~/Downloads/LakeSupRegBouguer.gxf","Projection"]

GeoProjectionData["UTMZone16"]

Extract the ranges of projected coordinates and construct a polygon:

gridrange=Import["~/Downloads/LakeSupRegBouguer.gxf","SpatialRange"]

{{xmin,xmax},{ymin,ymax}}=gridrange;

pol=Polygon[GeoGridPosition[{{xmin,ymin},{xmin,ymax},{xmax,ymax},{xmax,ymin}},proj]];

Overlay that image over a relief map of the same area:

GeoGraphics[{GeoStyling[{image,"Projection"->proj,"Bounds"->gridrange},Opacity[0.7]],pol},
GeoProjection->proj,GeoGridRange->gridrange,GeoGridRangePadding->50000,GeoBackground->"ReliefMap",ImageSize->Large]

3. Geo Computation

Unprojected polygon data has currently form Polygon[GeoPosition[...]] and expresses holes using FilledCurve:

world=Entity["GeographicRegion","World"]["Polygon"]

GIS` utility GIS`FilledCurveToPolygon:

GIS`FilledCurveToPolygon[world]

GeoGraphics[%]

Show[%,GeoRange->(Entity["Netherlands", Entity["Country","Netherlands"], "Entity[\"Country\", \"Netherlands\"]", "country"])]

4. New Satellite imagery

• Recall the correspondence between zoom and geo range at the equator:

Table[zoom->Quantity[40. 10^6,"Meters"]/2^zoom,{zoom,1,18}]

For each pixel:

Table[zoom->Quantity[40. 10^6,"Meters"]/256/2^zoom,{zoom,1,18}]

• For zooms 1 (whole world) to 8 (~100 miles range) we still use the same tiles:

GeoImage[]

Switch to the new data provider:

GeoImage[GeoServer->"Mapbox"]

GeoImage[Here]

• Satellite background in GeoGraphics:

GeoGraphics[Here,GeoBackground->"Satellite",ImageSize->Large]

5. Obliqueness

Spherical rotations are fundamental in the Astro project.

Example with an Earth polygon:

Manipulate[Graphics[GIS`EquirectangularView[Polygon[GeoPosition[]],{{lat,lon},alpha},All,0.1],ImageSize->1000],{{lat,0},-90,90},{{lon,0},-180,180},{{alpha,0},-180,180}]

Example with a coordinate grid (or “graticule”):

grid={
{Brown,Line[GeoPosition[Table[{lat,lon},{lat,-75,75,15},{lon,-180,180}]]]},
{ColorNegate[Brown],Line[GeoPosition[Table[{lat,lon},{lon,-180,180,15},{lat,If[Mod[lon,90]==0,-90,-75],If[Mod[lon,90]==0,90,75]}]]]},
{Red,Line[GeoPosition[Table[{0,lon},{lon,-180,180}]]]},
{ColorNegate[Red],Line[GeoPosition[Table[{lat,0},{lat,-90,90}]]]}
};
Manipulate[Graphics[{Thick,GIS`EquirectangularView[grid,{{lat,lon},alpha},All,0.1]},ImageSize->1000],{{lat,0},-90,90},{{lon,0},-180,180},{{alpha,0},-180,180}]

• We have some projections with built-in obliqueness:

GeoProjectionData["ObliqueMercator"]

GeoGraphics[GeoRange->"World",GeoProjection->{"ObliqueMercator","Centering"->{Here,90}},GeoGridLines->Automatic,PlotRange->Pi{{-1,1},{-1,1}}]

6. Reference frames and time systems

A position in space/time is always given wrt something else:

Physics: (reference) frame, time system. Datum in Geo.

Maths: coordinate system, calendar. Projection in Geo.

Precise results only possible using precise frames.

Astro: Spherical geometry. Geo: Ellipsoidal geometry.

Astro: Frames move; high speed (relativity). Geo: Tectonic drift.

6.1. Time Systems

Atomic time at the surface of the Earth.

Variants: “TAI”, “TT”, “UTC”.

Time determined by the rotation of the Earth.

Variants: “UT1”.

Time at places other than Earth surface.

Variants: “TCG”, “TCB”, “TDB”, ...

dates=DateRange[DateObject[{1960,1,1}],Today,"Month"];
dt=GeoOrientationData[dates,"TAIMinusUT1"];
leap=GeoOrientationData[dates,"TAIMinusUTC"];

Orange: TT - UT1 ± 0.9s,  Blue: TT - UTC

diff=Quantity[32.184,"Seconds"];
Show[DateListPlot[{Transpose[{dates,dt+diff-Quantity[0.9,"Seconds"]}],Transpose[{dates,dt+diff+Quantity[0.9,"Seconds"]}]},PlotStyle->Orange],DateListPlot[Transpose[{dates,leap+diff}]],ImageSize->Large,PlotRange->{32,72}]

Note: J2000 (or J2000.0) is

DateObject[{2000.},CalendarType->"JulianYear",TimeSystem->"TT"]

CalendarConvert[%]

TimeSystemConvert[%,"UTC"]

First[%]

New in 13.2:

RandomDate[]

RandomDate[DateObject[{2022,10}]]

RandomDate[Now+Quantity[{-10,10},"Days"]]

RandomDate[Quantity[-14,"Days"]]

RandomDate[Quantity[-14,"Days"],5]

RandomDate[Quantity[-14,"Days"],5,DateGranularity->"Day"]

6.1. Main Four Families of Frames

“Equatorial”: use Earth’s equatorial plane. Variants:

({{Frame, Default Origin, Rotation, X Axis}, {"ICRS", SSB, none, ~ mean equinox of J2000}, {"J2000", SSB, none, mean equinox of J2000}, {"GCRS", Earth, none, same as ICRS}, {"TETE", Earth, prec/nut, true equinox}, {"MEME", Earth, prec, mean equinox}, {"CIRS", Earth, prec/nut, CIO}, {"TIRS", Earth, day/prec/nut, TIO}, {"ITRS", Earth, polar/day/prec/nut, geo X axis}})

“Ecliptic”: use Earth’s orbital plane. Variants:

({{Frame, Default Origin, Rotation, X Axis}, {"TrueEcliptic", Earth, ecl prec, true equinox}, {"MeanEcliptic", SSB, ecl prec, mean equinox}})

“Galactic”: use Milky Way’s plane.

({{Frame, Default Origin, Rotation, X Axis}, {"Galactic", SSB, none, galaxy center (~1950)}})

“Horizon” frames: use local horizontal plane.

({{Frame, Default Origin, Rotation, X Axis}, {"Horizon", Here, all, local North or South}})

Conclusions

The Geo functionality continues advancing in version 13.2.

New geo formats, new satellite data.

Strong interaction between the Geo and Astro projects.

Next steps:

Complete the transition to vector maps.

Stronger geo region computation.

Many more new geo formats.

Thank you!