R package leafem has recently gained a suite of functions to visualise cloud optimised vector data. Therse are:
addPMPolygons(), addPMPoints(), addPMPolylines()addFgb() (has been introduced some time ago and is now the default mechanism used in mapview)For details on how to install and use software needed to create cloud optimised data that can be visualised with these function, please see this tutorial.
In this tutorial we will introduce these functions in detail and highlight the different options we can use to control appearance of the visualisations.
I’ve prepared an AWS S3 bucket with some vector tile data. Let’s see what is in the bucket.
library(paws)
svc = s3(config = list(region = "eu-central-1"))
obj = svc$list_objects(Bucket = "vector-tiles-data")
sapply(obj$Contents, '[[', "Key")[1] "depoints.pmtiles" "nz-building-outlines.fgb"
[3] "nz-building-outlines.pmtiles" "rivers_africa.fgb"
[5] "rivers_africa.pmtiles" So, there’s a couple of FlatGeobuf files (.fgb) and a few PMTiles files (.pmtiles). Both of these formats are capable of streaming http range requests so they both classify as cloud optimised data formats.
FlatGeobuf is a performant binary encoding format for geographic vector data. Please refer to the Flatgeobuf github repository for details about this format.
PMTiles is a single-file archive format for tiled data that is also able to hold raster data (but we will use it for vector data only in this tutorial). Again, please refer to the PMTiles github repository for details about this format.
Let’s start by visualising one of the .fgb files. When visualising large FlatGeobuf data sets we need to be carful to not overload our browsers capacity because in the standard function definition, the function addFgb() will attempt to fetch & visualise all data from the data set. However, there are two arguments called minZoom and maxZoom that can be used to stream only those parts of the data that are in the current view of the map.
library(leaflet)
library(leafem)
nzb_fgb_url = "https://vector-tiles-data.s3.eu-central-1.amazonaws.com/rivers_africa.fgb"
leaflet(height = "800px", width = "100%") |>
addTiles(group = "osm") |>
addProviderTiles(
"Esri.WorldImagery"
, group = "esri"
) |>
addFgb(
url = nzb_fgb_url
, layerId = "rivers_fgb"
, group = "rivers_fgb"
, popup = TRUE
, minZoom = 6
, maxZoom = 10
, weight = 2
) |>
addMouseCoordinates() |>
addLayersControl(
baseGroups = c("osm", "esri")
, overlayGroups = c("rivers_fgb")
) |>
setView(0, 0, 2)So, if we want to also see our data at lower zoom levels, we can use .pmtiles data. Let’s investigate the rivers_africa.pmtiles.
pmtiles-show /home/tim/tappelhans/privat/data/rivers_africa_37333/rivers_africa.pmtilesspec version: 2
metadata:
name = rivers_africa.mbtiles
description = rivers_africa.mbtiles
version = 2
minzoom = 0
maxzoom = 5
center = 16.875000,16.560724,5
bounds = -17.239583,-34.764583,54.350000,37.218750
type = overlay
format = pbf
generator = tippecanoe v2.1.0
generator_options = tippecanoe -zg '--projection=EPSG:4326' --no-tile-compression --no-feature-limit --no-tile-size-limit -o rivers_africa.mbtiles -l rivers_africa rivers_africa.fgb
json = {"vector_layers": [ { "id": "rivers_africa", "description": "", "minzoom": 0, "maxzoom": 5, "fields": {"ARCID": "Mixed", "A_Strahler": "Mixed", "FID_af_str": "Mixed", "FID_sub_ba": "Mixed", "FROM_NODE": "Mixed", "LEGEND": "Mixed", "MAJ_AREA": "Mixed", "MAJ_BAS": "Mixed", "MAJ_NAME": "String", "RASTERVA_1": "Mixed", "RASTERVA_2": "Number", "Regime": "String", "SUBBAS_ID": "Mixed", "SUB_BAS": "Mixed", "SUB_NAME": "String", "Strahler": "Mixed", "TOBAS_ID": "Mixed", "TO_NODE": "Mixed"} } ],"tilestats": {"layerCount": 1,"layers": [{"layer": "rivers_africa","count": 185476,"geometry": "LineString","attributeCount": 18,"attributes": [{"attribute": "ARCID","count": 1000,"type": "mixed","values": ["1","10","100","1000","10000","100000","100001","100002","100003","100004","100005","100006","100007","100008","100009","10001","100010","100011","100012","100013","100014","100015","100016","100017","100018","100019","10002","100020","100021","100022","100023","100024","100025","100026","100027","100028","100029","10003","100030","100031","100032","100033","100034","100035","100036","100037","100038","100039","10004","100040","100041","100042","100043","100044","100045","100046","100047","100048","100049","10005","100050","100051","100052","100053","100054","100055","100056","100057","100058","100059","10006","100060","100061","100062","100063","100064","100065","100066","100067","100068","100069","10007","100070","100071","100072","100073","100074","100075","100076","100077","100078","100079","10008","100080","100081","100082","100083","100084","100085","100086"]},{"attribute": "A_Strahler","count": 8,"type": "mixed","values": ["1","2","3","4","5","6","7","8"]},{"attribute": "FID_af_str","count": 1000,"type": "mixed","values": ["0","1","10","100","1000","10000","100000","100001","100002","100003","100004","100005","100006","100007","100008","100009","10001","100010","100011","100012","100013","100014","100015","100016","100017","100018","100019","10002","100020","100021","100022","100023","100024","100025","100026","100027","100028","100029","10003","100030","100031","100032","100033","100034","100035","100036","100037","100038","100039","10004","100040","100041","100042","100043","100044","100045","100046","100047","100048","100049","10005","100050","100051","100052","100053","100054","100055","100056","100057","100058","100059","10006","100060","100061","100062","100063","100064","100065","100066","100067","100068","100069","10007","100070","100071","100072","100073","100074","100075","100076","100077","100078","100079","10008","100080","100081","100082","100083","100084","100085"]},{"attribute": "FID_sub_ba","count": 603,"type": "mixed","values": ["32334","32335","32336","32337","32338","32339","32340","32341","32342","32343","32344","32345","32346","32347","32348","32349","32350","32351","32352","32353","32354","32355","32356","32357","32358","32359","32360","32361","32362","32363","32364","32365","32366","32367","32368","32369","32370","32371","32372","32373","32374","32375","32376","32377","32378","32379","32380","32381","32382","32383","32384","32385","32386","32387","32388","32389","32390","32391","32392","32393","32394","32395","32396","32397","32398","32399","32400","32401","32402","32403","32404","32405","32406","32407","32408","32409","32410","32411","32412","32413","32414","32415","32416","32417","32418","32419","32420","32421","32422","32423","32424","32425","32426","32427","32428","32429","32430","32431","32432","32433"]},{"attribute": "FROM_NODE","count": 1000,"type": "mixed","values": ["10","1000","10000","100000","100001","100002","100003","100004","100005","100006","100007","100008","100009","10001","100010","100011","100012","100013","100014","100015","100016","100017","100018","100019","10002","100020","100021","100022","100023","100024","100025","100026","100027","100028","100029","10003","100030","100031","100032","100033","100034","100035","100036","100037","100038","100039","10004","100040","100041","100042","100043","100044","100045","100046","100047","100048","100049","10005","100050","100051","100053","100054","100055","100056","100057","100058","100059","10006","100060","100061","100062","100063","100064","100065","100066","100067","100068","100069","10007","100070","100071","100072","100073","100074","100075","100076","100077","100078","100079","10008","100080","100081","100082","100083","100084","100085","100086","100087","100088","100089"]},{"attribute": "LEGEND","count": 25,"type": "mixed","values": ["1","10","11","14","15","16","17","18","19","2","20","21","22","23","24","25","26","27","3","4","5","6","7","8","9"]},{"attribute": "MAJ_AREA","count": 25,"type": "mixed","values": ["1010044","102100","1041192","1373296","2136941","244531","2461890","260457","3074955","3696670","403126","411058","411553","477345","499542","5531340","558292","596220","638878","699755","772710","796599","809724","863869","984867"]},{"attribute": "MAJ_BAS","count": 25,"type": "mixed","values": ["7001","7002","7003","7004","7005","7006","7007","7008","7009","7010","7011","7014","7015","7016","7017","7018","7019","7020","7021","7022","7023","7024","7025","7026","7027"]},{"attribute": "MAJ_NAME","count": 25,"type": "string","values": ["Africa, East Central Coast","Africa, Indian Ocean Coast","Africa, North Interior","Africa, North West Coast","Africa, Red Sea - Gulf of Aden Coast","Africa, South Interior","Africa, West Coast","Angola, Coast","Congo","Gulf of Guinea","Lake Chad","Limpopo","Madasgacar","Mediterranean South Coast","Namibia, Coast","Niger","Nile","Orange","Rift Valley","Senegal","Shebelli - Juba","South Africa, South Coast","South Africa, West Coast","Volta","Zambezi"]},{"attribute": "RASTERVA_1","count": 1000,"type": "mixed","values": ["10000","100000","100002","1000062","10001","1000104","100011","1000146","100016","100028","10003","1000312","1000360","10004","100040","1000420","100044","100054","1000560","10006","100060","10006017","100062","10006380","100065","10006689","10007","100073","100079","1000797","10008","100086","1000860","100089424","10010","100103984","100107","100111","100113","100120","10012120","100128","10013","100130","100134","100140","1001466","100147","100149","1001532","10015680","1001580","100164840","1001658","1001679","100168","10017","100170","100191","100191376","1001970","10020","100200","100206","10021","100212","10022","100220","100225","10023","100233","10023489","10023531","1002372","100240","100244","10024460","100249","100251936","1002535","10026","100260","1002603","100263","1002645","10027","100275136","10028","1002834","100296","1002978","10030","100300","10030160","1003032","1003107","1003160","100317","10032","100320"]},{"attribute": "RASTERVA_2","count": 1000,"type": "number","values": [0,1,10,100,1000,1001.00006104,1002,1003,1003.99993896,1005,1006.00006104,1007,1008,1008.99993896,101,1010,1011.00006104,1011.99993896,1013,1014.00006104,1015,1016.00006104,1016.99993896,1018,1019.00006104,102,1020,1021,1021.99993896,1023,1024,1025,1026,1027,1028,1029,103,1030,1031,1032,1033,1034,1035,1036,1037,1038,1039,104,1040,1041,1042,1043,1044,1045,1046,1047,1048,1049,105,1050,1051,1052,1053,1054,1055,1056,1057,1058,1059,106,1060,1061,1062,1063,1064,1065,1066,1067,1068,1069,107,1070,1071,1072,1073,1074,1075,1076,1077,1078,1079,108,1080,1081,1082,1083,1084,1085,1086,1087],"min": 0,"max": 3058},{"attribute": "Regime","count": 2,"type": "string","values": ["I","P"]},{"attribute": "SUBBAS_ID","count": 597,"type": "mixed","values": ["7010978","7010981","7010982","7010983","7010984","7010985","7010986","7010987","7010988","7010989","7019831","7020081","7020086","7020812","7020813","7020814","7020816","7020819","7020821","7020822","7020823","7020824","7020825","7020826","7020827","7020828","7020829","7020831","7020832","7020834","7020835","7020836","7020837","7020838","7020839","7020841","7020842","7020843","7020844","7020845","7020846","7020848","7020849","7020851","7020852","7020853","7020854","7020856","7020858","7020859","7020871","7020872","7020873","7020874","7020875","7020877","7020878","7020879","7020881","7020882","7020883","7020884","7020885","7020886","7020887","7020888","7020889","7020891","7020892","7020893","7020894","7020895","7020896","7020897","7020898","7020899","7028332","7028911","7030003","7030029","7030198","7030210","7030211","7030212","7030213","7030214","7030215","7030216","7030217","7030218","7030219","7030220","7030221","7030222","7030223","7030227","7030228","7030229","7030234","7030239"]},{"attribute": "SUB_BAS","count": 597,"type": "mixed","values": ["100030","100302","100306","100309","100324","100327","100329","100370","103190","103192","103271","103701","103702","103703","103704","103705","103706","103707","10978","10981","10982","10983","10984","10985","10986","10987","10988","10989","110570","115701","115702","115703","115704","115705","115709","149900","149902","149903","149904","149905","149906","149907","149908","149909","149910","150395","150396","150398","153971","153972","153973","153974","153975","153976","160552","160554","160555","160556","165530","165531","165532","165533","165534","165535","165537","165538","165539","165570","170512","170513","170514","170515","170518","170521","170523","170524","170525","170526","170527","170528","170529","170532","170539","180371","180372","180374","180376","180377","180378","180379","180381","180382","180383","180384","180386","180387","180391","180392","180393","180394"]},{"attribute": "SUB_NAME","count": 592,"type": "string","values": ["Abu Hut","Adrar Sotuf","Afram","Afrera Ye Chew Hayk","Agneby","Akoba","Al Ghallah","Al Hawad","Alal","Algerian Atlas","Algerian east coast","Algerian west coast","Alibori","Alima","Anambra","Andokat","Angwa","Ankobra","Aruwimi","Asa","Auob","Awash Wenz 1","Awash Wenz 2","Az Zayn","Ba","Bafing","Bagoe 1","Bagoe 2","Bahi swamps","Bahr Azum","Bahr Keita","Bahr al Arab","Bakoy","Bandah","Bandama 1","Bandama 2","Bandama 3","Bandama Blanc 1","Bandama Blanc 2","Bani 1","Bani 2","Bani 3","Baoule","Baro Wenz","Bay al Kebir","Beli","Bengo","Benito","Benue 1","Benue 2","Benue 3","Benue 4","Benue 5","Betsiboka","Bia","Bibamba","Black Umbuluzi","Blue Nile 1","Blue Nile 2","Blue Nile 3","Blue Nile 4","Blue Nile 5","Bombo","Borkou","Botletli","Bou","Bou Regreg","Bouluo","Boumba","Brak","Brass River/Niger Delta","Breede","Buhayrat Abyad","Bunsuru","Busira 1","Busira 2","Buzi","Cabinda Coast","Caledon","Capoche","Carnarvon Leegte","Casamance","Catumbela","Cavalla","Central Mauritania","Cestas","Chanchaga","Changane","Chari 1","Chari 2","Chebeika","Chelif","Chinko","Chire","Chissonga","Chongwe","Chott Chergui","Chott Hodna","Coast south of Congo","Comoe 1"]},{"attribute": "Strahler","count": 8,"type": "mixed","values": ["1","2","3","4","5","6","7","8"]},{"attribute": "TOBAS_ID","count": 211,"type": "mixed","values": ["-888","-999","7010981","7010982","7010983","7010985","7010987","7020081","7020819","7020821","7020823","7020825","7020827","7020831","7020834","7020835","7020837","7020839","7020841","7020843","7020845","7020848","7020851","7020853","7020859","7020871","7020873","7020875","7020877","7020879","7020881","7020883","7020885","7020887","7020893","7020895","7020897","7028911","7030003","7030029","7030211","7030213","7030215","7030217","7030218","7030219","7030220","7030221","7030223","7030227","7030239","7030241","7030242","7030243","7030245","7030247","7030251","7030253","7030255","7030256","7030262","7030267","7030269","7030272","7030281","7030286","7030288","7030295","7030297","7032961","7040361","7040362","7040364","7040367","7040369","7050069","7050611","7050613","7050614","7050615","7050617","7050618","7050619","7050623","7050624","7050625","7050627","7050641","7050643","7050645","7050655","7050661","7050663","7050665","7050667","7050670","7050682","7050683","7050685","7050687"]},{"attribute": "TO_NODE","count": 1000,"type": "mixed","values": ["1","100","1000","100002","100007","100008","100009","10001","100012","100013","100015","100016","100017","100018","10002","100022","100023","100024","100025","100026","100027","100028","100030","100031","100034","100035","100037","100039","100043","100045","100046","100047","100049","100050","100051","100052","100055","100056","100057","100060","100061","100062","100064","100067","100068","100069","100070","100071","100072","100073","100074","100076","100077","10008","100081","100082","100084","100087","100088","10009","100090","100093","100094","100097","10010","100100","100101","100102","100103","100105","100107","100108","100110","100111","100112","100118","100120","100121","100122","100124","100126","100127","100129","100134","100136","100140","100142","100143","100146","100148","100151","100152","100154","100155","100156","100158","100159","10016","100162","100163"]}]}]}}
compression = gzip
root dir tiles: 63
leaf directories: 0So let’s see
library(leaflet)
library(leafem)
nzb_pmt_url = "https://vector-tiles-data.s3.eu-central-1.amazonaws.com/rivers_africa.pmtiles"
leaflet() |>
addTiles(group = "osm") |>
addProviderTiles(
"Esri.WorldImagery"
, group = "esri"
) |>
addPMPolylines(
url = nzb_pmt_url
, layerId = "rivers_pmt"
, group = "rivers_pmt"
, style = paintRules(
layer = "rivers_africa"
, color = "blue"
, width = 2
)
) |>
addMouseCoordinates() |>
addLayersControl(
baseGroups = c("osm", "esri")
, overlayGroups = c("rivers_fgb")
) |>
setView(0, 0, 2)And now together
leaflet(height = "800px", width = "100%") |>
addTiles(group = "osm") |>
addProviderTiles(
"Esri.WorldImagery"
, group = "esri"
) |>
addPMPolylines(
url = nzb_pmt_url
, layerId = "rivers_pmt2"
, group = "rivers_pmt2"
, style = paintRules(
layer = "rivers_africa"
, color = "blue"
, width = 2
)
) |>
addFgb(
url = nzb_fgb_url
, layerId = "rivers_fgb2"
, group = "rivers_fgb2"
, color = "black"
, popup = TRUE
, minZoom = 8
, weight = 4
) |>
addMouseCoordinates() |>
addLayersControl(
baseGroups = c("osm", "esri")
, overlayGroups = c("rivers_pmt2", "rivers_fgb2")
) |>
setView(0, 0, 2)So, what else is in the bucket? There’s a files called depoints.pmtiles which is all points contained in the Geofabrik OSM snapshot that I downloaded in mid August 2022. This files was created with thefollowing commands:
ogr2ogr \
-f FlatGeoBuf \ # output file format
depoints.fgb \ # output file name
germany-latest.osm.pbf \ # input file name
points # layer name
tippecanoe \
-z 14 \ # maximum zoom level
--projection=EPSG:4326 \ # crs definition
-x other_tags \ # drop feature column called "other tags"
--cluster-densest-as-needed \ # if too many points in tile, cluster
--extend-zooms-if-still-dropping \ # if still clustering, keep increasing max zoom
-o depoints.mbtiles \ # output file name
-l depoints \ # layer name to be created
depoints.fgb # input file name
pmtiles-convert
--gzip \ # compress the tiles using gzip
depoints.mbtiles \ # input file name
depoints.pmtiles # output file nameThe important bits here are in the tippecanoe call:
-x flag we drop one of the feature columns called “other_tags”--cluster-densest-as-needed to make sure our tiles won’t get too big at low zoom levels--extend-zoom-levels-if-still-dropping (actually still-clustering in this case) to make sure we provide all data at lowest zoom levelsThis setup should allow us to make sure we have performant tiles at all zoom levels. Clustered, i.e. thinned out, at low zoom levels and all points at high zoom levels.
So how many points are in the file?
ogrinfo -so -al /home/tim/tappelhans/privat/data/osm_germany/depoints.fgbINFO: Open of `/home/tim/tappelhans/privat/data/osm_germany/depoints.fgb'
using driver `FlatGeobuf' successful.
Layer name: points
Geometry: Point
Feature Count: 15372518
Extent: (1.404256, 47.144467) - (22.950335, 60.436330)
Layer SRS WKT:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
osm_id: String (0.0)
name: String (0.0)
barrier: String (0.0)
highway: String (0.0)
ref: String (0.0)
address: String (0.0)
is_in: String (0.0)
place: String (0.0)
man_made: String (0.0)
other_tags: String (0.0)So let’s see how performant this visualisation will be.
url_depoints = "https://vector-tiles-data.s3.eu-central-1.amazonaws.com/depoints.pmtiles"
leaflet(height = "800px", width = "100%") |>
addProviderTiles("CartoDB.Positron", group = "CartoDB.Positron") |>
addProviderTiles("Esri.WorldImagery", group = "Esri.WorldImagery") |>
addPMPoints(
url = url_depoints
, layerId = "depoints"
, group = "depoints"
, style = paintRules(
layer = "depoints"
, fillColor = "#295fcc"
, stroke = "black"
, radius = 2
)
) |>
addMouseCoordinates() |>
setView(10, 50, zoom = 4) |>
addLayersControl(
baseGroups = c(
"CartoDB.Positron"
, "Esri.WorldImagery"
)
, overlayGroups = c("depoints")
)Given that we outsource data storage to the cloud and have performant data access when rendering, we can create really massive data visualisations now. Here’s all the vector data we have, all put together on one map.
Yes, we can also mix raster and vector data.
natearth_url = "https://raster-tiles-data.s3.eu-central-1.amazonaws.com/natearth_3857_cog_defl.tif"
leaflet() |>
addTiles(group = "osm") |>
addCOG(
url = natearth_url
, group = "natearth"
, layerId = "natearth"
, opacity = 1
, resolution = 256
, autozoom = FALSE
) |>
addPMPoints(
url = url_depoints
, layerId = "depoints"
, group = "depoints"
, style = paintRules(
layer = "depoints"
, fillColor = "black"
, stroke = "transparent"
, radius = 2
)
) |>
addMouseCoordinates() |>
addLayersControl(
baseGroups = c("osm")
, overlayGroups = c("natearth", "depoints")
) |>
setView(0, 0, 2)