bin/area.py

   1 #!/usr/bin/env python
   2 # -*- encoding: utf-8
   3 """
   4 Library for areas.
   5 Basic usage is checking whether a point is inside an area.
   6 Supports GoogleEarth KML coordinates.
   7 FIXME: It should works using polar coordinated, but now works in 2D.
   8 """
   9
  10 from __future__ import division
  11
  12 __all__ = [ 'Area', 'load_area_from_kml_polygon', 'list_areas' ]
  13
  14 import os
  15 import re
  16
  17 AREA_DIR = u'/var/lib/ais/areas/'
  18
  19 class Area:
  20     """
  21     That class defines an area (on the Earth)
  22     It provides testing whether a point is inside or not
  23     """
  24     def __init__(self, points=()):
  25         self.points = []
  26         for point in points:
  27             self.addpoint(point)
  28
  29     def addpoint(self, point):
  30         """
  31         Add a (x,y,z) point.
  32         """
  33         self.points.append(point)
  34         # min/max doesn't work around the change date line...
  35         #if len(self.points)==1:
  36         #    self.min = point
  37         #    self.max = point
  38         #    return
  39         #if point[0] < self.min[0]:
  40         #    self.min = (point[0], self.min[1])
  41         #elif point[0] > self.max[0]:
  42         #    self.max = (point[0], self.max[1])
  43         #if point[1] < self.min[1]:
  44         #    self.min = (self.min[0], point[1])
  45         #elif point[1] > self.max[1]:
  46         #    self.max = (self.max[0], point[1])
  47
  48     def reverse(self):
  49         '''
  50         Invert the area
  51         From clockwise to counter-clockwise
  52         Or from counter-clockwise to clockwise
  53         '''
  54         self.points.reverse()
  55
  56     def check(self):
  57         """
  58         Area library self-test:
  59         We only support counter-clockwise and convex areas.
  60         """
  61         for point in self.points:
  62             if not self.contains(point):
  63                 return False
  64         return True
  65
  66     def contains(self, point):
  67         """
  68         Tests whether a point is in self.
  69         """
  70         if not self.points:
  71             return False
  72         # first test the bounding box
  73         #if point[0] < self.min[0] \
  74         # or point[0] > self.max[0] \
  75         # or point[1] < self.min[1] \
  76         # or point[1] > self.max[1] :
  77         #    return False
  78         for i in range(len(self.points)):
  79             p1 = self.points[i]
  80             x1, y1 = p1
  81             p2 = self.points[(i+1)%len(self.points)]
  82             x2, y2 = p2
  83             xa = point[0] - x1
  84             ya = point[1] - y1
  85             xb = x2 - x1
  86             yb = y2 - y1
  87             if xa * yb < xb * ya:
  88                 return False
  89         return True
  90
  91 def load_area_from_kml_polygon(filename, reverse=False):
  92     """
  93     Loads a kml file into an Area structure.
  94     The kml must contains exacly one set of <coordinates>.
  95     The first and last points must be the same.
  96     It must also be counter-clockwise and convex.
  97     Actually, it may be clockwise, but then you need reverse=True.
  98     """
  99     kml = open(filename).read()
 100     coordinates_lines = re.findall('<coordinates>(.*)</coordinates>', kml, re.IGNORECASE|re.DOTALL)
 101     assert len(coordinates_lines) == 1, \
 102         'There should be exactly one set of <coordinates> %s' % filename
 103     coordinates = coordinates_lines[0].replace('\n', ' ').replace('\r', ' ').replace('\t', ' ')
 104     coordinates = [ xyz for xyz in coordinates.split(' ') if xyz ]
 105     assert coordinates[0] == coordinates[-1], \
 106         'First and last coordinates of %s should be the same: %s, %s' % \
 107         (filename, coordinates[0], coordinates[-1])
 108     assert len(coordinates)>3, 'polygon should have 3 edges minimum'
 109
 110     area = Area()
 111     for xyz in coordinates[:-1]:
 112         x, y, z = xyz.split(',')
 113         area.addpoint((float(y), float(x)))
 114     if reverse:
 115         area.reverse()
 116     assert area.check(), 'Polygon should be counter-clockwise and convex.'
 117     return area
 118
 119
 120 def list_areas():
 121     """
 122     return a list of areas as tupples:
 123     [(nicename, fullpath), ...]
 124     """
 125     results = []
 126     for filename in os.listdir(AREA_DIR):
 127         if not filename.endswith(u'.kml'):
 128             continue # ignore non-kml files
 129         results.append((filename[:-4], AREA_DIR+filename))
 130     return results
 131
 132 if __name__ == '__main__':
 133     print list_areas()
 134     pelagos = load_area_from_kml_polygon('/var/lib/ais/areas/pelagos.kml')
 135     # counter clock-wise : Positive
 136     #pelagos = Area([
 137     #    (42.91, 12.5),
 138     #    (45.3612930132714, 10.01843703552244),
 139     #    (43.6,5.5),
 140     #    (40.57,8.6)
 141     #    ])
 142     for p in [
 143         (42,9),
 144         (41,5),
 145         (40,12),
 146         (45,13),
 147         (45,7),
 148         ]:
 149         print "testing", p
 150         if pelagos.contains(p):
 151             print "inside"
 152         else:
 153             print"outside"
 154