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)
24 DB_NS="http://docbook.org/ns/docbook"
26 XI_NS="http://www.w3.org/2001/XInclude"
28 XLINK_NS="http://www.w3.org/1999/xlink"
29 XLINK = "{%s}" % XLINK_NS
30 HTML_NS="http://www.w3.org/1999/xhtml"
31 HTML = "{%s}" % HTML_NS
32 NSMAP = {None : DB_NS,
39 def __init__(self, lat, lon):
40 self.latitude = float(lat)
41 self.longitude = float(lon)
45 return "http://www.openstreetmap.org/?mlat=%s&mlon=%s&zoom=18&layers=M"\
46 % (self.latitude,self.longitude)
51 mnt,sec = divmod(ns*3600,60)
52 deg,mnt = divmod(mnt,60)
53 out = u'''%d°%2d'%5.2f"%s''' % ( deg,mnt,sec,'N')
54 mnt,sec = divmod(ew*3600,60)
55 deg,mnt = divmod(mnt,60)
56 out += u''' %d°%2d'%05.2f"%s''' % ( deg,mnt,sec,'E')
59 def lontile(lon, zoom):
60 tile = ((lon + 180) / 360) * (2**zoom)
63 def lattile(lat, zoom):
64 tile = (1 - log(tan(lat * pi/180) + 1 / cos(lat* pi/180)) / pi) /2 * 2**zoom
65 # tile = (1-log(tan(radians(lat))+1/cos(radians(lat)))/pi)/2*2**zoom
68 def coordtile(coord, zoom):
69 x = lontile(coord[1],zoom)
70 y = lattile(coord[0],zoom)
74 return 'http://tile.openstreetmap.org/%s/%s/%s.png' % (z,x,y)
77 (x, y) = coordtile(coord,zoom)
82 def offset(coord, zoom):
83 (x, y) = coordtile(coord,zoom)
84 xo = int(floor((x-floor(x))*TS))
85 yo = int(floor((y-floor(y))*TS))
88 def distance(self, point):
89 res = Geodesic.WGS84.Inverse(self.latitude, self.longitude,
90 point.latitude, point.longitude)
93 def png(self,zoom=15,size=(TS,TS)):
94 filename = encode(self.latitude, self.longitude)+'.png'
95 if path.isfile(filename):
96 if path.getctime(filename) > time.time() - 60*60*24*2:
98 im = Image.new("RGB", size, None)
99 center = (size[0]/2, size[1]/2)
100 nst = ceil(center[0]-TS/2)/TS
101 ewt = ceil(center[1]-TS/2)/TS
102 x, y = offset(co,zoom)
103 (ns, ew) = coordtile(co,zoom)
116 ul = (ns-nst-n, ew-ewt-e)
117 lr = (ns+nst+s, ew+ewt+w)
118 lattiles = range(int(ul[0]),int(lr[0])+1)
119 lontiles = range(int(ul[1]),int(lr[1])+1)
121 size = (len(lontiles)*TS,len(lattiles)*TS)
122 grid = Image.new("RGB", size, None)
123 for yi, y in enumerate(lattiles):
124 for xi, x in enumerate(lontiles):
127 request, content = h.request(url)
128 img = Image.open(StringIO(content))
131 t = coordtile(co,zoom)
133 yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
134 xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
136 gridc = grid.crop((xp-TS/2,yp-TS/2,xp+TS/2,yp+TS/2))
140 uri = etree.Element(DB+'uri',nsmap=NSMAP)
141 ln = etree.SubElement(uri, DB+'link')
142 ln.set(XLINK+'href',self.osmlink())
143 imo = etree.SubElement(ln, DB+'inlinemediaobject')
144 io = etree.SubElement(imo, DB+'imageobject', condition="web")
145 idat = etree.SubElement(io , DB+'imagedata',
146 fileref=encode(self.latitude, self.longitude)+'.png',
148 to = etree.SubElement(imo, DB+'textobject')
149 ph = etree.SubElement(to, DB+'phrase')
150 ph.text = "geo:"+str(self.latitude)+","+str(self.longitude)
151 para = etree.SubElement(ln, DB+'para')
152 para.text = self.dms()
155 class Address(object):
156 """Address object to contain everything known about an address"""
157 def __init__(self,address):
158 self._address_string = address
159 self._root = etree.Element(DB+'address',nsmap=NSMAP)
162 def geocode(self,country=None):
163 params = { 'q': self._address_string,
170 params['countrycodes'] = country
172 base_url = 'http://nominatim.openstreetmap.org/search?%s'
173 url = base_url % urllib.urlencode(params)
174 resp, content = h.request(url)
175 root = etree.fromstring(content)
176 place = root.find("place")
177 if place is not None:
178 print (etree.tostring(root, pretty_print=True))
179 self._coord=Coord(place.get("lat"),place.get("lon"))
187 return self._coord.db_xml()
190 def distance(p1, p2):
191 res = Geodesic.WGS84.Inverse(p1[0], p1[1], p2[0], p2[1])
194 def mark(image, coord):
195 draw = ImageDraw.Draw(image)
198 bbox = (x-r, y-r, x+r, y+r)
199 draw.ellipse(bbox, outline="red")
201 def mapimages(coords, zoom=15,size=(TS,TS)):
207 minlat = min(minlat,c[0])
208 maxlat = max(maxlat,c[0])
209 minlon = min(minlon,c[1])
210 maxlon = max(maxlon,c[1])
211 # Find minimal bounding box and expand it 5%
212 hyp = distance((maxlat,minlon),(minlat,maxlon))
214 uld = Geodesic.WGS84.Direct(maxlat, minlon, 315, hyp)
215 urd = Geodesic.WGS84.Direct(maxlat, maxlon, 45, hyp)
216 lld = Geodesic.WGS84.Direct(minlat, minlon, 225, hyp)
217 lrd = Geodesic.WGS84.Direct(minlat, maxlon, 135, hyp)
219 ul = (uld['lat2'],uld['lon2'])
220 ur = (urd['lat2'],urd['lon2'])
221 ll = (lld['lat2'],lld['lon2'])
222 lr = (lrd['lat2'],lrd['lon2'])
223 top = distance(ul,ur)
224 bottom = distance(ll,lr)
225 left = distance(ul,ll)
226 right = distance(ur,lr)
227 # m_per_pix = EC*abs(cos(lat))/2**(zoomlevel+8)
228 # m_per_pix = ER*2*pi*abs(cos(lat))/2**(zoomlevel+8)
229 for zoom in range(1,18):
230 t = 2**(zoom)*TS/(ER*2*pi*abs(cos(ul[0])))
233 b = 2**(zoom)*TS/(ER*2*pi*abs(cos(ll[0])))
236 # print "Zoom: %s : %s %s %s %s" % (zoom, toppix, bottompix, leftpix, rightpix)
237 if max(toppix,bottompix,leftpix,rightpix) > TS*4:
239 # print "Zoom to use : %s" % (zoom)
240 ult = coordtile(ul,zoom)
241 lrt = coordtile(lr,zoom)
242 lattiles = range(int(ult[0]),int(lrt[0])+1)
243 lontiles = range(int(ult[1]),int(lrt[1])+1)
244 size = (len(lontiles)*TS,len(lattiles)*TS)
245 grid = Image.new("RGB", size, None)
247 for yi, y in enumerate(lattiles):
248 for xi, x in enumerate(lontiles):
252 request, content = h.request(url)
253 img = Image.open(StringIO(content))
257 t = coordtile(c,zoom)
259 yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
260 xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
263 t = coordtile(ul,zoom)
265 yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
266 xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
267 t = coordtile(lr,zoom)
269 ypl = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
270 xpl = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
271 gridc = grid.crop((xp,yp,xpl,ypl))
273 # gridc.save("cap-un.png")