xinclude: added pic and latests plugins
[treecutter.git] / xinclude / address.py
1 #!/usr/bin/python
2 # -*- coding: utf-8 -*-
3
4 import time
5 import sys
6 from os import path
7 from httplib2 import Http
8 from urllib import urlencode
9 from math import log,tan,pi,cos,radians,ceil,floor
10 from lxml import etree
11 from lxml.builder import ElementMaker
12 from PIL import Image, ImageDraw
13 from StringIO import StringIO
14 from geohash import encode
15 from geographiclib.geodesic import Geodesic
16 from treecutter import constants as const
17
18 # EC Equator lenght
19 EC = 40075016.686 # meter
20 # ER Earth radius
21 ER = 6372798.2 # meter
22 # Availible zoom levels in OSM
23 ZOOMRANGE = range(1, 18)
24 # tile size
25 TS = 256
26
27 h = Http(".cache")
28
29 class Coord(object):
30     def __init__(self, lat, lon):
31         self.latitude  = float(lat)
32         self.longitude = float(lon)
33         self.image = None
34
35     def osmlink(self):
36         return "http://www.openstreetmap.org/?mlat=%s&mlon=%s&zoom=18&layers=M"\
37             % (self.latitude,self.longitude)
38
39     def dms(self):
40         ns = self.latitude
41         ew = self.longitude
42         mnt,sec = divmod(ns*3600,60)
43         deg,mnt = divmod(mnt,60)
44         out = u'''%d°%2d'%5.2f"%s''' % ( deg,mnt,sec,'N')
45         mnt,sec = divmod(ew*3600,60)
46         deg,mnt = divmod(mnt,60)
47         out +=  u''' %d°%2d'%05.2f"%s''' % ( deg,mnt,sec,'E')
48         return out
49
50     def lontile(self, zoom):
51         tile = ((self.longitude + 180) / 360) * (2**zoom)
52         return tile
53
54     def lattile(self, zoom):
55         rad = radians(self.latitude)
56         tile = (1-log(tan(rad)+1/cos(rad))/pi)/2*2**zoom
57         return tile
58
59     def offset(self,zoom):
60         x = self.lontile(zoom)
61         y = self.lattile(zoom)
62         xo = int(floor((x-floor(x))*TS))
63         yo = int(floor((y-floor(y))*TS))
64         return (xo, yo)
65
66     def distance(self, point):
67         res = Geodesic.WGS84.Inverse(self.latitude, self.longitude,
68                                      point.latitude, point.longitude)
69         return res['s12']
70
71     def direct(self, direction, lenght):
72         point = Geodesic.WGS84.Direct(self.latitude, self.longitude,
73                                       direction, length)
74         return self.__class__(point['lat2'],point['lon2'])
75
76     def png(self,zoom=15,size=(400,150)):
77         filename = encode(self.latitude, self.longitude)+'.png'
78 #        if path.isfile(filename):
79 #            if path.getctime(filename) > time.time() - 60*60*24*2:
80 #                return
81         im = Image.new("RGB", size, None)
82
83         ew = int(self.lontile(zoom))
84         ns = int(self.lattile(zoom))
85
86         (xo, yo) = self.offset(zoom)
87         et = int(floor((xo - ceil(size[0]/2))/TS))
88         nt = int(floor((yo - ceil(size[1]/2))/TS))
89         wt = int(floor((xo + ceil(size[0]/2))/TS))
90         st = int(floor((yo + ceil(size[1]/2))/TS))
91
92         lontiles =  range(ew+et,ew+wt+1)
93         lattiles = range(ns+nt,ns+st+1)
94         imsize = (len(lontiles)*TS,len(lattiles)*TS)
95         grid = Image.new("RGB", imsize, None)
96         for yi, y in enumerate(lattiles):
97             for xi, x in enumerate(lontiles):
98                 url = 'http://tile.openstreetmap.org/%s/%s/%s.png' % (zoom,x,y)
99                 request, content = h.request(url)
100                 img = Image.open(StringIO(content))
101 #                dr = ImageDraw.Draw(img)
102 #                dr.rectangle([0,0,TS,TS], outline=0)
103                 box = (xi*TS, yi*TS)
104                 grid.paste(img, box)
105
106         yp = [i for i,j in enumerate(lattiles) if j == int(ns)][0]*TS+yo
107         xp = [i for i,j in enumerate(lontiles) if j == int(ew)][0]*TS+xo
108         mark(grid, (xp,yp))
109         xc = int(ceil(size[0]/2))
110         yc = int(ceil(size[1]/2))
111
112 #        draw = ImageDraw.Draw(grid)
113 #        draw.rectangle([xp-xc,yp-yc,xp+xc,yp+yc], outline="red")
114         gridc = grid.crop((xp-xc,yp-yc,xp+xc,yp+yc))
115         gridc.save(filename)
116
117     def db_xml(self):
118         self.png()
119         img = encode(self.latitude, self.longitude)+'.png'
120         phr = "geo:"+str(self.latitude)+","+str(self.longitude)
121
122         db = ElementMaker(namespace=const.DB_NS, nsmap=const.NSMAP)
123         uri = db.uri(db.link(
124                 db.inlinemediaobject(
125                     db.imageobject(db.imagedata(
126                             fileref=img,
127                             format='PNG')),
128                     db.textobject(db.phrase(phr))
129                     ),
130                 self.dms(),
131                 **{const.XLINK+"href": self.osmlink()}))
132         return uri
133
134 def c(s):
135     return s or ''
136
137 def s(s,n):
138     if n is not None:
139         return s or n.text
140     return s
141
142
143 class Address(object):
144     """Address object to contain everything known about an address"""
145     def __init__(self,street=None,postcode=None,city=None,country=None):
146         self.street = street
147         self.postcode = postcode
148         self.city = city
149         self.country = country
150         self.country_code = None
151         self.phone = None
152         self.name = None
153         self.coord = None
154
155
156     def geocode(self, language='en'):
157         base_url = 'http://nominatim.openstreetmap.org/search?%s'
158         params = { 'addressdetails': 1,
159                    'limit': 1,
160                    'format': 'xml',
161                    'polygon': 0,
162                    'accept-language': language}
163
164         if self.country:
165             t = etree.parse('/usr/share/xml/iso-codes/iso_3166.xml')
166             r = t.xpath('//iso_3166_entry[@name="'+self.country+'"]')
167             if len(r)==1:
168                 self.country_code = r[0].get("alpha_2_code")
169         if self.country_code:
170             params['countrycodes'] = self.country_code
171
172         addrlist=[]
173         addrlist.append(u''+c(self.name)+', '+c(self.street)+', '+c(self.city))
174         addrlist.append(u''+c(self.street)+', '+c(self.postcode)+', '+c(self.city))
175         addrlist.append(u''+c(self.street)+', '+c(self.city))
176
177         for addr in addrlist:
178             params['q'] = addr.encode('utf-8')
179             url = base_url % urlencode(params)
180             time.sleep(1)
181             resp, content = h.request(url)
182             root = etree.fromstring(content)
183             places = int(root.xpath('count(//place[@place_id])'))
184             if places == 1:
185                 place = root.find("place")
186 #                print etree.tostring(place,encoding='UTF-8',pretty_print=True)
187                 self.postcode = s(self.postcode,place.find("postcode"))
188                 self.city = s(self.city,place.find("city"))
189                 self.country = s(self.country,place.find("country"))
190                 self.country_code = s(self.country_code,place.find("country_code"))
191                 self.coord=Coord(place.get("lat"),place.get("lon"))
192                 return
193
194     def add_phone(self, phone):
195         self.phone = phone
196
197     def set_name(self, name):
198         self.name = name
199
200     def db_xml(self):
201         db = ElementMaker(namespace=const.DB_NS, nsmap=const.NSMAP)
202         co = ''
203         if self.coord:
204             co = self.coord.db_xml()
205         adr = db.address(db.street(c(self.street)),
206                          db.postcode(c(self.postcode)),
207                          db.city(c(self.city)),
208                          db.country(c(self.country)),
209                          db.phone(c(self.phone)),
210                          co)
211 #                       type=self.type,
212         return adr
213
214
215 def distance(p1, p2):
216     res = Geodesic.WGS84.Inverse(p1[0], p1[1], p2[0], p2[1])
217     return res['s12']
218
219 def mark(image, coord):
220     draw = ImageDraw.Draw(image)
221     x, y = coord
222     r = 5
223     bbox = (x-r, y-r, x+r, y+r)
224     draw.ellipse(bbox, outline="red")
225
226 def box(image,box):
227     draw = ImageDraw.Draw(image)
228     draw.rectangle(box, outline="red")
229
230 def mapimages(coords, zoom=15,size=(TS,TS)):
231
232         minlat = min(coords,key=attrgetter('latitude'))
233         maxlat = max(coords,key=attrgetter('latitude'))
234         minlon = min(coords,key=attrgetter('longitude'))
235         maxlon = max(coords,key=attrgetter('longitude'))
236
237         # Find minimal bounding box and expand it 5%
238         hyp = distance((maxlat,minlon),(minlat,maxlon))
239         hyp = hyp*0.05
240         uld =  Geodesic.WGS84.Direct(maxlat, minlon, 315, hyp)
241         urd =  Geodesic.WGS84.Direct(maxlat, maxlon,  45, hyp)
242         lld =  Geodesic.WGS84.Direct(minlat, minlon, 225, hyp)
243         lrd =  Geodesic.WGS84.Direct(minlat, maxlon, 135, hyp)
244
245         ul = Coord(maxlat,minlon).direct(315, hyp)
246
247         ul = (uld['lat2'],uld['lon2'])
248         ur = (urd['lat2'],urd['lon2'])
249         ll = (lld['lat2'],lld['lon2'])
250         lr = (lrd['lat2'],lrd['lon2'])
251         top = distance(ul,ur)
252         bottom = distance(ll,lr)
253         left = distance(ul,ll)
254         right = distance(ur,lr)
255 #    m_per_pix = EC*abs(cos(lat))/2**(zoomlevel+8)
256 #    m_per_pix = ER*2*pi*abs(cos(lat))/2**(zoomlevel+8)
257         for zoom in range(1,18):
258             t = 2**(zoom)*TS/(ER*2*pi*abs(cos(ul[0])))
259             toppix = t*top
260             leftpix = t*left
261             b = 2**(zoom)*TS/(ER*2*pi*abs(cos(ll[0])))
262             bottompix = b*bottom
263             rightpix  = b*right
264 #            print "Zoom: %s : %s %s %s %s" % (zoom, toppix, bottompix, leftpix, rightpix)
265             if max(toppix,bottompix,leftpix,rightpix) > TS*4:
266                 break
267 #        print "Zoom to use : %s" % (zoom)
268         ult = coordtile(ul,zoom)
269         lrt = coordtile(lr,zoom)
270         lattiles = range(int(ult[0]),int(lrt[0])+1)
271         lontiles =  range(int(ult[1]),int(lrt[1])+1)
272         size = (len(lontiles)*TS,len(lattiles)*TS)
273         grid = Image.new("RGB", size, None)
274         img = []
275         for yi, y in enumerate(lattiles):
276             for xi, x in enumerate(lontiles):
277                 url = tile(x,y,zoom)
278 #                print url
279                 time.sleep(1)
280                 request, content = h.request(url)
281                 img = Image.open(StringIO(content))
282                 box = (xi*TS, yi*TS)
283                 grid.paste(img, box)
284         for c in coords:
285             t = coordtile(c,zoom)
286             o = offset(c,zoom)
287             yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
288             xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
289             mark(grid, (xp,yp))
290
291         t = coordtile(ul,zoom)
292         o = offset(ul,zoom)
293         yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
294         xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
295         t = coordtile(lr,zoom)
296         o = offset(lr,zoom)
297         ypl = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
298         xpl = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
299         gridc = grid.crop((xp,yp,xpl,ypl))
300         gridc.show()
301 #        gridc.save("cap-un.png")
302
303 if __name__ == "__main__":
304     for arg in sys.argv[1:]:
305         al = arg.split("=")
306         if al[0] == "lang":
307             lang = al[1]
308         if al[0] == "xptr":
309             argument = al[1].decode('utf-8')
310
311     addrlist = argument.split(',')
312     addrfmt = "street,postcode,city,country"
313     adict = addrfmt.split(',')
314     argdict = dict(zip(adict,addrlist))
315     addr = Address(**argdict)
316     addr.geocode()
317     axml = addr.db_xml()
318 #    clean_db(axml)
319
320     #print(etree.tostring(cxml, pretty_print=True))
321     #sys.stdout.write(out.encode('utf-8'))
322     sys.stdout.write(etree.tostring(axml,encoding='UTF-8',pretty_print=False))