Local Projections in a World of Spherical Mercator

02 Jul 2012

This text was originally posted on my previous employer's blog (Kartena), but since they've closed it down, I'm publishing it here as well.

Update 2014-02-24: A more thorough and updated version on the topic is The Hitchhiker's Guide to Tiled Maps. The premise: you want to use the new and cool open source mapping solutions like Leaflet, TileStache and TileMill. While you love projects like OpenStreetMap, you are in a position where you cannot throw out your old map tiles, that are rendered in some local projection, instead of the the one projection to rule them all, spherical mercator. You need to coerce these tools to use your projection and your tiles. Fortunately, these tools do support local projections, or at least were designed with this support in mind. Yay! Less fortunately, it seems very few actually use that functionality, and if they do, they are not bloggers, since we at Kartena couldn't find any useful information on how to go about configuring local projections. That's why we're here! We want to share the wonderful world of local projections with you.

What are projections?

A projection (in cartography) is a way of transforming a location on the earth, typically described by a latitude and longitude, to a flat surface (a map) coordinate, like x and y. Projections come in a myriad and shapes and flavors, and different projections have different advantages and disadvantages (you can read up on that on Wikipedia, for example).

What is tiling?

Tiling means that we chop up a large (sometimes really, truly, huge) bitmap into smaller more manageable bitmaps that can later be put together to appear as one seamless image again. Tiling is used since it would not be possible to render those huge bitmaps (too large to keep in memory), and because time to transfer the bitmap to a client would be  unacceptable. Also, the client rarely looks at the entire map area, and only need small parts of it.

A contract for tiling

To make tiling possible, we need to set a standard for how to reference a certain part of the larger bitmap - if the client needs to show a certain geographic area, it needs to know:

1. Tile size and scales: how many tiles fit into that area
2. Tile origin: what the (projected) coordinates of those tiles' corners are
3. Grid alignment: given the tile grid origin, in what directions do the row and column coordinates of the grid grow?
4. Tile naming: how to request a tile from the server, given its row and column coordinate in the grid

These four points constitute a contract about how tiling is done, and server and client need to agree on this. For the spherical Mercator projection, Google Maps has set up what has ended being a de facto standard for tiling in this projection. This standard makes it easy to use spherical mercator with many different combinations of map clients and tile servers, since the contract is well defined, well documented and implemented almost everywhere (some variations exist, more on that later). For other projections, this represents more of a challenge: all of the points above can be specified in a number of ways. We need to break down the points to something we can express in code.

Tile size and scales

We first have to agree on how large each tile is, both in pixels and in projected coordinates. Deciding the tile size in pixels is fairly easy. Google Maps uses 256x256 pixels, and if you don't have any particular reason to do otherwise, so should probably you. The size of a tile in projected coordinates is somewhat trickier. Another way of formulating this is: what zoom/scale levels do you need. Each level will result in a different tile size in projected coordinates. For example, let's say that you have a level with the scale of one pixel per 100 projected coordinates (scale 1:100). If your tile size, in pixels, is 256x256, the tile size in projected coordinates will be 256 / (1/100) = 256x100 = 25 600 projected coordinate units. The other way around, you can take projected coordinate units and multiply them by scale to get pixels: 25 600 projected coordinate units = 25 600x(1/100) = 256 pixels. To sum up, we need to specify these two things:

• The tile size in pixels
• The scale for each zoom/scale level.

Tile origin

Tiles are organized in a grid. The rows and columns of the grid are numbered, and to make this numbering unambiguous, we need to define where in the projected coordinate space the grids origin (0,0) is. Each projection transforms a certain area of the earth into a certain coordinate space. This coordinate space varies for each projection. Take for example EPSG:2400 (Swedish coordinate system RT90). It is well defined for latitudes 55.2 to 69.1 degrees north, and longitudes 10.57 to 24.18 degrees east. This is projected to the x-coordinates 1166653.6161 to 2032341.6763 and y-coordinates 6131388.6471 to 7690505.5552. In this case, it might be suitable to choose an origin at top left (1166653.6161, 7690505.5552) or bottom left (1166653.6161, 6131388.6471). Or you might even define the origin as (0,0), or any other coordinate that suits you. The important thing is that client and server agrees on the origin, or they will actually use two different grids.

Grid alignment

To further define the grid, we need to decide how rows and columns in the grid are numbered. We know the coordinate of the origin of (0,0), but what does it mean, in terms of projected coordinates, to increment the row or column coordinate? In the world of spherical Mercator, there are to approaches here:

• Google Maps and OpenStreetMap set the grid's origin at the top left and row coordinate grows downward, such that the row coordinate grows the further south you go. Columns grow from west to east (left to right).
• TMS (Tile Map Service) set the origin at bottom left, such that the row coordinate grows the further north you go. Columns grow from west to east (left to right), like Google Maps.

Tile naming

Given the definitions above, we can unambiguously calculate the grid row and column for a projected coordinate. What remains is how to request this tile from the server. This basically boils down to deciding on a naming standard for the tiles. One of the most popular, employed by both Google Maps/OSM and TMS is often summarized as xyz. This standard calls the row y, the column x and the zoom level z, where zoom levels are numbered from 0 (most zoomed out, largest scale) and more zoomed in with increasing numbers. An example of a tile can look like this:

http://mytileserver.com/sweden/3/127/32.png

A great number of naming strategies exist, instead of using the rows and columns, the projected coordinate of one of the tile's corners can for example be used. Zoom level can be denoted by the actual scale used, or resolution (that is, effectively the inverse of the scale).

Getting down to details

This is all fine and well, but also very theoretical. How do we actually implement this, if we want to use our own definition instead of spherical Mercator? A widely used and pretty flexible way to break this problem down is to use these four components:

• Projection - projects latitude and longitude to x/y
• Scale definitions - defines what zoom levels are available and the scale for each level
• Coordinate to grid transformation matrix - describes how a projected coordinate is transformed to a row and column in the grid
• Tile naming template - given zoom level, row and column forms the filename, URL, et c. for the tile

Projection

Projection can be done in a lot of ways. One of the most flexible is to use an implementation of Proj.4 for the language you're working in. At Kartena, we use Proj4Js in Javascript and pyproj for Python. spatialreference.org has Proj.4 definitions for all frequently used coordinate reference systems.

Scale definitions

Depending on your requirements, this can either be a list of your zoom levels' scales, or a function, if zoom levels follow such a pattern.

Coordinate to grid transformation matrix

Getting projected coordinates to grid coordinates involves first applying the grid origin and the grid alignment, and then the scale to convert the projected coordinate to pixel coordinates. These are all so called linear transformations, which can be expressed by these formulas:

grid_col = int(scale * (a * proj_x + b))
grid_row = int(scale * (c * proj_y + d))

where scale is the current scale level and a, b, c and d are constants; these constants can be expressed as a transformation matrix. For spherical Mercator, this matrix looks like this: (0.5 / π, 0.5, -0.5 / π, 0.5). This can be interpreted as the origin is at (b,d) = (0.5, 0.5). Both x and y axis is scaled by a factor of 0.5 / π, but the y axis is inverted (grows from north to south) since c has a negative sign. At first, this transformation might be the hardest one to grasp (especially if you're not a fan of math). To make it easier, here are two variants of the transformation matrix that generally works well for many local projections.

Local projection with google maps/OSM tiling scheme

For a local projection where you want to use the Google Maps/OpenStreetMap tiling scheme, where rows grow from north to south (and your projected coordinates go the other way, so that y grows further north), you want to do this:

• Set origo at the top left of the covered area (the point furthest north and west). Call this (origo_x, origo_y)
• The transformation matrix is: (1, origo_x, -1, origo_y)

LOCAL PROJECTION WITH TMS TILING SCHEME

For a local projection where you want to use the TMS tiling scheme, where rows grow from south to north (and your projected coordinates go the same way, so that y grows further north), you want to do this:

• Set origo at the bottom left of the covered area (the point furthest south and west). Call this (origo_x, origo_y)
• The transformation matrix is: (1, origo_x, 1, -origo_y)

Tile Naming Template

This is usually one of the easiest steps, usually just a set of string replacements.

Practical example

Wow, you've almost reached the end! Let's take a quick look how this might look in practice. At Kartena, our favourite map client is currently Leaflet. Our only real problem with it, was that out of the box it doesn't support any local projections. To solve this, we have written a small project called Proj4Leaflet, that acts as glue between Leaflet and Proj4Js. With this project, and the information above, you can easily have tiling maps with any projection supported by Proj.4. Using what we've talked about above, we can put this all together in a few lines of code:

var resolutions = [8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1, 0.5]
,crs = L.CRS.proj4js('EPSG:2400'
    ,'+lon_0=15.808277777799999 +lat_0=0.0 +k=1.0 +x_0=1500000.0 +y_0=0.0 +proj=tmerc +ellps=bessel +units=m +towgs84=414.1,41.3,603.1,-0.855,2.141,-7.023,0 +no_defs'
    ,new L.Transformation(1, 0, -1, 0))
    ,scale: function(zoom) {
      return 1 / resolutions[zoom];
    }
    [...]
})
,mapUrl = 'http://api.geosition.com/tile/lmv/{z}/{x}/{y}.png'

First, we define the resolutions that we are going to use for the different zoom levels. Note that resolution (meters/pixel) is the inverse of scale (pixels/meter). This is reflected by the use of our custom scale function, which returns the inverted value.

We then set up the CRS (Coordinate Reference System) with its name, "EPSG:2400", and the Proj.4 definition (the really long string).

Next, we set up the transformation matrix, exactly as described above. In this case, we have put origo at (0,0) and use the Google Maps/OSM tiling scheme, with rows growing from north to south. (Note that this results in row coordinates being negative; this is something we use for the "cleanness" of using a simple origo as (0,0)).

The last line is the string template for getting the actual tile URL, following the xyz scheme mentioned above.