Mapping Large Data

Jared P. Lander

Lander Analytics

What is Big?

Files

locations <- readr::read_csv(here::here('data', 'locations_full.csv')) |> 
    sf::st_as_sf(coords=c('longitude', 'latitude'), remove=FALSE, crs = sf::st_crs(4326))
locations
## Simple feature collection with 325600 features and 8 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -19.17664 ymin: 23.94552 xmax: 46.06759 ymax: 59.4917
## Geodetic CRS:  WGS 84
## # A tibble: 325,600 × 9
##    longitude latitude timestamp           id     value1 value2 size   count
##  *     <dbl>    <dbl> <dttm>              <chr>   <dbl>  <dbl> <chr>  <dbl>
##  1   -0.0558     51.3 2022-09-16 23:57:43 0a65x6   20.6  1892. medium  6578
##  2   -1.32       51.8 2022-10-03 19:57:26 bs7c6b   28.5  1677. medium  5182
##  3   -1.35       51.5 2022-12-05 02:45:18 o8s1vn   27.9  1818. small   3843
##  4   -0.248      52.7 2023-01-08 20:35:19 j9ge45   27.6  1884. small   6909
##  5   -1.47       51.2 2023-01-12 04:56:20 3qj7wc   28.4  1447. large   3996
##  6    2.28       50.7 2022-08-06 21:48:26 gszw57   27.3  1535. large   4651
##  7   -0.838      50.6 2022-10-09 18:55:31 cnb96e   27.4  1544. large   5798
##  8    1.88       51.3 2022-06-14 23:54:40 v87j5b   25.3  1936. small   7620
##  9    0.148      51.7 2022-07-24 13:06:09 qrq4sr   25.9  1844. large   1832
## 10   -0.279      50.9 2022-12-08 18:11:21 5bzw3p   21.8  1935. large   3171
## # ℹ 325,590 more rows
## # ℹ 1 more variable: geometry <POINT [°]>
locations |> nrow() |> scales::comma()
## [1] "325,600"

longitude,latitude,timestamp,id,value1,value2,size,count\n-0.05577217978412706,51.260731483483646,2022-09-16T23:57:43Z,0a65x6,20.62621257932002,1891.5277283524965,medium,6578\n-1.3193385232774468,51.79962207361285,2022-10-03T19:57:26Z,bs7c6b,28.471566960996277,1676.5703173405125,medium,5182\n-1.3539010673951777,51.52585002509342,2022-12-05T02:45:18Z,o8s1vn,27.88808188675451,1817.9992027257351,small,3843\n-0.24789018310844563,52.73734346769109,2023-01-08T20:35:19Z,j9ge45,27.567753958960875,1883.6389906329923,small,6909\n-1.4727693445858168,51.163188062999986,2023-01-12T04:56:20Z,3qj7wc,28.39833723924046,1446.5678160872817,large,3996

{ "type": "FeatureCollection", "name": "locations", "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } }, "features": [ { "type": "Feature", "properties": { "longitude": -0.05577217978412706, "latitude": 51.260731483483646, "timestamp": "2022-09-16T23:57:43Z", "id": "0a65x6", "value1": 20.626212579320018, "value2": 1891.5277283524963, "size": "medium", "count": 6578.0 }, "geometry": { "type": "Point", "coordinates": [ -0.055772179784127, 51.260731483483646 ] } }, { "type": "Feature", "properties": { "longitude": -1.3193385232774468, "latitude": 51.799622073612852, "timestamp": "2022-10-03T19:57:26Z", "id": "bs7c6b", "value1": 28.471566960996277, "value2": 1676.5703173405125, "size": "medium", "count": 5182.0 }, "geometry": { "type": "Point", "coordinates": [ -1.319338523277447, 51.799622073612852 ] } }, { "type": "Feature", "properties": { "longitude": -1.3539010673951777, "latitude": 51.525850025093419, "timestamp": "2022-12-05T02:45:18Z", "id": "o8s1vn", "value1": 27.888081886754509, "value2": 1817.9992027257351, "size": "small", "count": 3843.0 }, "geometry": { "type": "Point", "coordinates": [ -1.353901067395178, 51.525850025093419 ] } }, { "type": "Feature", "properties": { "longitude": -0.24789018310844563, "latitude": 52.73734346769109, "timestamp": "2023-01-08T20:35:19Z", "id": "j9ge45", "value1": 27.567753958960875, "value2": 1883.6389906329923, "size": "small", "count": 6909.0 }, "geometry": { "type": "Point", "coordinates": [ -0.247890183108446, 52.73734346769109 ] } }, { "type": "Feature", "properties": { "longitude": -1.4727693445858168, "latitude": 51.163188062999986, "timestamp": "2023-01-12T04:56:20Z", "id": "3qj7wc", "value1": 28.398337239240458, "value2": 1446.5678160872817, "size": "large", "count": 3996.0 }, "geometry": { "type": "Point", "coordinates": [ -1.472769344585817, 51.163188062999986 ] } }, { "type": "Feature", "properties": { "longitude": 2.2792216108954779, "latitude": 50.666979435294792, "timestamp": "2022-08-06T21:48:26Z", "id": "gszw57", "value1": 27.264727010027567, "value2": 1534.6726730410919, "size": "large", "count": 4651.0 }, "geometry": { "type": "Point", "coordinates": [ 2.279221610895478, 50.666979435294792 ] } }, { "type": "Feature", "properties": { "longitude": -0.8375930520472068, "latitude": 50.552178039176447, "timestamp": "2022-10-09T18:55:31Z", "id": "cnb96e", "value1": 27.402930506049032, "value2": 1544.3676924122717, "size": "large", "count": 5798.0 }, "geometry": { "type": "Point", "coordinates": [ -0.837593052047207, 50.552178039176447 ] } },

<?xml version="1.0" encoding="utf-8" ?> <kml xmlns="http://www.opengis.net/kml/2.2"> <Document id="root_doc"> <Schema name="locations_full" id="locations_full">    <SimpleField name="longitude" type="float"></SimpleField>   <SimpleField name="latitude" type="float"></SimpleField>    <SimpleField name="timestamp" type="string"></SimpleField>  <SimpleField name="id" type="string"></SimpleField>     <SimpleField name="value1" type="float"></SimpleField>  <SimpleField name="value2" type="float"></SimpleField>  <SimpleField name="size" type="string"></SimpleField>   <SimpleField name="count" type="float"></SimpleField>

fs::file_size(here::here("data", glue::glue("locations_full.{c('csv', 'geojson', 'kml')}")))
##  35.3M 107M 144.1M

Geoparquet

fs::file_size(here::here("data", glue::glue("locations_full.{c('csv', 'geojson', 'kml', 'parquet')}")))
##  35.3M 107M 144.1M  16M

Create the Files

locations |> sf::write_sf('data/locs.csv')
locations |> sf::write_sf('data/locs.geojson')
locations |> sf::write_sf('data/locs.kml')
locations |> geoarrow::write_geoparquet('data/locs.paqquet')

Points on a Map

library(leaflet)
set.seed(2861)

leaflet(locations |> dplyr::slice_sample(n=4000), elementId='points-sample') |> 
    addTiles() |> 
    addCircles()

Polygons

countries <- rnaturalearth::countries110 |> 
    sf::st_as_sf() |> 
    tibble::as_tibble() |> 
    sf::st_as_sf()

head(countries)
## Simple feature collection with 6 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 169
##   featurecla    scalerank LABELRANK SOVEREIGNT SOV_A3 ADM0_DIF LEVEL TYPE  TLC  
##   <chr>             <int>     <int> <chr>      <chr>     <int> <int> <chr> <chr>
## 1 Admin-0 coun…         1         6 Fiji       FJI           0     2 Sove… 1    
## 2 Admin-0 coun…         1         3 United Re… TZA           0     2 Sove… 1    
## 3 Admin-0 coun…         1         7 Western S… SAH           0     2 Inde… 1    
## 4 Admin-0 coun…         1         2 Canada     CAN           0     2 Sove… 1    
## 5 Admin-0 coun…         1         2 United St… US1           1     2 Coun… 1    
## 6 Admin-0 coun…         1         3 Kazakhstan KA1           1     1 Sove… 1    
## # ℹ 160 more variables: ADMIN <chr>, ADM0_A3 <chr>, GEOU_DIF <int>,
## #   GEOUNIT <chr>, GU_A3 <chr>, SU_DIF <int>, SUBUNIT <chr>, SU_A3 <chr>,
## #   BRK_DIFF <int>, NAME <chr>, NAME_LONG <chr>, BRK_A3 <chr>, BRK_NAME <chr>,
## #   BRK_GROUP <chr>, ABBREV <chr>, POSTAL <chr>, FORMAL_EN <chr>,
## #   FORMAL_FR <chr>, NAME_CIAWF <chr>, NOTE_ADM0 <chr>, NOTE_BRK <chr>,
## #   NAME_SORT <chr>, NAME_ALT <chr>, MAPCOLOR7 <int>, MAPCOLOR8 <int>,
## #   MAPCOLOR9 <int>, MAPCOLOR13 <int>, POP_EST <dbl>, POP_RANK <int>, …

leaflet(elementId='countries') |> 
    addTiles() |> 
    addPolygons(data=countries |> dplyr::slice_head(n=50)) |> 
    setView(0, 0, 2)

countries |> 
    sf::st_cast("POINT") |> 
    nrow()
## [1] 10654

{leafgl}

library(leafgl)

leaflet(elementId='points-gl') |> 
    addTiles() |> 
    addGlPoints(locations |> dplyr::slice_sample(n=250000))

leaflet(elementId='countries-gl') |> 
    addTiles() |> 
    addGlPolygons(data=countries |> sf::st_cast("POLYGON")) |> 
    leaflet::setView(0, 0, 2)

leaflet(elementId='countries-gl-color') |> 
    addTiles() |> 
    addGlPolygons(data=countries |> sf::st_cast("POLYGON"), fillColor='BRK_A3') |> 
    setView(0, 0, 2)

Hex Bins

library(h3jsr)

binned <- locations |> 
    point_to_cell(res=5, simple=FALSE) |> 
    dplyr::count(hex=h3_resolution_5, sort=TRUE) |> 
    cell_to_polygon(simple=FALSE)

binned
## Simple feature collection with 17630 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -19.27516 ymin: 23.80341 xmax: 46.1625 ymax: 59.59969
## Geodetic CRS:  WGS 84
## First 10 features:
##       n      h3_address                       geometry
## 1  1772 85182c17fffffff POLYGON ((-6.270301 53.3549...
## 2  1377 85182c07fffffff POLYGON ((-6.342615 53.2170...
## 3  1334 85182c8bfffffff POLYGON ((-6.368699 53.6089...
## 4  1313 85182c13fffffff POLYGON ((-6.441145 53.4714...
## 5  1023 85182c8ffffffff POLYGON ((-6.197661 53.4924...
## 6   667 85182c03fffffff POLYGON ((-6.513264 53.3335...
## 7   659 85182c9bfffffff POLYGON ((-6.295924 53.7461...
## 8   559 85182c3bfffffff POLYGON ((-6.414605 53.0788...
## 9   317 85182c83fffffff POLYGON ((-6.124691 53.6297...
## 10  306 85182c33fffffff POLYGON ((-6.17286 53.10026...

leaflet(elementId='binned') |> 
    addTiles() |> 
    addGlPolygons(data=binned, fillColor=~log(n))

Different Layers

map_layers <- leaflet(elementId='binned') |> 
    addTiles() |> 
    addGlPolygons(data=binned, fillColor=~log(n), group='Bins') |> 
    addGlPoints(locations |> dplyr::slice_sample(n=250000), group='Points') |> 
    groupOptions("Bins", zoomLevels = 1:3) |> 
    groupOptions("Points", zoomLevels = 4:15) |> 
    leafem::addMouseCoordinates()



PMTiles

tippecanoe \
    -z 15 \
    --projection=EPSG:4326 \
    --cluster-densest-as-needed --extend-zooms-if-still-dropping \
    -o locations.pmtiles \
    --force -L spots:data/locations_full.geojson

From Protomaps Documentation:

Nginx supports HTTP Range Requests. CORS headers should be set by configuration files.

library(leafem)

leaflet() |> addTiles() |> 
    addPMPoints(
        url="https://vector-tiles-data.s3.eu-central-1.amazonaws.com/depoints.pmtiles"
        , layerId="depoints", group="depoints"
        , style=paintRules(layer="depoints", fillColor="black", stroke="white", radius=4)
    ) |> 
    setView(10, 51, 6)

Wrapping Up

File Types

{leafgl}

Binning

PMTiles

Thank You