From 1ed4abf91e059a83c82b07be9931d5b605085677 Mon Sep 17 00:00:00 2001 From: Fredrik Unger Date: Wed, 16 May 2012 12:44:12 +0200 Subject: [PATCH] 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. --- xinclude/addr.py | 236 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 236 insertions(+) create mode 100755 xinclude/addr.py diff --git a/xinclude/addr.py b/xinclude/addr.py new file mode 100755 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") -- 2.30.2