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()}))
141 class Address(object):
142 """Address object to contain everything known about an address"""
143 def __init__(self,street=None,postcode=None,city=None,country=None):
145 self.postcode = postcode
147 self.country = country
148 self.country_code = None
154 base_url = 'http://nominatim.openstreetmap.org/search?%s'
155 params = { 'addressdetails': 1,
161 t = etree.parse('/usr/share/xml/iso-codes/iso_3166.xml')
162 r = t.xpath('//iso_3166_entry[@name="'+self.country+'"]')
164 self.country_code = r[0].get("alpha_2_code")
165 if self.country_code:
166 params['countrycodes'] = self.country_code
169 addrlist.append(u''+c(self.name)+', '+c(self.street)+', '+c(self.city))
170 addrlist.append(u''+c(self.street)+', '+c(self.postcode)+', '+c(self.city))
171 addrlist.append(u''+c(self.street)+', '+c(self.city))
173 for addr in addrlist:
174 params['q'] = addr.encode('utf-8')
175 url = base_url % urlencode(params)
177 resp, content = h.request(url)
178 root = etree.fromstring(content)
179 places = int(root.xpath('count(//place[@place_id])'))
181 place = root.find("place")
182 self.postcode = s(self.postcode,place.find("boundary").text)
183 self.city = s(self.city,place.find("city").text)
184 self.country = s(self.country,place.find("country").text)
185 self.country_code = s(self.country_code,place.find("country_code").text)
186 self.coord=Coord(place.get("lat"),place.get("lon"))
189 def add_phone(self, phone):
192 def set_name(self, name):
196 db = ElementMaker(namespace=const.DB_NS, nsmap=const.NSMAP)
199 co = self.coord.db_xml()
200 adr = db.address(db.street(c(self.street)),
201 db.postcode(c(self.postcode)),
202 db.city(c(self.city)),
203 db.country(c(self.country)),
204 db.phone(c(self.phone)),
210 def distance(p1, p2):
211 res = Geodesic.WGS84.Inverse(p1[0], p1[1], p2[0], p2[1])
214 def mark(image, coord):
215 draw = ImageDraw.Draw(image)
218 bbox = (x-r, y-r, x+r, y+r)
219 draw.ellipse(bbox, outline="red")
221 def mapimages(coords, zoom=15,size=(TS,TS)):
227 minlat = min(minlat,c[0])
228 maxlat = max(maxlat,c[0])
229 minlon = min(minlon,c[1])
230 maxlon = max(maxlon,c[1])
231 # Find minimal bounding box and expand it 5%
232 hyp = distance((maxlat,minlon),(minlat,maxlon))
234 uld = Geodesic.WGS84.Direct(maxlat, minlon, 315, hyp)
235 urd = Geodesic.WGS84.Direct(maxlat, maxlon, 45, hyp)
236 lld = Geodesic.WGS84.Direct(minlat, minlon, 225, hyp)
237 lrd = Geodesic.WGS84.Direct(minlat, maxlon, 135, hyp)
239 ul = (uld['lat2'],uld['lon2'])
240 ur = (urd['lat2'],urd['lon2'])
241 ll = (lld['lat2'],lld['lon2'])
242 lr = (lrd['lat2'],lrd['lon2'])
243 top = distance(ul,ur)
244 bottom = distance(ll,lr)
245 left = distance(ul,ll)
246 right = distance(ur,lr)
247 # m_per_pix = EC*abs(cos(lat))/2**(zoomlevel+8)
248 # m_per_pix = ER*2*pi*abs(cos(lat))/2**(zoomlevel+8)
249 for zoom in range(1,18):
250 t = 2**(zoom)*TS/(ER*2*pi*abs(cos(ul[0])))
253 b = 2**(zoom)*TS/(ER*2*pi*abs(cos(ll[0])))
256 # print "Zoom: %s : %s %s %s %s" % (zoom, toppix, bottompix, leftpix, rightpix)
257 if max(toppix,bottompix,leftpix,rightpix) > TS*4:
259 # print "Zoom to use : %s" % (zoom)
260 ult = coordtile(ul,zoom)
261 lrt = coordtile(lr,zoom)
262 lattiles = range(int(ult[0]),int(lrt[0])+1)
263 lontiles = range(int(ult[1]),int(lrt[1])+1)
264 size = (len(lontiles)*TS,len(lattiles)*TS)
265 grid = Image.new("RGB", size, None)
267 for yi, y in enumerate(lattiles):
268 for xi, x in enumerate(lontiles):
272 request, content = h.request(url)
273 img = Image.open(StringIO(content))
277 t = coordtile(c,zoom)
279 yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
280 xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
283 t = coordtile(ul,zoom)
285 yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
286 xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
287 t = coordtile(lr,zoom)
289 ypl = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
290 xpl = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
291 gridc = grid.crop((xp,yp,xpl,ypl))
293 # gridc.save("cap-un.png")
295 if __name__ == "__main__":
296 for arg in sys.argv[1:]:
301 argument = al[1].decode('utf-8')
303 addrlist = argument.split(',')
304 addrfmt = "street,postcode,city,country"
305 adict = addrfmt.split(',')
306 argdict = dict(zip(adict,addrlist))
307 addr = Address(**argdict)
312 #print(etree.tostring(cxml, pretty_print=True))
313 #sys.stdout.write(out.encode('utf-8'))
314 sys.stdout.write(etree.tostring(axml,encoding='UTF-8',pretty_print=False))