Stitching MapProxy tiles to receive one georeferenced image?

by Bernd V.   Last Updated April 15, 2019 23:22 PM

In my struggle to get a decent aerial image for playing around with as a layer in Biosphere3D (and for other purposes as well), I managed to install Mapproxy and then cached the area in my desired zoom level from a 2m-Orthophoto WMS through the proxy within QGIS..

Well, now I am in possession of tiles worth 60MB, but i see no way how to stitch them together to one image.

The tiles are stored in a folder-structure like this: enter image description here

The following sets of tiles are in 32/000/781/000/25, 32/000/781/000/25 etc. containing 408.png to 431.png

I have no clue how to interpret this structure.

I wonder if there is any way (script) to stitch those tiles to receive one georeferenced image.

The tiles were captured in WGS84.

Grateful for any hint on how to proceed.



Answers 3


By default MapProxy uses a TileCache compatible directory layout (zz/xxx/xxx/xxx/yyy/yyy/yyy.format), but you can use directory_layout for switch to TMS compatible directories (zz/xxxx/yyyy.format).

I've wrote little python script for stitching tiles from TileCache compatible directories into one big image. PIL library should be installed:

import os
import Image

z_level = 0
x_level = 3
y_level = 6
tile_size = 256

# Example of tile path: /home/denis/tilecache2tms/mymapproxy/cache_data/osm_cache_EPSG900913/05/000/000/023/000/000/016.png

base = "/home/denis/tilecache2tms/mymapproxy/cache_data/osm_cache_EPSG900913/"
path = "05"
tiles_paths = []

merged_path = os.path.join(base, path, 'merged.png') 
if os.path.isfile(merged_path):
    os.remove(merged_path)

for root, dirs, files in os.walk(os.path.join(base, path)):
    for f in files:
        fname = f.split('.')[0]
        tiles_paths.append(os.path.join(root, fname).split('/')[-7:])

x_num_tiles = len(os.listdir(os.path.join(base, "/".join(tiles_paths[0][:x_level]))))
y_num_tiles = max([len(os.listdir(os.path.join(base, "/".join(tiles_paths[i][:y_level])))) for i in range(len(tiles_paths))])

width = tile_size*x_num_tiles
height = tile_size*y_num_tiles

print "Image size: %s x %s" % (str(width), str(height))

max_x_tile = max([int(d) for d in os.listdir(os.path.join(base, "/".join(tiles_paths[0][:x_level])))])
max_y_tile = max([int(tiles_paths[i][y_level]) for i in range(len(tiles_paths))])

print "Xmax: %s, Ymax: %s" % (str(max_x_tile), str(max_y_tile))

img = Image.new('RGBA', (width, height))

for tile in tiles_paths:
    x_tile = int("".join(tile[1:x_level+1]))
    y_tile = int("".join(tile[4:y_level+1]))
    tile_img = Image.open("%s.%s" % (os.path.join(base, "/".join(tile)), 'png'))
    img.paste(tile_img, (width - (max_x_tile-x_tile+1)*tile_size, (max_y_tile-y_tile)*tile_size))

img.save(merged_path)

print "Image saved into: %s" % merged_path

Now you can use any appropriate software for georeferencing your image by coordinates.

UPD:

For checking PIL works correctly try to create an empty image using PIL and paste not transparent tile into it (use your own paths):

import Image
img = Image.new('RGBA', (256, 256))
tile = Image.open('/home/denis/tilecache2tms/mymapproxy/cache_data/osm_cache_EPSG900913/06/000/000/042/000/000/035.png')
img.paste(tile, (0,0))
img.save('/home/denis/test.png')
drnextgis
drnextgis
September 09, 2013 08:55 AM

If you are running on Windows, the taho.exe tool might be helpful for you:

http://www.dimitri-junker.de/eng/html/openstreetmap.html

Although written to receive online tiles, it can be easily configured to use locally available tiles from file:/// or localhost

The output can be set to larger tiles, or one complete image.

For Linux, you would have to extend the taho.pl by yourself.

AndreJ
AndreJ
September 10, 2013 06:14 AM

Thanks for the excellent piece of code! I made a small bugfix for y positions of tiles:

import os
# import Image from Pillow
from PIL import Image

z_level = 0
x_level = 3
y_level = 6
tile_size = 256

# Example of tile path: /home/denis/tilecache2tms/mymapproxy/cache_data/osm_cache_EPSG900913/05/000/000/023/000/000/016.png

base = "/home/denis/tilecache2tms/mymapproxy/cache_data/osm_cache_EPSG900913/"

path = "05"    
tiles_paths = []

merged_path = os.path.join(base, path, 'merged.png')
if os.path.isfile(merged_path):
    os.remove(merged_path)

for root, dirs, files in os.walk(os.path.join(base, path)):
    for f in files:
        fname = f.split('.')[0]
        tiles_paths.append(os.path.join(root, fname).split('/')[-7:])

x_num_tiles = len(os.listdir(os.path.join(base, "/".join(tiles_paths[0][:x_level]))))
y_num_tiles = max([len(os.listdir(os.path.join(base, "/".join(tiles_paths[i][:y_level])))) for i in range(len(tiles_paths))])

width = tile_size*x_num_tiles
height = tile_size*y_num_tiles

print("Image tiles: %s x %s" % (str(x_num_tiles), str(y_num_tiles)))

print("Image size: %s x %s" % (str(width), str(height)))

min_x_tile = min([int(d) for d in os.listdir(os.path.join(base, "/".join(tiles_paths[0][:x_level])))])
min_y_tile = min([int(tiles_paths[i][y_level]) for i in range(len(tiles_paths))])

print("Xmin: %s, Ymin: %s" % (str(min_x_tile), str(min_y_tile)))

max_x_tile = max([int(d) for d in os.listdir(os.path.join(base, "/".join(tiles_paths[0][:x_level])))])
max_y_tile = max([int(tiles_paths[i][y_level]) for i in range(len(tiles_paths))])

print("Xmax: %s, Ymax: %s" % (str(max_x_tile), str(max_y_tile)))

img = Image.new('RGBA', (width, height))

for tile in tiles_paths:
    x_tile = int("".join(tile[1:x_level+1]))
    y_tile = int("".join(tile[4:y_level+1]))

    tile_img = Image.open("%s.%s" % (os.path.join(base, "/".join(tile)), 'png'))

    x_tile_position = width - (max_x_tile-x_tile+1)*tile_size
    y_tile_position = height - (max_y_tile-y_tile)*tile_size

    print('assemble tile:',x_tile,',',y_tile,'->',x_tile_position,',',y_tile_position)

    img.paste(tile_img, (x_tile_position, y_tile_position))

img.save(merged_path)

print("Image saved into: %s" % merged_path)
user3788120
user3788120
April 15, 2019 22:50 PM

Related Questions


Updated July 21, 2015 14:09 PM

Updated February 18, 2019 12:22 PM

Updated July 26, 2017 17:22 PM

Updated May 28, 2018 05:22 AM

Updated May 02, 2016 08:09 AM