gdal - Converting Ordnance Survey raster maps to WGS84 webmap tiles

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

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.