Regression in xy_bbox when tile matrices bounds differ from the boundingBox?

[samn]

We've been using morecantile to figure out the tile coordinates covering a polygon in a WMTS' tile matrix set. We map fields in WMTSCapabilities.xml to fields in morecantile models to construct a custom TileMatrixSet, and compute the coordinates of tiles in the TMS for a given lat lng bbox.

Many WMTSCapabilities return a different bounding box in the layer's ows:BoundingBox than in the TileMatrixSet. A recent change to the implementation of TileMatrixSet#xy_bbox caused incorrect tile coordinates to be generated in this case, and I'm not totally sure if it's just that the WMTS' tile matrices are malformed, or if it's something morecantile should be able to handle.

Here's an example WMTS with this problem: https://geodata.md.gov/imap/rest/services/Imagery/MD_SixInchImagery/MapServer/WMTS/1.0.0/WMTSCapabilities.xml

The example script below shows that the ows:BoundingBox and the bbox that morecantile computes from the tile matrices differ. The previous implementation of xy_bbox would use the boundingBox (instead of computing one from the tile matrices), and tile coordinates ended up getting computed correctly.

On morecantile >= 3.4.0 the tile coordinates are calculated incorrectly (25,318 total, instead of 4 as expected):

Example Script (morecantile 5)

import morecantile

# https://geodata.md.gov/imap/rest/services/Imagery/MD_SixInchImagery/MapServer/WMTS/1.0.0/WMTSCapabilities.xml

min_zoom = 0
max_zoom = 14

# The ScaleDenominator, MatrixWidth, and MatrixHeight for the TileMatrix with Identifier 0
scale_denominator = 8735660.375447182
matrix_width = 10
matrix_height = 21

tile_matrices = []
for zoom in range(min_zoom, max_zoom + 1):
  tile_matrices.append(morecantile.models.TileMatrix(
    id=str(zoom),
    scaleDenominator=scale_denominator,
    cellSize=matrix_height / 2.0**zoom,
    tileWidth=512,
    tileHeight=512,
    matrixWidth=matrix_width,
    matrixHeight=matrix_height,
    cornerOfOrigin="topLeft",
    # All the TileMatrixes in the WMTS have the same TopLeftCorner 
    pointOfOrigin=(-2.00377E7, 3.02411E7),
  ))
  scale_denominator /= 2
  matrix_width *=  2
  matrix_height *= 2



crs = 'http://www.opengis.net/def/crs/EPSG/0/3857'

# Taken from the `Imagery_MD_SixInchImagery` Layer ows:BoundingBox
tms_bbox = morecantile.models.TMSBoundingBox(crs=crs, lowerLeft=(-8852727.0526, 4559519.201499999), upperRight=(-8349977.3074, 4830929.361399999))

tms = morecantile.models.TileMatrixSet(
  crs=crs,
  boundingBox=tms_bbox,
  tileMatrices=tile_matrices
)

z = 14
bbox_lat_lng = (-76.19430541992188, 38.575547964508196, -76.19396209716797, 38.575816369153586)

tiles = list(tms.tiles(*bbox_lat_lng, z))
print(len(tiles))
assert len(tiles) == 4

But on morecantile==3.3.0 the correct number of tiles are computed:

Example Script (morecantile 3.3.0)

import morecantile

# https://geodata.md.gov/imap/rest/services/Imagery/MD_SixInchImagery/MapServer/WMTS/1.0.0/WMTSCapabilities.xml

min_zoom = 0
max_zoom = 14

scale_denominator = 8735660.375447182
matrix_width = 10
matrix_height = 21

tile_matrices = []
for zoom in range(min_zoom, max_zoom + 1):
  tile_matrices.append(morecantile.models.TileMatrix(
    identifier=str(zoom),
    scaleDenominator=scale_denominator,
    tileWidth=512,
    tileHeight=512,
    matrixWidth=matrix_width,
    matrixHeight=matrix_height,
    topLeftCorner=(-2.00377E7, 3.02411E7),
  ))
  scale_denominator /= 2
  matrix_width *=  2
  matrix_height *= 2



crs = 'http://www.opengis.net/def/crs/EPSG/0/3857'

tms_bbox = morecantile.models.TMSBoundingBox(crs=crs, lowerCorner=(-8852727.0526, 4559519.201499999), upperCorner=(-8349977.3074, 4830929.361399999))

tms = morecantile.models.TileMatrixSet(
  identifier='test',
  title='test',
  supportedCRS=crs,
  boundingBox=tms_bbox,
  tileMatrix=tile_matrices
)

z = 14
bbox_lat_lng = (-76.19430541992188, 38.575547964508196, -76.19396209716797, 38.575816369153586)

tiles = list(tms.tiles(*bbox_lat_lng, z))
print(len(tiles))
assert len(tiles) == 4

I realize the root cause here might be the inconsistency between the overall bounds, and the bounds on the tile matrices, but would appreciate any guidance on how to handle this. Here's another example WMTS with this problem. It looks like both are from ArcGIS.

[vincentsarago]

👋 @samn, I'm really sorry about this. We decided to stop using the top level bounds because:

• it's optional
• it was sometime wrong
• I'm pretty sure the specification says that the TMS first level should describe the full bbox

In the same CRS as the TileMatrixSet, boundingBox should be considered informative about the area covered by this TileMatrixSet. It SHOULD NOT be used to calculate the position of the tiles in the CRS space. Instead use the cornerOfOrigin and the pointOfOrigin of the corresponding TileMatrix. If data is not available for the entire tiled space, TileMatrixSetLimits will declare what tiles have data (see Clause 8.2.1, Requirements class 7).

from https://docs.ogc.org/is/17-083r4/17-083r4.html#toc16

ref #106

[samn]

Ahh thanks for explaining, @vincentsarago

do you think this because the WMTS and TMS specs are not exactly aligned with how they treat the bounding box? I was surprised to find so many examples of WMTS GetCapabilities with this problem.