2 # -*- coding: utf-8 -*-
7 from httplib2 import Http
8 from urllib import urlencode
9 from math import log,tan,pi,cos,radians,ceil,floor
10 from lxml import etree
11 from lxml.builder import ElementMaker
12 from PIL import Image, ImageDraw
13 from StringIO import StringIO
14 from geohash import encode
15 from geographiclib.geodesic import Geodesic
16 from treecutter import constants as const
19 EC = 40075016.686 # meter
21 ER = 6372798.2 # meter
22 # Availible zoom levels in OSM
23 ZOOMRANGE = range(1, 18)
30 def __init__(self, lat, lon):
31 self.latitude = float(lat)
32 self.longitude = float(lon)
36 return "http://www.openstreetmap.org/?mlat=%s&mlon=%s&zoom=18&layers=M"\
37 % (self.latitude,self.longitude)
42 mnt,sec = divmod(ns*3600,60)
43 deg,mnt = divmod(mnt,60)
44 out = u'''%d°%2d'%5.2f"%s''' % ( deg,mnt,sec,'N')
45 mnt,sec = divmod(ew*3600,60)
46 deg,mnt = divmod(mnt,60)
47 out += u''' %d°%2d'%05.2f"%s''' % ( deg,mnt,sec,'E')
50 def lontile(self, zoom):
51 tile = ((self.longitude + 180) / 360) * (2**zoom)
54 def lattile(self, zoom):
55 rad = radians(self.latitude)
56 tile = (1-log(tan(rad)+1/cos(rad))/pi)/2*2**zoom
59 def offset(self,zoom):
60 x = self.lontile(zoom)
61 y = self.lattile(zoom)
62 xo = int(floor((x-floor(x))*TS))
63 yo = int(floor((y-floor(y))*TS))
66 def distance(self, point):
67 res = Geodesic.WGS84.Inverse(self.latitude, self.longitude,
68 point.latitude, point.longitude)
71 def direct(self, direction, lenght):
72 point = Geodesic.WGS84.Direct(self.latitude, self.longitude,
74 return self.__class__(point['lat2'],point['lon2'])
76 def png(self,zoom=15,size=(400,150)):
77 filename = encode(self.latitude, self.longitude)+'.png'
78 # if path.isfile(filename):
79 # if path.getctime(filename) > time.time() - 60*60*24*2:
81 im = Image.new("RGB", size, None)
83 ew = int(self.lontile(zoom))
84 ns = int(self.lattile(zoom))
86 (xo, yo) = self.offset(zoom)
87 et = int(floor((xo - ceil(size[0]/2))/TS))
88 nt = int(floor((yo - ceil(size[1]/2))/TS))
89 wt = int(floor((xo + ceil(size[0]/2))/TS))
90 st = int(floor((yo + ceil(size[1]/2))/TS))
92 lontiles = range(ew+et,ew+wt+1)
93 lattiles = range(ns+nt,ns+st+1)
94 imsize = (len(lontiles)*TS,len(lattiles)*TS)
95 grid = Image.new("RGB", imsize, None)
96 for yi, y in enumerate(lattiles):
97 for xi, x in enumerate(lontiles):
98 url = 'http://tile.openstreetmap.org/%s/%s/%s.png' % (zoom,x,y)
99 request, content = h.request(url)
100 img = Image.open(StringIO(content))
101 # dr = ImageDraw.Draw(img)
102 # dr.rectangle([0,0,TS,TS], outline=0)
106 yp = [i for i,j in enumerate(lattiles) if j == int(ns)][0]*TS+yo
107 xp = [i for i,j in enumerate(lontiles) if j == int(ew)][0]*TS+xo
109 xc = int(ceil(size[0]/2))
110 yc = int(ceil(size[1]/2))
112 # draw = ImageDraw.Draw(grid)
113 # draw.rectangle([xp-xc,yp-yc,xp+xc,yp+yc], outline="red")
114 gridc = grid.crop((xp-xc,yp-yc,xp+xc,yp+yc))
119 img = encode(self.latitude, self.longitude)+'.png'
120 phr = "geo:"+str(self.latitude)+","+str(self.longitude)
122 db = ElementMaker(namespace=const.DB_NS, nsmap=const.NSMAP)
123 uri = db.uri(db.link(
124 db.inlinemediaobject(
125 db.imageobject(db.imagedata(
128 db.textobject(db.phrase(phr))
131 **{const.XLINK+"href": self.osmlink()}))
135 class Address(object):
136 """Address object to contain everything known about an address"""
137 def __init__(self,street=None,postcode=None,city=None,country=None):
140 self.postcode = postcode
142 self.country = country
147 base_url = 'http://nominatim.openstreetmap.org/search?%s'
148 params = { 'addressdetails': 1,
154 t = etree.parse('/usr/share/xml/iso-codes/iso_3166.xml')
155 r = t.xpath('//iso_3166_entry[@name="'+self.country+'"]')
157 params['countrycodes'] = r[0].get("alpha_2_code")
160 if self.name and len(self.name)>0:
161 addrlist.append(u''+self.name+', '+self.street+', '+self.city)
162 addrlist.append(u''+self.street+', '+self.postcode+', '+self.city)
163 addrlist.append(u''+self.street+', '+self.city)
165 for addr in addrlist:
166 params['q'] = addr.encode('utf-8')
167 url = base_url % urlencode(params)
169 resp, content = h.request(url)
170 root = etree.fromstring(content)
171 places = int(root.xpath('count(//place[@place_id])'))
173 place = root.find("place")
174 self.coord=Coord(place.get("lat"),place.get("lon"))
176 sys.stderr.write(u'FAILURE: Did not find:\n')
177 sys.stderr.write(addrlist[0].encode('utf-8'))
178 sys.stderr.write(url)
180 def add_phone(self, phone):
183 def set_name(self, name):
187 db = ElementMaker(namespace=const.DB_NS, nsmap=const.NSMAP)
188 adr = db.address(db.street(self.street),
189 db.postcode(self.postcode),
191 db.country(self.country),
192 db.phone(self.phone),
198 def distance(p1, p2):
199 res = Geodesic.WGS84.Inverse(p1[0], p1[1], p2[0], p2[1])
202 def mark(image, coord):
203 draw = ImageDraw.Draw(image)
206 bbox = (x-r, y-r, x+r, y+r)
207 draw.ellipse(bbox, outline="red")
209 def mapimages(coords, zoom=15,size=(TS,TS)):
215 minlat = min(minlat,c[0])
216 maxlat = max(maxlat,c[0])
217 minlon = min(minlon,c[1])
218 maxlon = max(maxlon,c[1])
219 # Find minimal bounding box and expand it 5%
220 hyp = distance((maxlat,minlon),(minlat,maxlon))
222 uld = Geodesic.WGS84.Direct(maxlat, minlon, 315, hyp)
223 urd = Geodesic.WGS84.Direct(maxlat, maxlon, 45, hyp)
224 lld = Geodesic.WGS84.Direct(minlat, minlon, 225, hyp)
225 lrd = Geodesic.WGS84.Direct(minlat, maxlon, 135, hyp)
227 ul = (uld['lat2'],uld['lon2'])
228 ur = (urd['lat2'],urd['lon2'])
229 ll = (lld['lat2'],lld['lon2'])
230 lr = (lrd['lat2'],lrd['lon2'])
231 top = distance(ul,ur)
232 bottom = distance(ll,lr)
233 left = distance(ul,ll)
234 right = distance(ur,lr)
235 # m_per_pix = EC*abs(cos(lat))/2**(zoomlevel+8)
236 # m_per_pix = ER*2*pi*abs(cos(lat))/2**(zoomlevel+8)
237 for zoom in range(1,18):
238 t = 2**(zoom)*TS/(ER*2*pi*abs(cos(ul[0])))
241 b = 2**(zoom)*TS/(ER*2*pi*abs(cos(ll[0])))
244 # print "Zoom: %s : %s %s %s %s" % (zoom, toppix, bottompix, leftpix, rightpix)
245 if max(toppix,bottompix,leftpix,rightpix) > TS*4:
247 # print "Zoom to use : %s" % (zoom)
248 ult = coordtile(ul,zoom)
249 lrt = coordtile(lr,zoom)
250 lattiles = range(int(ult[0]),int(lrt[0])+1)
251 lontiles = range(int(ult[1]),int(lrt[1])+1)
252 size = (len(lontiles)*TS,len(lattiles)*TS)
253 grid = Image.new("RGB", size, None)
255 for yi, y in enumerate(lattiles):
256 for xi, x in enumerate(lontiles):
260 request, content = h.request(url)
261 img = Image.open(StringIO(content))
265 t = coordtile(c,zoom)
267 yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
268 xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
271 t = coordtile(ul,zoom)
273 yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
274 xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
275 t = coordtile(lr,zoom)
277 ypl = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
278 xpl = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
279 gridc = grid.crop((xp,yp,xpl,ypl))
281 # gridc.save("cap-un.png")
283 if __name__ == "__main__":
284 for arg in sys.argv[1:]:
289 argument = al[1].decode('utf-8')
291 addrlist = argument.split(',')
292 addrfmt = "street,postcode,city,country"
293 adict = addrfmt.split(',')
294 argdict = dict(zip(adict,addrlist))
295 addr = Address(**argdict)
300 #print(etree.tostring(cxml, pretty_print=True))
301 #sys.stdout.write(out.encode('utf-8'))
302 sys.stdout.write(etree.tostring(axml,encoding='UTF-8',pretty_print=False))