Codebase list twms / b8f28d2
New upstream version 0.07z+git20200829+cb7d39a Andrej Shadura 3 years ago
19 changed file(s) with 1542 addition(s) and 1210 deletion(s). Raw diff Collapse all Expand all
77 from twms import twms
88
99
10 if __name__ != '__main__':
11 try:
12 from mod_python import apache, util
13 import datetime
14 except ImportError:
15 pass
10 if __name__ != "__main__":
11 try:
12 from mod_python import apache, util
13 import datetime
14 except ImportError:
15 pass
1616
17
17
1818 def handler(req):
1919 """
2020 A handler for mod_python.
99
1010 __name__ = "twms"
1111
12
1213 def read(fname):
13 return open(os.path.join(os.path.dirname(__file__), fname), 'rb').read().decode('utf-8')
14 return (
15 open(os.path.join(os.path.dirname(__file__), fname), 'rb')
16 .read()
17 .decode('utf-8')
18 )
19
1420
1521 def glob(fname):
1622 return abs_glob(os.path.join(os.path.dirname(__file__), fname))
23
1724
1825 def man_path(fname):
1926 category = fname.rsplit('.', 1)[1]
2229 def man_files(pattern):
2330 return list(map(man_path, glob(pattern)))
2431
32
2533 def config_files():
2634 if not is_windows:
2735 return [(os.path.join('/etc', __name__), [os.path.join('twms', 'twms.conf')])]
2836 else:
2937 return []
3038
39
3140 # monkey patch setuptools to use distutils owner/group functionality
3241 from setuptools.command import sdist
42
3343 sdist_org = sdist.sdist
44
45
3446 class sdist_new(sdist_org):
3547 def initialize_options(self):
3648 sdist_org.initialize_options(self)
3749 self.owner = self.group = 'root'
50
51
3852 sdist.sdist = sdist_new
3953
4054 setup(
55 import sys
66 from lxml import etree
77
8 tiles_cache = "/var/www/latlon/wms/cache/"
9 layers = ["irs", "yhsat", "DGsat", "yasat","SAT"]
8 tiles_cache = "/var/www/latlon/wms/cache/"
9 layers = ["irs", "yhsat", "DGsat", "yasat", "SAT"]
1010 default_user = "Komzpa"
1111 user = default_user
1212 ts = ""
1515 nodes = {}
1616 curway = []
1717 for lay in layers:
18 print(lay)
19 src = lay+".osm"
20 if os.path.exists(src):
21 corrfile = tiles_cache + lay+ "/rectify.txt"
22 corr = open(corrfile, "w")
23 context = etree.iterparse(open(src))
24 for action, elem in context:
25 items = dict(elem.items())
26 ts = items.get("timestamp",ts)
27
28 # print items
29 if elem.tag == "node":
30 nodes[int(items["id"])] = (float(items["lon"]), float(items["lat"]))
31 elif elem.tag == "nd":
32 curway.append(nodes[int(items["ref"])])
33 elif elem.tag == "tag":
34 if items["k"] == "user":
35 user = items["v"]
36 if items["k"] == "timestamp":
37 ts = items["v"]
38 elif elem.tag == "way":
39 ts = items.get("timestamp",ts)
40 print("%s %s %s %s %s %s"% (curway[0][0],curway[0][1],curway[1][0],curway[1][1], user, ts), file=corr)
41 curway = []
42 user = default_user
43 ts = ""
18 print(lay)
19 src = lay + ".osm"
20 if os.path.exists(src):
21 corrfile = tiles_cache + lay + "/rectify.txt"
22 corr = open(corrfile, "w")
23 context = etree.iterparse(open(src))
24 for action, elem in context:
25 items = dict(elem.items())
26 ts = items.get("timestamp", ts)
27
28 # print items
29 if elem.tag == "node":
30 nodes[int(items["id"])] = (float(items["lon"]), float(items["lat"]))
31 elif elem.tag == "nd":
32 curway.append(nodes[int(items["ref"])])
33 elif elem.tag == "tag":
34 if items["k"] == "user":
35 user = items["v"]
36 if items["k"] == "timestamp":
37 ts = items["v"]
38 elif elem.tag == "way":
39 ts = items.get("timestamp", ts)
40 print(
41 "%s %s %s %s %s %s"
42 % (
43 curway[0][0],
44 curway[0][1],
45 curway[1][0],
46 curway[1][1],
47 user,
48 ts,
49 ),
50 file=corr,
51 )
52 curway = []
53 user = default_user
54 ts = ""
1010
1111 for lay in layers:
1212 nid, wid = 1, 1
13 out = open(lay+".osm","w")
14 corrfile = tiles_cache + lay+ "/rectify.txt"
13 out = open(lay + ".osm", "w")
14 corrfile = tiles_cache + lay + "/rectify.txt"
1515 corr = open(corrfile, "r")
1616 print('<osm version="0.6">', file=out)
1717 for line in corr:
18 d,c,b,a,user,ts = line.split()
19 d,c,b,a = (float(d),float(c),float(b),float(a))
20 print('<node id="-%s" lon="%s" lat="%s" version="1" />'%(nid, d, c), file=out)
21 nid += 1
22 print('<node id="-%s" lon="%s" lat="%s" version="1" />'%(nid, b, a), file=out)
23 nid += 1
24 print('<way id="-%s" version="1">'% wid, file=out)
25 print(' <nd ref="-%s" />'%(nid-2), file=out)
26 print(' <nd ref="-%s" />'%(nid-1), file=out)
27 print(' <tag k="%s" v="%s" />"'%("user", user), file=out)
28 print(' <tag k="%s" v="%s" />"'%("timestamp", ts), file=out)
29 print("</way>", file=out)
30 wid += 1
31 print('</osm>', file=out)
18 d, c, b, a, user, ts = line.split()
19 d, c, b, a = (float(d), float(c), float(b), float(a))
20 print('<node id="-%s" lon="%s" lat="%s" version="1" />' % (nid, d, c), file=out)
21 nid += 1
22 print('<node id="-%s" lon="%s" lat="%s" version="1" />' % (nid, b, a), file=out)
23 nid += 1
24 print('<way id="-%s" version="1">' % wid, file=out)
25 print(' <nd ref="-%s" />' % (nid - 2), file=out)
26 print(' <nd ref="-%s" />' % (nid - 1), file=out)
27 print(' <tag k="%s" v="%s" />"' % ("user", user), file=out)
28 print(' <tag k="%s" v="%s" />"' % ("timestamp", ts), file=out)
29 print("</way>", file=out)
30 wid += 1
31 print("</osm>", file=out)
00 # -*- coding: utf-8 -*-
11
2 __all__ = ["twms", "bbox", "canvas", "correctify", "drawing", "projections", "reproject"]
2 __all__ = [
3 "twms",
4 "bbox",
5 "canvas",
6 "correctify",
7 "drawing",
8 "projections",
9 "reproject",
10 ]
99
1010
1111 def point_is_in(bbox, point):
12 """
13 Checks whether EPSG:4326 point is in bbox
14 """
15 #bbox = normalize(bbox)[0]
16
17 return point[0]>=bbox[0] and point[0]<=bbox[2] and point[1]>=bbox[1] and point[1]<=bbox[3]
12 """
13 Check whether EPSG:4326 point is in bbox
14 """
15 # bbox = normalize(bbox)[0]
1816
19 def bbox_is_in(bbox_outer, bbox_to_check, fully = True):
20 """
21 Checks whether EPSG:4326 bbox is inside outer
22 """
23 bo = normalize(bbox_outer)[0]
24 bc = normalize(bbox_to_check)[0]
25 if fully:
26 return (bo[0]<=bc[0] and bo[2]>=bc[2]) and (bo[1]<=bc[1] and bo[3]>=bc[3])
27 else:
28 if bo[0] > bc[0]:
29 bo, bc = bc, bo
30 if bc[0] <= bo[2]:
31 if bo[1] > bc[1]:
32 bo, bc = bc, bo
33 return bc[1] <= bo[3]
34 return False
17 return (
18 point[0] >= bbox[0]
19 and point[0] <= bbox[2]
20 and point[1] >= bbox[1]
21 and point[1] <= bbox[3]
22 )
3523
3624
25 def bbox_is_in(bbox_outer, bbox_to_check, fully=True):
26 """
27 Check whether EPSG:4326 bbox is inside outer
28 """
29 bo = normalize(bbox_outer)[0]
30 bc = normalize(bbox_to_check)[0]
31 if fully:
32 return (bo[0] <= bc[0] and bo[2] >= bc[2]) and (
33 bo[1] <= bc[1] and bo[3] >= bc[3]
34 )
35 else:
36 if bo[0] > bc[0]:
37 bo, bc = bc, bo
38 if bc[0] <= bo[2]:
39 if bo[1] > bc[1]:
40 bo, bc = bc, bo
41 return bc[1] <= bo[3]
42 return False
3743
38 return ((bo[0]<=bc[0] and bo[2]>=bc[0]) or (bo[0]<=bc[2] and bo[2]>=bc[2])) and ((bo[1]<=bc[1] and bo[3]>=bc[1]) or (bo[1]<=bc[3] and bo[3]>=bc[3])) or ((bc[0]<=bo[0] and bc[2]>=bo[0]) or (bc[0]<=bo[2] and bc[2]>=bo[2])) and ((bc[1]<=bo[1] and bc[3]>=bo[1]) or (bc[1]<=bo[3] and bc[3]>=bo[3]))
44 return (
45 ((bo[0] <= bc[0] and bo[2] >= bc[0]) or (bo[0] <= bc[2] and bo[2] >= bc[2]))
46 and (
47 (bo[1] <= bc[1] and bo[3] >= bc[1])
48 or (bo[1] <= bc[3] and bo[3] >= bc[3])
49 )
50 or (
51 (bc[0] <= bo[0] and bc[2] >= bo[0])
52 or (bc[0] <= bo[2] and bc[2] >= bo[2])
53 )
54 and (
55 (bc[1] <= bo[1] and bc[3] >= bo[1])
56 or (bc[1] <= bo[3] and bc[3] >= bo[3])
57 )
58 )
59
3960
4061 def add(b1, b2):
41 """
42 Returns bbox that contains two bboxes.
43 """
44 return (min(b1[0],b2[0]),min(b1[1],b2[1]),max(b1[2],b2[2]),max(b1[3],b2[3]))
62 """
63 Return bbox containing two bboxes.
64 """
65 return (min(b1[0], b2[0]), min(b1[1], b2[1]), max(b1[2], b2[2]), max(b1[3], b2[3]))
66
4567
4668 def expand_to_point(b1, p1):
47 """
48 Expands bbox b1 to contain p1: [(x,y),(x,y)]
49 """
50 for p in p1:
51 b1 = add(b1, (p[0],p[1],p[0],p[1]))
52 return b1
53 def normalize (bbox):
54 """
55 Normalizes EPSG:4326 bbox order. Returns normalized bbox, and whether it was flipped on horizontal axis.
56 """
57
58 flip_h = False
59 bbox = list(bbox)
60 while bbox[0] < -180.:
61 bbox[0] += 360.
62 bbox[2] += 360.
63 if bbox[0] > bbox[2]:
64 bbox = (bbox[0],bbox[1],bbox[2]+360,bbox[3])
65 #bbox = (bbox[2],bbox[1],bbox[0],bbox[3])
66 if bbox[1] > bbox[3]:
67 flip_h = True
68 bbox = (bbox[0],bbox[3],bbox[2],bbox[1])
69
70
71
72 return bbox, flip_h
69 """
70 Expand bbox b1 to contain p1: [(x,y),(x,y)]
71 """
72 for p in p1:
73 b1 = add(b1, (p[0], p[1], p[0], p[1]))
74 return b1
7375
7476
75 def zoom_for_bbox (bbox, size, layer, min_zoom = 1, max_zoom = 18,max_size = (10000,10000)):
76 """
77 Calculate a best-fit zoom level
78 """
77 def normalize(bbox):
78 """
79 Normalise EPSG:4326 bbox order. Returns normalized bbox, and whether it was flipped on horizontal axis.
80 """
7981
80 h,w = size
81 for i in range (min_zoom,max_zoom):
82 cx1, cy1, cx2, cy2 = projections.tile_by_bbox (bbox, i, layer["proj"])
83 if w is not 0:
84 if (cx2-cx1)*256 >= w*0.9 :
85 return i
86 if h is not 0:
87 if (cy1-cy2)*256 >= h*0.9:
88 return i
89 if (cy1-cy2)*256 >= max_size[0]/2:
90 return i
91 if (cx2-cx1)*256 >= max_size[1]/2:
92 return i
93 return max_zoom
82 flip_h = False
83 bbox = list(bbox)
84 while bbox[0] < -180.0:
85 bbox[0] += 360.0
86 bbox[2] += 360.0
87 if bbox[0] > bbox[2]:
88 bbox = (bbox[0], bbox[1], bbox[2] + 360, bbox[3])
89 # bbox = (bbox[2],bbox[1],bbox[0],bbox[3])
90 if bbox[1] > bbox[3]:
91 flip_h = True
92 bbox = (bbox[0], bbox[3], bbox[2], bbox[1])
93
94 return bbox, flip_h
95
96
97 def zoom_for_bbox(bbox, size, layer, min_zoom=1, max_zoom=18, max_size=(10000, 10000)):
98 """
99 Calculate a best-fit zoom level
100 """
101
102 h, w = size
103 for i in range(min_zoom, max_zoom):
104 cx1, cy1, cx2, cy2 = projections.tile_by_bbox(bbox, i, layer["proj"])
105 if w != 0:
106 if (cx2 - cx1) * 256 >= w * 0.9:
107 return i
108 if h != 0:
109 if (cy1 - cy2) * 256 >= h * 0.9:
110 return i
111 if (cy1 - cy2) * 256 >= max_size[0] / 2:
112 return i
113 if (cx2 - cx1) * 256 >= max_size[1] / 2:
114 return i
115 return max_zoom
1111 ##
1212
1313
14
1514 import projections
1615
1716 try:
1918 except ImportError:
2019 import Image, ImageFilter
2120
22 import urllib2
21 import urllib
2322 from io import BytesIO
2423 import datetime
2524 import sys
2625 import threading
2726
27
2828 def debug(st):
29 sys.stderr.write(str(st)+"\n")
29 sys.stderr.write(str(st) + "\n")
30
3031
3132 class WmsCanvas:
33 def __init__(
34 self,
35 wms_url=None,
36 proj="EPSG:4326",
37 zoom=4,
38 tile_size=None,
39 mode="RGBA",
40 tile_mode="WMS",
41 ):
42 self.wms_url = wms_url
43 self.zoom = zoom
44 self.proj = proj
45 self.mode = mode
46 self.tile_mode = tile_mode
47 self.tile_height = 256
48 self.tile_width = 256
3249
50 if tile_size:
51 self.tile_width, self.tile_height = tile_size
52 self.tiles = {}
3353
34 def __init__(self, wms_url = None, proj = "EPSG:4326", zoom = 4, tile_size = None, mode = "RGBA", tile_mode = "WMS"):
35 self.wms_url = wms_url
36 self.zoom = zoom
37 self.proj = proj
38 self.mode = mode
39 self.tile_mode = tile_mode
40 self.tile_height = 256
41 self.tile_width = 256
54 def __setitem__(self, x, v):
55 x, y = x
56 tile_x = int(x / self.tile_height)
57 x = x % self.tile_height
58 tile_y = int(y / self.tile_width)
59 y = y % self.tile_width
60 try:
61 self.tiles[(tile_x, tile_y)]["pix"][x, y] = v
62 except KeyError:
63 self.FetchTile(tile_x, tile_y)
64 self.tiles[(tile_x, tile_y)]["pix"][x, y] = v
4265
43 if tile_size:
44 self.tile_width, self.tile_height = tile_size
45 self.tiles = {}
46 def __setitem__ (self, x, v):
47 x, y = x
48 tile_x = int(x/self.tile_height)
49 x = x % self.tile_height
50 tile_y = int(y/self.tile_width)
51 y = y % self.tile_width
52 try:
53 self.tiles[(tile_x, tile_y)]["pix"][x,y] = v
54 except KeyError:
55 self.FetchTile(tile_x, tile_y)
56 self.tiles[(tile_x, tile_y)]["pix"][x,y] = v
57
58 def __getitem__ (self, x):
59 x, y = x
60 tile_x = int(x/self.tile_height)
61 x = x % self.tile_height
62 tile_y = int(y/self.tile_width)
63 y = y % self.tile_width
64 try:
65 return self.tiles[(tile_x, tile_y)]["pix"][x,y]
66 except KeyError:
67 self.FetchTile(tile_x, tile_y)
68 return self.tiles[(tile_x, tile_y)]["pix"][x,y]
69
70 def ConstructTileUrl (self, x, y):
71 if self.tile_mode == "WMS":
72 a,b,c,d = projections.from4326(projections.bbox_by_tile(self.zoom, x, y, self.proj), self.proj)
73 return self.wms_url + "width=%s&height=%s&srs=%s&bbox=%s,%s,%s,%s" % (self.tile_width, self.tile_height, self.proj, a,b,c,d)
74 else:
75 return self.wms_url + "z=%s&x=%s&y=%s&width=%s&height=%s" % (self.zoom-1, x, y, self.tile_width, self.tile_height)
66 def __getitem__(self, x):
67 x, y = x
68 tile_x = int(x / self.tile_height)
69 x = x % self.tile_height
70 tile_y = int(y / self.tile_width)
71 y = y % self.tile_width
72 try:
73 return self.tiles[(tile_x, tile_y)]["pix"][x, y]
74 except KeyError:
75 self.FetchTile(tile_x, tile_y)
76 return self.tiles[(tile_x, tile_y)]["pix"][x, y]
7677
78 def ConstructTileUrl(self, x, y):
79 if self.tile_mode == "WMS":
80 a, b, c, d = projections.from4326(
81 projections.bbox_by_tile(self.zoom, x, y, self.proj), self.proj
82 )
83 return self.wms_url + "width=%s&height=%s&srs=%s&bbox=%s,%s,%s,%s" % (
84 self.tile_width,
85 self.tile_height,
86 self.proj,
87 a,
88 b,
89 c,
90 d,
91 )
92 else:
93 return self.wms_url + "z=%s&x=%s&y=%s&width=%s&height=%s" % (
94 self.zoom - 1,
95 x,
96 y,
97 self.tile_width,
98 self.tile_height,
99 )
77100
78 def FetchTile(self, x, y):
79 if (x,y) in self.tiles:
80 if self.tiles[(x, y)]["status"] == "DL":
81 self.tiles[(x, y)]["thread"].join()
82 else:
83 im = ""
84 if self.wms_url:
85 remote = self.ConstructTileUrl (x, y)
86 debug(remote)
87 ttz = datetime.datetime.now()
88 contents = urllib2.urlopen(remote).read()
89 debug("Download took %s" % str(datetime.datetime.now() - ttz))
90 im = Image.open(BytesIO(contents))
91 if im.mode != self.mode:
92 im = im.convert(self.mode)
93 else:
94 im = Image.new(self.mode, (self.tile_width,self.tile_height))
95 debug("blanking tile")
96 self.tiles[(x,y)] = {}
97 self.tiles[(x,y)]["im"] = im
98 self.tiles[(x,y)]["pix"] = im.load()
99 self.tiles[(x, y)]["status"] = "RD"
100
101 def PreparePixel(self, x, y):
102 tile_x = int(x/self.tile_height)
103 x = x % self.tile_height
104 tile_y = int(y/self.tile_width)
105 y = y % self.tile_width
106 if not (tile_x, tile_y):
107 self.tiles[(tile_x, tile_y)] = {}
108 self.tiles[(tile_x, tile_y)]["status"] = "PP"
109 if not "pix" in self.tiles[(tile_x, tile_y)]:
110 if self.tiles[(tile_x, tile_y)]["status"] == "PP":
111 self.tiles[(tile_x, tile_y)]["status"] = "DL"
112 self.tiles[(tile_x, tile_y)]["thread"] = threading.Thread(group=None, target=self.FetchTile, name=None, args=(self, tile_x, tile_y), kwargs={})
113 self.tiles[(tile_x, tile_y)]["thread"].start()
114 return self.tiles[(tile_x, tile_y)]["status"] == "RD"
101 def FetchTile(self, x, y):
102 if (x, y) in self.tiles:
103 if self.tiles[(x, y)]["status"] == "DL":
104 self.tiles[(x, y)]["thread"].join()
105 else:
106 im = ""
107 if self.wms_url:
108 remote = self.ConstructTileUrl(x, y)
109 debug(remote)
110 ttz = datetime.datetime.now()
111 contents = urllib.urlopen(remote).read()
112 debug("Download took %s" % str(datetime.datetime.now() - ttz))
113 im = Image.open(BytesIO(contents))
114 if im.mode != self.mode:
115 im = im.convert(self.mode)
116 else:
117 im = Image.new(self.mode, (self.tile_width, self.tile_height))
118 debug("blanking tile")
119 self.tiles[(x, y)] = {}
120 self.tiles[(x, y)]["im"] = im
121 self.tiles[(x, y)]["pix"] = im.load()
122 self.tiles[(x, y)]["status"] = "RD"
115123
124 def PreparePixel(self, x, y):
125 tile_x = int(x / self.tile_height)
126 x = x % self.tile_height
127 tile_y = int(y / self.tile_width)
128 y = y % self.tile_width
129 if not (tile_x, tile_y):
130 self.tiles[(tile_x, tile_y)] = {}
131 self.tiles[(tile_x, tile_y)]["status"] = "PP"
132 if not "pix" in self.tiles[(tile_x, tile_y)]:
133 if self.tiles[(tile_x, tile_y)]["status"] == "PP":
134 self.tiles[(tile_x, tile_y)]["status"] = "DL"
135 self.tiles[(tile_x, tile_y)]["thread"] = threading.Thread(
136 group=None,
137 target=self.FetchTile,
138 name=None,
139 args=(self, tile_x, tile_y),
140 kwargs={},
141 )
142 self.tiles[(tile_x, tile_y)]["thread"].start()
143 return self.tiles[(tile_x, tile_y)]["status"] == "RD"
116144
117 def PixelAs4326(self,x,y):
118 return projections.coords_by_tile(self.zoom, 1.*x/self.tile_width, 1.*y/self.tile_height, self.proj)
119
120
121 def PixelFrom4326(self,lon,lat):
122 a,b = projections.tile_by_coords((lon, lat), self.zoom, self.proj)
123 return a*self.tile_width, b*self.tile_height
145 def PixelAs4326(self, x, y):
146 return projections.coords_by_tile(
147 self.zoom, 1.0 * x / self.tile_width, 1.0 * y / self.tile_height, self.proj
148 )
124149
125 def MaxFilter(self, size = 5):
126 tiles = self.tiles.keys()
127 for tile in tiles:
128 self.tiles[tile]["im"] = self.tiles[tile]["im"].filter(ImageFilter.MedianFilter(size))
129 self.tiles[tile]["pix"] = self.tiles[tile]["im"].load()
150 def PixelFrom4326(self, lon, lat):
151 a, b = projections.tile_by_coords((lon, lat), self.zoom, self.proj)
152 return a * self.tile_width, b * self.tile_height
153
154 def MaxFilter(self, size=5):
155 tiles = self.tiles.keys()
156 for tile in tiles:
157 self.tiles[tile]["im"] = self.tiles[tile]["im"].filter(
158 ImageFilter.MedianFilter(size)
159 )
160 self.tiles[tile]["pix"] = self.tiles[tile]["im"].load()
77 import config
88 import projections
99
10
1011 def get(version, ref):
11 content_type = "text/xml"
12
13 if version == "1.0.0":
14 req = """
12 content_type = "text/xml"
13
14 if version == "1.0.0":
15 req = (
16 """
1517 <?xml version="1.0" standalone="no"?>
1618 <!-- The DTD (Document Type Definition) given here must correspond to the version number declared in the WMT_MS_Capabilities element below. -->
1719 <!DOCTYPE WMT_MS_Capabilities SYSTEM "http://www2.demis.nl/WMS/capabilities_1_0_0.dtd"
2830
2931 <!-- end of DOCTYPE declaration -->
3032 <!-- The version number listed in the WMT_MS_Capabilities element here must correspond to the DTD declared above. See the WMT specification document for how to respond when a client requests a version number not implemented by the server. -->
31 <WMT_MS_Capabilities version=\"""" +str(version)+ """">
33 <WMT_MS_Capabilities version=\""""
34 + str(version)
35 + """">
3236 <Service>
3337 <!-- The WMT-defined name for this type of service -->
3438 <Name>GetMap</Name>
3539 <!-- Human-readable title for pick lists -->
36 <Title>""" + config.wms_name + """</Title>
40 <Title>"""
41 + config.wms_name
42 + """</Title>
3743 <!-- Narrative description providing additional information -->
3844
3945 <Abstract>None</Abstract>
4046 <Keywords></Keywords>
4147 <!-- Top-level address of service or service provider. See also onlineResource attributes of <DCPType> children. -->
42 <OnlineResource>"""+ ref +"""</OnlineResource>
48 <OnlineResource>"""
49 + ref
50 + """</OnlineResource>
4351 <!-- Fees or access constraints imposed. -->
4452 <Fees>none</Fees>
4553 <AccessConstraints>none</AccessConstraints>
5866 <DCPType>
5967 <HTTP>
6068 <!-- The URL here for HTTP GET requests includes only the prefix before the query string.-->
61 <Get onlineResource=\"""" + ref + """?"/>
69 <Get onlineResource=\""""
70 + ref
71 + """?"/>
6272 </HTTP>
6373 </DCPType>
6474 </Map>
7080 <DCPType>
7181 <HTTP>
7282 <!-- The URL here for HTTP GET requests includes only the prefix before the query string.-->
73 <Get onlineResource=\"""" + ref + """?"/>
83 <Get onlineResource=\""""
84 + ref
85 + """?"/>
7486 </HTTP>
7587 </DCPType>
7688
8597 </Format>
8698 </Exception>
8799 <Layer>
88 <Title>""" + config.wms_name + """</Title>
100 <Title>"""
101 + config.wms_name
102 + """</Title>
89103 <Abstract/>"""
90 pset = set(projections.projs.keys())
91 pset = pset.union(set(projections.proj_alias.keys()))
92 for proj in pset:
93 req += "<SRS>%s</SRS>" % proj
94 req += """<LatLonBoundingBox minx="-180" miny="-85.0511287798" maxx="180" maxy="85.0511287798"/>
104 )
105 pset = set(projections.projs.keys())
106 pset = pset.union(set(projections.proj_alias.keys()))
107 for proj in pset:
108 req += "<SRS>%s</SRS>" % proj
109 req += """<LatLonBoundingBox minx="-180" miny="-85.0511287798" maxx="180" maxy="85.0511287798"/>
95110 <BoundingBox SRS="EPSG:4326" minx="-184" miny="85.0511287798" maxx="180" maxy="85.0511287798"/>
96111 """
97
98 lala = """<Layer queryable="1">
112
113 lala = """<Layer queryable="1">
99114 <Name>%s</Name>
100115 <Title>%s</Title>
101116 <BoundingBox SRS="EPSG:4326" minx="%s" miny="%s" maxx="%s" maxy="%s"/>
102117 <ScaleHint min="0" max="124000"/>
103118 </Layer>"""
104 for i in config.layers.keys():
105 b = config.layers[i].get("bbox", config.default_bbox)
106 req += lala%(i,config.layers[i]["name"],b[0],b[1],b[2],b[3])
107
108 req += """ </Layer>
119 for i in config.layers.keys():
120 b = config.layers[i].get("bbox", config.default_bbox)
121 req += lala % (i, config.layers[i]["name"], b[0], b[1], b[2], b[3])
122
123 req += """ </Layer>
109124 </Capability>
110125 </WMT_MS_Capabilities>"""
111126
112
113
114
115
116 else:
117 content_type = "application/vnd.ogc.wms_xml"
118 req = """<?xml version="1.0"?>
127 else:
128 content_type = "application/vnd.ogc.wms_xml"
129 req = (
130 """<?xml version="1.0"?>
119131 <!DOCTYPE WMT_MS_Capabilities SYSTEM "http://www2.demis.nl/WMS/capabilities_1_1_1.dtd" [
120132 <!-- Vendor-specific elements are defined here if needed. -->
121133 <!-- If not needed, just leave this EMPTY declaration. Do not
122134 delete the declaration entirely. -->
123135 <!ELEMENT VendorSpecificCapabilities EMPTY>
124136 ]>
125 <WMT_MS_Capabilities version=\""""+ str(version) +"""">
137 <WMT_MS_Capabilities version=\""""
138 + str(version)
139 + """">
126140 <!-- Service Metadata -->
127141 <Service>
128142 <!-- The WMT-defined name for this type of service -->
129143 <Name>twms</Name>
130144 <!-- Human-readable title for pick lists -->
131 <Title>""" + config.wms_name + """</Title>
145 <Title>"""
146 + config.wms_name
147 + """</Title>
132148 <!-- Narrative description providing additional information -->
133149 <Abstract>None</Abstract>
134150 <!-- Top-level web address of service or service provider. See also OnlineResource
135151 elements under <DCPType>. -->
136 <OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href=\"""" + ref+ """"/>
152 <OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href=\""""
153 + ref
154 + """"/>
137155 <!-- Contact information -->
138156 <ContactInformation>
139157 <ContactPersonPrimary>
140 <ContactPerson>"""+config.contact_person["real_name"]+"""</ContactPerson>
141 <ContactOrganization>"""+config.contact_person["organization"]+"""</ContactOrganization>
158 <ContactPerson>"""
159 + config.contact_person["real_name"]
160 + """</ContactPerson>
161 <ContactOrganization>"""
162 + config.contact_person["organization"]
163 + """</ContactOrganization>
142164 </ContactPersonPrimary>
143 <ContactElectronicMailAddress>"""+config.contact_person["mail"]+"""</ContactElectronicMailAddress>
165 <ContactElectronicMailAddress>"""
166 + config.contact_person["mail"]
167 + """</ContactElectronicMailAddress>
144168 </ContactInformation>
145169 <!-- Fees or access constraints imposed. -->
146170 <Fees>none</Fees>
155179 <Get>
156180 <!-- The URL here for invoking GetCapabilities using HTTP GET
157181 is only a prefix to which a query string is appended. -->
158 <OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href=\"""" + ref + """?"/>
182 <OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href=\""""
183 + ref
184 + """?"/>
159185 </Get>
160186 </HTTP>
161187 </DCPType>
170196 <Get>
171197 <!-- The URL here for invoking GetCapabilities using HTTP GET
172198 is only a prefix to which a query string is appended. -->
173 <OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href=\"""" + ref + """?"/>
199 <OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href=\""""
200 + ref
201 + """?"/>
174202 </Get>
175203 </HTTP>
176204 </DCPType>
186214 <VendorSpecificCapabilities/>
187215 <Layer>
188216 <Title>World Map</Title>"""
189 pset = set(projections.projs.keys())
190 pset = pset.union(set(projections.proj_alias.keys()))
191 for proj in pset:
192 req += "<SRS>%s</SRS>" % proj
193 req += """
217 )
218 pset = set(projections.projs.keys())
219 pset = pset.union(set(projections.proj_alias.keys()))
220 for proj in pset:
221 req += "<SRS>%s</SRS>" % proj
222 req += """
194223 <LatLonBoundingBox minx="-180" miny="-85.0511287798" maxx="180" maxy="85.0511287798"/>
195224 <BoundingBox SRS="EPSG:4326" minx="-180" miny="-85.0511287798" maxx="180" maxy="85.0511287798"/>
196225 """
197 lala = """
226 lala = """
198227 <Layer queryable="0" opaque="1">
199228 <Name>%s</Name>
200229 <Title>%s</Title>
202231 <ScaleHint min="0" max="124000"/>
203232 </Layer>
204233 """
205 for i in config.layers.keys():
206 b = config.layers[i].get("bbox", config.default_bbox)
207 req += lala%(i,config.layers[i]["name"],b[0],b[1],b[2],b[3])
208
209 req += """ </Layer>
234 for i in config.layers.keys():
235 b = config.layers[i].get("bbox", config.default_bbox)
236 req += lala % (i, config.layers[i]["name"], b[0], b[1], b[2], b[3])
237
238 req += """ </Layer>
210239 </Capability>
211240 </WMT_MS_Capabilities>"""
212241
213 return content_type, req
242 return content_type, req
44 # the extent permitted by applicable law. You can redistribute it
55 # and/or modify it under the terms specified in COPYING.
66
7 #import sys
87 import projections
98 import os
109 import config
1110
1211
13 distance = lambda z,x,y,g: ((z-y)**2+(x-g)**2)**(0.5)
12 distance = lambda z, x, y, g: ((z - y) ** 2 + (x - g) ** 2) ** (0.5)
13
1414
1515 def has_corrections(layer):
16 corrfile = config.tiles_cache + layer.get("prefix", "")+ "/rectify.txt"
17 return os.path.exists(corrfile)
16 corrfile = config.tiles_cache + layer.get("prefix", "") + "/rectify.txt"
17 return os.path.exists(corrfile)
18
1819
1920 def corr_wkt(layer):
20 corrfile = config.tiles_cache + layer.get("prefix", "")+ "/rectify.txt"
21 corr = open(corrfile, "r")
22 wkt = ""
23 for line in corr:
24 d,c,b,a,user,ts = line.split()
25 d,c,b,a = (float(d),float(c),float(b),float(a))
26 wkt +="POINT(%s %s),LINESTRING(%s %s,%s %s),"%(d,c,d,c,b,a)
27 return wkt[:-1]
28
21 corrfile = config.tiles_cache + layer.get("prefix", "") + "/rectify.txt"
22 corr = open(corrfile, "r")
23 wkt = ""
24 for line in corr:
25 d, c, b, a, user, ts = line.split()
26 d, c, b, a = (float(d), float(c), float(b), float(a))
27 wkt += "POINT(%s %s),LINESTRING(%s %s,%s %s)," % (d, c, d, c, b, a)
28 return wkt[:-1]
29
2930
3031 def rectify(layer, point):
31 corrfile = config.tiles_cache + layer.get("prefix", "")+ "/rectify.txt"
32 corrfile = config.tiles_cache + layer.get("prefix", "") + "/rectify.txt"
3233 srs = layer["proj"]
3334 if not os.path.exists(corrfile):
34 return point
35 return point
3536 corr = open(corrfile, "r")
3637 lons, lats = point
37 loni, lati, lona, lata = projections.projs[projections.proj_alias.get(srs,srs)]["bounds"]
38 loni, lati, lona, lata = projections.projs[projections.proj_alias.get(srs, srs)][
39 "bounds"
40 ]
3841 if (lons is loni and lats is lati) or (lons is lona and lats is lata):
39 return point
40 #print(pickle.dumps(coefs[layer]), file=sys.stderr)
41 # sys.stderr.flush()
42 return point
43 # print(pickle.dumps(coefs[layer]), file=sys.stderr)
44 # sys.stderr.flush()
4245 lonaz, loniz, lataz, latiz = lona, loni, lata, lati
43 maxdist = (1.80)
46 maxdist = 1.80
4447 for line in corr:
45 d,c,b,a,user,ts = line.split()
46 d,c,b,a = (float(d),float(c),float(b),float(a))
47 #for d,c,b,a in coefs[layer]:
48 # print(a,b, distance(lons, lats, b, a), file=sys.stderr)
49 if distance(b,a, lons, lats) < maxdist:
50 if a > lats:
51 if distance(a,b,lats,lons) <= distance(lata,lona,lats,lons):
52 lata = a
53 lataz = c
54 if a < lats:
55 if distance(a,b,lats,lons) <= distance(lati,loni,lats,lons):
56 lati = a
57 latiz = c
58 if b > lons:
59 if distance(a,b,lats,lons) <= distance(lata,lona,lats,lons):
60 lona = b
61 lonaz = d
62 if b < lons:
63 if distance(a,b,lats,lons) <= distance(lati,loni,lats,lons):
64 loni = b
65 loniz = d
66 # print(loni, lati, lona, lata, distance(loni, lati, lona, lata), file=sys.stderr)
67 # print("clat:", (lata-lati)/(lataz-latiz), (lona-loni)/(lonaz-loniz), file=sys.stderr)
68 # sys.stderr.flush()
48 d, c, b, a, user, ts = line.split()
49 d, c, b, a = (float(d), float(c), float(b), float(a))
50 # for d,c,b,a in coefs[layer]:
51 # print(a,b, distance(lons, lats, b, a), file=sys.stderr)
52 if distance(b, a, lons, lats) < maxdist:
53 if a > lats:
54 if distance(a, b, lats, lons) <= distance(lata, lona, lats, lons):
55 lata = a
56 lataz = c
57 if a < lats:
58 if distance(a, b, lats, lons) <= distance(lati, loni, lats, lons):
59 lati = a
60 latiz = c
61 if b > lons:
62 if distance(a, b, lats, lons) <= distance(lata, lona, lats, lons):
63 lona = b
64 lonaz = d
65 if b < lons:
66 if distance(a, b, lats, lons) <= distance(lati, loni, lats, lons):
67 loni = b
68 loniz = d
69 # print(loni, lati, lona, lata, distance(loni, lati, lona, lata), file=sys.stderr)
70 # print("clat:", (lata-lati)/(lataz-latiz), (lona-loni)/(lonaz-loniz), file=sys.stderr)
71 # sys.stderr.flush()
6972
7073 lons, lats = projections.from4326(point, srs)
71 lona, lata = projections.from4326((lona,lata), srs)
72 loni, lati = projections.from4326((loni,lati), srs)
73 lonaz, lataz = projections.from4326((lonaz,lataz), srs)
74 loniz, latiz = projections.from4326((loniz,latiz), srs)
74 lona, lata = projections.from4326((lona, lata), srs)
75 loni, lati = projections.from4326((loni, lati), srs)
76 lonaz, lataz = projections.from4326((lonaz, lataz), srs)
77 loniz, latiz = projections.from4326((loniz, latiz), srs)
78
79 latn = (lats - lati) / (lata - lati)
80 latn = (latn * (lataz - latiz)) + latiz
81 lonn = (lons - loni) / (lona - loni)
82 lonn = (lonn * (lonaz - loniz)) + loniz
83
84 return projections.to4326((lonn, latn), srs)
7585
7686
77 latn = (lats-lati)/(lata-lati)
78 latn = (latn * (lataz-latiz))+latiz
79 lonn = (lons-loni)/(lona-loni)
80 lonn = (lonn * (lonaz-loniz))+loniz
81
82 return projections.to4326((lonn,latn), srs)
83
84 #print(rectify("yasat", (27.679068, 53.885122), ""))
87 # print(rectify("yasat", (27.679068, 53.885122), ""))
8588 def r_bbox(layer, bbox):
86 corrfile = config.tiles_cache + layer.get("prefix", "")+ "/rectify.txt"
89 corrfile = config.tiles_cache + layer.get("prefix", "") + "/rectify.txt"
8790 srs = layer["proj"]
8891 if not os.path.exists(corrfile):
89 return bbox
90 a,b,c,d = projections.from4326(bbox,srs)
91 cx, cy = (a+c)/2, (b+d)/2
92 cx1,cy1 = projections.from4326(rectify(layer,projections.to4326((cx,cy), srs)),srs)
93 a1,b1 = projections.from4326(rectify(layer,projections.to4326((a,b), srs)),srs)
94 c1,d1 = projections.from4326(rectify(layer,projections.to4326((c,d), srs)),srs)
95
96 dx,dy = ((cx1-cx)+(a1-a)+(c1-c))/3, ((cy1-cy)+(b1-b)+(d1-d))/3
97 # print(dx,dy, file=sys.stderr)
98 # sys.stderr.flush()
99 return projections.to4326((a+dx,b+dy,c+dx,d+dy),srs)
100
92 return bbox
93 a, b, c, d = projections.from4326(bbox, srs)
94 cx, cy = (a + c) / 2, (b + d) / 2
95 cx1, cy1 = projections.from4326(
96 rectify(layer, projections.to4326((cx, cy), srs)), srs
97 )
98 a1, b1 = projections.from4326(rectify(layer, projections.to4326((a, b), srs)), srs)
99 c1, d1 = projections.from4326(rectify(layer, projections.to4326((c, d), srs)), srs)
100
101 dx, dy = (
102 ((cx1 - cx) + (a1 - a) + (c1 - c)) / 3,
103 ((cy1 - cy) + (b1 - b) + (d1 - d)) / 3,
104 )
105 # print(dx,dy, file=sys.stderr)
106 # sys.stderr.flush()
107 return projections.to4326((a + dx, b + dy, c + dx, d + dy), srs)
1414 import sys, socket
1515
1616 try:
17 import psyco
18 psyco.full()
17 import psyco
18
19 psyco.full()
1920 except ImportError:
20 pass
21 pass
2122
2223 OK = 200
2324 ERROR = 500
2829 A handler for web.py.
2930 """
3031 resp, ctype, content = twms_main(data)
31 web.header('Content-Type', ctype)
32 web.header("Content-Type", ctype)
3233 return content
3334
3435
36 urls = (
37 "/(.*)/([0-9]+)/([0-9]+)/([0-9]+)(\.[a-zA-Z]+)?(.*)",
38 "tilehandler",
39 "/(.*)",
40 "mainhandler",
41 )
3542
36 urls = (
37 '/(.*)/([0-9]+)/([0-9]+)/([0-9]+)(\.[a-zA-Z]+)?(.*)', 'tilehandler',
38 '/(.*)', 'mainhandler'
39 )
4043
4144 class tilehandler:
4245 def GET(self, layers, z, x, y, format, rest):
5053 "format": format.strip("."),
5154 "z": z,
5255 "x": x,
53 "y": y
56 "y": y,
5457 }
5558 return handler(data)
5659
6063 data = web.input()
6164 data = dict((k.lower(), data[k]) for k in iter(data))
6265 if "ref" not in data:
63 if web.ctx.env['HTTP_HOST']:
64 data["ref"] = web.ctx.env['wsgi.url_scheme'] + "://" + web.ctx.env['HTTP_HOST'] + "/"
66 if web.ctx.env["HTTP_HOST"]:
67 data["ref"] = (
68 web.ctx.env["wsgi.url_scheme"]
69 + "://"
70 + web.ctx.env["HTTP_HOST"]
71 + "/"
72 )
6573 return handler(data)
66
6774
6875
6976 def main():
7077 try:
71 if sys.argv[1] == "josm": # josm mode
72 import cgi
73 url, params = sys.argv[2].split("/?", 1)
74 data = cgi.parse_qs(params)
75 for t in data.keys():
76 data[t] = data[t][0]
77 resp, ctype, content = twms_main(data)
78 print(content)
79 exit()
78 if sys.argv[1] == "josm": # josm mode
79 import cgi
80
81 url, params = sys.argv[2].split("/?", 1)
82 data = cgi.parse_qs(params)
83 for t in data.keys():
84 data[t] = data[t][0]
85 resp, ctype, content = twms_main(data)
86 print(content)
87 exit()
8088 except IndexError:
81 pass
89 pass
8290
8391 try:
84 app = web.application(urls, globals())
85 app.run() # standalone run
92 app = web.application(urls, globals())
93 app.run() # standalone run
8694 except socket.error:
87 print("Can't open socket. Abort.", file=sys.stderr)
88 sys.exit(1)
95 print("Can't open socket. Abort.", file=sys.stderr)
96 sys.exit(1)
8997
9098
91 application = web.application(urls, globals()).wsgifunc() # mod_wsgi
99 application = web.application(urls, globals()).wsgifunc() # mod_wsgi
92100
93101 if __name__ == "__main__":
94102 main()
2121
2222 HAVE_CAIRO = True
2323 try:
24 import cairo
25 import pango
26 import pangocairo
24 import cairo
25 import pango
26 import pangocairo
2727 except ImportError:
28 HAVE_CAIRO = False
28 HAVE_CAIRO = False
2929
30 def wkt(wkt, img, bbox, srs, color, blend = 0.5):
30
31 def wkt(wkt, img, bbox, srs, color, blend=0.5):
3132 """
3233 Simple WKT renderer
3334 """
34 wkt = wkt.replace("((","(")
35 wkt = wkt.replace("((", "(")
3536 for obj in wkt.split("),"):
3637
37 canvas = img.copy()
38 name, coords = obj.split("(")
39 coords = coords.replace(")","")
40 text = ""
41 if name == "TEXT":
42 text, coords = coords.split(",", 1)
43 coords = coords.split(",")
44 coords = [ [float(t) for t in x.split(" ")] for x in coords]
45 canvas = render_vector(name, canvas, bbox, coords, srs, color, text=text)
46 img = Image.blend(img, canvas, blend)
38 canvas = img.copy()
39 name, coords = obj.split("(")
40 coords = coords.replace(")", "")
41 text = ""
42 if name == "TEXT":
43 text, coords = coords.split(",", 1)
44 coords = coords.split(",")
45 coords = [[float(t) for t in x.split(" ")] for x in coords]
46 canvas = render_vector(name, canvas, bbox, coords, srs, color, text=text)
47 img = Image.blend(img, canvas, blend)
4748 return img
4849
49 def gpx(track, img, bbox, srs, color, blend = 0.5):
50
51 def gpx(track, img, bbox, srs, color, blend=0.5):
5052 """
5153 Simple GPX renderer
5254 """
5355 for i in track.tracks.keys():
54 coords = track.getTrack(i)
55 canvas = img.copy()
56 canvas = render_vector("LINESTRING", canvas, bbox, coords, srs, color)
57 img = Image.blend(img, canvas, blend)
56 coords = track.getTrack(i)
57 canvas = img.copy()
58 canvas = render_vector("LINESTRING", canvas, bbox, coords, srs, color)
59 img = Image.blend(img, canvas, blend)
5860 return img
5961
6062
61 def render_vector(geometry, img, bbox, coords, srs, color=None, renderer=None,
62 linestring_width=None, text=None, cairo_font=None, pil_font=None):
63 def render_vector(
64 geometry,
65 img,
66 bbox,
67 coords,
68 srs,
69 color=None,
70 renderer=None,
71 linestring_width=None,
72 text=None,
73 cairo_font=None,
74 pil_font=None,
75 ):
6376 """
6477 Renders a vector geometry on image.
6578 """
6679 if not color:
67 color = config.geometry_color[geometry]
80 color = config.geometry_color[geometry]
6881 if not renderer:
69 renderer = config.default_vector_renderer
82 renderer = config.default_vector_renderer
7083 if not linestring_width:
71 if 'linestring_width' in dir(config):
72 linestring_width = config.linestring_width
73 else:
74 linestring_width = 3
84 if "linestring_width" in dir(config):
85 linestring_width = config.linestring_width
86 else:
87 linestring_width = 3
7588 if not cairo_font:
76 if 'cairo_font' in dir(config):
77 cairo_font = config.cairo_font
89 if "cairo_font" in dir(config):
90 cairo_font = config.cairo_font
7891 if not pil_font:
79 if 'pil_font' in dir(config):
80 pil_font = config.pil_font
92 if "pil_font" in dir(config):
93 pil_font = config.pil_font
8194 bbox = projections.from4326(bbox, srs)
8295 lo1, la1, lo2, la2 = bbox
8396 coords = projections.from4326(coords, srs)
84 W,H = img.size
97 W, H = img.size
8598 prevcoord = False
86 coords = [((coord[0]-lo1)*(W-1)/abs(lo2-lo1), (la2-coord[1])*(H-1)/(la2-la1)) for coord in coords]
99 coords = [
100 (
101 (coord[0] - lo1) * (W - 1) / abs(lo2 - lo1),
102 (la2 - coord[1]) * (H - 1) / (la2 - la1),
103 )
104 for coord in coords
105 ]
87106
88107 if renderer == "cairo" and HAVE_CAIRO:
89 "rendering as cairo"
90 imgd = img.tostring()
91 a = array.array('B',imgd)
92 surface = cairo.ImageSurface.create_for_data (a, cairo.FORMAT_ARGB32, W, H, W*4)
93 cr = cairo.Context(surface)
94 cr.move_to(*coords[0])
95 color = ImageColor.getrgb(color)
96 cr.set_source_rgba(color[2]/255.0, color[1]/255.0, color[0]/255.0, 1)
97 if geometry == "LINESTRING" or geometry == "POLYGON":
98 for k in coords:
99 cr.line_to(*k)
100 if geometry == "LINESTRING":
101 if linestring_width is not None:
102 cr.set_line_width(linestring_width)
103 cr.stroke()
104 elif geometry == "POLYGON":
105 cr.fill()
106 elif geometry == "POINT":
107 cr.arc(coords[0][0],coords[0][1],6,0,2*math.pi)
108 cr.fill()
109 elif geometry == "TEXT" and text and cairo_font:
110 pcr = pangocairo.CairoContext(cr)
111 pcr.set_antialias(cairo.ANTIALIAS_SUBPIXEL)
112 layout = pcr.create_layout()
113 font = pango.FontDescription(cairo_font)
114 layout.set_font_description(font)
115 layout.set_text(text)
116 pcr.update_layout(layout)
117 pcr.show_layout(layout)
108 "rendering as cairo"
109 imgd = img.tostring()
110 a = array.array("B", imgd)
111 surface = cairo.ImageSurface.create_for_data(
112 a, cairo.FORMAT_ARGB32, W, H, W * 4
113 )
114 cr = cairo.Context(surface)
115 cr.move_to(*coords[0])
116 color = ImageColor.getrgb(color)
117 cr.set_source_rgba(color[2] / 255.0, color[1] / 255.0, color[0] / 255.0, 1)
118 if geometry == "LINESTRING" or geometry == "POLYGON":
119 for k in coords:
120 cr.line_to(*k)
121 if geometry == "LINESTRING":
122 if linestring_width is not None:
123 cr.set_line_width(linestring_width)
124 cr.stroke()
125 elif geometry == "POLYGON":
126 cr.fill()
127 elif geometry == "POINT":
128 cr.arc(coords[0][0], coords[0][1], 6, 0, 2 * math.pi)
129 cr.fill()
130 elif geometry == "TEXT" and text and cairo_font:
131 pcr = pangocairo.CairoContext(cr)
132 pcr.set_antialias(cairo.ANTIALIAS_SUBPIXEL)
133 layout = pcr.create_layout()
134 font = pango.FontDescription(cairo_font)
135 layout.set_font_description(font)
136 layout.set_text(text)
137 pcr.update_layout(layout)
138 pcr.show_layout(layout)
118139
119 img = Image.frombuffer("RGBA",( W,H ),surface.get_data(),"raw","RGBA",0,1)
120
140 img = Image.frombuffer("RGBA", (W, H), surface.get_data(), "raw", "RGBA", 0, 1)
141
121142 else:
122 "falling back to PIL"
123 coord = [(int(coord[0]),int(coord[1])) for coord in coords] # PIL dislikes subpixels
124 draw = ImageDraw.Draw(img)
125 if geometry == "LINESTRING":
126 draw.line (coords, fill=color, width=linestring_width)
127 elif geometry == "POINT":
128 draw.ellipse((coords[0][0]-3,coords[0][1]-3,coords[0][0]+3,coords[0][1]+3),fill=color,outline=color)
129 elif geometry == "POLYGON":
130 draw.polygon(coords, fill=color, outline=color)
131 elif geometry == "TEXT" and text and pil_font:
132 font = ImageFont.truetype(pil_font[0], pil_font[1])
133 draw.text(coords[0], text, color, font=font)
143 "falling back to PIL"
144 coord = [
145 (int(coord[0]), int(coord[1])) for coord in coords
146 ] # PIL dislikes subpixels
147 draw = ImageDraw.Draw(img)
148 if geometry == "LINESTRING":
149 draw.line(coords, fill=color, width=linestring_width)
150 elif geometry == "POINT":
151 draw.ellipse(
152 (
153 coords[0][0] - 3,
154 coords[0][1] - 3,
155 coords[0][0] + 3,
156 coords[0][1] + 3,
157 ),
158 fill=color,
159 outline=color,
160 )
161 elif geometry == "POLYGON":
162 draw.polygon(coords, fill=color, outline=color)
163 elif geometry == "TEXT" and text and pil_font:
164 font = ImageFont.truetype(pil_font[0], pil_font[1])
165 draw.text(coords[0], text, color, font=font)
134166 return img
44 # the extent permitted by applicable law. You can redistribute it
55 # and/or modify it under the terms specified in COPYING.
66
7 try:
8 from urllib.request import urlopen
9 except ImportError:
10 from urllib2 import urlopen
7 from urllib.request import urlopen
118 import filecmp
129 import time
1310 import os
1411 import math
1512 import sys
1613 from io import BytesIO
14
1715 try:
1816 from PIL import Image
1917 except ImportError:
3028 thread_responses = {}
3129 zhash_lock = {}
3230
31
3332 def fetch(z, x, y, this_layer):
34 zhash = repr((z,x,y,this_layer))
35 try:
36 zhash_lock[zhash] += 1
37 except KeyError:
38 zhash_lock[zhash] = 1
39 if zhash not in fetching_now:
40 atomthread = threading.Thread(None, threadwrapper, None, (z, x, y, this_layer, zhash))
41 atomthread.start()
42 fetching_now[zhash] = atomthread
43 if fetching_now[zhash].is_alive():
44 fetching_now[zhash].join()
45 resp = thread_responses[zhash]
46 zhash_lock[zhash] -= 1
47 if not zhash_lock[zhash]:
48 del thread_responses[zhash]
49 del fetching_now[zhash]
50 del zhash_lock[zhash]
51 return resp
33 zhash = repr((z, x, y, this_layer))
34 try:
35 zhash_lock[zhash] += 1
36 except KeyError:
37 zhash_lock[zhash] = 1
38 if zhash not in fetching_now:
39 atomthread = threading.Thread(
40 None, threadwrapper, None, (z, x, y, this_layer, zhash)
41 )
42 atomthread.start()
43 fetching_now[zhash] = atomthread
44 if fetching_now[zhash].is_alive():
45 fetching_now[zhash].join()
46 resp = thread_responses[zhash]
47 zhash_lock[zhash] -= 1
48 if not zhash_lock[zhash]:
49 del thread_responses[zhash]
50 del fetching_now[zhash]
51 del zhash_lock[zhash]
52 return resp
5253
53 def threadwrapper(z,x,y,this_layer, zhash):
54 try:
55 thread_responses[zhash] = this_layer["fetch"](z,x,y,this_layer)
56 except OSError:
57 for i in range(20):
58 time.sleep(0.1)
59 try:
60 thread_responses[zhash] = this_layer["fetch"](z,x,y,this_layer)
61 return
62 except OSError:
63 continue
64 thread_responses[zhash] = None
6554
66 def WMS (z, x, y, this_layer):
67 if "max_zoom" in this_layer:
68 if z >= this_layer["max_zoom"]:
69 return None
70 wms = this_layer["remote_url"]
71 req_proj = this_layer.get("wms_proj", this_layer["proj"])
72 width = 384 # using larger source size to rescale better in python
73 height = 384
74 local = config.tiles_cache + this_layer["prefix"] + "/z%s/%s/x%s/%s/y%s."%(z, x/1024, x, y/1024,y)
75 tile_bbox = "bbox=%s,%s,%s,%s" % tuple( projections.from4326(projections.bbox_by_tile(z,x,y,req_proj),req_proj))
55 def threadwrapper(z, x, y, this_layer, zhash):
56 try:
57 thread_responses[zhash] = this_layer["fetch"](z, x, y, this_layer)
58 except OSError:
59 for i in range(20):
60 time.sleep(0.1)
61 try:
62 thread_responses[zhash] = this_layer["fetch"](z, x, y, this_layer)
63 return
64 except OSError:
65 continue
66 thread_responses[zhash] = None
7667
77 wms += tile_bbox + "&width=%s&height=%s&srs=%s"%(width, height, req_proj)
78 if this_layer.get("cached", True):
79 if not os.path.exists("/".join(local.split("/")[:-1])):
80 os.makedirs("/".join(local.split("/")[:-1]))
68
69 def WMS(z, x, y, this_layer):
70 if "max_zoom" in this_layer:
71 if z >= this_layer["max_zoom"]:
72 return None
73 wms = this_layer["remote_url"]
74 req_proj = this_layer.get("wms_proj", this_layer["proj"])
75 width = 384 # using larger source size to rescale better in python
76 height = 384
77 local = (
78 config.tiles_cache
79 + this_layer["prefix"]
80 + "/z%s/%s/x%s/%s/y%s." % (z, x / 1024, x, y / 1024, y)
81 )
82 tile_bbox = "bbox=%s,%s,%s,%s" % tuple(
83 projections.from4326(projections.bbox_by_tile(z, x, y, req_proj), req_proj)
84 )
85
86 wms += tile_bbox + "&width=%s&height=%s&srs=%s" % (width, height, req_proj)
87 if this_layer.get("cached", True):
88 if not os.path.exists("/".join(local.split("/")[:-1])):
89 os.makedirs("/".join(local.split("/")[:-1]))
90 try:
91 os.mkdir(local + "lock")
92 except OSError:
93 for i in range(20):
94 time.sleep(0.1)
95 try:
96 if not os.path.exists(local + "lock"):
97 im = Image.open(local + this_layer["ext"])
98 return im
99 except (IOError, OSError):
100 return None
101 im = Image.open(BytesIO(urlopen(wms).read()))
102 if width != 256 and height != 256:
103 im = im.resize((256, 256), Image.ANTIALIAS)
104 im = im.convert("RGBA")
105
106 if this_layer.get("cached", True):
107 ic = Image.new(
108 "RGBA", (256, 256), this_layer.get("empty_color", config.default_background)
109 )
110 if im.histogram() == ic.histogram():
111 tne = open(local + "tne", "wb")
112 when = time.localtime()
113 tne.write(
114 "%02d.%02d.%04d %02d:%02d:%02d"
115 % (when[2], when[1], when[0], when[3], when[4], when[5])
116 )
117 tne.close()
118 return False
119 im.save(local + this_layer["ext"])
120 os.rmdir(local + "lock")
121 return im
122
123
124 def Tile(z, x, y, this_layer):
125 global OSError, IOError
126 d_tuple = z, x, y
127 if "max_zoom" in this_layer:
128 if z >= this_layer["max_zoom"]:
129 return None
130 if "transform_tile_number" in this_layer:
131 d_tuple = this_layer["transform_tile_number"](z, x, y)
132
133 remote = this_layer["remote_url"] % d_tuple
134 if this_layer.get("cached", True):
135 local = (
136 config.tiles_cache
137 + this_layer["prefix"]
138 + "/z%s/%s/x%s/%s/y%s." % (z, x / 1024, x, y / 1024, y)
139 )
140 if not os.path.exists("/".join(local.split("/")[:-1])):
141 os.makedirs("/".join(local.split("/")[:-1]))
142 try:
143 os.mkdir(local + "lock")
144 except OSError:
145 for i in range(20):
146 time.sleep(0.1)
147 try:
148 if not os.path.exists(local + "lock"):
149 im = Image.open(local + this_layer["ext"])
150 return im
151 except (IOError, OSError):
152 return None
81153 try:
82 os.mkdir(local+"lock")
83 except OSError:
84 for i in range(20):
85 time.sleep(0.1)
154 contents = urlopen(remote).read()
155 im = Image.open(BytesIO(contents))
156 except IOError:
157 if this_layer.get("cached", True):
158 os.rmdir(local + "lock")
159 return False
160 if this_layer.get("cached", True):
161 os.rmdir(local + "lock")
162 open(local + this_layer["ext"], "wb").write(contents)
163 if "dead_tile" in this_layer:
86164 try:
87 if not os.path.exists(local+"lock"):
88 im = Image.open(local + this_layer["ext"])
89 return im
90 except (IOError, OSError):
91 return None
92 im = Image.open(BytesIO(urlopen(wms).read()))
93 if width is not 256 and height is not 256:
94 im = im.resize((256,256),Image.ANTIALIAS)
95 im = im.convert("RGBA")
96
97 if this_layer.get("cached", True):
98 ic = Image.new("RGBA", (256, 256), this_layer.get("empty_color", config.default_background) )
99 if im.histogram() == ic.histogram():
100 tne = open (local+"tne", "wb")
101 when = time.localtime()
102 tne.write("%02d.%02d.%04d %02d:%02d:%02d"%(when[2],when[1],when[0],when[3],when[4],when[5]))
103 tne.close()
104 return False
105 im.save(local+this_layer["ext"])
106 os.rmdir(local+"lock")
107 return im
108
109 def Tile (z, x, y, this_layer):
110 global OSError, IOError
111 d_tuple = z,x,y
112 if "max_zoom" in this_layer:
113 if z >= this_layer["max_zoom"]:
114 return None
115 if "transform_tile_number" in this_layer:
116 d_tuple = this_layer["transform_tile_number"](z,x,y)
117
118 remote = this_layer["remote_url"] % d_tuple
119 if this_layer.get("cached", True):
120 local = config.tiles_cache + this_layer["prefix"] + "/z%s/%s/x%s/%s/y%s."%(z, x/1024, x, y/1024,y)
121 if not os.path.exists("/".join(local.split("/")[:-1])):
122 os.makedirs("/".join(local.split("/")[:-1]))
123 try:
124 os.mkdir(local+"lock")
125 except OSError:
126 for i in range(20):
127 time.sleep(0.1)
128 try:
129 if not os.path.exists(local+"lock"):
130 im = Image.open(local + this_layer["ext"])
131 return im
132 except (IOError, OSError):
133 return None
134 try:
135 contents = urlopen(remote).read()
136 im = Image.open(BytesIO(contents))
137 except IOError:
138 if this_layer.get("cached", True):
139 os.rmdir(local+"lock")
140 return False
141 if this_layer.get("cached", True):
142 os.rmdir(local+"lock")
143 open(local+ this_layer["ext"], "wb").write(contents)
144 if "dead_tile" in this_layer:
145 try:
146 dt = open(this_layer["dead_tile"],"rb").read()
147 if contents == dt:
148 if this_layer.get("cached", True):
149 tne = open (local + "tne", "wb")
150 when = time.localtime()
151 tne.write("%02d.%02d.%04d %02d:%02d:%02d"%(when[2],when[1],when[0],when[3],when[4],when[5]))
152 tne.close()
153 os.remove(local + this_layer["ext"])
154 return False
155 except IOError:
156 pass
157 return im
165 dt = open(this_layer["dead_tile"], "rb").read()
166 if contents == dt:
167 if this_layer.get("cached", True):
168 tne = open(local + "tne", "wb")
169 when = time.localtime()
170 tne.write(
171 "%02d.%02d.%04d %02d:%02d:%02d"
172 % (when[2], when[1], when[0], when[3], when[4], when[5])
173 )
174 tne.close()
175 os.remove(local + this_layer["ext"])
176 return False
177 except IOError:
178 pass
179 return im
1010 import Image, ImageFilter, ImageEnhance, ImageOps
1111
1212 try:
13 import numpy
14 NUMPY_AVAILABLE = True
13 import numpy
14
15 NUMPY_AVAILABLE = True
1516 except ImportError:
16 NUMPY_AVAILABLE = False
17 NUMPY_AVAILABLE = False
1718 import datetime
1819 from twms import getimg
1920
2021 try:
21 import config
22
22 import config
2323 except:
24 pass
24 pass
2525
2626
27 def raster(result_img, filt, bbox = (-999,-999,999,9999), srs="EPSG:3857"):
27 def raster(result_img, filt, bbox=(-999, -999, 999, 9999), srs="EPSG:3857"):
2828 """
2929 Applies various filters to image.
3030 """
3131 for ff in filt:
32 if ff.split(":") == [ff,]:
33 if ff == "bw":
34 result_img = result_img.convert("L")
35 result_img = result_img.convert("RGBA")
32 if ff.split(":") == [ff]:
33 if ff == "bw":
34 result_img = result_img.convert("L")
35 result_img = result_img.convert("RGBA")
3636
37 if ff == "contour":
38 result_img = result_img.filter(ImageFilter.CONTOUR)
39 if ff == "median":
40 result_img = result_img.filter(ImageFilter.MedianFilter(5))
41 if ff == "blur":
42 result_img = result_img.filter(ImageFilter.BLUR)
43 if ff == "edge":
44 result_img = result_img.filter(ImageFilter.EDGE_ENHANCE)
45 if ff == "equalize":
46 result_img = result_img.convert("RGB")
47 result_img = ImageOps.equalize(result_img)
48 result_img = result_img.convert("RGBA")
49 if ff == "autocontrast":
50 result_img = result_img.convert("RGB")
51 result_img = ImageOps.autocontrast(result_img)
52 result_img = result_img.convert("RGBA")
53 if ff == "swaprb":
54 ii = result_img.split()
55 result_img = Image.merge("RGBA", (ii[2],ii[1],ii[0],ii[3]))
37 if ff == "contour":
38 result_img = result_img.filter(ImageFilter.CONTOUR)
39 if ff == "median":
40 result_img = result_img.filter(ImageFilter.MedianFilter(5))
41 if ff == "blur":
42 result_img = result_img.filter(ImageFilter.BLUR)
43 if ff == "edge":
44 result_img = result_img.filter(ImageFilter.EDGE_ENHANCE)
45 if ff == "equalize":
46 result_img = result_img.convert("RGB")
47 result_img = ImageOps.equalize(result_img)
48 result_img = result_img.convert("RGBA")
49 if ff == "autocontrast":
50 result_img = result_img.convert("RGB")
51 result_img = ImageOps.autocontrast(result_img)
52 result_img = result_img.convert("RGBA")
53 if ff == "swaprb":
54 ii = result_img.split()
55 result_img = Image.merge("RGBA", (ii[2], ii[1], ii[0], ii[3]))
5656
57 else:
58 ff, tts = ff.split(":")
59 try:
60 tt = float(tts)
61 except ValueError:
62 tt = 1
63 pass
64 if ff == "brightness":
65 enhancer = ImageEnhance.Brightness(result_img)
66 result_img = enhancer.enhance(tt)
67 if ff == "contrast":
68 enhancer = ImageEnhance.Contrast(result_img)
69 result_img = enhancer.enhance(tt)
70 if ff == "sharpness":
71 enhancer = ImageEnhance.Sharpness(result_img)
72 result_img = enhancer.enhance(tt)
73 if ff == "autocontrast":
74 result_img = result_img.convert("RGB")
75 result_img = ImageOps.autocontrast(result_img, tt)
76 result_img = result_img.convert("RGBA")
77 if ff == "fusion" and NUMPY_AVAILABLE:
78 pix = numpy.array(result_img, dtype=int)
79 a,b = result_img.size
80 pan_img = getimg (bbox, srs, [b, a], config.layers[tts], datetime.datetime.now(), [])
81 pan_img = pan_img.convert("L")
82 print(pix.dtype)
83 print(pix[...,1].shape)
84
85 pan = numpy.array(pan_img)
86
87 pix[...,0] = pix[...,0]*pan/(pix[...,0] + pix[...,1] + pix[...,2])
88 pix[...,1] = pix[...,1]*pan/(pix[...,0] + pix[...,1] + pix[...,2])
89 pix[...,2] = pix[...,2]*pan/(pix[...,0] + pix[...,1] + pix[...,2])
57 else:
58 ff, tts = ff.split(":")
59 try:
60 tt = float(tts)
61 except ValueError:
62 tt = 1
63 pass
64 if ff == "brightness":
65 enhancer = ImageEnhance.Brightness(result_img)
66 result_img = enhancer.enhance(tt)
67 if ff == "contrast":
68 enhancer = ImageEnhance.Contrast(result_img)
69 result_img = enhancer.enhance(tt)
70 if ff == "sharpness":
71 enhancer = ImageEnhance.Sharpness(result_img)
72 result_img = enhancer.enhance(tt)
73 if ff == "autocontrast":
74 result_img = result_img.convert("RGB")
75 result_img = ImageOps.autocontrast(result_img, tt)
76 result_img = result_img.convert("RGBA")
77 if ff == "fusion" and NUMPY_AVAILABLE:
78 pix = numpy.array(result_img, dtype=int)
79 a, b = result_img.size
80 pan_img = getimg(
81 bbox, srs, [b, a], config.layers[tts], datetime.datetime.now(), []
82 )
83 pan_img = pan_img.convert("L")
84 print(pix.dtype)
85 print(pix[..., 1].shape)
9086
91 print(pix.shape)
92 result_img = Image.fromarray(numpy.uint8(pix))
93
94
95 result_img = result_img.convert("RGBA")
87 pan = numpy.array(pan_img)
88
89 pix[..., 0] = (
90 pix[..., 0] * pan / (pix[..., 0] + pix[..., 1] + pix[..., 2])
91 )
92 pix[..., 1] = (
93 pix[..., 1] * pan / (pix[..., 0] + pix[..., 1] + pix[..., 2])
94 )
95 pix[..., 2] = (
96 pix[..., 2] * pan / (pix[..., 0] + pix[..., 1] + pix[..., 2])
97 )
98
99 print(pix.shape)
100 result_img = Image.fromarray(numpy.uint8(pix))
101
102 result_img = result_img.convert("RGBA")
96103 return result_img
97
98
77 import sys, string, bz2, gzip, os
88 from xml.dom import minidom, Node
99
10
1011 class GPXParser:
11 def __init__(self, filename):
12 self.tracks = {}
13 self.pointnum = 0
14 self.trknum = 0
15 self.bbox = (999,999,-999,-999)
16 try:
17 file = open(filename)
18 signature = file.read(2)
19 file.close()
20 file = {
21 "BZ": lambda f: bz2.BZ2File(f),
22 "\x1f\x8b": lambda f: gzip.GzipFile(f),
23 "<?": lambda f: open(f)
24 }[signature](filename)
25 except (OSError, IOError):
26 return
27 try:
28 doc = minidom.parse(file)
29 doc.normalize()
30 except (KeyboardInterrupt, SystemExit):
31 raise
32 except:
33 return # handle this properly later
34 gpx = doc.documentElement
35 for node in gpx.getElementsByTagName('trk'):
36 self.parseTrack(node)
12 def __init__(self, filename):
13 self.tracks = {}
14 self.pointnum = 0
15 self.trknum = 0
16 self.bbox = (999, 999, -999, -999)
17 try:
18 file = open(filename)
19 signature = file.read(2)
20 file.close()
21 file = {
22 "BZ": lambda f: bz2.BZ2File(f),
23 "\x1f\x8b": lambda f: gzip.GzipFile(f),
24 "<?": lambda f: open(f),
25 }[signature](filename)
26 except (OSError, IOError):
27 return
28 try:
29 doc = minidom.parse(file)
30 doc.normalize()
31 except (KeyboardInterrupt, SystemExit):
32 raise
33 except:
34 return # handle this properly later
35 gpx = doc.documentElement
36 for node in gpx.getElementsByTagName("trk"):
37 self.parseTrack(node)
3738
38 def parseTrack(self, trk):
39 #name = trk.getElementsByTagName('name')[0].firstChild.data
40 name = self.trknum
41 self.trknum += 1
42 if not name in self.tracks:
43 self.tracks[name] = {}
44 minlat, minlon, maxlat, maxlon = self.bbox
45 for trkseg in trk.getElementsByTagName('trkseg'):
46 for trkpt in trkseg.getElementsByTagName('trkpt'):
47 lat = float(trkpt.getAttribute('lat'))
48 lon = float(trkpt.getAttribute('lon'))
49 if lat > maxlat:
50 maxlat = lat
51 if lat < minlat:
52 minlat = lat
53 if lon > maxlon:
54 maxlon = lon
55 if lon < minlon:
56 minlon = lon
57 # ele = float(trkpt.getElementsByTagName('ele')[0].firstChild.data)
58 rfc3339 = trkpt.getElementsByTagName('time')[0].firstChild.data
59 self.pointnum += 1
60 self.tracks[name][self.pointnum]={'lat':lat,'lon':lon}
61 self.bbox = (minlon, minlat, maxlon, maxlat)
39 def parseTrack(self, trk):
40 # name = trk.getElementsByTagName('name')[0].firstChild.data
41 name = self.trknum
42 self.trknum += 1
43 if not name in self.tracks:
44 self.tracks[name] = {}
45 minlat, minlon, maxlat, maxlon = self.bbox
46 for trkseg in trk.getElementsByTagName("trkseg"):
47 for trkpt in trkseg.getElementsByTagName("trkpt"):
48 lat = float(trkpt.getAttribute("lat"))
49 lon = float(trkpt.getAttribute("lon"))
50 if lat > maxlat:
51 maxlat = lat
52 if lat < minlat:
53 minlat = lat
54 if lon > maxlon:
55 maxlon = lon
56 if lon < minlon:
57 minlon = lon
58 # ele = float(trkpt.getElementsByTagName('ele')[0].firstChild.data)
59 rfc3339 = trkpt.getElementsByTagName("time")[0].firstChild.data
60 self.pointnum += 1
61 self.tracks[name][self.pointnum] = {"lat": lat, "lon": lon}
62 self.bbox = (minlon, minlat, maxlon, maxlat)
6263
63 def getTrack(self, name):
64 times = self.tracks[name].keys()
65 times.sort()
66 points = [self.tracks[name][time] for time in times]
67 return [(point['lon'],point['lat']) for point in points]
64 def getTrack(self, name):
65 times = self.tracks[name].keys()
66 times.sort()
67 points = [self.tracks[name][time] for time in times]
68 return [(point["lon"], point["lat"]) for point in points]
66
77 from config import *
88 import projections
9
910
1011 def html(ref):
1112 """
1819 resp += wms_name
1920 resp += "</h2><table>"
2021 for i in layers:
21 bbox = layers[i].get("data_bounding_box", projections.projs[layers[i]["proj"]]["bounds"] )
22 resp += "<tr><td><img src=\""
23 resp += ref + "?layers=" + i+ "&amp;bbox=%s,%s,%s,%s" % bbox + "&amp;width=200&amp;format=image/png\" width=\"200\" /></td><td><h3>"
24 resp += layers[i]["name"]
25 resp += "</h3><b>Bounding box:</b> "+ str(bbox) +" (show on <a href=\"http://openstreetmap.org/?minlon=%s&amp;minlat=%s&amp;maxlon=%s&amp;maxlat=%s&amp;box=yes\">OSM</a>" % bbox +")<br />"
26 resp += "<b>Projection:</b> "+ layers[i]["proj"] +"<br />"
27 resp += "<b>WMS half-link:</b> "+ ref + "?layers=" + i + "&amp;<br />"
28 resp += "<b>Tiles URL:</b> "+ ref + "" + i + "/!/!/!." + layers[i].get("ext", "jpg") + "<br />"
29 resp += "</td></tr>"
22 bbox = layers[i].get(
23 "data_bounding_box", projections.projs[layers[i]["proj"]]["bounds"]
24 )
25 resp += '<tr><td><img src="'
26 resp += (
27 ref
28 + "?layers="
29 + i
30 + "&amp;bbox=%s,%s,%s,%s" % bbox
31 + '&amp;width=200&amp;format=image/png" width="200" /></td><td><h3>'
32 )
33 resp += layers[i]["name"]
34 resp += (
35 "</h3><b>Bounding box:</b> "
36 + str(bbox)
37 + ' (show on <a href="http://openstreetmap.org/?minlon=%s&amp;minlat=%s&amp;maxlon=%s&amp;maxlat=%s&amp;box=yes">OSM</a>'
38 % bbox
39 + ")<br />"
40 )
41 resp += "<b>Projection:</b> " + layers[i]["proj"] + "<br />"
42 resp += "<b>WMS half-link:</b> " + ref + "?layers=" + i + "&amp;<br />"
43 resp += (
44 "<b>Tiles URL:</b> "
45 + ref
46 + ""
47 + i
48 + "/!/!/!."
49 + layers[i].get("ext", "jpg")
50 + "<br />"
51 )
52 resp += "</td></tr>"
3053 resp += "</table></body></html>"
3154 return resp
99 import pyproj
1010 except ImportError:
1111 class pyproj:
12
1312 class Proj:
14
1513 def __init__(self, pstring):
1614 self.pstring = pstring
1715
2018 return c1, c2
2119 else:
2220 raise NotImplementedError(
23 "Pyproj is not installed - can't convert between projectios. Install pyproj please.")
21 "Pyproj is not installed - can't convert between projectios. Install pyproj please."
22 )
2423
2524
2625 projs = {
2726 "EPSG:4326": {
28 "proj":
29 pyproj.Proj(
30 "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"),
27 "proj": pyproj.Proj("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"),
3128 "bounds": (-180.0, -90.0, 180.0, 90.0),
3229 },
3330 "NASA:4326": {
34 "proj":
35 pyproj.Proj(
36 "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"),
31 "proj": pyproj.Proj("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"),
3732 "bounds": (-180.0, -166.0, 332.0, 346.0),
3833 },
3934 "EPSG:3395": {
40 "proj":
41 pyproj.Proj(
42 "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"),
35 "proj": pyproj.Proj(
36 "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
37 ),
4338 "bounds": (-180.0, -85.0840591556, 180.0, 85.0840590501),
4439 },
4540 "EPSG:3857": {
46 "proj":
47 pyproj.Proj(
48 "+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs"),
41 "proj": pyproj.Proj(
42 "+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs"
43 ),
4944 "bounds": (-180.0, -85.0511287798, 180.0, 85.0511287798),
5045 },
5146 "EPSG:32635": {
52 "proj":
53 pyproj.Proj(
54 "+proj=utm +zone=35 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"),
47 "proj": pyproj.Proj(
48 "+proj=utm +zone=35 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
49 ),
5550 "bounds": (-180.0, -90.0, 180.0, 90.0),
5651 },
5752 "EPSG:32636": {
58 "proj":
59 pyproj.Proj(
60 "+proj=utm +zone=36 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"),
53 "proj": pyproj.Proj(
54 "+proj=utm +zone=36 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
55 ),
6156 "bounds": (-180.0, -90.0, 180.0, 90.0),
6257 },
6358 "EPSG:32637": {
64 "proj":
65 pyproj.Proj(
66 "+proj=utm +zone=37 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"),
59 "proj": pyproj.Proj(
60 "+proj=utm +zone=37 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
61 ),
6762 "bounds": (-180.0, -90.0, 180.0, 90.0),
6863 },
6964 "EPSG:32638": {
70 "proj":
71 pyproj.Proj(
72 "+proj=utm +zone=38 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"),
65 "proj": pyproj.Proj(
66 "+proj=utm +zone=38 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
67 ),
7368 "bounds": (-180.0, -90.0, 180.0, 90.0),
7469 },
7570 "EPSG:32639": {
76 "proj":
77 pyproj.Proj(
78 "+proj=utm +zone=39 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"),
71 "proj": pyproj.Proj(
72 "+proj=utm +zone=39 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
73 ),
7974 "bounds": (-180.0, -90.0, 180.0, 90.0),
8075 },
8176 "EPSG:32640": {
82 "proj":
83 pyproj.Proj(
84 "+proj=utm +zone=40 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"),
77 "proj": pyproj.Proj(
78 "+proj=utm +zone=40 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
79 ),
8580 "bounds": (-180.0, -90.0, 180.0, 90.0),
8681 },
8782 "EPSG:32641": {
88 "proj":
89 pyproj.Proj(
90 "+proj=utm +zone=41 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"),
83 "proj": pyproj.Proj(
84 "+proj=utm +zone=41 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
85 ),
9186 "bounds": (-180.0, -90.0, 180.0, 90.0),
9287 },
9388 "EPSG:32642": {
94 "proj":
95 pyproj.Proj(
96 "+proj=utm +zone=42 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"),
89 "proj": pyproj.Proj(
90 "+proj=utm +zone=42 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
91 ),
9792 "bounds": (-180.0, -90.0, 180.0, 90.0),
9893 },
9994 }
100 proj_alias = {
101 "EPSG:900913": "EPSG:3857",
102 "EPSG:3785": "EPSG:3857",
103 }
95 proj_alias = {"EPSG:900913": "EPSG:3857", "EPSG:3785": "EPSG:3857"}
10496
10597
10698 def _c4326t3857(t1, t2, lon, lat):
109101 """
110102 lat_rad = math.radians(lat)
111103 xtile = lon * 111319.49079327358
112 ytile = math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / \
113 math.pi * 20037508.342789244
114 return(xtile, ytile)
104 ytile = (
105 math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad)))
106 / math.pi
107 * 20037508.342789244
108 )
109 return (xtile, ytile)
115110
116111
117112 def _c3857t4326(t1, t2, lon, lat):
119114 Pure python 3857 -> 4326 transform. About 12x faster than pyproj.
120115 """
121116 xtile = lon / 111319.49079327358
122 ytile = math.degrees(
123 math.asin(math.tanh(lat / 20037508.342789244 * math.pi)))
124 return(xtile, ytile)
117 ytile = math.degrees(math.asin(math.tanh(lat / 20037508.342789244 * math.pi)))
118 return (xtile, ytile)
125119
126120
127121 def _c4326t3395(t1, t2, lon, lat):
133127 F = 53.5865938
134128 tmp = math.tan(0.78539816339744830962 + math.radians(lat) / 2.0)
135129 pow_tmp = math.pow(
136 math.tan(0.78539816339744830962 +
137 math.asin(E * math.sin(math.radians(lat))) / 2.0),
138 E)
130 math.tan(
131 0.78539816339744830962 + math.asin(E * math.sin(math.radians(lat))) / 2.0
132 ),
133 E,
134 )
139135 x = lon * 111319.49079327358
140136 y = 6378137.0 * math.log(tmp / pow_tmp)
141137 return (x, y)
157153 TOL = 1e-7
158154 i = N_ITER
159155 dphi = 0.1
160 while ((abs(dphi) > TOL) and (i > 0)):
156 while (abs(dphi) > TOL) and (i > 0):
161157 i -= 1
162158 con = eccent * math.sin(Phi)
163 dphi = HALFPI - 2.0 * \
164 math.atan(ts * math.pow((1.0 - con) / (1.0 + con), eccnth)) - Phi
159 dphi = (
160 HALFPI
161 - 2.0 * math.atan(ts * math.pow((1.0 - con) / (1.0 + con), eccnth))
162 - Phi
163 )
165164 Phi += dphi
166165
167166 x = lon / 111319.49079327358
201200 Converts (z,x,y) to coordinates of corner of srs-projected tile
202201 """
203202 z -= 1
204 normalized_tile = (x / (2. ** z), 1. - (y / (2. ** z)))
203 normalized_tile = (x / (2.0 ** z), 1.0 - (y / (2.0 ** z)))
205204 projected_bounds = from4326(projs[proj_alias.get(srs, srs)]["bounds"], srs)
206 maxp = [projected_bounds[2] - projected_bounds[0],
207 projected_bounds[3] - projected_bounds[1]]
208 projected_coords = [(
209 normalized_tile[0] * maxp[0]) + projected_bounds[0],
210 (normalized_tile[1] * maxp[1]) + projected_bounds[1]]
205 maxp = [
206 projected_bounds[2] - projected_bounds[0],
207 projected_bounds[3] - projected_bounds[1],
208 ]
209 projected_coords = [
210 (normalized_tile[0] * maxp[0]) + projected_bounds[0],
211 (normalized_tile[1] * maxp[1]) + projected_bounds[1],
212 ]
211213 return to4326(projected_coords, srs)
212214
213215
223225 projected_bounds = from4326(projs[proj_alias.get(srs, srs)]["bounds"], srs)
224226 point = from4326((lon, lat), srs)
225227 point = [point[0] - projected_bounds[0], point[1] - projected_bounds[1]]
226 # shifting (0,0)
227 maxp = [projected_bounds[2] - projected_bounds[0],
228 projected_bounds[3] - projected_bounds[1]]
229 point = [1. * point[0] / maxp[0], 1. * point[1] / maxp[1]]
230 # normalizing
228 # shifting (0,0)
229 maxp = [
230 projected_bounds[2] - projected_bounds[0],
231 projected_bounds[3] - projected_bounds[1],
232 ]
233 point = [1.0 * point[0] / maxp[0], 1.0 * point[1] / maxp[1]]
234 # normalizing
231235 return point[0] * (2 ** zoom), (1 - point[1]) * (2 ** zoom)
232236
233237
291295
292296 if __name__ == "__main__":
293297 import debug
298
294299 print(_c4326t3857(1, 2, 27.6, 53.2))
295300 print(from4326((27.6, 53.2), "EPSG:3857"))
296301 a = _c4326t3857(1, 2, 27.6, 53.2)
1212 import projections
1313 import sys
1414
15
1516 def reproject(image, bbox, srs_from, srs_to):
1617 out = image.copy()
1718 op = out.load()
1920 bbox_from = projections.from4326(bbox, srs_from)
2021 bbox_to = projections.from4326(bbox, srs_to)
2122 width, height = image.size
22
23 for x in range(0,width):
24 for y in range(0,height):
25 dx = ((1.*x/width)*(bbox_to[2]-bbox_to[0]))+bbox_to[0]
26 dy = ((1.*y/height)*(bbox_to[3]-bbox_to[1]))+bbox_to[1]
27 #print(x, dx, y, dy, file=sys.stderr)
28 dx,dy = tuple(projections.transform([dx,dy], srs_to, srs_from))
29
30 dx = width*(dx - bbox_from[0])/(bbox_from[2]-bbox_from[0])
31 dy = height*(dy - bbox_from[1])/(bbox_from[3]-bbox_from[1])
32 op[x,y] = il[int(dx),int(dy)]
3323
34 #sys.stderr.flush()
24 for x in range(0, width):
25 for y in range(0, height):
26 dx = ((1.0 * x / width) * (bbox_to[2] - bbox_to[0])) + bbox_to[0]
27 dy = ((1.0 * y / height) * (bbox_to[3] - bbox_to[1])) + bbox_to[1]
28 # print(x, dx, y, dy, file=sys.stderr)
29 dx, dy = tuple(projections.transform([dx, dy], srs_to, srs_from))
30
31 dx = width * (dx - bbox_from[0]) / (bbox_from[2] - bbox_from[0])
32 dy = height * (dy - bbox_from[1]) / (bbox_from[3] - bbox_from[1])
33 op[x, y] = il[int(dx), int(dy)]
34
35 # sys.stderr.flush()
3536 return out
36
37
38
99 string = "abcdefghijklmnopqrstuvwxyz012345ABCDEFGHIJKLMNOPQRSTUVWXYZ6789{}"
1010
1111
12 def decode(bbox, sketch):
13 version, sketch = sketch.split(";", 1)
1214
13 def decode(bbox, sketch):
14 version, sketch = sketch.split(";",1)
15
16 def encode_point(bbox, point, length, length_out = None):
17 if not length_out:
18 length_out = length
19 code = "."
20 if not point_is_in(bbox, point):
21 bbox = (-180,-90,180,90)
22 code += "@"
23 length = length_out - 1
24 lon,lat = point
25 lon = (lon-bbox[0])/(bbox[2]-bbox[0]) #normalizing points to bbox
26 lat = (lat-bbox[1])/(bbox[3]-bbox[1])
27 lats, lons = [], []
28
29 for i in range(0,length):
30 latt = int(lat*8)
31 lont = int(lon*8)
32 lat = lat*8 - int(lat*8)
33 lon = lon*8 - int(lon*8)
34 code += string[lont*8+latt]
35 return code
3615
16 def encode_point(bbox, point, length, length_out=None):
17 if not length_out:
18 length_out = length
19 code = "."
20 if not point_is_in(bbox, point):
21 bbox = (-180, -90, 180, 90)
22 code += "@"
23 length = length_out - 1
24 lon, lat = point
25 lon = (lon - bbox[0]) / (bbox[2] - bbox[0]) # normalizing points to bbox
26 lat = (lat - bbox[1]) / (bbox[3] - bbox[1])
27 lats, lons = [], []
28
29 for i in range(0, length):
30 latt = int(lat * 8)
31 lont = int(lon * 8)
32 lat = lat * 8 - int(lat * 8)
33 lon = lon * 8 - int(lon * 8)
34 code += string[lont * 8 + latt]
35 return code
3736
3837
3938 def decode_point(bbox, code):
40 lat,lon = (0,0)
41 if code[0] == ".":
42 code = code[1:]
43 if code[0] == "@":
44 code = code[1:]
45 bbox = (-180,-90,180,90)
46 c = ""
47 code = " " + code #reverse
48 for a in range(0,len(code)-1):
49 c += code[-1]
50 code = code[:-1]
39 lat, lon = (0, 0)
40 if code[0] == ".":
41 code = code[1:]
42 if code[0] == "@":
43 code = code[1:]
44 bbox = (-180, -90, 180, 90)
45 c = ""
46 code = " " + code # reverse
47 for a in range(0, len(code) - 1):
48 c += code[-1]
49 code = code[:-1]
5150
52 code = c
51 code = c
52
53 for t in code:
54 z = string.find(t)
55 lont = int(z / 8.0)
56 latt = (z / 8.0 - int(z / 8.0)) * 8
57 lat += latt
58 lat /= 8.0
59 lon += lont
60 lon /= 8.0
61 lat = lat * (bbox[3] - bbox[1]) + bbox[1]
62 lon = lon * (bbox[2] - bbox[0]) + bbox[0]
63 return lon, lat
5364
5465
55 for t in code:
56 z = string.find(t)
57 lont = int(z/8.)
58 latt = (z/8. - int(z/8.))*8
59 lat += latt
60 lat /= 8.
61 lon += lont
62 lon /= 8.
63 lat = lat*(bbox[3]-bbox[1])+bbox[1]
64 lon = lon*(bbox[2]-bbox[0])+bbox[0]
65 return lon, lat
66
67 #Window.alert(code_point((0,0,0,0), 53.11, 27.3434))
68 #print(decode_point((0,0,0,0), ".@aaa"))
66 # Window.alert(code_point((0,0,0,0), 53.11, 27.3434))
67 # print(decode_point((0,0,0,0), ".@aaa"))
2525
2626 config_path = "/etc/twms/twms.conf"
2727 if os.path.exists(config_path):
28 try:
29 config = imp.load_source("twms.config", config_path)
30 except:
31 config = imp.load_source("config", config_path)
28 try:
29 config = imp.load_source("twms.config", config_path)
30 except:
31 config = imp.load_source("config", config_path)
3232 else:
33 try:
34 config_path = os.path.join(os.path.dirname(__file__), "twms.conf")
35 config = imp.load_source("twms.config", config_path)
36 except:
37 config_path = os.path.join(os.path.realpath(sys.path[0]), "twms.conf")
38 config = imp.load_source("config", os.path.join(os.path.realpath(sys.path[0]), "twms.conf"))
39 sys.stderr.write("Configuration file not found, using defaults from %s\n" % config_path)
40 sys.stderr.flush()
33 try:
34 config_path = os.path.join(os.path.dirname(__file__), "twms.conf")
35 config = imp.load_source("twms.config", config_path)
36 except:
37 config_path = os.path.join(os.path.realpath(sys.path[0]), "twms.conf")
38 config = imp.load_source(
39 "config", os.path.join(os.path.realpath(sys.path[0]), "twms.conf")
40 )
41 sys.stderr.write(
42 "Configuration file not found, using defaults from %s\n" % config_path
43 )
44 sys.stderr.flush()
4145
4246
4347 import correctify
4448 import capabilities
4549 import fetchers
46 #import config
50
51 # import config
4752 import bbox
4853 import bbox as bbox_utils
4954 import projections
5358 from reproject import reproject
5459
5560 try:
56 import psyco
57 psyco.full()
61 import psyco
62
63 psyco.full()
5864 except ImportError:
59 pass
60
61
65 pass
6266
6367
6468 OK = 200
6569 ERROR = 500
6670
67 cached_objs = {} # a dict. (layer, z, x, y): PIL image
71 cached_objs = {} # a dict. (layer, z, x, y): PIL image
6872 cached_hist_list = []
6973
7074 formats = {
7781
7882 mimetypes = dict(zip(formats.values(), formats.keys()))
7983
84
8085 def twms_main(data):
8186 """
8287 Do main TWMS work.
8994 content_type = "text/html"
9095 resp = ""
9196 srs = data.get("srs", "EPSG:4326")
92 gpx = data.get("gpx","").split(",")
93 if gpx == ['']:
94 gpx = []
95 wkt = data.get("wkt","")
96 trackblend = float(data.get("trackblend","0.5"))
97 color = data.get("color",data.get("colour","")).split(",")
97 gpx = data.get("gpx", "").split(",")
98 if gpx == [""]:
99 gpx = []
100 wkt = data.get("wkt", "")
101 trackblend = float(data.get("trackblend", "0.5"))
102 color = data.get("color", data.get("colour", "")).split(",")
98103 track = False
99104 tracks = []
100105 if len(gpx) == 0:
101 req_bbox = projections.from4326((27.6518898,53.8683186,27.6581944,53.8720359), srs)
106 req_bbox = projections.from4326(
107 (27.6518898, 53.8683186, 27.6581944, 53.8720359), srs
108 )
102109 else:
103 for g in gpx:
104 local_gpx = config.gpx_cache + "%s.gpx" % g
105 if not os.path.exists (config.gpx_cache):
106 os.makedirs(config.gpx_cache)
107 if not os.path.exists (local_gpx):
108 urllib.urlretrieve ("http://www.openstreetmap.org/trace/%s/data" % g, local_gpx)
109 if not track:
110 track = GPXParser(local_gpx)
111 req_bbox = projections.from4326(track.bbox, srs)
112 else:
113 track = GPXParser(local_gpx)
114 req_bbox = bbox.add(req_bbox, projections.from4326(track.bbox, srs))
115 tracks.append(track)
116
117 req_type = data.get("request","GetMap")
118 version = data.get("version","1.1.1")
119 ref = data.get("ref",config.service_url)
110 for g in gpx:
111 local_gpx = config.gpx_cache + "%s.gpx" % g
112 if not os.path.exists(config.gpx_cache):
113 os.makedirs(config.gpx_cache)
114 if not os.path.exists(local_gpx):
115 urllib.urlretrieve(
116 "http://www.openstreetmap.org/trace/%s/data" % g, local_gpx
117 )
118 if not track:
119 track = GPXParser(local_gpx)
120 req_bbox = projections.from4326(track.bbox, srs)
121 else:
122 track = GPXParser(local_gpx)
123 req_bbox = bbox.add(req_bbox, projections.from4326(track.bbox, srs))
124 tracks.append(track)
125
126 req_type = data.get("request", "GetMap")
127 version = data.get("version", "1.1.1")
128 ref = data.get("ref", config.service_url)
120129 if req_type == "GetCapabilities":
121 content_type, resp = capabilities.get(version, ref)
122 return (OK, content_type, resp)
123
124 layer = data.get("layers",config.default_layers).split(",")
130 content_type, resp = capabilities.get(version, ref)
131 return (OK, content_type, resp)
132
133 layer = data.get("layers", config.default_layers).split(",")
125134 if ("layers" in data) and not layer[0]:
126 layer = ["transparent"]
135 layer = ["transparent"]
127136
128137 if req_type == "GetCorrections":
129 points = data.get("points",data.get("POINTS", "")).split("=")
130 resp = ""
131 points = [a.split(",") for a in points]
132 points = [(float(a[0]), float(a[1])) for a in points]
133
134 req.content_type = "text/plain"
135 for lay in layer:
136 for point in points:
137 resp += "%s,%s;"% tuple(correctify.rectify(config.layers[lay], point))
138 resp += "\n"
139 return (OK, content_type, resp)
140
138 points = data.get("points", data.get("POINTS", "")).split("=")
139 resp = ""
140 points = [a.split(",") for a in points]
141 points = [(float(a[0]), float(a[1])) for a in points]
142
143 req.content_type = "text/plain"
144 for lay in layer:
145 for point in points:
146 resp += "%s,%s;" % tuple(correctify.rectify(config.layers[lay], point))
147 resp += "\n"
148 return (OK, content_type, resp)
141149
142150 force = data.get("force", "")
143151 if force != "":
149157 filt = filt.split(",")
150158 filt = tuple(filt)
151159
152 if layer == [""] :
153 content_type = "text/html"
154 resp = overview.html(ref)
155 return (OK, content_type, resp)
160 if layer == [""]:
161 content_type = "text/html"
162 resp = overview.html(ref)
163 return (OK, content_type, resp)
156164
157165 format = data.get("format", config.default_format).lower()
158166 format = formats.get("image/" + format, format)
159167 format = formats.get(format, format)
160168 if format not in formats.values():
161 return (ERROR, content_type, "Invalid format")
169 return (ERROR, content_type, "Invalid format")
162170 content_type = mimetypes[format]
163171
164 width=0
165 height=0
166 resp_cache_path, resp_ext = "",""
172 width = 0
173 height = 0
174 resp_cache_path, resp_ext = "", ""
167175 if req_type == "GetTile":
168 width=256
169 height=256
170 height = int(data.get("height",height))
171 width = int(data.get("width",width))
172 srs = data.get("srs", "EPSG:3857")
173 x = int(data.get("x",0))
174 y = int(data.get("y",0))
175 z = int(data.get("z",1)) + 1
176 if "cache_tile_responses" in dir(config) and not wkt and (len(gpx) == 0):
177 if (srs, tuple(layer), filt, width, height, force, format) in config.cache_tile_responses:
178
179 resp_cache_path, resp_ext = config.cache_tile_responses[(srs, tuple(layer), filt, width, height, force, format)]
180 resp_cache_path = resp_cache_path+"/%s/%s/%s.%s"%(z-1,x,y,resp_ext)
181 if os.path.exists(resp_cache_path):
182 return (OK, content_type, open(resp_cache_path, "r").read())
183 if len(layer) == 1:
184 if layer[0] in config.layers:
185 if config.layers[layer[0]]["proj"] == srs and width is 256 and height is 256 and not filt and not force and not correctify.has_corrections(config.layers[layer[0]]):
186 local = config.tiles_cache + config.layers[layer[0]]["prefix"] + "/z%s/%s/x%s/%s/y%s."%(z, x/1024, x, y/1024,y)
187 ext = config.layers[layer]["ext"]
188 adds = ["","ups."]
189 for add in adds:
190 if os.path.exists(local+add+ext):
191 tile_file = open(local+add+ext, "r")
192 resp = tile_file.read()
193 return (OK, content_type, resp)
194 req_bbox = projections.from4326(projections.bbox_by_tile(z,x,y,srs),srs)
195
196 if data.get("bbox",None):
197 req_bbox = tuple(map(float,data.get("bbox",req_bbox).split(",")))
176 width = 256
177 height = 256
178 height = int(data.get("height", height))
179 width = int(data.get("width", width))
180 srs = data.get("srs", "EPSG:3857")
181 x = int(data.get("x", 0))
182 y = int(data.get("y", 0))
183 z = int(data.get("z", 1)) + 1
184 if "cache_tile_responses" in dir(config) and not wkt and (len(gpx) == 0):
185 if (
186 srs,
187 tuple(layer),
188 filt,
189 width,
190 height,
191 force,
192 format,
193 ) in config.cache_tile_responses:
194
195 resp_cache_path, resp_ext = config.cache_tile_responses[
196 (srs, tuple(layer), filt, width, height, force, format)
197 ]
198 resp_cache_path = resp_cache_path + "/%s/%s/%s.%s" % (
199 z - 1,
200 x,
201 y,
202 resp_ext,
203 )
204 if os.path.exists(resp_cache_path):
205 return (OK, content_type, open(resp_cache_path, "r").read())
206 if len(layer) == 1:
207 if layer[0] in config.layers:
208 if (
209 config.layers[layer[0]]["proj"] == srs
210 and width is 256
211 and height is 256
212 and not filt
213 and not force
214 and not correctify.has_corrections(config.layers[layer[0]])
215 ):
216 local = (
217 config.tiles_cache
218 + config.layers[layer[0]]["prefix"]
219 + "/z%s/%s/x%s/%s/y%s." % (z, x / 1024, x, y / 1024, y)
220 )
221 ext = config.layers[layer]["ext"]
222 adds = ["", "ups."]
223 for add in adds:
224 if os.path.exists(local + add + ext):
225 tile_file = open(local + add + ext, "r")
226 resp = tile_file.read()
227 return (OK, content_type, resp)
228 req_bbox = projections.from4326(projections.bbox_by_tile(z, x, y, srs), srs)
229
230 if data.get("bbox", None):
231 req_bbox = tuple(map(float, data.get("bbox", req_bbox).split(",")))
198232
199233 req_bbox = projections.to4326(req_bbox, srs)
200234
201235 req_bbox, flip_h = bbox.normalize(req_bbox)
202236 box = req_bbox
203 #print(req_bbox, file=sys.stderr)
204 #sys.stderr.flush()
205
206 height = int(data.get("height",height))
207 width = int(data.get("width",width))
237 # print(req_bbox, file=sys.stderr)
238 # sys.stderr.flush()
239
240 height = int(data.get("height", height))
241 width = int(data.get("width", width))
208242 width = min(width, config.max_width)
209243 height = min(height, config.max_height)
210244 if (width == 0) and (height == 0):
211 width = 350
212
213
214 # layer = layer.split(",")
215
216 imgs = 1.
245 width = 350
246
247 # layer = layer.split(",")
248
249 imgs = 1.0
217250 ll = layer.pop(0)
218251 if ll[-2:] == "!c":
219 ll = ll[:-2]
220 if wkt:
221 wkt = ","+wkt
222 wkt = correctify.corr_wkt(config.layers[ll]) + wkt
223 srs = config.layers[ll]["proj"]
252 ll = ll[:-2]
253 if wkt:
254 wkt = "," + wkt
255 wkt = correctify.corr_wkt(config.layers[ll]) + wkt
256 srs = config.layers[ll]["proj"]
224257 try:
225 result_img = getimg(box,srs, (height, width), config.layers[ll], start_time, force)
258 result_img = getimg(
259 box, srs, (height, width), config.layers[ll], start_time, force
260 )
226261 except KeyError:
227 result_img = Image.new("RGBA", (width,height))
228
229
230 #width, height = result_img.size
262 result_img = Image.new("RGBA", (width, height))
263
264 # width, height = result_img.size
231265 for ll in layer:
232 if ll[-2:] == "!c":
233 ll = ll[:-2]
234 if wkt:
235 wkt = ","+wkt
236 wkt = correctify.corr_wkt(config.layers[ll]) + wkt
237 srs = config.layers[ll]["proj"]
238
239 im2 = getimg(box, srs,(height, width), config.layers[ll], start_time, force)
240
241 if "empty_color" in config.layers[ll]:
242 ec = ImageColor.getcolor(config.layers[ll]["empty_color"], "RGBA")
243 sec = set(ec)
244 if "empty_color_delta" in config.layers[ll]:
245 delta = config.layers[ll]["empty_color_delta"]
246 for tr in range(-delta, delta):
247 for tg in range(-delta, delta):
248 for tb in range(-delta, delta):
249 if (ec[0]+tr) >= 0 and (ec[0]+tr) < 256 and (ec[1]+tr) >= 0 and (ec[1]+tr) < 256 and(ec[2]+tr) >= 0 and (ec[2]+tr) < 256:
250 sec.add((ec[0]+tr,ec[1]+tg,ec[2]+tb,ec[3]))
251 i2l = im2.load()
252 for x in range(0,im2.size[0]):
253 for y in range(0,im2.size[1]):
254 t = i2l[x,y]
255 if t in sec:
256 i2l[x,y] = (t[0],t[1],t[2],0)
257 if not im2.size == result_img.size:
258 im2 = im2.resize(result_img.size, Image.ANTIALIAS)
259 im2 = Image.composite(im2,result_img, im2.split()[3]) # imgs/(imgs+1.))
260
261 if "noblend" in force:
262 result_img = im2
263 else:
264 result_img = Image.blend(im2, result_img, 0.5)
265 imgs += 1.
266 if ll[-2:] == "!c":
267 ll = ll[:-2]
268 if wkt:
269 wkt = "," + wkt
270 wkt = correctify.corr_wkt(config.layers[ll]) + wkt
271 srs = config.layers[ll]["proj"]
272
273 im2 = getimg(box, srs, (height, width), config.layers[ll], start_time, force)
274
275 if "empty_color" in config.layers[ll]:
276 ec = ImageColor.getcolor(config.layers[ll]["empty_color"], "RGBA")
277 sec = set(ec)
278 if "empty_color_delta" in config.layers[ll]:
279 delta = config.layers[ll]["empty_color_delta"]
280 for tr in range(-delta, delta):
281 for tg in range(-delta, delta):
282 for tb in range(-delta, delta):
283 if (
284 (ec[0] + tr) >= 0
285 and (ec[0] + tr) < 256
286 and (ec[1] + tr) >= 0
287 and (ec[1] + tr) < 256
288 and (ec[2] + tr) >= 0
289 and (ec[2] + tr) < 256
290 ):
291 sec.add((ec[0] + tr, ec[1] + tg, ec[2] + tb, ec[3]))
292 i2l = im2.load()
293 for x in range(0, im2.size[0]):
294 for y in range(0, im2.size[1]):
295 t = i2l[x, y]
296 if t in sec:
297 i2l[x, y] = (t[0], t[1], t[2], 0)
298 if not im2.size == result_img.size:
299 im2 = im2.resize(result_img.size, Image.ANTIALIAS)
300 im2 = Image.composite(im2, result_img, im2.split()[3]) # imgs/(imgs+1.))
301
302 if "noblend" in force:
303 result_img = im2
304 else:
305 result_img = Image.blend(im2, result_img, 0.5)
306 imgs += 1.0
266307
267308 ##Applying filters
268309 result_img = filter.raster(result_img, filt, req_bbox, srs)
269310
270 #print(wkt, file=sys.stderr)
271 #sys.stderr.flush()
311 # print(wkt, file=sys.stderr)
312 # sys.stderr.flush()
272313 if wkt:
273 result_img = drawing.wkt(wkt, result_img, req_bbox, srs, color if len(color) > 0 else None, trackblend)
314 result_img = drawing.wkt(
315 wkt,
316 result_img,
317 req_bbox,
318 srs,
319 color if len(color) > 0 else None,
320 trackblend,
321 )
274322 if len(gpx) > 0:
275 last_color = None
276 c = iter(color)
277 for track in tracks:
323 last_color = None
324 c = iter(color)
325 for track in tracks:
326 try:
327 last_color = c.next()
328 except StopIteration:
329 pass
330 result_img = drawing.gpx(
331 track, result_img, req_bbox, srs, last_color, trackblend
332 )
333
334 if flip_h:
335 result_img = ImageOps.flip(result_img)
336 image_content = BytesIO()
337
338 if format == "JPEG":
278339 try:
279 last_color = c.next();
280 except StopIteration:
281 pass
282 result_img = drawing.gpx(track, result_img, req_bbox, srs, last_color, trackblend)
283
284 if flip_h:
285 result_img = ImageOps.flip(result_img)
286 image_content = BytesIO()
287
288 if format == "JPEG":
289 try:
290 result_img.save(image_content, format, quality=config.output_quality, progressive=config.output_progressive)
291 except IOError:
292 result_img.save(image_content, format, quality=config.output_quality)
340 result_img.save(
341 image_content,
342 format,
343 quality=config.output_quality,
344 progressive=config.output_progressive,
345 )
346 except IOError:
347 result_img.save(image_content, format, quality=config.output_quality)
293348 elif format == "PNG":
294 result_img.save(image_content, format, progressive=config.output_progressive, optimize =config.output_optimize)
349 result_img.save(
350 image_content,
351 format,
352 progressive=config.output_progressive,
353 optimize=config.output_optimize,
354 )
295355 elif format == "GIF":
296 result_img.save(image_content, format, quality=config.output_quality, progressive=config.output_progressive)
297 else: ## workaround for GIF
298 result_img = result_img.convert("RGB")
299 result_img.save(image_content, format, quality=config.output_quality, progressive=config.output_progressive)
356 result_img.save(
357 image_content,
358 format,
359 quality=config.output_quality,
360 progressive=config.output_progressive,
361 )
362 else: ## workaround for GIF
363 result_img = result_img.convert("RGB")
364 result_img.save(
365 image_content,
366 format,
367 quality=config.output_quality,
368 progressive=config.output_progressive,
369 )
300370 resp = image_content.getvalue()
301371 if resp_cache_path:
302 try:
303 "trying to create local cache directory, if it doesn't exist"
304 os.makedirs("/".join(resp_cache_path.split("/")[:-1]))
305 except OSError:
306 pass
307 try:
308 a = open(resp_cache_path, "w")
309 a.write(resp)
310 a.close()
311 except (OSError, IOError):
312 print("error saving response answer to file %s." % (resp_cache_path), file=sys.stderr)
313 sys.stderr.flush()
372 try:
373 "trying to create local cache directory, if it doesn't exist"
374 os.makedirs("/".join(resp_cache_path.split("/")[:-1]))
375 except OSError:
376 pass
377 try:
378 a = open(resp_cache_path, "w")
379 a.write(resp)
380 a.close()
381 except (OSError, IOError):
382 print(
383 "error saving response answer to file %s." % (resp_cache_path),
384 file=sys.stderr,
385 )
386 sys.stderr.flush()
314387
315388 return (OK, content_type, resp)
316389
317390
318 def tile_image (layer, z, x, y, start_time, again=False, trybetter = True, real = False):
319 """
391 def tile_image(layer, z, x, y, start_time, again=False, trybetter=True, real=False):
392 """
320393 Returns asked image.
321394 again - is this a second pass on this tile?
322395 trybetter - should we try to combine this tile from better ones?
323396 real - should we return the tile even in not good quality?
324397 """
325 x = x % (2 ** (z-1))
326 if y<0 or y >= (2 ** (z-1)):
327 return None
328 if not bbox.bbox_is_in(projections.bbox_by_tile(z,x,y,layer["proj"]), layer.get("data_bounding_box",config.default_bbox), fully=False):
329 return None
330 global cached_objs, cached_hist_list
331 if "prefix" in layer:
332 if (layer["prefix"], z, x, y) in cached_objs:
333 return cached_objs[(layer["prefix"], z, x, y)]
334 if layer.get("cached", True):
335 local = config.tiles_cache + layer["prefix"] + "/z%s/%s/x%s/%s/y%s."%(z, x/1024, x, y/1024,y)
336 ext = layer["ext"]
337 if "cache_ttl" in layer:
338 for ex in [ext, "dsc."+ext, "ups."+ext, "tne"]:
339 f = local+ex
340 if os.path.exists(f):
341 if (os.stat(f).st_mtime < (time.time()-layer["cache_ttl"])):
342 os.remove(f)
343
344 gpt_image = False
345 try:
346 "trying to create local cache directory, if it doesn't exist"
347 os.makedirs("/".join(local.split("/")[:-1]))
348 except OSError:
349 pass
350 if not os.path.exists(local+"tne") and not os.path.exists(local+"lock"):
351 if os.path.exists(local+ext): # First, look for tile in cache
398 x = x % (2 ** (z - 1))
399 if y < 0 or y >= (2 ** (z - 1)):
400 return None
401 if not bbox.bbox_is_in(
402 projections.bbox_by_tile(z, x, y, layer["proj"]),
403 layer.get("data_bounding_box", config.default_bbox),
404 fully=False,
405 ):
406 return None
407 global cached_objs, cached_hist_list
408 if "prefix" in layer:
409 if (layer["prefix"], z, x, y) in cached_objs:
410 return cached_objs[(layer["prefix"], z, x, y)]
411 if layer.get("cached", True):
412 local = (
413 config.tiles_cache
414 + layer["prefix"]
415 + "/z%s/%s/x%s/%s/y%s." % (z, x / 1024, x, y / 1024, y)
416 )
417 ext = layer["ext"]
418 if "cache_ttl" in layer:
419 for ex in [ext, "dsc." + ext, "ups." + ext, "tne"]:
420 f = local + ex
421 if os.path.exists(f):
422 if os.stat(f).st_mtime < (time.time() - layer["cache_ttl"]):
423 os.remove(f)
424
425 gpt_image = False
352426 try:
353 im1 = Image.open(local+ext)
354 im1.is_ok = True
355 return im1
356 except IOError:
357 if os.path.exists(local+"lock"):
358 return None
359 else:
360 os.remove(local+ext) # # Cached tile is broken - remove it
361
362
363 if layer["scalable"] and (z<layer.get("max_zoom", config.default_max_zoom)) and trybetter: # Second, try to glue image of better ones
364 if os.path.exists(local+"ups."+ext):
365 try:
366 im = Image.open(local+"ups."+ext)
367 im.is_ok = True
368 return im
369 except IOError:
370 pass
371 ec = ImageColor.getcolor(layer.get("empty_color", config.default_background), "RGBA")
372 ec = (ec[0],ec[1],ec[2],0)
373 im = Image.new("RGBA", (512, 512), ec)
374 im1 = tile_image(layer, z+1,x*2,y*2, start_time)
375 if im1:
376 im2 = tile_image(layer, z+1,x*2+1,y*2, start_time)
377 if im2:
378 im3 = tile_image(layer, z+1,x*2,y*2+1, start_time)
379 if im3:
380 im4 = tile_image(layer, z+1,x*2+1,y*2+1, start_time)
381 if im4:
382 im.paste(im1,(0,0))
383 im.paste(im2,(256,0))
384 im.paste(im3,(0,256))
385 im.paste(im4,(256,256))
386 im = im.resize((256,256),Image.ANTIALIAS)
387 if layer.get("cached", True):
388 try:
389 im.save(local+"ups."+ext)
390 except IOError:
391 pass
392 im.is_ok = True
427 "trying to create local cache directory, if it doesn't exist"
428 os.makedirs("/".join(local.split("/")[:-1]))
429 except OSError:
430 pass
431 if not os.path.exists(local + "tne") and not os.path.exists(local + "lock"):
432 if os.path.exists(local + ext): # First, look for tile in cache
433 try:
434 im1 = Image.open(local + ext)
435 im1.is_ok = True
436 return im1
437 except IOError:
438 if os.path.exists(local + "lock"):
439 return None
440 else:
441 os.remove(local + ext) # # Cached tile is broken - remove it
442
443 if (
444 layer["scalable"]
445 and (z < layer.get("max_zoom", config.default_max_zoom))
446 and trybetter
447 ): # Second, try to glue image of better ones
448 if os.path.exists(local + "ups." + ext):
449 try:
450 im = Image.open(local + "ups." + ext)
451 im.is_ok = True
452 return im
453 except IOError:
454 pass
455 ec = ImageColor.getcolor(
456 layer.get("empty_color", config.default_background), "RGBA"
457 )
458 ec = (ec[0], ec[1], ec[2], 0)
459 im = Image.new("RGBA", (512, 512), ec)
460 im1 = tile_image(layer, z + 1, x * 2, y * 2, start_time)
461 if im1:
462 im2 = tile_image(layer, z + 1, x * 2 + 1, y * 2, start_time)
463 if im2:
464 im3 = tile_image(layer, z + 1, x * 2, y * 2 + 1, start_time)
465 if im3:
466 im4 = tile_image(
467 layer, z + 1, x * 2 + 1, y * 2 + 1, start_time
468 )
469 if im4:
470 im.paste(im1, (0, 0))
471 im.paste(im2, (256, 0))
472 im.paste(im3, (0, 256))
473 im.paste(im4, (256, 256))
474 im = im.resize((256, 256), Image.ANTIALIAS)
475 if layer.get("cached", True):
476 try:
477 im.save(local + "ups." + ext)
478 except IOError:
479 pass
480 im.is_ok = True
481 return im
482 if not again:
483 if "fetch" in layer:
484 delta = datetime.datetime.now() - start_time
485 delta = delta.seconds + delta.microseconds / 1000000.0
486 if (config.deadline > delta) or (z < 4):
487 im = fetchers.fetch(z, x, y, layer) # Try fetching from outside
488 if im:
489 im.is_ok = True
490 return im
491 if real and (z > 1):
492 im = tile_image(
493 layer,
494 z - 1,
495 int(x / 2),
496 int(y / 2),
497 start_time,
498 again=False,
499 trybetter=False,
500 real=True,
501 )
502 if im:
503 im = im.crop(
504 (
505 128 * (x % 2),
506 128 * (y % 2),
507 128 * (x % 2) + 128,
508 128 * (y % 2) + 128,
509 )
510 )
511 im = im.resize((256, 256), Image.BILINEAR)
512 im.is_ok = False
393513 return im
394 if not again:
514 else:
395515 if "fetch" in layer:
396 delta = (datetime.datetime.now() - start_time)
397 delta = delta.seconds + delta.microseconds/1000000.
398 if (config.deadline > delta) or (z < 4):
399 im = fetchers.fetch(z,x,y,layer) # Try fetching from outside
400 if im:
401 im.is_ok = True
402 return im
403 if real and (z>1):
404 im = tile_image(layer, z-1, int(x/2), int(y/2), start_time, again=False, trybetter=False, real=True)
405 if im:
406 im = im.crop((128 * (x % 2), 128 * (y % 2), 128 * (x % 2) + 128, 128 * (y % 2) + 128))
407 im = im.resize((256,256), Image.BILINEAR)
408 im.is_ok = False
409 return im
410 else:
411 if "fetch" in layer:
412 delta = (datetime.datetime.now() - start_time)
413 delta = delta.seconds + delta.microseconds/1000000.
414 if (config.deadline > delta) or (z < 4):
415 im = fetchers.fetch(z,x,y,layer) # Try fetching from outside
416 if im:
417 im.is_ok = True
418 return im
419
420 def getimg (bbox, request_proj, size, layer, start_time, force):
421 orig_bbox = bbox
422 ## Making 4-corner maximal bbox
423 bbox_p = projections.from4326(bbox, request_proj)
424 bbox_p = projections.to4326((bbox_p[2],bbox_p[1],bbox_p[0],bbox_p[3]), request_proj)
425
426 bbox_4 = ( (bbox_p[2],bbox_p[3]),(bbox[0], bbox[1]),(bbox_p[0],bbox_p[1]),(bbox[2],bbox[3]) )
427 if "nocorrect" not in force:
428 bb4 = []
429 for point in bbox_4:
430 bb4.append(correctify.rectify(layer, point))
431 bbox_4 = bb4
432 bbox = bbox_utils.expand_to_point(bbox, bbox_4)
433 #print(bbox)
434 #print(orig_bbox)
435
436
437 global cached_objs
438 H,W = size
439
440 max_zoom = layer.get("max_zoom",config.default_max_zoom)
441 min_zoom = layer.get("min_zoom",1)
442
443 zoom = bbox_utils.zoom_for_bbox (bbox,size,layer, min_zoom, max_zoom, (config.max_height,config.max_width))
444 lo1, la1, lo2, la2 = bbox
445 from_tile_x, from_tile_y, to_tile_x, to_tile_y = projections.tile_by_bbox(bbox, zoom, layer["proj"])
446 cut_from_x = int(256*(from_tile_x - int(from_tile_x)))
447 cut_from_y = int(256*(from_tile_y - int(from_tile_y)))
448 cut_to_x = int(256*(to_tile_x - int(to_tile_x)))
449 cut_to_y = int(256*(to_tile_y - int(to_tile_y)))
450
451 from_tile_x, from_tile_y = int(from_tile_x), int(from_tile_y)
452 to_tile_x, to_tile_y = int(to_tile_x), int(to_tile_y)
453 bbox_im = (cut_from_x, cut_to_y, 256*(to_tile_x-from_tile_x)+cut_to_x, 256*(from_tile_y-to_tile_y)+cut_from_y )
454 x = 256*(to_tile_x-from_tile_x+1)
455 y = 256*(from_tile_y-to_tile_y+1)
456 #print(x, y, file=sys.stderr)
457 #sys.stderr.flush()
458 out = Image.new("RGBA", (x, y))
459 for x in range (from_tile_x, to_tile_x+1):
460 for y in range (to_tile_y, from_tile_y+1):
461 got_image = False
462 im1 = tile_image (layer,zoom,x,y, start_time, real = True)
463 if im1:
464 if "prefix" in layer:
465 if (layer["prefix"], zoom, x, y) not in cached_objs:
466 if im1.is_ok:
467 cached_objs[(layer["prefix"], zoom, x, y)] = im1
468 cached_hist_list.append((layer["prefix"], zoom, x, y))
469 #print((layer["prefix"], zoom, x, y), cached_objs[(layer["prefix"], zoom, x, y)], file=sys.stderr)
470 #sys.stderr.flush()
471 if len(cached_objs) >= config.max_ram_cached_tiles:
472 del cached_objs[cached_hist_list.pop(0)]
473 #print("Removed tile from cache", file=sys.stderr)
474 #sys.stderr.flush()
475 else:
476 ec = ImageColor.getcolor(layer.get("empty_color", config.default_background), "RGBA")
477 #ec = (ec[0],ec[1],ec[2],0)
478 im1 = Image.new("RGBA", (256, 256), ec)
479
480
481 out.paste(im1,((x - from_tile_x)*256, (-to_tile_y + y )*256,))
482 if "filter" in layer:
483 out = filter.raster(out, layer["filter"], orig_bbox, request_proj)
484
485 ## TODO: Here's a room for improvement. we could drop this crop in case user doesn't need it.
486 out = out.crop(bbox_im)
487 if "noresize" not in force:
488 if (H == W) and (H == 0):
489 W, H = out.size
490 if H == 0:
491 H = out.size[1]*W // out.size[0]
492 if W == 0:
493 W = out.size[0]*H // out.size[1]
494 #bbox = orig_bbox
495 quad = []
496 trans_needed = False
497 for point in bbox_4:
498 x = (point[0]-bbox[0])/(bbox[2]-bbox[0])*(out.size[0])
499 y = (1-(point[1]-bbox[1])/(bbox[3]-bbox[1]))*(out.size[1])
500 x = int(round(x))
501 y = int(round(y))
502 if (x is not 0 and x is not out.size[0]) or (y is not 0 and y is not out.size[1]):
503 trans_needed = True
504 quad.append(x)
505 quad.append(y)
506
507 if trans_needed:
508 quad = tuple(quad)
509 out = out.transform((W,H), Image.QUAD, quad, Image.BICUBIC)
510 elif (W != out.size[0]) or (H != out.size[1]):
511 "just resize"
512 out = out.resize((W,H), Image.ANTIALIAS)
513 # out = reproject(out, bbox, layer["proj"], request_proj)
514 return out
515
516
517
518 import filter
516 delta = datetime.datetime.now() - start_time
517 delta = delta.seconds + delta.microseconds / 1000000.0
518 if (config.deadline > delta) or (z < 4):
519 im = fetchers.fetch(z, x, y, layer) # Try fetching from outside
520 if im:
521 im.is_ok = True
522 return im
523
524
525 def getimg(bbox, request_proj, size, layer, start_time, force):
526 orig_bbox = bbox
527 ## Making 4-corner maximal bbox
528 bbox_p = projections.from4326(bbox, request_proj)
529 bbox_p = projections.to4326(
530 (bbox_p[2], bbox_p[1], bbox_p[0], bbox_p[3]), request_proj
531 )
532
533 bbox_4 = (
534 (bbox_p[2], bbox_p[3]),
535 (bbox[0], bbox[1]),
536 (bbox_p[0], bbox_p[1]),
537 (bbox[2], bbox[3]),
538 )
539 if "nocorrect" not in force:
540 bb4 = []
541 for point in bbox_4:
542 bb4.append(correctify.rectify(layer, point))
543 bbox_4 = bb4
544 bbox = bbox_utils.expand_to_point(bbox, bbox_4)
545 # print(bbox)
546 # print(orig_bbox)
547
548 global cached_objs
549 H, W = size
550
551 max_zoom = layer.get("max_zoom", config.default_max_zoom)
552 min_zoom = layer.get("min_zoom", 1)
553
554 zoom = bbox_utils.zoom_for_bbox(
555 bbox, size, layer, min_zoom, max_zoom, (config.max_height, config.max_width)
556 )
557 lo1, la1, lo2, la2 = bbox
558 from_tile_x, from_tile_y, to_tile_x, to_tile_y = projections.tile_by_bbox(
559 bbox, zoom, layer["proj"]
560 )
561 cut_from_x = int(256 * (from_tile_x - int(from_tile_x)))
562 cut_from_y = int(256 * (from_tile_y - int(from_tile_y)))
563 cut_to_x = int(256 * (to_tile_x - int(to_tile_x)))
564 cut_to_y = int(256 * (to_tile_y - int(to_tile_y)))
565
566 from_tile_x, from_tile_y = int(from_tile_x), int(from_tile_y)
567 to_tile_x, to_tile_y = int(to_tile_x), int(to_tile_y)
568 bbox_im = (
569 cut_from_x,
570 cut_to_y,
571 256 * (to_tile_x - from_tile_x) + cut_to_x,
572 256 * (from_tile_y - to_tile_y) + cut_from_y,
573 )
574 x = 256 * (to_tile_x - from_tile_x + 1)
575 y = 256 * (from_tile_y - to_tile_y + 1)
576 # print(x, y, file=sys.stderr)
577 # sys.stderr.flush()
578 out = Image.new("RGBA", (x, y))
579 for x in range(from_tile_x, to_tile_x + 1):
580 for y in range(to_tile_y, from_tile_y + 1):
581 got_image = False
582 im1 = tile_image(layer, zoom, x, y, start_time, real=True)
583 if im1:
584 if "prefix" in layer:
585 if (layer["prefix"], zoom, x, y) not in cached_objs:
586 if im1.is_ok:
587 cached_objs[(layer["prefix"], zoom, x, y)] = im1
588 cached_hist_list.append((layer["prefix"], zoom, x, y))
589 # print((layer["prefix"], zoom, x, y), cached_objs[(layer["prefix"], zoom, x, y)], file=sys.stderr)
590 # sys.stderr.flush()
591 if len(cached_objs) >= config.max_ram_cached_tiles:
592 del cached_objs[cached_hist_list.pop(0)]
593 # print("Removed tile from cache", file=sys.stderr)
594 # sys.stderr.flush()
595 else:
596 ec = ImageColor.getcolor(
597 layer.get("empty_color", config.default_background), "RGBA"
598 )
599 # ec = (ec[0],ec[1],ec[2],0)
600 im1 = Image.new("RGBA", (256, 256), ec)
601
602 out.paste(im1, ((x - from_tile_x) * 256, (-to_tile_y + y) * 256))
603 if "filter" in layer:
604 out = filter.raster(out, layer["filter"], orig_bbox, request_proj)
605
606 ## TODO: Here's a room for improvement. we could drop this crop in case user doesn't need it.
607 out = out.crop(bbox_im)
608 if "noresize" not in force:
609 if (H == W) and (H == 0):
610 W, H = out.size
611 if H == 0:
612 H = out.size[1] * W // out.size[0]
613 if W == 0:
614 W = out.size[0] * H // out.size[1]
615 # bbox = orig_bbox
616 quad = []
617 trans_needed = False
618 for point in bbox_4:
619 x = (point[0] - bbox[0]) / (bbox[2] - bbox[0]) * (out.size[0])
620 y = (1 - (point[1] - bbox[1]) / (bbox[3] - bbox[1])) * (out.size[1])
621 x = int(round(x))
622 y = int(round(y))
623 if (x != 0 and x != out.size[0]) or (y != 0 and y != out.size[1]):
624 trans_needed = True
625 quad.append(x)
626 quad.append(y)
627
628 if trans_needed:
629 quad = tuple(quad)
630 out = out.transform((W, H), Image.QUAD, quad, Image.BICUBIC)
631 elif (W != out.size[0]) or (H != out.size[1]):
632 "just resize"
633 out = out.resize((W, H), Image.ANTIALIAS)
634 # out = reproject(out, bbox, layer["proj"], request_proj)
635 return out