+++ /dev/null
-#!/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")
--- /dev/null
+#!/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")
projects / treecutter.git / commitdiff