xinclude/address.py

   1 #!/usr/bin/python
   2 # -*- coding: utf-8 -*-
   3
   4 import time
   5 from os import path
   6 from httplib2 import Http
   7 import urllib
   8 from math import *
   9 from lxml import etree
  10 from PIL import Image, ImageDraw
  11 from StringIO import StringIO
  12 from geohash import encode
  13 from geographiclib.geodesic import Geodesic
  14
  15 # EC Equator lenght
  16 EC = 40075016.686 # meter
  17 # ER Earth radius
  18 ER = 6372798.2 # meter
  19 # Availible zoom levels in OSM
  20 ZOOMRANGE = range(1, 18)
  21 # tile size
  22 TS = 256
  23
  24 DB_NS="http://docbook.org/ns/docbook"
  25 DB = "{%s}" % DB_NS
  26 XI_NS="http://www.w3.org/2001/XInclude"
  27 XI = "{%s}" % XI_NS
  28 XLINK_NS="http://www.w3.org/1999/xlink"
  29 XLINK = "{%s}" % XLINK_NS
  30 HTML_NS="http://www.w3.org/1999/xhtml"
  31 HTML = "{%s}" % HTML_NS
  32 NSMAP = {None : DB_NS,
  33          'xlink' : XLINK_NS}
  34
  35
  36 h = Http(".cache")
  37
  38 class Coord(object):
  39     def __init__(self, lat, lon):
  40         self.latitude  = float(lat)
  41         self.longitude = float(lon)
  42         self.image = None
  43
  44     def osmlink(self):
  45         return "http://www.openstreetmap.org/?mlat=%s&mlon=%s&zoom=18&layers=M"\
  46             % (self.latitude,self.longitude)
  47
  48     def dms(self):
  49         ns = self.latitude
  50         ew = self.longitude
  51         mnt,sec = divmod(ns*3600,60)
  52         deg,mnt = divmod(mnt,60)
  53         out = u'''%d°%2d'%5.2f"%s''' % ( deg,mnt,sec,'N')
  54         mnt,sec = divmod(ew*3600,60)
  55         deg,mnt = divmod(mnt,60)
  56         out +=  u''' %d°%2d'%05.2f"%s''' % ( deg,mnt,sec,'E')
  57         return out
  58
  59     def lontile(lon, zoom):
  60         tile = ((lon + 180) / 360) * (2**zoom)
  61         return tile
  62
  63     def lattile(lat, zoom):
  64         tile = (1 - log(tan(lat * pi/180) + 1 / cos(lat* pi/180)) / pi) /2 * 2**zoom
  65     #    tile = (1-log(tan(radians(lat))+1/cos(radians(lat)))/pi)/2*2**zoom
  66         return tile
  67
  68     def coordtile(coord, zoom):
  69         x = lontile(coord[1],zoom)
  70         y = lattile(coord[0],zoom)
  71         return (y,x)
  72
  73     def tile(x,y,z):
  74         return 'http://tile.openstreetmap.org/%s/%s/%s.png' % (z,x,y)
  75
  76     def link(coord,zoom):
  77         (x, y) = coordtile(coord,zoom)
  78         x = int(floor(x))
  79         y = int(floor(y))
  80         return tile(x,y,zoom)
  81
  82     def offset(coord, zoom):
  83         (x, y) = coordtile(coord,zoom)
  84         xo = int(floor((x-floor(x))*TS))
  85         yo = int(floor((y-floor(y))*TS))
  86         return (xo, yo)
  87
  88     def distance(self, point):
  89         res = Geodesic.WGS84.Inverse(self.latitude, self.longitude,
  90                                      point.latitude, point.longitude)
  91         return res['s12']
  92
  93     def png(self,zoom=15,size=(TS,TS)):
  94         filename = encode(self.latitude, self.longitude)+'.png'
  95         if path.isfile(filename):
  96             if path.getctime(filename) > time.time() - 60*60*24*2:
  97                 return
  98         im = Image.new("RGB", size, None)
  99         center = (size[0]/2, size[1]/2)
 100         nst = ceil(center[0]-TS/2)/TS
 101         ewt = ceil(center[1]-TS/2)/TS
 102         x, y = offset(co,zoom)
 103         (ns, ew) = coordtile(co,zoom)
 104         if x < TS/2:
 105             e = 1
 106             w = 0
 107         if x > TS/2:
 108             e = 0
 109             w = 1
 110         if y < TS/2:
 111             n = 1
 112             s = 0
 113         if y > TS/2:
 114             n = 0
 115             s = 1
 116         ul = (ns-nst-n, ew-ewt-e)
 117         lr = (ns+nst+s, ew+ewt+w)
 118         lattiles = range(int(ul[0]),int(lr[0])+1)
 119         lontiles =  range(int(ul[1]),int(lr[1])+1)
 120 #        url = link(c,zoom)
 121         size = (len(lontiles)*TS,len(lattiles)*TS)
 122         grid = Image.new("RGB", size, None)
 123         for yi, y in enumerate(lattiles):
 124             for xi, x in enumerate(lontiles):
 125                 url = tile(x,y,zoom)
 126                 time.sleep(1)
 127                 request, content = h.request(url)
 128                 img = Image.open(StringIO(content))
 129                 box = (xi*TS, yi*TS)
 130                 grid.paste(img, box)
 131         t = coordtile(co,zoom)
 132         o = offset(co,zoom)
 133         yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
 134         xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
 135         mark(grid, (xp,yp))
 136         gridc = grid.crop((xp-TS/2,yp-TS/2,xp+TS/2,yp+TS/2))
 137         gridc.save(filename)
 138
 139     def db_xml(self):
 140         uri = etree.Element(DB+'uri',nsmap=NSMAP)
 141         ln  = etree.SubElement(uri,  DB+'link')
 142         ln.set(XLINK+'href',self.osmlink())
 143         imo  = etree.SubElement(ln,  DB+'inlinemediaobject')
 144         io   = etree.SubElement(imo, DB+'imageobject', condition="web")
 145         idat = etree.SubElement(io , DB+'imagedata',
 146                                 fileref=encode(self.latitude, self.longitude)+'.png',
 147                                 format='PNG')
 148         to   = etree.SubElement(imo, DB+'textobject')
 149         ph   = etree.SubElement(to,  DB+'phrase')
 150         ph.text = "geo:"+str(self.latitude)+","+str(self.longitude)
 151         para = etree.SubElement(ln,  DB+'para')
 152         para.text = self.dms()
 153         return uri
 154
 155 class Address(object):
 156     """Address object to contain everything known about an address"""
 157     def __init__(self,address):
 158         self._address_string = address
 159         self._root = etree.Element(DB+'address',nsmap=NSMAP)
 160         self._coord = None
 161
 162     def geocode(self,country=None):
 163         params = { 'q': self._address_string,
 164                    'addressdetails': 1,
 165                    'limit': 1,
 166                    'format': 'xml',
 167                    'polygon': 0 }
 168
 169         if country:
 170             params['countrycodes'] = country
 171
 172         base_url = 'http://nominatim.openstreetmap.org/search?%s'
 173         url = base_url % urllib.urlencode(params)
 174         resp, content = h.request(url)
 175         root = etree.fromstring(content)
 176         place = root.find("place")
 177         if place is not None:
 178             print (etree.tostring(root, pretty_print=True))
 179             self._coord=Coord(place.get("lat"),place.get("lon"))
 180             return 1
 181         else:
 182             print resp
 183             print content
 184             return 0
 185
 186     def db_xml(self):
 187         return self._coord.db_xml()
 188
 189
 190 def distance(p1, p2):
 191     res = Geodesic.WGS84.Inverse(p1[0], p1[1], p2[0], p2[1])
 192     return res['s12']
 193
 194 def mark(image, coord):
 195     draw = ImageDraw.Draw(image)
 196     x, y = coord
 197     r = 5
 198     bbox = (x-r, y-r, x+r, y+r)
 199     draw.ellipse(bbox, outline="red")
 200
 201 def mapimages(coords, zoom=15,size=(TS,TS)):
 202         minlat = 1000
 203         maxlat = 0
 204         minlon = 1000
 205         maxlon = 0
 206         for c in coords:
 207             minlat = min(minlat,c[0])
 208             maxlat = max(maxlat,c[0])
 209             minlon = min(minlon,c[1])
 210             maxlon = max(maxlon,c[1])
 211         # Find minimal bounding box and expand it 5%
 212         hyp = distance((maxlat,minlon),(minlat,maxlon))
 213         hyp = hyp*0.05
 214         uld =  Geodesic.WGS84.Direct(maxlat, minlon, 315, hyp)
 215         urd =  Geodesic.WGS84.Direct(maxlat, maxlon,  45, hyp)
 216         lld =  Geodesic.WGS84.Direct(minlat, minlon, 225, hyp)
 217         lrd =  Geodesic.WGS84.Direct(minlat, maxlon, 135, hyp)
 218
 219         ul = (uld['lat2'],uld['lon2'])
 220         ur = (urd['lat2'],urd['lon2'])
 221         ll = (lld['lat2'],lld['lon2'])
 222         lr = (lrd['lat2'],lrd['lon2'])
 223         top = distance(ul,ur)
 224         bottom = distance(ll,lr)
 225         left = distance(ul,ll)
 226         right = distance(ur,lr)
 227 #    m_per_pix = EC*abs(cos(lat))/2**(zoomlevel+8)
 228 #    m_per_pix = ER*2*pi*abs(cos(lat))/2**(zoomlevel+8)
 229         for zoom in range(1,18):
 230             t = 2**(zoom)*TS/(ER*2*pi*abs(cos(ul[0])))
 231             toppix = t*top
 232             leftpix = t*left
 233             b = 2**(zoom)*TS/(ER*2*pi*abs(cos(ll[0])))
 234             bottompix = b*bottom
 235             rightpix  = b*right
 236 #            print "Zoom: %s : %s %s %s %s" % (zoom, toppix, bottompix, leftpix, rightpix)
 237             if max(toppix,bottompix,leftpix,rightpix) > TS*4:
 238                 break
 239 #        print "Zoom to use : %s" % (zoom)
 240         ult = coordtile(ul,zoom)
 241         lrt = coordtile(lr,zoom)
 242         lattiles = range(int(ult[0]),int(lrt[0])+1)
 243         lontiles =  range(int(ult[1]),int(lrt[1])+1)
 244         size = (len(lontiles)*TS,len(lattiles)*TS)
 245         grid = Image.new("RGB", size, None)
 246         img = []
 247         for yi, y in enumerate(lattiles):
 248             for xi, x in enumerate(lontiles):
 249                 url = tile(x,y,zoom)
 250 #                print url
 251                 time.sleep(1)
 252                 request, content = h.request(url)
 253                 img = Image.open(StringIO(content))
 254                 box = (xi*TS, yi*TS)
 255                 grid.paste(img, box)
 256         for c in coords:
 257             t = coordtile(c,zoom)
 258             o = offset(c,zoom)
 259             yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
 260             xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
 261             mark(grid, (xp,yp))
 262
 263         t = coordtile(ul,zoom)
 264         o = offset(ul,zoom)
 265         yp = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
 266         xp = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
 267         t = coordtile(lr,zoom)
 268         o = offset(lr,zoom)
 269         ypl = [i for i,x in enumerate(lattiles) if x == int(t[0])][0]*TS+o[1]
 270         xpl = [i for i,x in enumerate(lontiles) if x == int(t[1])][0]*TS+o[0]
 271         gridc = grid.crop((xp,yp,xpl,ypl))
 272         gridc.show()
 273 #        gridc.save("cap-un.png")