diff options
| -rw-r--r-- | get-wms.py | 112 |
1 files changed, 112 insertions, 0 deletions
diff --git a/get-wms.py b/get-wms.py new file mode 100644 index 0000000..e199d89 --- /dev/null +++ b/get-wms.py @@ -0,0 +1,112 @@ +from math import ceil +import os, sys + +from owslib.wms import WebMapService +from pyproj import Proj, transform +from PIL import Image + +wms = WebMapService(sys.argv[1], version='1.1.1') + +if len(sys.argv) < 3: + first_layer = None + for layer in wms.contents: + if first_layer is None: + first_layer = layer + print layer + print wms[first_layer].boundingBox + print wms[first_layer].boundingBoxWGS84 + sys.exit() +else: + layers = sys.argv[8:] + bbox = [float(b) for b in sys.argv[2:6]] + coord_sys = wms[layers[0]].boundingBox[4] + wgs = Proj(init='epsg:4326') + out_proj = Proj(init=coord_sys) + result_bbox = transform(wgs, out_proj, bbox[0], bbox[1]) + transform(wgs, out_proj, bbox[2], bbox[3]) + + image_size = [int(sys.argv[6]), int(sys.argv[7])] + bbox_size = [abs(result_bbox[2] - result_bbox[0]), abs(result_bbox[3] - result_bbox[1])] + bbox_ratio = bbox_size[0] / bbox_size[1] + image_size_corrected = [ + int(max(image_size[0], image_size[1]*bbox_ratio)), + int(max(image_size[1], image_size[0]/bbox_ratio))] + print 'Correcting image size to %dx%d due to coordinates box ratio' % ( + image_size_corrected[0], image_size_corrected[1]) + + image_size = [int(ceil(i/400.0))*400 for i in image_size_corrected] + print 'Correcting image size to %dx%d due to tile size (400x400)' % ( + image_size[0], image_size[1]) + + min_long = min(result_bbox[0], result_bbox[2]) + max_long = max(result_bbox[0], result_bbox[2]) + long_span = max_long - min_long + min_lat = min(result_bbox[1], result_bbox[3]) + max_lat = max(result_bbox[1], result_bbox[3]) + lat_span = max_lat - min_lat + tile_long_span = long_span * 400.0 / image_size_corrected[0] + tile_lat_span = lat_span * 400.0 / image_size_corrected[1] + + tiles = [] + current_lat = min_lat + current_long = min_long + while current_lat < max_lat: + current_row = [] + while current_long < max_long: + current_row.append([ + current_long, current_lat, + current_long+tile_long_span, current_lat+tile_lat_span]) + current_long += tile_long_span + if len(current_row): + tiles.append(current_row) + current_lat += tile_lat_span + current_long = min_long + + files = [] + x = 0 + y = 0 + for row in tiles: + filerow = [] + for tile in row: + print 'Fetching tile %dx%d: %s' % (x, y, tile) + img = wms.getmap(layers=layers, + srs=coord_sys, + bbox=tile, + size=(400,400), + format='image/png', + transparent=True) + filename = 'out-%d-%d.png' % (x, y) + out = open(filename, 'wb') + out.write(img.read()) + out.close() + filerow.append(filename) + y += 1 + files.append(filerow) + x += 1 + y = 0 + + strips = [] + r = 0 + for row in files: + tiles = [Image.open(f) for f in row] + row_image = Image.new('RGBA', (400*len(tiles), 400)) + offset = 0 + for tile in tiles: + row_image.paste(tile, (offset, 0)) + offset += 400 + strip_filename = 'out-%d.png' % (r) + row_image.save(strip_filename) + strips.append(strip_filename) + for f in row: + os.remove(f) + r += 1 + + strips.reverse() + strip_images = [Image.open(f) for f in strips] + final_image = Image.new('RGBA', (strip_images[0].width, 400*len(strip_images))) + offset = 0 + for strip_image in strip_images: + final_image.paste(strip_image, (0, offset)) + offset += 400 + final_image.save('out.png') + for s in strips: + os.remove(s) |