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 polylines.
   7 FIXME: It should works using polar coordinated, but now works in 2D.
   8 """
   9
  10 __all__ = [ 'Area', 'load_area_from_kml_polygon' ]
  11
  12 class Area:
  13     """
  14     That class defines an area (on the Earth)
  15     It provides testing whether a point is inside or not
  16     """
  17     def __init__(self, points=()):
  18         self.points = []
  19         for point in points:
  20             self.addpoint(point)
  21
  22     def addpoint(self, point):
  23         """
  24         Add a (x,y,z) point.
  25         """
  26         self.points.append(point)
  27         # min/max doesn't work around the change date line...
  28         #if len(self.points)==1:
  29         #    self.min = point
  30         #    self.max = point
  31         #    return
  32         #if point[0] < self.min[0]:
  33         #    self.min = (point[0], self.min[1])
  34         #elif point[0] > self.max[0]:
  35         #    self.max = (point[0], self.max[1])
  36         #if point[1] < self.min[1]:
  37         #    self.min = (self.min[0], point[1])
  38         #elif point[1] > self.max[1]:
  39         #    self.max = (self.max[0], point[1])
  40
  41     def check(self):
  42         """
  43         Area library self-test:
  44         We only support counter-clockwise and convex areas.
  45         """
  46         for point in self.points:
  47             if not self.contains(point):
  48                 return False
  49         return True
  50
  51     def contains(self, point):
  52         """
  53         Tests whether a point is in self.
  54         """
  55         if not self.points:
  56             return False
  57         # first test the bounding box
  58         #if point[0] < self.min[0] \
  59         # or point[0] > self.max[0] \
  60         # or point[1] < self.min[1] \
  61         # or point[1] > self.max[1] :
  62         #    return False
  63         for i in range(len(self.points)):
  64             p1 = self.points[i]
  65             x1, y1 = p1
  66             p2 = self.points[(i+1)%len(self.points)]
  67             x2, y2 = p2
  68             xa = point[0] - x1
  69             ya = point[1] - y1
  70             xb = x2 - x1
  71             yb = y2 - y1
  72             if xa * yb < xb * ya:
  73                 return False
  74         return True
  75
  76 def load_area_from_kml_polygon(filename):
  77     """
  78     Loads a kml file into an Area structure.
  79     The kml must contains exacly 1 polyline structure.
  80     The first and last points must be the same.
  81     It must also be counter-clockwise and convex.
  82
  83     FIXME: This makes a lot of assumption about the way GoogleEarth output the
  84     XML file.
  85     """
  86     kmlfile = open(filename)
  87     coordinates_lines = [ line for line in kmlfile.readlines() if '</coordinates>' in line ]
  88     assert len(coordinates_lines) == 1, \
  89         'There should be exactly one line with coordinates in %s' % filename
  90     coordinates = coordinates_lines[0].replace('</coordinates>', '').replace('\n', '').replace('\r', '')
  91     coordinates = [ xyz for xyz in coordinates.split(' ') if xyz ]
  92     assert coordinates[0] == coordinates[-1], \
  93         'First and last coordinates of %s should be the same: %s, %s' % \
  94         (filename, coordinates[0], coordinates[-1])
  95     assert len(coordinates)>3, 'polygon should have 3 edges minimum'
  96
  97     area = Area()
  98     for xyz in coordinates[0:-1]:
  99         x, y, z = xyz.split(',')
 100         area.addpoint((float(y), float(x)))
 101     assert area.check(), 'Polygon should be counter-clockwise and convex.'
 102     return area
 103
 104
 105 #if __name__ == '__main__':
 106 # counter clock-wise : Positive
 107 #pelagos = Area([
 108 #    (42.91, 12.5),
 109 #    (45.3612930132714, 10.01843703552244),
 110 #    (43.6,5.5),
 111 #    (40.57,8.6)
 112 #    ])
 113 #for p in [
 114 #    (42,9),
 115 #    (41,5),
 116 #    (40,12),
 117 #    (45,13),
 118 #    (45,7),
 119 #    ]:
 120 #    print "testing", p
 121 #    if pelagos.contains(p):
 122 #        print "inside"
 123 #    else:
 124 #        print"outside"
 125
 126