4

I'm trying to generate some map tiles based on the Ordnance Survey Open Data (available from here). I don't have much experience with this, but I've read a few tutorials online and have managed to produce some tiles.

The process I'm using is:

# Convert to expanded RGBA
for i in *.tif
do
        gdal_translate -of vrt -expand rgba $i $(basename $i .tif).vrt
done
gdal_merge.py -o merged.tif *.vrt
gdal2tiles.py -r METHOD -z 12 -e -w leaflet -t "Ordnance Survey Raster 250k" merged.tif output_folder

I've tried with METHOD as antialias and as lanczos, although I've done smaller experiments with all the other options.

It works and produces tiles, but they look terrible, so I assume I must be doing something wrong.

The first image below is the result of the lanczos generation, with the tiles combined for a small area around Bath (Zoom: 12, X: 2020-2022, Y: 1363-1365 in google format or 2730-2732 in TMS format). With antialias as the method, the image looks almost identical.

a

By contrast, below is the same area that I downloaded a few years ago from (I think) cz.tileserver.com/osnew, although that doesn't seem to exist anymore, so I can't check whether that's right. As you can see, the image quality is much, much better than what I'm able to produce.

Can anyone point me in the direction of what I need to be doing differently? Or is it just a case of the warping being a fundamental problem and maybe the tiles I've seen before were un-warped and hence not quite in the right place geographically?

enter image description here

Edit: as noted in a comment below, it's more important to me that the map looks good than that it's exactly right in terms of projection, so I'd be happy to find a method that generates tiles and puts them in as close a place as possible if that's always going to look better than the warped tiles. I know what I'm trying to achieve is possible as the lower of the two images above was created from tiles downloaded from a map server. That map server had achieved what I want to do, I just want to be able to do the same myself from an OS Open Data download.

  • With profile=raster the quality loss should be minimal. I tried to test but generated OpenLayers client did not work for me. You can check the quality of the tiles from the tile directories. I made the .vrt from one .TIF only and my command was gdal2tiles -r bilinear -s epsg:27700 -e -w all -p raster -t "Ordnance Survey Raster 250k" temp.vrt bilinear_raster_test. – user30184 Jun 18 at 15:16
  • @user30184 I've tried this now, but the problem with the -p raster option is that it completely disconnects the resulting map from the grid that is used in a map viewer. The tile server from which the bottom image was downloaded shared tiles in a format that could be loaded by Leaflet, Google Maps or (when put in a GEMF file) Locus Pro. The tiles that are produced with -p raster are at a different scale (even with -z 12 added to the command line) and with a different coordinate datum, so can't be used in the same way. – DrAl Jun 19 at 14:15
  • Certainly yes, profile=raster creates tiles to suit a local grid that matches the SRS and extent of the source layer. If you want to use the common Google maps grid then you must re-project and re-scale the tiles which can lead to poorer quality. There may be a place to improve gdal2tiles with this scenario. – user30184 Jun 19 at 21:19
  • @user30184 Thanks. I was coming to that conclusion: that gdal2tiles can't do this for me. I'd just love to know how the site from which I downloaded the tiles that make up the second image managed to do it without the distortion. It must be possible as they did it. Unfortunately the site doesn't seem to exist anymore so I can't even contact the webmaster and ask them. – DrAl Jun 20 at 6:59
1

Reprojecting rasters of this type will always look "odd", if you need a projected GB base map I would recommend using the OS Zoomstack which is a nicely styled vector data set that you can reproject and then tile so the pixels don't get "squashed".

  • Thanks Ian. Part of what I was hoping for was the style of map which I find very easy to read and I was hoping to make a GEMF map for off-line use. As an alternative, is there a good way of making tiles out of the map without reprojecting, but with the least bad positional error when displayed with a WGS84 based mapping tool? – DrAl Jun 13 at 15:25
0

Having got no answers for how to do this the "proper" way that the web service I'd previously found had managed, I ended up doing this in a slightly iffy way, but it seems to work okay-ish as far as I can tell.

I wrote a simple script that works out where the centre of each tile is and creates the tile at that point from the geotiff data. It doesn't try to distort it in any way, so inevitably there will be some discontinuities spread around the map. However, these are relatively infrequent and in general the map is usable and a lot better than the version produced by using gdal2tiles.

Here's my script for any future seekers who are interested:

import os
import glob
import sys
import subprocess

def run_subprocess(args):
    print("\n => Running '%s'" % " ".join(args))
    subprocess.check_call(args)

sys.path.append('/usr/bin')

import PIL

from osgeo import gdal
from gdal2tiles import GlobalMercator

USE_CENTRAL = True
SKIP_IF_PRESENT = True

PIL.Image.MAX_IMAGE_PIXELS = 2000000000

input_data = '/zfs.mount/notbackedup_or_encrypted/ras250_gb_2019-06.zip'
workdir = '/zfs.mount/notbackedup_or_encrypted/raster_working_central'

tile_root = os.path.join(workdir, 'tiles', 'OS-250k-Raster')

if not os.path.exists(workdir):
    os.makedirs(workdir)
os.chdir(workdir)

data_dir = os.path.join(workdir, 'ras250_gb/data')
input_geotiff = os.path.join(workdir, 'merged.tif')
png_file = input_geotiff.replace('.tif', '.png')

if not os.path.exists(data_dir):
    print("Unzipping data")
    run_subprocess(['unzip', input_data])
input_files = glob.glob(os.path.join(data_dir, '*.tif'))
# Speed up...
# subset_tiles = [i for i in input_files if i in ['SO.tif', 'ST.tif', 'SP.tif', 'SU.tif']]
# subset_tiles = [i for i in input_files if i in ['SQ.tif', 'SR.tif', 'SV.tif', 'SW.tif']]
preprocessed_files = [i.replace('.tif', '_expanded.tif') for i in input_files]

if not SKIP_IF_PRESENT or not os.path.exists(input_geotiff):
    print("Expanding RGBA")
    for fi, fo in zip(input_files, preprocessed_files):
        print("Expanding " + fi)
        run_subprocess(['gdal_translate', '-of', 'GTiff', '-expand', 'rgba', fi, fo])
    print("Merging expanded files")
    run_subprocess(['gdal_merge.py', '-o', input_geotiff] + preprocessed_files)
if not SKIP_IF_PRESENT or not os.path.exists(png_file):
    print("Making PNG")
    # Consider "-expand rgb" and then no convert, or "-expand rgba" and later convert
    run_subprocess(['gdal_translate', '-of', 'PNG', input_geotiff, png_file])

mercator = GlobalMercator()
src = gdal.osr.SpatialReference()
src.SetWellKnownGeogCS("WGS84")

print("Loading GeoTiff")
dataset = gdal.Open(input_geotiff)
projection = dataset.GetProjection()
dst = gdal.osr.SpatialReference(projection)
transform = dataset.GetGeoTransform()
ct = gdal.osr.CoordinateTransformation(src,dst)

print("Loading Image")
image = PIL.Image.open(png_file).convert('RGB')
imw, imh = image.size

# Need to define tiles as array of x,y,z
for z in [12]:
    for x in range(1949, 2070):
        print("X=%d" % x)
        dest_path = os.path.join(tile_root, str(z), str(x))
        if not os.path.exists(dest_path):
            os.makedirs(dest_path)
        for y in range(1168, 1393):
            llbounds = mercator.TileLatLonBounds(x, y, z)

            lat1 = -llbounds[0]
            lon1 = llbounds[1]
            lat2 = -llbounds[2]
            lon2 = llbounds[3]

            if USE_CENTRAL:
                lat = (lat1+lat2)/2.0
                lon = (lon1+lon2)/2.0
                xy = ct.TransformPoint(lon, lat)
                xp = int(((xy[0] - transform[0]) / transform[1]))
                yp = int(((xy[1] - transform[3]) / transform[5]))
                xp1 = xp - 121
                xp2 = xp + 122
                yp1 = yp - 121
                yp2 = yp + 122
            else:
                xy1 = ct.TransformPoint(lon1, lat1)
                xy2 = ct.TransformPoint(lon2, lat2)
                xp1 = int(((xy1[0] - transform[0]) / transform[1]))
                xp2 = int(((xy2[0] - transform[0]) / transform[1]))
                yp1 = int(((xy1[1] - transform[3]) / transform[5]))
                yp2 = int(((xy2[1] - transform[3]) / transform[5]))

            if xp1 > imw or xp2 > imw:
                # Tile beyond image boundary
                continue
            elif yp1 > imh or yp2 > imh:
                # Tile beyond image boundary
                continue
            elif xp1 < 0 or xp2 < 0:
                # Tile beyond image boundary
                continue
            elif yp1 < 0 or yp2 < 0:
                # Tile beyond image boundary
                continue
            else:
                tile_img = image.crop((xp1, yp1, xp2, yp2))
                extrema = tile_img.convert("L").getextrema()
                if extrema[0] == extrema[1]:
                    # All pixels the same
                    continue
                tile_img = tile_img.resize((256,256), PIL.Image.LANCZOS)
                tile_img.save(os.path.join(dest_path, '%d.jpg'  % y))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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