Creating cloud optimised data

Intro

The future of data analysis likely lies in the cloud. Meaning, most data and even maybe analysis tools will not live on laptops or desktops anymore, but will reside in centralised data stores and analysis computing infrastructure. As such, it becomes more and more important to enable remote access to the data and have tools that can visualise that data locally without having to copy the whole data first.

In the last decade or so, there has been much development in creating tools that optimise data structures so that they can be streamed from a remote location rather than having to download the whole thing. So called https range requests have become the de facto standard for partial data streaming. This is also the technology used when we watch Netflix or listen to Spotify, we can start watching/listening instantly as we don’t need to download everything first, but get served what we are seeing/hearing chunk after chunk. Geospatial cloud optimised data gets served a the similar manner. Only parts of the data that we are currently interested in are being served from the remote location. However, while movie or music data ususally get streamed linearly, i.e. from beginning to end, geospatial data will be served with regard to the spatial extent we are interested in. These geospatial chunks are generally referred to as tiles. Hence, cloud optimised data needs some sort of spatial index so that spatial queries are fast enough and the data lookup does not hinder streaming. Spatial indexes are a topic of their own and will not be covered in this document.

This document serves as a compilation of possible approaches to create and store cloud optimised geospatial data.

Software

In order to create cloud optimised data sets, we will need a few different tools (R is actually not one of them). As always, we can (and should) make a distinction between raster data and vector data.

Raster data

In order to create cloud optimised tiled raster data sets, the ever-powerful GDAL software is what we will use. There are many installation instructions on the wild wild web for GDAL, but the R inclined user is referred to the sf package README. GDAL has a very powerful command line interface (CLI) that we can use to verify that we have a working GDAL installation on our system.

gdalinfo --version

output

GDAL 3.4.3, released 2022/04/22

Vector data

Creating vector tiles is a bit more work as the structure of the data is not quite as organised as for regular raster data. Chunking raster data to created tiles is easy as it is made up of squares (pixels) so cutting out some of them is straight forward. With vector data we generally follow the same approach for creating vector tiles. In addition to GDAL we need some specialised tiling and indexing software to create cloud optimised vector data sets, namely tippecanoe and PMTiles. tippecanoe was created by MapBox, has now found a new home at @protomaps. PMTiles is also part of the protomaps suite of software.

Tippencanoe

To install tippecanoe we can follow these installation instructions

git clone https://github.com/protomaps/tippecanoe.git
cd tippecanoe
make -j
make install

To verify we have a working installation of tippecanoe

tippecanoe --version

output

tippecanoe v2.1.0

PMTiles

Installation instructions for PMTiles can be found here. I’ve chosen the python version of the program, as I had a working python installation, so installing that was easiest:

pip install pmtiles

To verify the installation was successful:

pmtiles-convert -h

output

usage: pmtiles-convert [-h] [--maxzoom MAXZOOM] [--gzip] [--overwrite]
                       input output

Convert between PMTiles and other archive formats.

positional arguments:
  input              Input .mbtiles or .pmtiles
  output             Output .mbtiles, .pmtiles, or directory

optional arguments:
  -h, --help         show this help message and exit
  --maxzoom MAXZOOM  the maximum zoom level to include in the output.
  --gzip             The output should be gzip-compressed.
  --overwrite        Overwrite the existing output.

That’s all the software we need to create our cloud optimised data. So let’s see how.

Creation

Cloud optimised raster data

GDAL has a dedicated driver to create cloud optimised geotiffs (COGs). But to get a working COG, there’s one more step we need. Assuming we have a geotiff that we want to make cloud optimised, i.e. turn into a COG, it is not enough to simply gdal_translate it. As a first step we need to create overviews of the data using gdaladdo.

gdaladdo \
  -r average \         # resampling method
  -b 1 -b 2 -b 3 \     # bands to cerate overviews for
  natearth_3857.tif    # file name

Here, we create overviews for bands 1, 2 & 3 using resampling method average. This is what the resulting file looks like (see the Band information towards the end of the output):

gdalinfo /home/tim/tappelhans/privat/data/HYP_HR_SR_OB_DR/natearth_3857.tif

output

Driver: GTiff/GeoTIFF
Files: /home/tim/tappelhans/privat/data/HYP_HR_SR_OB_DR/natearth_3857.tif
       /home/tim/tappelhans/privat/data/HYP_HR_SR_OB_DR/natearth_3857.aux
Size is 4007, 48507
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["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]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06°S and 85.06°N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-20037508.342789243906736,242528680.943742722272873)
Pixel Size = (10000.000000000000000,-10000.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2012:07:16 09:21:17
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS5 Macintosh
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-20037508.343,242528680.944) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-20037508.343,-242541319.056) (180d 0' 0.00"W, 90d 0' 0.00"S)
Upper Right (20032491.657,242528680.944) (179d57'17.76"E, 90d 0' 0.00"N)
Lower Right (20032491.657,-242541319.056) (179d57'17.76"E, 90d 0' 0.00"S)
Center      (   -2508.343,   -6319.056) (  0d 1'21.12"W,  0d 3'24.35"S)
Band 1 Block=4007x1 Type=Byte, ColorInterp=Red
  Description = Layer_1
  Overviews: 2004x24254, 1002x12127, 668x8085, 501x6064, 251x3032, 126x1516, 63x758, 32x379, 16x190
  Metadata:
    LAYER_TYPE=athematic

Band 2 Block=4007x1 Type=Byte, ColorInterp=Green
  Description = Layer_2
  Overviews: 2004x24254, 1002x12127, 668x8085, 501x6064, 251x3032, 126x1516, 63x758, 32x379, 16x190
  Metadata:
    LAYER_TYPE=athematic

Band 3 Block=4007x1 Type=Byte, ColorInterp=Blue
  Description = Layer_3
  Overviews: 2004x24254, 1002x12127, 668x8085, 501x6064, 251x3032, 126x1516, 63x758, 32x379, 16x190
  Metadata:
    LAYER_TYPE=athematic

Now that we have thos overviews, we can gdal_translate this file to a COG. Note that the file ending of COGs is also .tif.

gdal_translate \
  natearth_3857.tif \            # input file
  natearth_3857_cog_defl.tif \   # output file
  -of COG \                      # output format (driver)
  -co COMPRESS=DEFLATE           # use DEFLATE compression

And the result:

gdalinfo /home/tim/tappelhans/privat/data/HYP_HR_SR_OB_DR/natearth_3857_cog_defl.tif

output

Driver: GTiff/GeoTIFF
Files: /home/tim/tappelhans/privat/data/HYP_HR_SR_OB_DR/natearth_3857_cog_defl.tif
Size is 4007, 48507
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["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]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06°S and 85.06°N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-20037508.342789243906736,242528680.943742722272873)
Pixel Size = (10000.000000000000000,-10000.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2012:07:16 09:21:17
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS5 Macintosh
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
  LAYOUT=COG
Corner Coordinates:
Upper Left  (-20037508.343,242528680.944) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-20037508.343,-242541319.056) (180d 0' 0.00"W, 90d 0' 0.00"S)
Upper Right (20032491.657,242528680.944) (179d57'17.76"E, 90d 0' 0.00"N)
Lower Right (20032491.657,-242541319.056) (179d57'17.76"E, 90d 0' 0.00"S)
Center      (   -2508.343,   -6319.056) (  0d 1'21.12"W,  0d 3'24.35"S)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  Description = Layer_1
  Overviews: 2004x24254, 1002x12127, 668x8085, 501x6064, 251x3032, 126x1516, 63x758, 32x379, 16x190
  Metadata:
    LAYER_TYPE=athematic
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  Description = Layer_2
  Overviews: 2004x24254, 1002x12127, 668x8085, 501x6064, 251x3032, 126x1516, 63x758, 32x379, 16x190
  Metadata:
    LAYER_TYPE=athematic
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  Description = Layer_3
  Overviews: 2004x24254, 1002x12127, 668x8085, 501x6064, 251x3032, 126x1516, 63x758, 32x379, 16x190
  Metadata:
    LAYER_TYPE=athematic

Not much different, but in the section called Image Structure Metadata: we see a Layout=COG. That is what we want and we can now use our COG. In the case of this file, it was uploaded to a public AWS S3 bucket located at https://raster-tiles-data.s3.eu-central-1.amazonaws.com/natearth_3857_cog_defl.tif that we can now visualise on a leaflet map using leafem::adCOG()

library(leaflet)
library(leafem)

url = "https://raster-tiles-data.s3.eu-central-1.amazonaws.com/natearth_3857_cog_defl.tif"

leaflet() |>
  addTiles() |>
  addCOG(
    url = url
    , group = "COG"
    , resolution = 192
    , autozoom = FALSE
  ) |>
  setView(0, 0, 2)

Cloud optimised vector tiles

In order to create vector tiles we need to ensure 2 things:

1. For use with leaflet, our data needs to be in geographic coordinates (EPSG:4326)
2. For tippecanoe we need our data to be either in .geojson or .fgb (FlatGeobuf) format. I recommend the latter.

The workflow of creating cloud optimised vector tiles will be shown using New Zealnd building footprint data sourced from the LINZ Data Service and licensed for reuse under CC BY 4.0. The data was downloaded from the link aboce as a GeoPackage (.gpkg).

gdalsrsinfo /home/tim/tappelhans/privat/data/lds-nz-building-outlines/nz-building-outlines.gpkg

output

PROJ.4 : +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

OGC WKT2:2018 :
PROJCRS["NZGD2000 / New Zealand Transverse Mercator 2000",
    BASEGEOGCRS["NZGD2000",
        DATUM["New Zealand Geodetic Datum 2000",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4167]],
    CONVERSION["New Zealand Transverse Mercator 2000",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",173,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",1600000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",10000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["New Zealand - North Island, South Island, Stewart Island - onshore."],
        BBOX[-47.33,166.37,-34.1,178.63]],
    ID["EPSG",2193]]

So, those data come projected in a CRS called “NZGD2000”. Hence, we need to transform them to EPSG:4326. In the same command we will also translate them to the FlatGeobuf format.

ogr2ogr \
  -t_srs EPSG:4326 \            # target coordinate reference system
  nz-building-outlines.fgb \    # output file
  nz-building-outlines.gpkg     # input file

Then, let’s use GDAL to get some information about that new .fgb file.

ogrinfo \
  -so \
  -al \
  /home/tim/tappelhans/privat/data/lds-nz-building-outlines/nz-building-outlines.fgb

output

INFO: Open of `/home/tim/tappelhans/privat/data/lds-nz-building-outlines/nz-building-outlines.fgb'
      using driver `FlatGeobuf' successful.

Layer name: nz_building_outlines
Geometry: Multi Polygon
Feature Count: 3320498
Extent: (167.398488, -46.897633) - (178.545566, -34.426026)
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
building_id: Integer (0.0)
name: String (250.0)
use: String (40.0)
suburb_locality: String (80.0)
town_city: String (80.0)
territorial_authority: String (100.0)
capture_method: String (40.0)
capture_source_group: String (80.0)
capture_source_id: Integer (0.0)
capture_source_name: String (100.0)
capture_source_from: DateTime (0.0)
capture_source_to: DateTime (0.0)
last_modified: DateTime (0.0)

So, now we have about 3.3M polygons projected in geographic coordinates.

In order to convert those to vector tiles, we’re going to use tippecanoe. There is a very informative blog post by Miles McBain on creating and using vector tiles in R with a very detailed tutorial on creating those tiles using tippecanoe. tippecanoe has a lot of options that can be set to create vector tiles. A detailed description of those can be found at the protomaps repository on github.

In the case of the New Zealand building footprints we used the following settings.

tippecanoe \
  -zg \                              # guess the zoom levels
  --projection=EPSG:4326 \           # the projection
  --no-tile-compression \            # do not compress the tiles
  --no-feature-limit \               # no feature limit for tiles
  --no-tile-size-limit \             # no size limit for tiles
  -l nz-building-outlines \          # layer name
  -o nz-building-outlines.mbtiles \  # output file
  nz-building-outlines.fgb           # input file

In theory, the created .mbtiles file can already be used for hosting the data in the cloud. Please refer the above mentioned blog post by Miles McBain to learn how to use {rdeck} package to visualise those. In our case, we will apply one more conversion to PMTiles which is optimised for rendering in the browser. Plus, it has the advantage that the data can be compressed, saving quite some space and hence money in our online bucket.

pmtiles-convert \
  --gzip \                          # gzip compress
  nz-building-outlines.mbtiles \    # input file
  nz-building-outlines.pmtiles      # output file

So now we have our final cloud optimised vector tile data that we can now upload to the cloud and render from there. Here’s what is in the data.

pmtiles-show /home/tim/tappelhans/privat/data/lds-nz-building-outlines/nz-building-outlines.pmtiles

output

spec version:  2
metadata:
name = nz-building-outlines.mbtiles
description = nz-building-outlines.mbtiles
version = 2
minzoom = 0
maxzoom = 14
center = 172.606201,-43.556510,14
bounds = 167.398488,-46.897633,178.545565,-34.426026
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 nz-building-outlines.mbtiles -l nz-building-outlines nz-building-outlines.fgb
json = {"vector_layers": [ { "id": "nz-building-outlines", "description": "", "minzoom": 0, "maxzoom": 14, "fields": {"building_id": "Mixed", "capture_method": "String", "capture_source_from": "String", "capture_source_group": "String", "capture_source_id": "Mixed", "capture_source_name": "String", "capture_source_to": "String", "last_modified": "String", "name": "String", "suburb_locality": "String", "territorial_authority": "String", "town_city": "String", "use": "String"} } ],"tilestats": {"layerCount": 1,"layers": [{"layer": "nz-building-outlines","count": 3320498,"geometry": "Polygon","attributeCount": 13,"attributes": [{"attribute": "building_id","count": 1000,"type": "mixed","values": ["1000000","1000001","1000002","1000003","1000004","1000005","1000006","1000007","1000008","1000009","1000010","1000011","1000014","1000015","1000016","1000017","1000018","1000019","1000020","1000021","1000022","1000023","1000024","1000026","1000027","1000028","1000029","1000030","1000031","1000032","1000033","1000034","1000035","1000036","1000037","1000038","1000039","1000040","1000041","1000043","1000044","1000045","1000046","1000047","1000048","1000049","1000050","1000051","1000052","1000053","1000054","1000055","1000058","1000059","1000060","1000061","1000062","1000063","1000064","1000065","1000066","1000067","1000068","1000069","1000071","1000072","1000073","1000074","1000075","1000076","1000077","1000078","1000079","1000080","1000081","1000082","1000083","1000084","1000085","1000086","1000088","1000090","1000091","1000092","1000093","1000094","1000095","1000097","1000098","1000099","1000100","1000101","1000102","1000103","1000104","1000105","1000106","1000107","1000109","1000110"]},{"attribute": "capture_method","count": 5,"type": "string","values": ["Feature Extraction","Trace Orthophotography","Trace Other Image","Trace Stereophotography","Unknown"]},{"attribute": "capture_source_from","count": 30,"type": "string","values": ["2010-01-01T00:00:00","2012-10-30T00:00:00","2013-10-29T00:00:00","2014-02-05T00:00:00","2014-11-09T00:00:00","2014-12-03T00:00:00","2015-01-02T00:00:00","2015-01-09T00:00:00","2015-02-26T00:00:00","2015-11-16T00:00:00","2015-12-26T00:00:00","2015-12-27T00:00:00","2016-01-07T00:00:00","2016-01-22T00:00:00","2016-03-03T00:00:00","2016-03-11T00:00:00","2016-09-10T00:00:00","2016-10-31T00:00:00","2016-11-09T00:00:00","2016-12-08T00:00:00","2017-01-06T00:00:00","2017-02-08T00:00:00","2017-02-14T00:00:00","2017-03-15T00:00:00","2017-12-02T00:00:00","2017-12-12T00:00:00","2017-12-28T00:00:00","2018-01-14T00:00:00","2019-01-02T00:00:00","2019-01-16T00:00:00"]},{"attribute": "capture_source_group","count": 2,"type": "string","values": ["NZ Aerial Imagery","Various NZ Imagery"]},{"attribute": "capture_source_id","count": 32,"type": "mixed","values": ["0","1000","1006","1013","1014","1016","1018","1028","1029","1032","1033","1034","1039","1042","1043","1052","1058","1073","1074","1075","1086","1091","1106","1110","1111","1112","1115","1145","1150","1151","1154","1158"]},{"attribute": "capture_source_name","count": 32,"type": "string","values": ["Auckland 0.075m Urban Aerial Photos (2017)","Bay of Plenty 0.125m Urban Aerial Photos (2014-2015)","Bay of Plenty 0.25m Rural Aerial Photos (2015-2017)","Bay of Plenty 0.3m Rural Aerial Photos (2016-2017)","Canterbury 0.3m Rural Aerial Photos (2014-2015)","Canterbury 0.3m Rural Aerial Photos (2015-2016)","Canterbury 0.3m Rural Aerial Photos (2017-2018)","Canterbury 0.4m Rural Aerial Photos (2012-2013)","Canterbury 0.4m Rural Aerial Photos (2013-2014)","Christchurch 0.075m Urban Aerial Photos (2015-2016)","Dunedin 0.1m Urban Aerial Photos (2018-2019)","Gisborne 0.3m Rural Aerial Photos (2017-2019)","Hawkes Bay 0.3m Rural Aerial Photos (2014-2015)","Manawatu Whanganui 0.3m Rural Aerial Photos (2015-2016)","Manawatu Whanganui 0.3m Rural Aerial Photos (2016-2017)","Marlborough 0.2m Rural Aerial Photos (2015-2016)","Marlborough 0.3m Rural Aerial Photos (2017-2018)","Nelson 0.3m Rural Aerial Photos (2018-2019)","Northland 0.4m Rural Aerial Photos (2014-2016)","Otago 0.3m Rural Aerial Photos (2017-2019)","Otago 0.4m Rural Aerial Photos (2013-2014)","Southland & Central Otago 0.4m Rural Aerial Photos (2013-2014)","Southland & Central Otago 0.4m Rural Aerial Photos (2015-2017)","Taranaki 0.3m Rural Aerial Photos (2016-2018)","Tasman 0.3m Rural Aerial Photos (2016-2017)","Tasman 0.3m Rural Aerial Photos (2018-2019)","Various NZ Imagery","Waikato 0.3m Rural Aerial Photos (2016-2019)","Waimakariri 0.075m Urban Aerial Photos (2015-2016)","Wellington 0.3m Rural Aerial Photos (2016-2017)","West Coast 0.3m Rural Aerial Photos (2015-2016)","West Coast 0.3m Rural Aerial Photos (2016-2017)"]},{"attribute": "capture_source_to","count": 30,"type": "string","values": ["2013-03-28T00:00:00","2014-03-13T00:00:00","2014-03-26T00:00:00","2014-03-27T00:00:00","2015-01-24T00:00:00","2015-03-01T00:00:00","2015-03-28T00:00:00","2016-01-22T00:00:00","2016-02-20T00:00:00","2016-03-07T00:00:00","2016-04-10T00:00:00","2016-04-21T00:00:00","2016-06-01T00:00:00","2016-07-01T00:00:00","2017-02-05T00:00:00","2017-02-10T00:00:00","2017-02-25T00:00:00","2017-03-08T00:00:00","2017-03-09T00:00:00","2017-03-17T00:00:00","2017-04-08T00:00:00","2017-05-06T00:00:00","2018-03-11T00:00:00","2018-03-30T00:00:00","2019-02-03T00:00:00","2019-03-04T00:00:00","2019-04-07T00:00:00","2019-04-15T00:00:00","2019-10-28T00:00:00","2020-04-24T00:00:00"]},{"attribute": "last_modified","count": 66,"type": "string","values": ["2018-10-17T00:00:00","2018-10-24T00:00:00","2018-10-26T00:00:00","2018-11-08T00:00:00","2018-11-13T00:00:00","2018-11-14T00:00:00","2018-11-22T00:00:00","2018-11-30T00:00:00","2018-12-06T00:00:00","2018-12-14T00:00:00","2019-01-04T00:00:00","2019-01-09T00:00:00","2019-01-15T00:00:00","2019-02-08T00:00:00","2019-02-11T00:00:00","2019-02-14T00:00:00","2019-02-18T00:00:00","2019-03-04T00:00:00","2019-03-15T00:00:00","2019-03-18T00:00:00","2019-03-19T00:00:00","2019-03-20T00:00:00","2019-03-21T00:00:00","2019-03-22T00:00:00","2019-03-25T00:00:00","2019-03-26T00:00:00","2019-03-27T00:00:00","2019-03-28T00:00:00","2019-03-29T00:00:00","2019-04-01T00:00:00","2019-04-02T00:00:00","2019-04-03T00:00:00","2019-04-04T00:00:00","2019-04-05T00:00:00","2019-04-10T00:00:00","2019-04-11T00:00:00","2019-04-17T00:00:00","2019-05-09T00:00:00","2019-08-20T00:00:00","2019-10-16T00:00:00","2019-11-13T00:00:00","2019-12-05T00:00:00","2019-12-20T00:00:00","2020-01-23T00:00:00","2020-02-24T00:00:00","2020-02-25T00:00:00","2020-04-08T00:00:00","2020-04-09T00:00:00","2020-04-10T00:00:00","2020-04-11T00:00:00","2020-04-14T00:00:00","2020-04-15T00:00:00","2020-04-16T00:00:00","2020-04-17T00:00:00","2020-07-01T00:00:00","2020-07-13T00:00:00","2020-07-22T00:00:00","2020-08-10T00:00:00","2020-08-21T00:00:00","2021-03-15T00:00:00","2021-03-30T00:00:00","2021-04-19T00:00:00","2021-05-07T00:00:00","2021-05-24T00:00:00","2021-05-25T00:00:00","2021-05-27T00:00:00"]},{"attribute": "name","count": 1000,"type": "string","values": ["","A1 Student Limited","ACG Parnell College","ACG Strathallan","ACG Sunderland","ACG Tauranga","ADDI Enrichment Academy","AGE School","Abbotsford School","Aberdeen School","Aberfeldy School","Addington Te Kura Taumatua","Adventure School","Ahipara School","Ahititi School","Ahuroa School","Aidanfield Christian School","Aka Aka School","Akaroa Area School","Akina School / Activity Centre","Ako Space","Al-Madinah School","Albany Junior High School","Albany School","Albany Senior High School","Albury School","Alexandra Hospital","Alexandra School","Alfredton School","Alfriston College","Alfriston School","Allandale School","Allenton Fresh","Allenton School","Allenvale Special School and Resource Centre","Amana Christian School","Amberley School","Ambury Park Centre","Amesbury School","Amisfield School","Amuri Area School","Anchorage Park School","Andersons Bay School","Andrew Street Foodcentre","Anglesea Hospital","Ao Tawhiti Unlimited Discovery","Aokautere School","Aoraki Mount Cook School","Aorangi School (Rotorua)","Aorere College","Aotea College","Apanui School","Aparima College","Apiti School","Appleby School","Aquinas College","Arahoe School","Arahunga School","Arakura School","Aranga School","Aranui School (Wanganui)","Arapohue School","Arataki School","Ardgowan School","Ardmore School","Argyll East School","Aria School","Arohanui Hospice","Arohanui Special School","Arohena School","Arowhenua Maori School","Arrowtown School","Arthur Miller School","Arthur Street School","Ascot Community School","Ascot Integrated Hospital","Ashbrook School","Ashburton Borough School","Ashburton Christian School","Ashburton College","Ashburton Hospital","Ashburton Intermediate","Ashburton Netherby School","Ashgrove School","Ashhurst School","Ashley Rakahuri School","Atea College","Auckland City Hospital","Auckland Girls' Grammar School","Auckland Grammar","Auckland International College","Auckland Normal Intermediate","Auckland Point School","Auckland Secondary Schools' Centre","Auckland Senior College","Auckland Seventh-Day Adventist High School","Auckland Spinal Rehabilitation","Auckland Surgical Centre","Auroa School","Aurora College"]},{"attribute": "suburb_locality","count": 1000,"type": "string","values": ["Abbey Caves","Abbotsford","Abel Tasman National Park","Acacia Bay","Acton","Adair","Addington","Admiral Hill","Admiralty Bay","Ahaura","Ahiaruhe","Ahikouka","Ahipara","Ahititi","Ahuriri","Ahuriri Conservation Park","Ahuriri Flat","Aidanfield","Airedale","Aka Aka","Akaroa","Akatarawa","Akatarawa Valley","Akatore","Akina","Akitio","Albany","Albany Heights","Albert Town","Albury","Alexandra","Alford Forest","Alfredton","Alfriston","Algies Bay","Alicetown","Allandale","Allanton","Allenton","Alma","Alton","Amberley","Amodeo Bay","Amuri Plain","Anakiwa","Anakoha","Anaura Bay","Anawhata","Anderson Park","Andersons Bay","Aniseed Valley","Annesbrook","Aokautere","Aongatete","Aoraki Mount Cook National Park","Aorangi","Aorangi Forest Park","Aotea","Aotuhia","Aparima","Apiti","Appleby","Arahura Valley","Aramoana","Aramoho","Aranga","Aranui","Arapawa Island","Arapohue","Arapuni","Ararata","Ararimu","Ararua","Aratapu","Ardgowan","Ardlussa","Ardmore","Argyle Corner","Argyle Hill","Argyll","Aria","Arkles Bay","Army Bay","Arnold Valley","Aro Valley","Arrow Junction","Arrowtown","Arthur's Pass","Arthur's Pass National Park","Arthurs Point","Arundel","Ascot","Ascot Park","Ashburton","Ashburton Forks","Ashburton Lakes","Ashers","Ashhurst","Ashley","Ashley Clinton"]},{"attribute": "territorial_authority","count": 67,"type": "string","values": ["Area Outside Territorial Authority","Ashburton District","Auckland","Buller District","Carterton District","Central Hawke's Bay District","Central Otago District","Christchurch City","Clutha District","Dunedin City","Far North District","Gisborne District","Gore District","Grey District","Hamilton City","Hastings District","Hauraki District","Horowhenua District","Hurunui District","Invercargill City","Kaikoura District","Kaipara District","Kapiti Coast District","Kawerau District","Lower Hutt City","Mackenzie District","Manawatu District","Marlborough District","Masterton District","Matamata-Piako District","Napier City","Nelson City","New Plymouth District","Opotiki District","Otorohanga District","Palmerston North City","Porirua City","Queenstown-Lakes District","Rangitikei District","Rotorua District","Ruapehu District","Selwyn District","South Taranaki District","South Waikato District","South Wairarapa District","Southland District","Stratford District","Tararua District","Tasman District","Taupo District","Tauranga City","Thames-Coromandel District","Timaru District","Upper Hutt City","Waikato District","Waimakariri District","Waimate District","Waipa District","Wairoa District","Waitaki District","Waitomo District","Wellington City","Western Bay of Plenty District","Westland District","Whakatane District","Whanganui District","Whangarei District"]},{"attribute": "town_city","count": 272,"type": "string","values": ["","Akaroa","Alexandra","Algies Bay","Amberley","Arrowtown","Ashburton","Ashhurst","Auckland","Balclutha","Big Omaha","Blenheim","Bluff","Brightwater","Bulls","Burnham","Cable Bay","Cable Bay Nelson District","Cambridge","Carterton","Cheviot","Christchurch","Clarks Beach","Clinton","Clive","Clyde","Coopers Beach","Coromandel","Cromwell","Dannevirke","Darfield","Dargaville","Doyleston","Drury","Dunedin","Edendale","Edgecumbe","Eketahuna","Eltham","Fairlie","Featherston","Feilding","Foxton","Geraldine","Gisborne","Gore","Granity","Greymouth","Greytown","Hamilton","Hampden","Hanmer Springs","Hastings","Haumoana","Havelock","Havelock North","Hawera","Hector","Helensville","Hikurangi","Hokitika","Huntly","Inglewood","Invercargill","Kaiapoi","Kaikohe","Kaikoura","Kairaki","Kaitaia","Kaitangata","Karangahake","Katikati","Kawakawa","Kawerau","Kerikeri","Kingston","Kirwee","Kumeu","Kurow","Lansdowne","Lawrence","Leeston","Levin","Lincoln","Lower Hutt","Lumsden","Lyttelton","Makara Beach","Mamaku","Manaia","Mangakino","Mangawhai","Mangonui","Mapua","Marlborough Sounds","Martinborough","Marton","Masterton","Matamata","Mataura"]},{"attribute": "use","count": 4,"type": "string","values": ["Hospital","School","Supermarket","Unknown"]}]}]}}
compression = gzip
root dir tiles: 6627
leaf directories: 3

This data can now be visualised using leafem::addPMPolygons

library(leaflet)
library(leafem)

url_nzb = "https://vector-tiles-data.s3.eu-central-1.amazonaws.com/nz-building-outlines.pmtiles"

leaflet() |>
  addTiles() |>
  addPMPolygons(
    url = url_nzb
    , layerId = "nzbuildings"
    , group = "nzbuildings"
    , style = paintRules(
      layer = "nz-building-outlines"
      , fillColor = "blue"
      , stroke = "black"
    )
    , attribution = '<a href="https://data.linz.govt.nz">LINZ Data Service licensed for reuse under CC BY 4.0</a>'
  ) |>
  setView(173.89, -40.65, zoom = 6)