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     return s or n
 139
 140
 141 class Address(object):
 142     """Address object to contain everything known about an address"""
 143     def __init__(self,street=None,postcode=None,city=None,country=None):
 144         self.street = street
 145         self.postcode = postcode
 146         self.city = city
 147         self.country = country
 148         self.country_code = None
 149         self.phone = None
 150         self.name = None
 151         self.coord = None
 152
 153
 154     def geocode(self, language='en'):
 155         base_url = 'http://nominatim.openstreetmap.org/search?%s'
 156         params = { 'addressdetails': 1,
 157                    'limit': 1,
 158                    'format': 'xml',
 159                    'polygon': 0,
 160                    'accept-language': language}
 161
 162         if self.country:
 163             t = etree.parse('/usr/share/xml/iso-codes/iso_3166.xml')
 164             r = t.xpath('//iso_3166_entry[@name="'+self.country+'"]')
 165             if len(r)==1:
 166                 self.country_code = r[0].get("alpha_2_code")
 167         if self.country_code:
 168             params['countrycodes'] = self.country_code
 169
 170         addrlist=[]
 171         addrlist.append(u''+c(self.name)+', '+c(self.street)+', '+c(self.city))
 172         addrlist.append(u''+c(self.street)+', '+c(self.postcode)+', '+c(self.city))
 173         addrlist.append(u''+c(self.street)+', '+c(self.city))
 174
 175         for addr in addrlist:
 176             params['q'] = addr.encode('utf-8')
 177             url = base_url % urlencode(params)
 178             time.sleep(1)
 179             resp, content = h.request(url)
 180             root = etree.fromstring(content)
 181             places = int(root.xpath('count(//place[@place_id])'))
 182             if places == 1:
 183                 place = root.find("place")
 184                 self.postcode = s(self.postcode,place.find("boundary").text)
 185                 self.city = s(self.city,place.find("city").text)
 186                 self.country = s(self.country,place.find("country").text)
 187                 self.country_code = s(self.country_code,place.find("country_code").text)
 188                 self.coord=Coord(place.get("lat"),place.get("lon"))
 189                 return
 190
 191     def add_phone(self, phone):
 192         self.phone = phone
 193
 194     def set_name(self, name):
 195         self.name = name
 196
 197     def db_xml(self):
 198         db = ElementMaker(namespace=const.DB_NS, nsmap=const.NSMAP)
 199         co = ''
 200         if self.coord:
 201             co = self.coord.db_xml()
 202         adr = db.address(db.street(c(self.street)),
 203                          db.postcode(c(self.postcode)),
 204                          db.city(c(self.city)),
 205                          db.country(c(self.country)),
 206                          db.phone(c(self.phone)),
 207                          co)
 208 #                       type=self.type,
 209         return adr
 210
 211
 212 def distance(p1, p2):
 213     res = Geodesic.WGS84.Inverse(p1[0], p1[1], p2[0], p2[1])
 214     return res['s12']
 215
 216 def mark(image, coord):
 217     draw = ImageDraw.Draw(image)
 218     x, y = coord
 219     r = 5
 220     bbox = (x-r, y-r, x+r, y+r)
 221     draw.ellipse(bbox, outline="red")
 222
 223 def box(image,box):
 224     draw = ImageDraw.Draw(image)
 225     draw.rectangle(box, outline="red")
 226
 227 def mapimages(coords, zoom=15,size=(TS,TS)):
 228
 229         minlat = min(coords,key=attrgetter('latitude'))
 230         maxlat = max(coords,key=attrgetter('latitude'))
 231         minlon = min(coords,key=attrgetter('longitude'))
 232         maxlon = max(coords,key=attrgetter('longitude'))
 233
 234         # Find minimal bounding box and expand it 5%
 235         hyp = distance((maxlat,minlon),(minlat,maxlon))
 236         hyp = hyp*0.05
 237         uld =  Geodesic.WGS84.Direct(maxlat, minlon, 315, hyp)
 238         urd =  Geodesic.WGS84.Direct(maxlat, maxlon,  45, hyp)
 239         lld =  Geodesic.WGS84.Direct(minlat, minlon, 225, hyp)
 240         lrd =  Geodesic.WGS84.Direct(minlat, maxlon, 135, hyp)
 241
 242         ul = Coord(maxlat,minlon).direct(315, hyp)
 243
 244         ul = (uld['lat2'],uld['lon2'])
 245         ur = (urd['lat2'],urd['lon2'])
 246         ll = (lld['lat2'],lld['lon2'])
 247         lr = (lrd['lat2'],lrd['lon2'])
 248         top = distance(ul,ur)
 249         bottom = distance(ll,lr)
 250         left = distance(ul,ll)
 251         right = distance(ur,lr)
 252 #    m_per_pix = EC*abs(cos(lat))/2**(zoomlevel+8)
 253 #    m_per_pix = ER*2*pi*abs(cos(lat))/2**(zoomlevel+8)
 254         for zoom in range(1,18):
 255             t = 2**(zoom)*TS/(ER*2*pi*abs(cos(ul[0])))
 256             toppix = t*top
 257             leftpix = t*left
 258             b = 2**(zoom)*TS/(ER*2*pi*abs(cos(ll[0])))
 259             bottompix = b*bottom
 260             rightpix  = b*right
 261 #            print "Zoom: %s : %s %s %s %s" % (zoom, toppix, bottompix, leftpix, rightpix)
 262             if max(toppix,bottompix,leftpix,rightpix) > TS*4:
 263                 break
 264 #        print "Zoom to use : %s" % (zoom)
 265         ult = coordtile(ul,zoom)
 266         lrt = coordtile(lr,zoom)
 267         lattiles = range(int(ult[0]),int(lrt[0])+1)
 268         lontiles =  range(int(ult[1]),int(lrt[1])+1)
 269         size = (len(lontiles)*TS,len(lattiles)*TS)
 270         grid = Image.new("RGB", size, None)
 271         img = []
 272         for yi, y in enumerate(lattiles):
 273             for xi, x in enumerate(lontiles):
 274                 url = tile(x,y,zoom)
 275 #                print url
 276                 time.sleep(1)
 277                 request, content = h.request(url)
 278                 img = Image.open(StringIO(content))
 279                 box = (xi*TS, yi*TS)
 280                 grid.paste(img, box)
 281         for c in coords:
 282             t = coordtile(c,zoom)
 283             o = offset(c,zoom)
 284             yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
 285             xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
 286             mark(grid, (xp,yp))
 287
 288         t = coordtile(ul,zoom)
 289         o = offset(ul,zoom)
 290         yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
 291         xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
 292         t = coordtile(lr,zoom)
 293         o = offset(lr,zoom)
 294         ypl = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
 295         xpl = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
 296         gridc = grid.crop((xp,yp,xpl,ypl))
 297         gridc.show()
 298 #        gridc.save("cap-un.png")
 299
 300 if __name__ == "__main__":
 301     for arg in sys.argv[1:]:
 302         al = arg.split("=")
 303         if al[0] == "lang":
 304             lang = al[1]
 305         if al[0] == "xptr":
 306             argument = al[1].decode('utf-8')
 307
 308     addrlist = argument.split(',')
 309     addrfmt = "street,postcode,city,country"
 310     adict = addrfmt.split(',')
 311     argdict = dict(zip(adict,addrlist))
 312     addr = Address(**argdict)
 313     addr.geocode()
 314     axml = addr.db_xml()
 315 #    clean_db(axml)
 316
 317     #print(etree.tostring(cxml, pretty_print=True))
 318     #sys.stdout.write(out.encode('utf-8'))
 319     sys.stdout.write(etree.tostring(axml,encoding='UTF-8',pretty_print=False))