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