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