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