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' ]
14 import re
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:
27         """
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 < self.min:
37         #    self.min = (point, self.min)
38         #elif point > self.max:
39         #    self.max = (point, self.max)
40         #if point < self.min:
41         #    self.min = (self.min, point)
42         #elif point > self.max:
43         #    self.max = (self.max, point)
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()
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
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 < self.min \
71         # or point > self.max \
72         # or point < self.min \
73         # or point > self.max :
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 - x1
81             ya = point - y1
82             xb = x2 - x1
83             yb = y2 - y1
84             if xa * yb < xb * ya:
85                 return False
86         return True
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     """
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.replace('\n', ' ').replace('\r', ' ').replace('\t', ' ')
101     coordinates = [ xyz for xyz in coordinates.split(' ') if xyz ]
102     assert coordinates == coordinates[-1], \
103         'First and last coordinates of %s should be the same: %s, %s' % \
104         (filename, coordinates, coordinates[-1])
105     assert len(coordinates)>3, 'polygon should have 3 edges minimum'
107     area = Area()
108     for xyz in coordinates[:-1]:
109         x, y, z = xyz.split(',')
111     if reverse:
112         area.reverse()
113     assert area.check(), 'Polygon should be counter-clockwise and convex.'
114     return area
116 if __name__ == '__main__':
117     # counter clock-wise : Positive
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"