Vector GeoGraphics

Tim Shedelbower, Francisco Rodriguez Arias,
 Jose Martin-Garcia,
 Wolfram

#WolframTechConf

Abstract

GeoGraphics in Wolfram Language 14.3 returns vector maps by default, and includes many new map styles. Version 15.0 extends the Geo functionality with new options and added flexibility.

Summary

1. New geo background styles.

2. New tile data protocol.

3. New spherical geometry engine. Additions in WL15.

1.1. New default map style

The most noticeable change is the terrain texture.
More subtle changes include label haloing and roads that appear gradually.

GeoGraphics[]

1.2. Five new GeoBackground styles

The new styles include:

Terrain

Plain

Basic

Monochrome

BlackWhite

A dark version of each style is automatically applied when dark mode is detected.

Terrain

Earth terrain texture

Dark mode

Plain

Default style for all eleven geo visualization functions

Dark mode

Basic

Flat color map

Dark mode

Monochrome

Gray level map

Dark mode

BlackWhite

Contrast map

Dark mode

1.3. Better geographic visualization

GeoSmoothHistogram

Less clutter for a simpler appearance.

Style[Grid[{Magnify[GeoSmoothHistogram[GeoPosition[],GeoBackground->#,ImageSize->{Automatic,300}],.5]&/@{"Plain","StreetMapVectorLabels"},{"new default","previous"}}],"Text",12]

GeoGraphValuePlot

A brighter, cleaner view to improve data focus.

Style[Grid[{Magnify[GeoGraphValuePlot[{{(Entity["Rome", Entity["City",{"Rome","Lazio","Italy"}], "Entity[\"City\", {\"Rome\", \"Lazio\", \"Italy\"}]", "city"]),(Entity["Milan", Entity["City",{"Milan","Lombardy","Italy"}], "Entity[\"City\", {\"Milan\", \"Lombardy\", \"Italy\"}]", "city"]),2.2},{(Entity["Milan", Entity["City",{"Milan","Lombardy","Italy"}], "Entity[\"City\", {\"Milan\", \"Lombardy\", \"Italy\"}]", "city"]),(Entity["Paris", Entity["City",{"Paris","IleDeFrance","France"}], "Entity[\"City\", {\"Paris\", \"IleDeFrance\", \"France\"}]", "city"]),2.4},{(Entity["Milan", Entity["City",{"Milan","Lombardy","Italy"}], "Entity[\"City\", {\"Milan\", \"Lombardy\", \"Italy\"}]", "city"]),(Entity["Frankfurt", Entity["City",{"Frankfurt","Hesse","Germany"}], "Entity[\"City\", {\"Frankfurt\", \"Hesse\", \"Germany\"}]", "city"]),5},{(Entity["Paris", Entity["City",{"Paris","IleDeFrance","France"}], "Entity[\"City\", {\"Paris\", \"IleDeFrance\", \"France\"}]", "city"]),(Entity["Frankfurt", Entity["City",{"Frankfurt","Hesse","Germany"}], "Entity[\"City\", {\"Frankfurt\", \"Hesse\", \"Germany\"}]", "city"]),3.2},{(Entity["Paris", Entity["City",{"Paris","IleDeFrance","France"}], "Entity[\"City\", {\"Paris\", \"IleDeFrance\", \"France\"}]", "city"]),(Entity["Madrid", Entity["City",{"Madrid","Madrid","Spain"}], "Entity[\"City\", {\"Madrid\", \"Madrid\", \"Spain\"}]", "city"]),4.3},{(Entity["Paris", Entity["City",{"Paris","IleDeFrance","France"}], "Entity[\"City\", {\"Paris\", \"IleDeFrance\", \"France\"}]", "city"]),(Entity["London", Entity["City",{"London","GreaterLondon","UnitedKingdom"}], "Entity[\"City\", {\"London\", \"GreaterLondon\", \"UnitedKingdom\"}]", "city"]),3.3},{(Entity["Madrid", Entity["City",{"Madrid","Madrid","Spain"}], "Entity[\"City\", {\"Madrid\", \"Madrid\", \"Spain\"}]", "city"]),(Entity["Seville", Entity["City",{"Seville","Seville","Spain"}], "Entity[\"City\", {\"Seville\", \"Seville\", \"Spain\"}]", "city"]),3.6},{(Entity["Madrid", Entity["City",{"Madrid","Madrid","Spain"}], "Entity[\"City\", {\"Madrid\", \"Madrid\", \"Spain\"}]", "city"]),(Entity["Barcelona", Entity["City",{"Barcelona","Barcelona","Spain"}], "Entity[\"City\", {\"Barcelona\", \"Barcelona\", \"Spain\"}]", "city"]),1.8}},GeoBackground->#,ImageSize->{Automatic,300}],.5]&/@{"Plain","StreetMapVectorLabels"},{"new default","previous"}}],"Text",12]

1.4. Fades

Terrain fade

The terrain fades to flat color as zoom level increases.

Road fade

Roads gradually appear as zoom level increases.

Initial road appearance now emerges at a higher zoom level.
The previous GeoBackground was ”StreetMapVectorLabels”.

2.1 GeoGraphics vector map protocol: Protomaps

GeoGraphics now uses Protomaps as main protocol to get the vector map data.

Protomaps format: PMTiles

PMTiles is the format specification used by Protomaps.

Very efficient way to store all vector tiles. It has a directory-like internal structure that allows saving space for empty tiles and duplicate tiles.

MBTiles is a well known container of tiles too, but the main difference is that PMTiles allows accessing the tiles using the HTTP-range specification directly from a remote server.

2.2. Protomaps structure

The PMTiles file specification divides the file in a header, root directory, metadata, leaf directores and tile data.

Header

The header contains important information about the internal structure of the content.

pmtilesgeo="http://geo.wolfram.com/vector/data/planet.pmtiles";
pmtiles=FileNameJoin[{GeoTools`GeoServerVectorResources`Private`$ResourcesPlanetDir,"planet-z0z4.pmtiles"}];
Dataset[Transpose@Association@{"zoom 0-4"->GeoTools`Protomaps`PMHead["Association",pmtiles],"full planet"-> GeoTools`Protomaps`PMHead["Association",pmtilesgeo]},MaxItems->11]

Metadata

Metadata contains the description of the tileset, like tile format, available zoom levels, and it includes the layers structure too.

Dataset[Transpose@Association@{"zoom 0-4"->GeoTools`ProtomapsMetadata[pmtiles],"full planet"-> GeoTools`ProtomapsMetadata[pmtilesgeo]}/.s_String/;StringLength[s]>20:>StringTake[s,10]<>"...",MaxItems->14]

Tile data

Individual tiles are stored as MVT (Mapbox Vector Tile format). Each tile is assigned an integer index which increases like a Hilbert curve corresponding to a Z/X/Y tile position.

tileCenter[{z_,x_,y_}]:=With[{semicirc=2^(z-1)},With[{proj={"Mercator","ReferenceModel"->semicirc/Pi,"GridOrigin"->semicirc {1,-1}}},GeoGridPosition[{x+0.5,-y-0.5},proj]]];
tileID[{z_,x_,y_}]:=Text[GeoTools`Protomaps`Private`zxyToTileID[{z,x,y}],tileCenter[{z,x,y}],Center];
tileIDTable[z_]:=Table[tileID[{z,x,y}],{x,0,2^z-1},{y,0,2^z-1}]
tileRegionBoundary2[{z_,x_,y_}]:=GeoBoundsRegionBoundary[GIS`TileGeoRange[{z,x,y},{"Wolfram","Projection"->"Mercator"}]];
tileTable[z_]:=Table[tileRegionBoundary2[{z,x,y}],{x,0,2^z-1},{y,0,2^z-1}];
tileRegion2[{z_,x_,y_}]:=GeoBoundsRegion[GIS`TileGeoRange[{z,x,y},{"Wolfram","Projection"->"Mercator"}]];
tileRegionTable[z_]:=With[{first=GeoTools`Protomaps`Private`zxyToTileID[{z,0,0}],n=4^z},Table[{GeoStyling[GrayLevel[(GeoTools`Protomaps`Private`zxyToTileID[{z,x,y}]-first)/n,0.8]],tileRegion2[{z,x,y}]},{x,0,2^z-1},{y,0,2^z-1}]]
rotate[n_,x_,y_,rx_,ry_]:=Module[{xx=x,yy=y},If[ry==0,If[rx!=0,xx=n-1-xx;
yy=n-1-yy;];
{yy,xx},{xx,yy}]];

IDToZxy[i_]:=Module[{z,acc,t,tx=0,ty=0,a,s,rx,ry},z=Floor[(BitLength[3 i+1]-1)/2];
acc=(2^(2 z)-1)/3;
t=i-acc;
For[a=0,a<z,a++,s=BitShiftLeft[1,a];
rx=BitAnd[1,BitShiftRight[t,1]];
ry=BitAnd[1,BitXor[t,rx]];
{tx,ty}=rotate[s,tx,ty,rx,ry];
tx+=BitShiftLeft[rx,a];
ty+=BitShiftLeft[ry,a];
t=BitShiftRight[t,2];];
{z,tx,ty}]
tileLine[z_]:=Line@With[{first=GeoTools`Protomaps`Private`zxyToTileID[{z,0,0}]},Table[tileCenter[IDToZxy[first+i]],{i,0,4^z-1}]]
GraphicsRow[Table[GeoGraphics[{{Gray,Thin,tileTable[z]},{Red,Thick,tileLine[z]},Haloing[GrayLevel[0.9,0.9],1,1],Darker@DarkRed,Style[tileIDTable[z],12-2^z,Bold]},GeoZoomLevel->Min[2,z],GeoRange->"World",GeoProjection->"Mercator",GeoBackground->"PlainNoLabels",ImageSize->Medium],{z,0,3}],ImageSize->Full]

This indexing means that neighbouring tiles are often found within the same block of data within the file.

This map of South America needs the 9 tiles: 31..33, 38..39, 49..52.

tileLine2[z_]:=GeoPath[With[{first=GeoTools`Protomaps`Private`zxyToTileID[{z,0,0}]},Table[tileCenter[IDToZxy[first+i]],{i,0,4^z-1}]],"Loxodrome"];
tileRegion2[id_Integer]:=GeoBoundsRegion[GIS`TileGeoRange[IDToZxy[id],{"Wolfram","Projection"->"Mercator"}]];
With[{z=3,h=200},{gs=GeoGraphics[{{Gray,Thin,tileTable[z]},{Red,Thick,tileLine2[z]},Haloing[GrayLevel[0.9,0.9],1,1],Darker@DarkRed,Style[tileIDTable[z],12-2^z,Bold]},GeoZoomLevel->Min[2,z],GeoRange->Entity["GeographicRegion","SouthAmerica"],GeoProjection->"Mercator",GeoBackground->"PlainNoLabels",ImageSize->{Automatic,h}]},GraphicsRow[{gs,GeoGraphics[{{Gray,Thin,tileTable[z]},{EdgeForm[{Black}],Gray,tileRegion2/@{31,32,33,38,39,49,50,51,52}},{EdgeForm[{DarkGray,Thick,Opacity[0.75]}],Red,GeoBoundsRegion[GeoRange/.Options[gs,GeoRange]]},{Red,Thick,tileLine2[z]},Haloing[GrayLevel[0.9,0.9],1,1],Darker@DarkRed,Style[tileIDTable[z],12-2^z,Bold]},GeoZoomLevel->Min[2,z],GeoRange->"World",GeoProjection->"Mercator",GeoBackground->"PlainNoLabels",ImageSize->{Automatic,h}]}]]

2.3. Advantages of using Protomaps

Access using any regular HTTP server or local file calls.

Faster tiles download by bulk downloading them with one HTTP call.

Hilbert-like indexing makes neighbouring tiles to be close within the file, allowing faster retrieval server-side.

Allows easier local permanent cache of tiles.

Reduced server storage needs: PMTiles is more compact.

Standard vector tile format.

Allows tiles with raster formats too.

3.1. New in WL14.3: New spherical engine

• Better support of very large primitives, like geo hemispheres, or larger.

• Current night hemisphere:

GeoGraphics[NightHemisphere[]]

• Twilight regions (new in 15.0):

nightHemispheres={
NightHemisphere[],
NightHemisphere[ReferenceAltitude->"Civil"],
NightHemisphere[ReferenceAltitude->"Nautical"],
NightHemisphere[ReferenceAltitude->"Astronomical"]
};

GeoGraphics[nightHemispheres]

Here we see more clearly the 6 degree bands:

Manipulate[GeoGraphics[nightHemispheres,GeoProjection->proj,GeoCenter->{lat,lon}],{{lat,0},-90,90},{{lon,0},-180,180},{proj,{"Orthographic","LambertAzimuthal","Bonne","Mollweide"}},SaveDefinitions->True]

3.2. New in WL15: GeoAxes

• Projected vs unprojected coordinates in the Mercator projection:

GeoGraphics["World",GeoBackground->"Basic",GeoProjection->"Mercator",Frame->True,GeoGridLines->Automatic,GeoAxes->True]

• Select where the geo axes will cross:

GeoGraphics["World",GeoBackground->"Basic",GeoProjection->"Mercator",GeoAxes->True,GeoAxesOrigin->Here]

• GeoAxes can be used with any projection:

GeoGraphics["World",GeoBackground->"Basic",GeoProjection->"PeirceQuincuncial",GeoGridLines->Quantity[10,"AngularDegrees"],GeoAxes->Automatic]

Manipulate[GeoGraphics[GeoRange->"World",GeoProjection->"LambertAzimuthal",GeoAxes->{<|"AxisStyle"->{Thick,Green}|>,<|"AxisStyle"->{Thick,Red}|>},GeoAxesOrigin->GeoPosition[{lat,lon}],GeoGridLines->True,ImageSize->360],{{lon,0},-180,180},{{lat,0},-90,90}]

• Use arbitrary AxisObject expressions:

cities={(Entity["Buenos Aires", Entity["City",{"BuenosAires","BuenosAires","Argentina"}], "Entity[\"City\", {\"BuenosAires\", \"BuenosAires\", \"Argentina\"}]", "city"]),(Entity["Tokyo", Entity["City",{"Tokyo","Tokyo","Japan"}], "Entity[\"City\", {\"Tokyo\", \"Tokyo\", \"Japan\"}]", "city"])};

dist=GeoDistance@@cities

ticks=Append[Range[(Quantity[0, "mi", miles, "Miles"]),dist,(Quantity[2000, "mi", miles, "Miles"])],dist]

GeoGraphics[AxisObject[GeoPath[cities],{0,QuantityMagnitude@dist},TickPositions->{{QuantityMagnitude@ticks}},TickLabels->{ticks}],
GeoGridLines->Automatic,GeoProjection->"Bonne",ImageSize->360,GeoBackground->"Plain",GeoRange->All]

3.3. New in WL15: NightHemisphere on other bodies

• Phase of the Moon next week, centering at {0, 0}:

date=Now+Quantity[7,"Days"];

GeoGraphics[{GeoStyling[Opacity[.7]],NightHemisphere[date]},GeoModel->"Moon",GeoProjection->"Orthographic",Background->Black,GeoGridRangePadding->0.1]

Full synodic month, as observed from the center of the Earth:

Manipulate[
GeoGraphics[{GeoStyling[Opacity[.7]],Black,NightHemisphere[date]},GeoRange->All,GeoModel->"Moon",GeoProjection->"Orthographic",Background->Black,GeoCenter->AstroSubpoint["Earth",{"Moon",date}],PlotRange->AstroDistance["Moon",Dated["Earth",date]]/Quantity[350000,"Kilometers"]],
{date,Now,Now+Quantity[1,"SynodicMonths"],Quantity[1,"Hours"]}
]

• Current phase of Venus, as observed from our location:

GeoGraphics[{GeoStyling[Opacity[.7]],NightHemisphere[date]},GeoCenter->AstroSubpoint[Here,{"Venus",date}],GeoModel->"Venus",GeoProjection->"Orthographic",Background->Black,GeoGridRangePadding->0.1]

Full period of 594 days:

Manipulate[GeoGraphics[{GeoStyling[Opacity[.7]],Black,NightHemisphere[date]},GeoRange->All,GeoModel->body,GeoProjection->"Orthographic",Background->Black,GeoCenter->AstroSubpoint["Earth",{body,date}],PlotRange->AstroDistance[body,Dated["Earth",date]]/mindist],{date,Now,Now+T,Quantity[1,"Hours"]},Initialization:>(body=(Entity["Venus", Entity["Planet","Venus"], "Entity[\"Planet\", \"Venus\"]", "planet"]);
T=Abs[1/(1/body["OrbitPeriod"]-1/(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"])["OrbitPeriod"])];
mindist=AstroDistance[body,Dated["Earth",FindAstroEvent[{"Distance",body,"Earth","GlobalMinimum"},DateInterval[{Now,Now+T}]]]];)]

• Reproduce view of Earth as seen from the Apollo 11 landing site in 1969:

https://upload.wikimedia.org/wikipedia/commons/8/83/AS11-40-5923_%2821669280451%29.jpg

Without rotation and clouds:

With[{date=DateObject[{1969,7,21,4,19,0},TimeZone->0]},
GeoGraphics[{{GeoStyling["SatelliteNight",GeoZoomLevel->2],NightHemisphere[date]}},GeoProjection->"Orthographic",GeoGridLines->Automatic,GeoCenter->AstroSubpoint[(Entity["Apollo 11", Entity["MannedSpaceMission","Apollo11"], "Entity[\"MannedSpaceMission\", \"Apollo11\"]", "manned space mission"]),date],Background->Black,GeoGridRangePadding->0.1]
]

3.4. New in 14.3: GeoReposition

• Move any polygon anywhere in the world (or even to other planets).

GeoGraphics[GeoReposition[(Entity["Antarctica", Entity["GeographicRegion","Antarctica"], "Entity[\"GeographicRegion\", \"Antarctica\"]", "geographic region"])["Polygon"],"SouthPole"->Here],GeoRange->"World",GeoBackground->"Plain"]

pol=(Entity["Australia", Entity["Country","Australia"], "Entity[\"Country\", \"Australia\"]", "country"])["Polygon"];
city=GeoPosition[(Entity["Melbourne", Entity["City",{"Melbourne","Victoria","Australia"}], "Entity[\"City\", {\"Melbourne\", \"Victoria\", \"Australia\"}]", "city"])];

Manipulate[GeoGraphics[{Green,Point[city],Red,Point[GeoPosition[{lat,lon}]],GeoReposition[pol,{city->GeoPosition[{lat,lon}],alpha}]},GeoRange->"World",ImageSize->300,GeoProjection->"Mercator",GeoBackground->"TerrainNoLabels"],{{lat,city["Lat"]},-90,90},{{lon,city["Lon"]},-180,180},{{alpha,0},0,360},SaveDefinitions->True]

Conclusions

• Wolfram GeoGraphics has been significantly upgraded in versions WL 14.3 and 15.0: vector maps, new styles, new data protocols, new options, and much more.

• Next year we add more geo computation functionality, more geo data, more geo formats, and also provide access to styling the many individual geometric elements of vector maps.

• Check out the extensive Geo documentation, including ref pages like GeoGraphics and GeoBackground, guides like GeoVisualization, and the Geography tutorial.