Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
    
I seem to have missed a useful post gis.stackexchange.com/questions/5101/using-gdal-to-cut-tiles that suggests using gdal2tiles.py. This correctly handles the projection, alignment to the webmap grid and zoom, and pyramids. It's said to be a bit on the slow side, and still requires a prior conversion to RGB. –  Pierre-Luc Paour Sep 19 '11 at 20:44
add comment

2 Answers

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

  • near near (11s)
  • bilinear bilinear (17s)
  • average average (27s, a bit soft)
  • cubic bicubic (30s, called cubic in GDAL, text thickness sometimes odd)
  • antialias antialias (45s, requires pil and numpy)
  • cubicspline cubicspline (1m53s, way too soft)
  • enter image description here 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.

share|improve this answer
add comment
up vote 2 down vote accepted

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.