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.
7 FIXME: It should works using polar coordinated, but now works in 2D.
8 """
10 from __future__ import division
12 __all__ = [ 'Area', 'load_area_from_kml_polygon', 'list_areas' ]
14 import os
15 import re
17 AREA_DIR = u'/var/lib/ais/areas/'
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:
30         """
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 < self.min:
40         #    self.min = (point, self.min)
41         #elif point > self.max:
42         #    self.max = (point, self.max)
43         #if point < self.min:
44         #    self.min = (self.min, point)
45         #elif point > self.max:
46         #    self.max = (self.max, point)
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()
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
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 < self.min \
74         # or point > self.max \
75         # or point < self.min \
76         # or point > self.max :
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 - x1
84             ya = point - y1
85             xb = x2 - x1
86             yb = y2 - y1
87             if xa * yb < xb * ya:
88                 return False
89         return True
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     """
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.replace('\n', ' ').replace('\r', ' ').replace('\t', ' ')
104     coordinates = [ xyz for xyz in coordinates.split(' ') if xyz ]
105     assert coordinates == coordinates[-1], \
106         'First and last coordinates of %s should be the same: %s, %s' % \
107         (filename, coordinates, coordinates[-1])
108     assert len(coordinates)>3, 'polygon should have 3 edges minimum'
110     area = Area()
111     for xyz in coordinates[:-1]:
112         x, y, z = xyz.split(',')
114     if reverse:
115         area.reverse()
116     assert area.check(), 'Polygon should be counter-clockwise and convex.'
117     return area
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     # sort by name
131     results.sort(cmp=lambda a1,a2: cmp(a1, a2))
132     return results
134 if __name__ == '__main__':
135     print list_areas()
137     # counter clock-wise : Positive
138     #pelagos = Area([
139     #    (42.91, 12.5),
140     #    (45.3612930132714, 10.01843703552244),
141     #    (43.6,5.5),
142     #    (40.57,8.6)
143     #    ])
144     for p in [
145         (42,9),
146         (41,5),
147         (40,12),
148         (45,13),
149         (45,7),
150         ]:
151         print "testing", p
152         if pelagos.contains(p):
153             print "inside"
154         else:
155             print"outside"