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