xinclude/address.py

   1 #!/usr/bin/python
   2 # -*- coding: utf-8 -*-
   3
   4 import time
   5 from os import path
   6 from httplib2 import Http
   7 import urllib
   8 from math import *
   9 from lxml import etree
  10 from PIL import Image, ImageDraw
  11 from StringIO import StringIO
  12 from geohash import encode
  13 from geographiclib.geodesic import Geodesic
  14
  15 # EC Equator lenght
  16 EC = 40075016.686 # meter
  17 # ER Earth radius
  18 ER = 6372798.2 # meter
  19 # Availible zoom levels in OSM
  20 ZOOMRANGE = range(1, 18)
  21 # tile size
  22 TS = 256
  23
  24
  25 h = Http(".cache")
  26
  27 def lontile(lon, zoom):
  28     tile = ((lon + 180) / 360) * (2**zoom)
  29     return tile
  30
  31 def lattile(lat, zoom):
  32     tile = (1 - log(tan(lat * pi/180) + 1 / cos(lat* pi/180)) / pi) /2 * 2**zoom
  33 #    tile = (1-log(tan(radians(lat))+1/cos(radians(lat)))/pi)/2*2**zoom
  34     return tile
  35
  36 def coordtile(coord, zoom):
  37     x = lontile(coord[1],zoom)
  38     y = lattile(coord[0],zoom)
  39     return (y,x)
  40
  41 def tile(x,y,z):
  42     return 'http://tile.openstreetmap.org/%s/%s/%s.png' % (z,x,y)
  43
  44 def link(coord,zoom):
  45     z = zoom
  46     x = int(floor(lontile(coord[1],z)))
  47     y = int(floor(lattile(coord[0],z)))
  48     return tile(x,y,z)
  49
  50 def offset(coord, zoom):
  51     z = zoom
  52     x = lontile(coord[1],z)
  53     y = lattile(coord[0],z)
  54     xo = int(floor((x-floor(x))*TS))
  55     yo = int(floor((y-floor(y))*TS))
  56     return (xo, yo)
  57
  58 def distance(p1, p2):
  59     res = Geodesic.WGS84.Inverse(p1[0], p1[1], p2[0], p2[1])
  60     return res['s12']
  61
  62 def geocode(address,country=None):
  63     params = { 'q': address,
  64                'addressdetails': 1,
  65                'limit': 1,
  66                'format': 'xml',
  67                'polygon': 0 }
  68     time.sleep(1)
  69     if country:
  70         params['countrycodes'] = country
  71
  72     base_url = 'http://nominatim.openstreetmap.org/search?%s'
  73     url = base_url % urllib.urlencode(params)
  74     resp, content = h.request(url)
  75 #    print resp
  76 #    print content
  77     root = etree.fromstring(content)
  78     etree.tostring(root)
  79 #    print "%s" % (address)
  80     lat = None
  81     lon = None
  82     for element in root.iter("place"):
  83         lat = element.get("lat")
  84         lon = element.get("lon")
  85 #        print(" : %s,%s" % (lat, lon))
  86         break
  87     if not lat or not lon:
  88         print resp
  89         print content
  90     return (float(lat),float(lon))
  91
  92 def mark(image, coord):
  93     draw = ImageDraw.Draw(image)
  94     x, y = coord
  95     r = 5
  96     bbox = (x-r, y-r, x+r, y+r)
  97     draw.ellipse(bbox, outline="red")
  98
  99 def mapimage(coords, zoom=15,size=(TS,TS)):
 100     if len(coords) == 1:
 101         co = coords[0]
 102         filename = encode(co[0],co[1])+'.png'
 103         if path.isfile(filename):
 104             if path.getctime(filename) > time.time() - 60*60*24*2:
 105                 return
 106         im = Image.new("RGB", size, None)
 107         center = (size[0]/2, size[1]/2)
 108         nst = ceil(center[0]-TS/2)/TS
 109         ewt = ceil(center[1]-TS/2)/TS
 110         x, y = offset(co,zoom)
 111         (ns, ew) = coordtile(co,zoom)
 112         if x < TS/2:
 113             e = 1
 114             w = 0
 115         if x > TS/2:
 116             e = 0
 117             w = 1
 118         if y < TS/2:
 119             n = 1
 120             s = 0
 121         if y > TS/2:
 122             n = 0
 123             s = 1
 124         ul = (ns-nst-n, ew-ewt-e)
 125         lr = (ns+nst+s, ew+ewt+w)
 126         lattiles = range(int(ul[0]),int(lr[0])+1)
 127         lontiles =  range(int(ul[1]),int(lr[1])+1)
 128 #        url = link(c,zoom)
 129         size = (len(lontiles)*TS,len(lattiles)*TS)
 130         grid = Image.new("RGB", size, None)
 131         img = []
 132         for yi, y in enumerate(lattiles):
 133             for xi, x in enumerate(lontiles):
 134                 url = tile(x,y,zoom)
 135 #                print url
 136                 time.sleep(1)
 137                 request, content = h.request(url)
 138                 img = Image.open(StringIO(content))
 139                 box = (xi*TS, yi*TS)
 140                 grid.paste(img, box)
 141         t = coordtile(co,zoom)
 142         o = offset(co,zoom)
 143         yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
 144         xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
 145         mark(grid, (xp,yp))
 146         gridc = grid.crop((xp-TS/2,yp-TS/2,xp+TS/2,yp+TS/2))
 147         gridc.save(filename)
 148     else:
 149         minlat = 1000
 150         maxlat = 0
 151         minlon = 1000
 152         maxlon = 0
 153         for c in coords:
 154             minlat = min(minlat,c[0])
 155             maxlat = max(maxlat,c[0])
 156             minlon = min(minlon,c[1])
 157             maxlon = max(maxlon,c[1])
 158         # Find minimal bounding box and expand it 5%
 159         hyp = distance((maxlat,minlon),(minlat,maxlon))
 160         hyp = hyp*0.05
 161         uld =  Geodesic.WGS84.Direct(maxlat, minlon, 315, hyp)
 162         urd =  Geodesic.WGS84.Direct(maxlat, maxlon,  45, hyp)
 163         lld =  Geodesic.WGS84.Direct(minlat, minlon, 225, hyp)
 164         lrd =  Geodesic.WGS84.Direct(minlat, maxlon, 135, hyp)
 165
 166         ul = (uld['lat2'],uld['lon2'])
 167         ur = (urd['lat2'],urd['lon2'])
 168         ll = (lld['lat2'],lld['lon2'])
 169         lr = (lrd['lat2'],lrd['lon2'])
 170         top = distance(ul,ur)
 171         bottom = distance(ll,lr)
 172         left = distance(ul,ll)
 173         right = distance(ur,lr)
 174 #    m_per_pix = EC*abs(cos(lat))/2**(zoomlevel+8)
 175 #    m_per_pix = ER*2*pi*abs(cos(lat))/2**(zoomlevel+8)
 176         for zoom in range(1,18):
 177             t = 2**(zoom)*TS/(ER*2*pi*abs(cos(ul[0])))
 178             toppix = t*top
 179             leftpix = t*left
 180             b = 2**(zoom)*TS/(ER*2*pi*abs(cos(ll[0])))
 181             bottompix = b*bottom
 182             rightpix  = b*right
 183 #            print "Zoom: %s : %s %s %s %s" % (zoom, toppix, bottompix, leftpix, rightpix)
 184             if max(toppix,bottompix,leftpix,rightpix) > TS*4:
 185                 break
 186 #        print "Zoom to use : %s" % (zoom)
 187         ult = coordtile(ul,zoom)
 188         lrt = coordtile(lr,zoom)
 189         lattiles = range(int(ult[0]),int(lrt[0])+1)
 190         lontiles =  range(int(ult[1]),int(lrt[1])+1)
 191         size = (len(lontiles)*TS,len(lattiles)*TS)
 192         grid = Image.new("RGB", size, None)
 193         img = []
 194         for yi, y in enumerate(lattiles):
 195             for xi, x in enumerate(lontiles):
 196                 url = tile(x,y,zoom)
 197 #                print url
 198                 time.sleep(1)
 199                 request, content = h.request(url)
 200                 img = Image.open(StringIO(content))
 201                 box = (xi*TS, yi*TS)
 202                 grid.paste(img, box)
 203         for c in coords:
 204             t = coordtile(c,zoom)
 205             o = offset(c,zoom)
 206             yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
 207             xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
 208             mark(grid, (xp,yp))
 209
 210         t = coordtile(ul,zoom)
 211         o = offset(ul,zoom)
 212         yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
 213         xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
 214         t = coordtile(lr,zoom)
 215         o = offset(lr,zoom)
 216         ypl = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
 217         xpl = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
 218         gridc = grid.crop((xp,yp,xpl,ypl))
 219         gridc.show()
 220 #        gridc.save("cap-un.png")