Fast aerial images part I - tiles

A WMTS (Web Map Tile Service) journey in the belgian coordinate system

•

We built several interactive visualizations of 3D models, like our Lidar demo. As our clients started using these visualisations more intensively (interactive or static screenshots) we ran into one particular issue: the ground picture is slow to load.

Loading ground image...

(a typical message on our visualization)

Why? Because we use a Web Map Service - WMS for short (in this case from Informatie Vlaanderen) to get the image corresponding to this ground square. WMS is an easy to use standard protocol, but it is computationally intensive: the WMS must crop and resize a tailor-made image for every single request, which may take several seconds.

The crucial part the URL of a WMS request looks like this:

.../wms?BBOX=173005,163450,173165,163610&CRS=EPSG:31370&WIDTH=1000&HEIGHT=1000...

which means: give me a 1000x1000px image filling the box bouded by xmin=173005, xmax=173165, ymin=163450 and ymax=163610 in the Belgian Lambert 72 coordinate system (EPSG:31370). Let’s time such a request:

Send WMS map request!

source: Informatie Vlaanderen

Tiles to the rescue

The most common workaround for WMS inconveniences is to provide imagery using tiles, and specifically through a Web Map Tile Service (WMTS): instead of allowing custom requests of any rectangle with any resolution, the web client is responsible for fetching 256x256px tiles along a fixed hierarchical grid. Instead of supplying a bounding box and a resolution like above, a WMTS url needs to specify a zoom level z and indices to identify the tile x and y.

source: qgis.org

Among the upsides of tiles are parallelization, pre-computability, and progressive access (depending on user interaction). The main downside of tiles is an additional complexity on the frontend: one needs to transform bounding boxes and desired resolution into a set of tiles indexed by x, y and z, and handle multiple tiles.

A special case: custom belgian coordinate system

Fortunately, this complexity of dealing with tiles is handled transparently by modern map frameworks (Leaflet, Open Layers, Google Maps, …) and tools (QGIS, …) However, we deviate from the standard toolset in two ways: we use the Belgian Lambert 72 coordinate system and we have a tree.js-based 3D visualization, which doesn’t handle WMTS out of the box. As a result, we rolled up our sleeves and delved into the WMTS specs to figure out how to go from a set of Lambert 72 coordinates to the correct tile indices.

The remainder of this section is fairly technical, in the hope it can save other people's time. Feel free to skip to the next section.

As required by the WMTS standard, the WMTS server publishes metadata, among which we find:

<TileMatrixSet>
    <ows:Identifier>BPL72VL</ows:Identifier>
    ...
    <ows:SupportedCRS>EPSG:31370</ows:SupportedCRS>
    <TileMatrix>
        ...
        <ScaleDenominator>3657142.85714285727590322495</ScaleDenominator>
        <TopLeftCorner>9928.000000 329072.000000</TopLeftCorner>
        ...
    </TileMatrix>
...
</TileMatrixSet>

Great, we have x_left=9928, y_top=329072 and a ScaleDenominator=3657142.85. We fed these in a fairly straightforward transformation

\[x_\text{tile} = \left \lfloor \frac{x_{\text{lambert}} - x_\text{left}}{\text{ScaleDenominator}} \cdot 2^\text{zoom} \right \rfloor\] \[y_\text{tile} = \left \lfloor \frac{y_\text{top} - y_\text{lambert}}{\text{ScaleDenominator}} \cdot 2^\text{zoom} \right \rfloor\]

but this doesn’t work.

Error

<Exception code="InvalidParameterValue" locator="tileCol" />

It seemed as if ScaleDenominator was off by a factor more than 10. After quite some digging, it turned out that ScaleDenominator is not expressed in meters (Lambert units), but in pixels at the highest zoom level. However, we did not have the information how big pixels were supposed to be.

Another look at the WMTS metadata revealed another piece of information:

<TileMatrixSet>
    <ows:Identifier>BPL72VL</ows:Identifier>
    <ows:BoundingBox crs="urn:ogc:def:crs:EPSG:6.3:31370">
        <ows:LowerCorner>9928.000000 66928.000000</ows:LowerCorner>
        <ows:UpperCorner>272072.000000 329072.000000</ows:UpperCorner>
    </ows:BoundingBox>
    <ows:SupportedCRS>EPSG:31370</ows:SupportedCRS>
    ...
</TileMatrixSet>

So x_left=9928, x_right=272072, y_bottom=66928, y_top=329072. This bounding box is square, so it looks like the coordinates of the zoom 0 tile. We revised our formula accordingly.

\[x_\text{tile} = \left \lfloor \frac{x_{\text{lambert}} - x_\text{left}}{x_\text{right} - x_\text{left}} \cdot 2^\text{zoom} \right \rfloor\] \[y_\text{tile} = \left \lfloor \frac{y_\text{top} - y_\text{lambert}}{y_\text{top} - y_\text{bottom}} \cdot 2^\text{zoom} \right \rfloor\]

and bingo we got WMS and WMTS to match!

Testing tile speed

Here’s a rendering of the same area as the WMS above using WMTS.

Send WMTS map request!

source: Informatie Vlaanderen

Conclusion

We managed to replace our usage of a WMS to that of a more efficient tile service for our custom 3D visualization. We found that for the large scale aerial imagery of Flanders, pixels at the highest zoom level (14) have a side lenght of \(\frac{y_\text{top} - y_\text{bottom}}{2^{14} \cdot \text{tile size}} = \frac{329072 - 66928}{2^{14} \cdot 256} = 1/16\text{m} =\)6.25cm.

Stay tuned for part II on speeding up imagery through Cloud Optimized Geotiffs.

As a slight nitpick, it seems ScaleDenominator is misreported in the metadata as it should be \(\frac{y_\text{top} - y_\text{bottom}}{16} =\) 4194304.