2

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.

4

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')
|improve this answer|||||
  • Many thanx for your script!! Unsuccessful so far to get it run properly. With the first attempts, it seemed to stitch the tiles, but the resulting png was empty/white (I think this might be the same PIL-issue I faced while installing MapProxy in a virtualenv, and resolved somehow). So my thought was, that it must be possible to run the script in this virtualenv of MapProxy. But from the moment I started MapProxy, the path in the script now appends the MapProxy-path to the script-path: No such file or directory: '/home/bernd/mymapproxy/cache_data/bydop_cache_EPSG4326/home/bernd/mymapproxy' – Bernd V. Sep 9 '13 at 10:02
  • Saw that you edited the script a bit. Now the path-issue doesn't occur anymore. The script says Image size: 11264 x 6144 Xmax: 823, Ymax: 431 But the resulting image is still blank. Do you think this is a PIL-issue or is it with the script? Very confusing subject ... Thanx for you help so far! (Am on Linux Mint 14 btw.) – Bernd V. Sep 9 '13 at 11:08
  • I've added additional information in answer how you can check your PIL works correctly. – drnextgis Sep 9 '13 at 11:44
  • Thank you very much! PIL worked. Spent some hours investigating the code till I found out, that i didn't set the grid in the mapproxy.yaml to geodetic. So all my tiles where stitched together somewhere in outer space. Strange the WMS proxy worked nevertheless. Phew! Took quite some time, but now I am looking into a brighter future ;) – Bernd V. Sep 9 '13 at 22:30
0

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.

|improve this answer|||||
  • I tried the program on my VirtualBox Win7. Does not work here (no reaction on any input). But thanx for the hint. – Bernd V. Sep 10 '13 at 10:31
  • It might be that a Visual C Runtime is needed to run it. I use it in real Windows without problems. – AndreJ Sep 10 '13 at 11:19
-1

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)
|improve this answer|||||

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.