Converting Ordnance Survey raster maps to WGS84 webmap tiles

Question

I would like to generate Google/Bing-like map tiles from Ordnance Survey's 1:25000 and 1:50000 raster maps, which are provided as 4000x4000 UK national grid tiles (GeoTIFF).

A sample tile can be found at http://www.ordnancesurvey.co.uk/docs/sample-data/25k-raster-sample-data.zip .

Thanks to this site and others (this is my first time using proper GIS tools), I have succeeded in using gdalwrap to reproject the tiles, and gdal_retile.py to cut the reprojected tile into smaller pieces.

I still have the following questions:

1. the source files use a color map, so based on my single-tile sample, I converted that to 24bit color in order to be able to use a better resampling algorithm than nearest neighbor, but I don't know how scalable that will be to the whole of the UK. Is there a way to get a quality resampling without first converting to 24bit color?

2. how do I get gdalwrap or gdal_retile.py (or another utility) to align the generated tiles to the webmap grid (the world divided into 2^zoom tiles of equal lat-lon delta)? Right now my tiles are the right pixel size, but most likely not aligned to the lat-lon grid nor the zoom level

3. the tiles generated by gdal_retile.py are named according to each one's location in the source file, but I need to translate that to the global tile number for my zoom level. Do the tools support doing that on the fly, or do I need to post-process?

4. using gdal_retile.py with the -levels option generates sub-directories but the tiles therein are all black; I was expecting the zoomed out pyramid (and fewer tiles, not more, than the default level), should I be using a different option?

5. even for this smallish sample, each command seems to take an inordinately long time, and this is on an octo-core MacPro and SSD. Is there a performance guide for GDAL?

Here's what I have been doing, starting with the sx99.tif sample file:

pct2rgb.py sx99.tif sx99_24b.tif

gdalwarp -t_srs EPSG:3857 -r lanczos -wm 150 -multi -of GTiff sx99_24b.tif sx99_24b_warped.tif

gdal_retile.py -of png -levels 1 -targetDir tiles sx99_24b_warped.tif 

---

Instead of gdal_retile.py, I've tried using gdal2tile.py on the pre-warped raster maps.

However, while the X coordinate for each of the generated tiles is correct, the Y coordinate seems to be at least double what I am expecting.

I've tried explicitly setting the source SRS and target profile, to no effect.

gdal2tiles.py sx99_24b_warped.tif

I can't understand what's going wrong. The warped intermediate image looks OK:

$ gdalinfo sx99_24b_warped.tif 
Driver: GTiff/GeoTIFF
Files: sx99_24b_warped.tif
Size is 4074, 4085
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-396511.819462091429159,6584391.706861065700650)
Pixel Size = (3.948971656093984,-3.948971656093984)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -396511.819, 6584391.707) (  3d33'42.93"W, 50d47'27.29"N)
Lower Left  ( -396511.819, 6568260.158) (  3d33'42.93"W, 50d41'57.18"N)
Upper Right ( -380423.709, 6584391.707) (  3d25' 2.66"W, 50d47'27.29"N)
Lower Right ( -380423.709, 6568260.158) (  3d25' 2.66"W, 50d41'57.18"N)
Center      ( -388467.764, 6576325.932) (  3d29'22.80"W, 50d44'42.32"N)

But the resulting files are in the range of (for zoom level 11): 1003 <= x <= 1004 and 1359 <= y <= 1360, where I was expecting y around 780...

Answer 1

Regarding execution times: Lanczos is extremely slow, but other resampling algorithms produce similar (perhaps even better) results:

• near (11s)
• bilinear (17s)
• average (27s, a bit soft)
• bicubic (30s, called cubic in GDAL, text thickness sometimes odd)
• antialias (45s, requires pil and numpy)
• cubicspline (1m53s, way too soft)
• lanczos (11m8s, extremely slow)

To my eye, the antialias result is the best compromise, and even without trying to use multiple cores, my desktop should be able to process the whole of the UK in about 48h.

Answer 2

gdal2tiles uses the bottom-left point of origin, not top-left, so the filename pattern is {zoom}/{x}/{tilenum-y}.png where tilenum is the number of vertical tiles for the zoom level.

It's also worth noting that gdal2tiles.py does not require the source data to be pre-warped, but it does require it to be RGB or RGBA.

Answer 3

Try MapTiler (http://www.maptiler.com/).

Installer for Mac OS X is available after you fill the form http://www.maptiler.com/demo/ - and you can drag&drop GeoTIFFs directly (do not make VRTs - they only slow down the tile rendering process).

Here is a video tutorial: https://youtu.be/WQbX4dqlHgY

See an example zoomable tile map of whole UK made of Ordnance Survey OpenData: http://os.tileserver.com/.

(BTW you can get from us the MBTiles from us too)