Lots of cleaning and documentation
[ais.git] / 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