2

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:

<ows:WGS84BoundingBox crs="urn:ogc:def:crs:OGC:2:84">
    <ows:LowerCorner>-7.793077 49.852539</ows:LowerCorner>
    <ows:UpperCorner>1.790425 60.894042</ows:UpperCorner>
</ows:WGS84BoundingBox>

and the WMTS TileMatrixSet has the tag:

<ows:Title>WGS84</ows:Title>
<ows:Abstract>WGS84 EPSG:4326</ows:Abstract>
<ows:Identifier>WGS84</ows:Identifier>
<ows:SupportedCRS>urn:ogc:def:crs:EPSG::4326</ows:SupportedCRS>

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?

  • 3
    The idea of tile map services is that there is no need to include georeferencing to the tiles because their location is known when the structure of the tile matrix is known as well as the zoom level-row-column combination of the tile. You probably need some library that can tell you the limits of each tile by its name like on the map in this web page maptiler.com/google-maps-coordinates-tile-bounds-projection. – user30184 May 20 at 21:02
3

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:

<TileMatrix>
  <ows:Identifier>GlobalCRS84Geometric:19</ows:Identifier>
  <ScaleDenominator>533.182395962446</ScaleDenominator>
  <TopLeftCorner>90.0 -180.0</TopLeftCorner>
  <TileWidth>256</TileWidth>
  <TileHeight>256</TileHeight>
  <MatrixWidth>1048576</MatrixWidth>
  <MatrixHeight>524288</MatrixHeight>
</TileMatrix>

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)
  • Thanks for your response. Not sure if I'm just being slow, but say I wanted to find the x of the top left corner of the 1000th row in your above example: -180 + (256*533.18...)*1000 = 136493900. Is this a WGS84 coordinate? – bm13563 May 21 at 12:56
  • sorry that should of course be /scaledenominator not *scaledenominator - that will teach me to go from memory with out looking at the code. – Ian Turton May 21 at 13:00
  • so the top left corner is at -43.5053066336, – Ian Turton May 21 at 13:02
  • I get that for x = TopLeftCorner.x + (TileWidth*ScaleDenominator)/columnCount. This also seems quite high, considering the top left corner is at -180 and this is only row 1000/104857, or have a I misunderstood? – bm13563 May 21 at 13:03
  • no - you are going rowcount rows down the grid and colcount across so they multiply the tile size (width/scaledenom) – Ian Turton May 21 at 13:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.