Initial code to build static images from addresses.

This code was developed over time, for the contact.py images.
Parts of it might be replaced by pythonstructures and it migth move
to a separate package if useful.

diff --git a/xinclude/addr.py b/xinclude/addr.py
new file mode 100755 (executable)
index 0000000..4ae59df
--- /dev/null
+++ b/xinclude/addr.py
@@ -0,0 +1,236 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+import time
+from os import path
+from httplib2 import Http
+import urllib
+from math import *
+from lxml import etree
+from PIL import Image, ImageDraw
+from StringIO import StringIO
+from geohash import encode
+from geographiclib.geodesic import Geodesic
+
+# EC Equator lenght
+EC = 40075016.686 # meter
+# ER Earth radius
+ER = 6372798.2 # meter
+# Availible zoom levels in OSM
+ZOOMRANGE = range(1, 18)
+# tile size
+TS = 256
+
+
+h = Http(".cache")
+
+class StaticOpenStreetMap(object):
+    """Setting up a static image from OSM with one or more markers"""
+    def __init__(self):
+        self._zoom = 0
+        self._width = 0
+        self._height = 0
+        self._markers = []
+        self._maptype = 'mapnik'
+
+    def add_marker(lon,lat):
+        self._markers.append((lon,lat))
+
+#    def construct_map(self):
+#        for coord in self._markers:
+#            print coord
+
+def lontile(lon, zoom):
+    tile = ((lon + 180) / 360) * (2**zoom)
+    return tile
+
+def lattile(lat, zoom):
+    tile = (1 - log(tan(lat * pi/180) + 1 / cos(lat* pi/180)) / pi) /2 * 2**zoom
+#    tile = (1-log(tan(radians(lat))+1/cos(radians(lat)))/pi)/2*2**zoom
+    return tile
+
+def coordtile(coord, zoom):
+    x = lontile(coord[1],zoom)
+    y = lattile(coord[0],zoom)
+    return (y,x)
+
+def tile(x,y,z):
+    return 'http://tile.openstreetmap.org/%s/%s/%s.png' % (z,x,y)
+
+def link(coord,zoom):
+    z = zoom
+    x = int(floor(lontile(coord[1],z)))
+    y = int(floor(lattile(coord[0],z)))
+    return tile(x,y,z)
+
+def offset(coord, zoom):
+    z = zoom
+    x = lontile(coord[1],z)
+    y = lattile(coord[0],z)
+    xo = int(floor((x-floor(x))*TS))
+    yo = int(floor((y-floor(y))*TS))
+    return (xo, yo)
+
+def distance(p1, p2):
+    res = Geodesic.WGS84.Inverse(p1[0], p1[1], p2[0], p2[1])
+    return res['s12']
+
+def geocode(address,country=None):
+    params = { 'q': address,
+               'addressdetails': 1,
+               'limit': 1,
+               'format': 'xml',
+               'polygon': 0 }
+    time.sleep(1)
+    if country:
+        params['countrycodes'] = country
+
+    base_url = 'http://nominatim.openstreetmap.org/search?%s'
+    url = base_url % urllib.urlencode(params)
+    resp, content = h.request(url)
+#    print resp
+#    print content
+    root = etree.fromstring(content)
+    etree.tostring(root)
+#    print "%s" % (address)
+    lat = None
+    lon = None
+    for element in root.iter("place"):
+        lat = element.get("lat")
+        lon = element.get("lon")
+#        print(" : %s,%s" % (lat, lon))
+        break
+    if not lat or not lon:
+        print resp
+        print content
+    return (float(lat),float(lon))
+
+def mark(image, coord):
+    draw = ImageDraw.Draw(image)
+    x, y = coord
+    r = 5
+    bbox = (x-r, y-r, x+r, y+r)
+    draw.ellipse(bbox, outline="red")
+
+def mapimage(coords, zoom=15,size=(TS,TS)):
+    if len(coords) == 1:
+        co = coords[0]
+        filename = encode(co[0],co[1])+'.png'
+        if path.isfile(filename):
+            if path.getctime(filename) > time.time() - 60*60*24*2:
+                return
+        im = Image.new("RGB", size, None)
+        center = (size[0]/2, size[1]/2)
+        nst = ceil(center[0]-TS/2)/TS
+        ewt = ceil(center[1]-TS/2)/TS
+        x, y = offset(co,zoom)
+        (ns, ew) = coordtile(co,zoom)
+        if x < TS/2:
+            e = 1
+            w = 0
+        if x > TS/2:
+            e = 0
+            w = 1
+        if y < TS/2:
+            n = 1
+            s = 0
+        if y > TS/2:
+            n = 0
+            s = 1
+        ul = (ns-nst-n, ew-ewt-e)
+        lr = (ns+nst+s, ew+ewt+w)
+        lattiles = range(int(ul[0]),int(lr[0])+1)
+        lontiles =  range(int(ul[1]),int(lr[1])+1)
+#        url = link(c,zoom)
+        size = (len(lontiles)*TS,len(lattiles)*TS)
+        grid = Image.new("RGB", size, None)
+        img = []
+        for yi, y in enumerate(lattiles):
+            for xi, x in enumerate(lontiles):
+                url = tile(x,y,zoom)
+#                print url
+                time.sleep(1)
+                request, content = h.request(url)
+                img = Image.open(StringIO(content))
+                box = (xi*TS, yi*TS)
+                grid.paste(img, box)
+        t = coordtile(co,zoom)
+        o = offset(co,zoom)
+        yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
+        xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
+        mark(grid, (xp,yp))
+        gridc = grid.crop((xp-TS/2,yp-TS/2,xp+TS/2,yp+TS/2))
+        gridc.save(filename)
+    else:
+        minlat = 1000
+        maxlat = 0
+        minlon = 1000
+        maxlon = 0
+        for c in coords:
+            minlat = min(minlat,c[0])
+            maxlat = max(maxlat,c[0])
+            minlon = min(minlon,c[1])
+            maxlon = max(maxlon,c[1])
+        # Find minimal bounding box and expand it 5%
+        hyp = distance((maxlat,minlon),(minlat,maxlon))
+        hyp = hyp*0.05
+        uld =  Geodesic.WGS84.Direct(maxlat, minlon, 315, hyp)
+        urd =  Geodesic.WGS84.Direct(maxlat, maxlon,  45, hyp)
+        lld =  Geodesic.WGS84.Direct(minlat, minlon, 225, hyp)
+        lrd =  Geodesic.WGS84.Direct(minlat, maxlon, 135, hyp)
+
+        ul = (uld['lat2'],uld['lon2'])
+        ur = (urd['lat2'],urd['lon2'])
+        ll = (lld['lat2'],lld['lon2'])
+        lr = (lrd['lat2'],lrd['lon2'])
+        top = distance(ul,ur)
+        bottom = distance(ll,lr)
+        left = distance(ul,ll)
+        right = distance(ur,lr)
+#    m_per_pix = EC*abs(cos(lat))/2**(zoomlevel+8)
+#    m_per_pix = ER*2*pi*abs(cos(lat))/2**(zoomlevel+8)
+        for zoom in range(1,18):
+            t = 2**(zoom)*TS/(ER*2*pi*abs(cos(ul[0])))
+            toppix = t*top
+            leftpix = t*left
+            b = 2**(zoom)*TS/(ER*2*pi*abs(cos(ll[0])))
+            bottompix = b*bottom
+            rightpix  = b*right
+#            print "Zoom: %s : %s %s %s %s" % (zoom, toppix, bottompix, leftpix, rightpix)
+            if max(toppix,bottompix,leftpix,rightpix) > TS*4:
+                break
+#        print "Zoom to use : %s" % (zoom)
+        ult = coordtile(ul,zoom)
+        lrt = coordtile(lr,zoom)
+        lattiles = range(int(ult[0]),int(lrt[0])+1)
+        lontiles =  range(int(ult[1]),int(lrt[1])+1)
+        size = (len(lontiles)*TS,len(lattiles)*TS)
+        grid = Image.new("RGB", size, None)
+        img = []
+        for yi, y in enumerate(lattiles):
+            for xi, x in enumerate(lontiles):
+                url = tile(x,y,zoom)
+#                print url
+                time.sleep(1)
+                request, content = h.request(url)
+                img = Image.open(StringIO(content))
+                box = (xi*TS, yi*TS)
+                grid.paste(img, box)
+        for c in coords:
+            t = coordtile(c,zoom)
+            o = offset(c,zoom)
+            yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
+            xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
+            mark(grid, (xp,yp))
+
+        t = coordtile(ul,zoom)
+        o = offset(ul,zoom)
+        yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
+        xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
+        t = coordtile(lr,zoom)
+        o = offset(lr,zoom)
+        ypl = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
+        xpl = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
+        gridc = grid.crop((xp,yp,xpl,ypl))
+        gridc.show()
+#        gridc.save("cap-un.png")