python - Getting spatially referenced tiles from WMTS

I have been able to download tiles from a WMTS request in .tif format using owslib in standalone python:

from owslib.wmts import WebMapTileService

wmts = WebMapTileService('wmts_url')
img = wmts.gettile(layer='v0Lw1953KMdH',

                   tilematrixset='WGS84',
                   tilematrix='16',
                   row='20480',
                   column='32001',
                   format='image/tif')

out = open('C:\\out\\test_out.tif', 'wb')
bytes_written = out.write(img.read())
out.close()

This .tif file is not spatially referenced. The WMTS Layer has the tag:

    -7.793077 49.852539
    1.790425 60.894042

and the WMTS TileMatrixSet has the tag:

WGS84
WGS84 EPSG:4326
WGS84

urn:ogc:def:crs:EPSG::4326

which to someone of my limited knoweledge indicates that there could be a way to obtain a spatially referenced .tif from these tiles. However, I haven't been able to obtain a spatially referenced .tif, despite trying to include the crs in gettile(). Where am I going wrong?

Failing this, are there any other methods that I could use to obtain spatially referenced tiles?

Answer

You should probably be using a WCS service end point to request georeferenced raster data, as @user30184 say's WMTS tiles are designed for display.

But if this is the only source of data you can easily calculate the corners of your tile using some simple maths. If you look in the getCapabilities document for your WMTS you will find the definition of the TileMatrix and for each zoom level it will include a statement like:

  GlobalCRS84Geometric:19
  533.182395962446

  90.0 -180.0
  256
  256
  1048576
  524288

In your request you know the zoomlevel (tilematrix) and row and column count, by combining them with this information you can work out the top left corner of your tile. You also need to know how many pixels there are per map unit, based on a pixel being .28cm (0.28e-3m). For degrees you have to divide by 111319 (see this which gives 60.1 NaMiles per Degree), for projected units just convert to metres and divide by it. See the full java code from GeoTools here.

x = TopLeftCorner.x + (pixelsperdegree * TileWidth * ScaleDenominator) * columnCount
y = TopLeftCorner.y - (pixelsperdegree * TileHeight * ScaleDenominator) * rowCount

As a worked example here is the bottom right tile:

x = -180 + ((256 * 0.28e-3/111319)*533.182395962446)  * 1048575
x = 180.001243876 (I think google calc has some rounding issues)
y = 90 - ((256 * 0.28e-3/111319)*533.182395962446)  * 524287
y = -90.0004 (again needs more precision)

and the bottom right is just

x2 = x +  (pixelsperdegree * TileWidth * ScaleDenominator)
y2 = y +  (pixelsperdegree * TileHeight * ScaleDenominator)