Projections

Web mapping tools using tiles use a variant of the Mercator Projection.

This can lead to some significant distortions: you can see this for yourself by zooming out in google maps, and scrolling around. The scale (bottom right) will change, depending on where you are on the globe. Locally, the distortion is not too bad, but tiles from different areas of the globe cannot be directly compared.

Internally, the projection maps the whole globe (between latitudes about -85 to 85 degrees) to a unit square:



In [1]:

    
import tilemapbase
tilemapbase.start_logging()
tilemapbase.tiles.build_OSM().get_tile(0,0,0)









    Out[1]:

Comparison to standards

The standards EPSG:3857 or EPSG:3785 are mentioned in relation to the "web mercator" projection.

The following (reproduced in tilemapbase) shows that the projections are the same, up to a rescaling, and flipping the y axis.

(It is curious to compare the text in the above two links, which to a causual reader might warn against using 3857 while 3785 is fine; they are both identical...)



In [2]:

    
import pyproj
p3857 = pyproj.Proj({"init": "epsg:3857"})
p3785 = pyproj.Proj({"init": "epsg:3785"})



In [3]:

    
scale = 20037508.342789244

def rescale(x, y):
    return ((x - 0.5) * 2 * scale, (0.5 - y) * 2 * scale)

import random, math

diffs = []
for _ in range(10000):
    x = random.random() * 360 - 180
    y = random.random() * 170 - 85

    assert p3857(x, y) == p3785(x, y)
    xx, yy = p3857(x, y)
    xxx, yyy = rescale(*tilemapbase.project(x, y))
    diffs.append( math.sqrt((xx - xxx)**2 + (yy - yyy)**2) )



In [4]:

    
max(diffs), sum(diffs) / 10000









    Out[4]:





(3.129243850708008e-07, 4.728993730872364e-09)

UK National Grid

To support the Ordance Survey (raster) tiles, we need to use the Ordnance Survey National Grid.

This is EPSG:7405
We need to be careful, as the British National Grid uses a different ellipsoid (approximation to the shape of the Earth) than the "usual" GPS based on.
So we need to use the transform method of pyproj. Credit

For example:

Corner of the Garstang Building at the university of Leeds has Lon/Lat coords -1.55532, 53.80474
See OS Maps Online
Claims to be SE 29383 34363
Land's end -5.71808, 50.06942
Is SW 34041 25435
John O'Groats -3.02516, 58.64389
Is ND 40594 73345

Further information: https://www.ordnancesurvey.co.uk/business-and-government/help-and-support/navigation-technology/os-net/formats-for-developers.html



In [5]:

    
import pyproj, functools

bng = pyproj.Proj(init="epsg:27700")
wgs84 = pyproj.Proj(init="epsg:4326")

def project(lon, lat):
    return pyproj.transform(wgs84, bng, lon, lat)

print(project(-1.55532, 53.80474))
print(project(-5.71808, 50.06942))
print(project(-3.02516, 58.64389))









    



(429383.15535284596, 434363.0962840691)
(134041.0757940985, 25435.90742218544)
(340594.4899132418, 973345.1181793772)



In [ ]: