2 # -*- coding: utf-8 -*-
6 from httplib2 import Http
10 from PIL import Image, ImageDraw
11 from StringIO import StringIO
12 from geohash import encode
13 from geographiclib.geodesic import Geodesic
16 EC = 40075016.686 # meter
18 ER = 6372798.2 # meter
19 # Availible zoom levels in OSM
20 ZOOMRANGE = range(1, 18)
27 class StaticOpenStreetMap(object):
28 """Setting up a static image from OSM with one or more markers"""
34 self._maptype = 'mapnik'
36 def add_marker(lon,lat):
37 self._markers.append((lon,lat))
39 # def construct_map(self):
40 # for coord in self._markers:
43 def lontile(lon, zoom):
44 tile = ((lon + 180) / 360) * (2**zoom)
47 def lattile(lat, zoom):
48 tile = (1 - log(tan(lat * pi/180) + 1 / cos(lat* pi/180)) / pi) /2 * 2**zoom
49 # tile = (1-log(tan(radians(lat))+1/cos(radians(lat)))/pi)/2*2**zoom
52 def coordtile(coord, zoom):
53 x = lontile(coord[1],zoom)
54 y = lattile(coord[0],zoom)
58 return 'http://tile.openstreetmap.org/%s/%s/%s.png' % (z,x,y)
62 x = int(floor(lontile(coord[1],z)))
63 y = int(floor(lattile(coord[0],z)))
66 def offset(coord, zoom):
68 x = lontile(coord[1],z)
69 y = lattile(coord[0],z)
70 xo = int(floor((x-floor(x))*TS))
71 yo = int(floor((y-floor(y))*TS))
75 res = Geodesic.WGS84.Inverse(p1[0], p1[1], p2[0], p2[1])
78 def geocode(address,country=None):
79 params = { 'q': address,
86 params['countrycodes'] = country
88 base_url = 'http://nominatim.openstreetmap.org/search?%s'
89 url = base_url % urllib.urlencode(params)
90 resp, content = h.request(url)
93 root = etree.fromstring(content)
95 # print "%s" % (address)
98 for element in root.iter("place"):
99 lat = element.get("lat")
100 lon = element.get("lon")
101 # print(" : %s,%s" % (lat, lon))
103 if not lat or not lon:
106 return (float(lat),float(lon))
108 def mark(image, coord):
109 draw = ImageDraw.Draw(image)
112 bbox = (x-r, y-r, x+r, y+r)
113 draw.ellipse(bbox, outline="red")
115 def mapimage(coords, zoom=15,size=(TS,TS)):
118 filename = encode(co[0],co[1])+'.png'
119 if path.isfile(filename):
120 if path.getctime(filename) > time.time() - 60*60*24*2:
122 im = Image.new("RGB", size, None)
123 center = (size[0]/2, size[1]/2)
124 nst = ceil(center[0]-TS/2)/TS
125 ewt = ceil(center[1]-TS/2)/TS
126 x, y = offset(co,zoom)
127 (ns, ew) = coordtile(co,zoom)
140 ul = (ns-nst-n, ew-ewt-e)
141 lr = (ns+nst+s, ew+ewt+w)
142 lattiles = range(int(ul[0]),int(lr[0])+1)
143 lontiles = range(int(ul[1]),int(lr[1])+1)
145 size = (len(lontiles)*TS,len(lattiles)*TS)
146 grid = Image.new("RGB", size, None)
148 for yi, y in enumerate(lattiles):
149 for xi, x in enumerate(lontiles):
153 request, content = h.request(url)
154 img = Image.open(StringIO(content))
157 t = coordtile(co,zoom)
159 yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
160 xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
162 gridc = grid.crop((xp-TS/2,yp-TS/2,xp+TS/2,yp+TS/2))
170 minlat = min(minlat,c[0])
171 maxlat = max(maxlat,c[0])
172 minlon = min(minlon,c[1])
173 maxlon = max(maxlon,c[1])
174 # Find minimal bounding box and expand it 5%
175 hyp = distance((maxlat,minlon),(minlat,maxlon))
177 uld = Geodesic.WGS84.Direct(maxlat, minlon, 315, hyp)
178 urd = Geodesic.WGS84.Direct(maxlat, maxlon, 45, hyp)
179 lld = Geodesic.WGS84.Direct(minlat, minlon, 225, hyp)
180 lrd = Geodesic.WGS84.Direct(minlat, maxlon, 135, hyp)
182 ul = (uld['lat2'],uld['lon2'])
183 ur = (urd['lat2'],urd['lon2'])
184 ll = (lld['lat2'],lld['lon2'])
185 lr = (lrd['lat2'],lrd['lon2'])
186 top = distance(ul,ur)
187 bottom = distance(ll,lr)
188 left = distance(ul,ll)
189 right = distance(ur,lr)
190 # m_per_pix = EC*abs(cos(lat))/2**(zoomlevel+8)
191 # m_per_pix = ER*2*pi*abs(cos(lat))/2**(zoomlevel+8)
192 for zoom in range(1,18):
193 t = 2**(zoom)*TS/(ER*2*pi*abs(cos(ul[0])))
196 b = 2**(zoom)*TS/(ER*2*pi*abs(cos(ll[0])))
199 # print "Zoom: %s : %s %s %s %s" % (zoom, toppix, bottompix, leftpix, rightpix)
200 if max(toppix,bottompix,leftpix,rightpix) > TS*4:
202 # print "Zoom to use : %s" % (zoom)
203 ult = coordtile(ul,zoom)
204 lrt = coordtile(lr,zoom)
205 lattiles = range(int(ult[0]),int(lrt[0])+1)
206 lontiles = range(int(ult[1]),int(lrt[1])+1)
207 size = (len(lontiles)*TS,len(lattiles)*TS)
208 grid = Image.new("RGB", size, None)
210 for yi, y in enumerate(lattiles):
211 for xi, x in enumerate(lontiles):
215 request, content = h.request(url)
216 img = Image.open(StringIO(content))
220 t = coordtile(c,zoom)
222 yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
223 xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
226 t = coordtile(ul,zoom)
228 yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
229 xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
230 t = coordtile(lr,zoom)
232 ypl = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
233 xpl = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
234 gridc = grid.crop((xp,yp,xpl,ypl))
236 # gridc.save("cap-un.png")