1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
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)
|
