How to convert UTM projected map to slippy map tilenames

Question

I work with a map of Paris in ETRS89 / UTM zone 31N projection (EPSG:25831) and need to write slippy map tiles in filenames like /z/x/y.png.

I made a custom python script to read the QGIS project file and write the tiles but I can't get the tile names right.

The extent of my map [minX, minY, maxX, maxY] is something like [230000, 6270000, 295000, 6225000]. So for example in zoom 10 they should be named something like:

/10/518/351.png
/10/519/351.png
/10/518/352.png
/10/519/352.png

How could I calculate the x and y tile names correctly? For now I try to calculate them starting with the transformation of the upper left corner from degrees to meters. The meter coordinates origin of slippy maps (85W 180N) in meters are:

originX = 471118.87
originY = 9446194.39

Now I calculate the tile size in meters using the resolution in meters (for example 76.43702828517624 for z=10) and tile size in pixels (for example 256):

tileSize = resolutionsInMeters[z] * tileSizeInPixels

With that numbers in mind I calculate the x and y tile names:

x = int( math.floor( ( minX-originX ) / tileSize ) )
y = int( math.floor( ( originY-maxY ) / tileSize ) )

The result definitely is wrong, any hints to get valid tile names?

Answer 1

The Slippy Map tile naming is fairly simple:

• Starting with zoom 0, each tile is divided into 4 children, with the origin in the top-left corner (increasing south- and eastwards), thus

  • for each zoom level, there are

    • n = 2^zoom tiles representing your extend

  • with your bounds, each tile represents (meter doesn´t matter for now)

    • ts_x = (maxX - minX) / n units
    • ts_y = (maxY - minY) / n units

  • a given coordinate pair for you rextend will then have an individual index of

    • x_idx = floor( (x - minX) / ts_x )
    • y_idx = floor( (maxY - y) / ts_y )
      
      

    • in Python

      n = 2**zoom
      x_idx = floor((n * (x - minX)) / (maxX - minX))
      y_idx = floor((n * (maxY - y)) / (maxY - minY))
      

---

Note that, since the whole naming scheme is used with our beloved Web Mercator projection, it assumes a square map as well as 'global' (as far as possible, at least) coverage; this means that

• this assumes you are creating a custom tile matrix covering your AOI

• your tiles will need to be of the same dimension as your [0, 0] tile on zoom 0, in pixels, for this to make sense. Or have your tile renderer create a square extend with blanks.

• since this will create a tile matrix based on your extend (i.e. origin at your extend's (minX, maxY)) and tile dimensions, e.g. Leaflet.js will not display your tiles correctly out-of-the-box. You'd need to set a custom tile origin and the custom dimension in the client viewer of your choice. Or, of course, implement your own tile map viewer.

• this is also a generic replacement for the usual LatLon index search used on Web Mercator maps in most practical environments. It works efficiently with the used CRS' units and would be needed if you want to stick to UTM.

---

Let me add:

You explicity asked for an index with the global Web Mercator reference...that is pretty much not going to work with the Slippy Map Tiles.

UTM doesn´t have a global extend that would be needed to create the index. You would need to create tiles of your area matching the Web Mercator tiles exactly to make this work. Even if you managed to do that, the tiles would, of course, not show the same area...

Projections are not meant to be interchangeable...

At this point, you should probably better reproject the whole thing and work with Web Mercator directly.

Answer 2

Yes, slippy map tile names assume a Web (spherical) Mercator projection so you need to reproject the map. Having done that, this code converts latitude and longitude to slippy map tile numbers. The results (x,y) are reals with the fractional part being the location within the tile.

But there are products - like MAPC2MAPC - that will both reproject and tile the result for you

Friend Sub ll2xyz(ByRef x As Double, ByRef y As Double, ByVal z As Integer, ByVal lat As Double, ByVal lon As Double)

  Dim w1, te, tn As Double

  w1 = 2 ^ z
  te = lon * kdeg2rad
  tn = System.Math.Log((1 + System.Math.Sin(lat * kdeg2rad)) / (1 - System.Math.Sin(lat * kdeg2rad))) / 2
  x = (te + kpi) * (w1 / (kpi * 2))
  y = (kpi - tn) * (w1 / (kpi * 2))
End Sub