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))