Getting tile bounds for EPSG:25833 WMTS tile using Python

Question

The Norwegian Polar Institute provides this WMTS endpoint. I can successfully open the tileset in QGIS, and observe the queries made, using the Debugging/Development Tools (F12):

One such query is this one:

https://geodata.npolar.no/arcgis/rest/services/Basisdata/NP_Ortofoto_Svalbard_WMTS_25833/MapServer/WMTS/tile/1.0.0/Basisdata_NP_Ortofoto_Svalbard_WMTS_25833/default/default028mm/17/31066/133128

Response:

Based on the URL, this tile has tile id 17/31066/133128. The endpoint defines the projection used: EPSG:25833.

How can I get the bounding box of a particular tile knowing its WMTS tile ID and the used projection using Python?

---

To help verify a possible solution, I found the tile on QGIS canvas and manually drew an approximate bounding box around it. The approximate bounding box (in EPSG:25833):

POLYGON((514863 8682927, 514905 8682927, 514905 8682970, 514863 8682970, 514863 8682927))

---

This is a similar question, but it seems to be R specific. Furthermore, EPSG:3857 is used (as explained in this answer), whereas I am interested in the more general case.

While this question is similar, the accepted answer assumes we know how to get the span of each tile (unlike in this thread), and the other answer is using 111319, stating:

111319 is a magic number for the number of metres in 60 Nautical Miles because we are using degrees here, for a projected CRS it would be the number of metres per unit of the projection.

The example above is clearly a "projected CRS", so it needs a "number of metres per unit of the projection" value. (I don't see how this value could be a constant, but that's another topic.)

Answer 1

Method I: Tedious calculation by hand

Start with the WMTS endpoint. Find <TileMatrix> for zoom level 17:

<TileMatrix>
 <ows:Identifier>17</ows:Identifier>
 <ScaleDenominator>590.588978796931</ScaleDenominator>
 <TopLeftCorner>-5120900.0 9998100.0</TopLeftCorner>
 <TileWidth>256</TileWidth>
 <TileHeight>256</TileHeight>
 <MatrixWidth>141182</MatrixWidth>
 <MatrixHeight>41237</MatrixHeight>
</TileMatrix>

Find ScaleDenominator value: 590.588978796931. Find TileWidth and TileHeight: 256. After reading What is WMTS' ScaleDenominator?, you'll know that

to calculate span of each tile in meters by multiplying TileWidth (or TileHeight) by ScaleDenominator * 0.00028.

Therefore, the edgelength of a zoom level 17 tile is: 256 * 590.588978796931 * 0.00028 = 42.3334180002 units.

Find <TopLeftCorner>: -5120900.0 9998100.0. To calculate the top left corner of 17/31066/133128, we need to add 133128 edgelengths to -5120900.0 and deduct 31066 edgelengths from 9998100.0. Ie:

x = -5120900.0 + 256 * 590.588978796931 * 0.00028 * 133128 # 514863.2715258347
y = 9998100.0 - 256 * 590.588978796931 * 0.00028 * 31066   # 8682970.036406904

To find the centre of this tile, we need to add and deduct half an edgelength more:

x = -5120900.0 + 256 * 590.588978796931 * 0.00028 * (133128 + 0.5) # 514884.438234834
y = 9998100.0 - 256 * 590.588978796931 * 0.00028 * (31066 + 0.5)   # 8682948.869697904

To create the tile squre, we can rely on the Creating square buffers around points using shapely? thread:

import shapely.geometry

p = shapely.geometry.Point(x, y).buffer(
    256 * 590.588978796931 * 0.00028 / 2, cap_style=3
)

where 256 * 590.588978796931 * 0.00028 / 2 is half an edgelength. Get WKT geometry via p.wkt:

POLYGON ((514905.60494383407 8682970.036406904, 514905.60494383407 8682927.702988904, 514863.27152583393 8682927.702988904, 514863.27152583393 8682970.036406904, 514905.60494383407 8682970.036406904))

Check in QGIS (using QuickWKT):

Cf query for 17/31066/133128. The result looks accurate.

---

Method II: just use a function

This method relies on OWSLib. pip3 install OWSLib. We can get tile bounds using tileid_to_wkt below:

from owslib.wmts import WebMapTileService
import shapely.geometry

url = "https://geodata.npolar.no/arcgis/rest/services/Basisdata/NP_Ortofoto_Svalbard_WMTS_25833/MapServer/WMTS/"

def tileid_to_wkt(z, x, y, url=url):
    wmts = WebMapTileService(url)
    m = wmts.tilematrixsets["default028mm"].tilematrix[str(z)]
    assert m.tileheight == m.tilewidth

    u = m.topleftcorner[0] + m.tileheight * m.scaledenominator * 0.00028 * (y + 0.5)
    v = m.topleftcorner[1] - m.tileheight * m.scaledenominator * 0.00028 * (x + 0.5)

    p = shapely.geometry.Point(u, v).buffer(
        m.tileheight * m.scaledenominator * 0.00028 / 2, cap_style=3
    )

    return p.wkt

tileid_to_wkt(17,31066,133128) will then return:

POLYGON ((514905.60494383407 8682970.036406904, 514905.60494383407 8682927.702988904, 514863.27152583393 8682927.702988904, 514863.27152583393 8682970.036406904, 514905.60494383407 8682970.036406904))

ie the same string we have arrived at using Method I.