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:
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?
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
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?
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?
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.
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...