4

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?

3
+25

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.

  • So do I understand you right that it's not possible to get a global valid index using my projection? Or could it be done reprojecting every coordinate to Web Mercator first? I want to avoid using setting a custom tile origin and want to make it usable out of the box, e. g. QGIS should be able to open it directly. The same NYC is doing here: maps.nyc.gov/tiles/#basemap-layers-xyz – geraldo Sep 6 '18 at 11:53
  • @geraldo yes, not with the Slippy Tile Index (or XYZ or others); this index is made for the Spherical (Web) Mercator projection, and if you manage to tile the UTM map according to that tile matrix (you could cut the UTM map into square tiles), it would not match any other native tiles (and would not fit the bounds etc.). a custom WMTS service could maybe work with UTM, but in general, reprojecting to Spherical (Web) Mercator prior to tiling would be the best to do here. – ThingumaBob Sep 8 '18 at 21:27
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
  • also (since this is .NET, no?) there are function definitions for both directions in a bunch of languages, including Python, on the linked OSM wiki page. – ThingumaBob Sep 4 '18 at 20:40
  • You should add that you are the author of the program. – AndreJ Sep 23 '18 at 11:39

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.