From: Fredrik Unger Date: Mon, 24 Sep 2012 18:17:41 +0000 (+0200) Subject: Renaming addr to a more proper name, will be used to do all address tasks. X-Git-Url: https://source.tree.se/git?p=treecutter.git;a=commitdiff_plain;h=0ca9b21003f541cb6c1007f6f2e4a51abb4e96f8 Renaming addr to a more proper name, will be used to do all address tasks. --- diff --git a/xinclude/addr.py b/xinclude/addr.py deleted file mode 100755 index 4ae59df..0000000 --- a/xinclude/addr.py +++ /dev/null @@ -1,236 +0,0 @@ -#!/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") diff --git a/xinclude/address.py b/xinclude/address.py new file mode 100755 index 0000000..4ae59df --- /dev/null +++ b/xinclude/address.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")