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