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 __all__ = [ 'Area', 'load_area_from_kml_polygon' ]
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:
23         """
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 < self.min:
33         #    self.min = (point, self.min)
34         #elif point > self.max:
35         #    self.max = (point, self.max)
36         #if point < self.min:
37         #    self.min = (self.min, point)
38         #elif point > self.max:
39         #    self.max = (self.max, point)
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
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 < self.min \
59         # or point > self.max \
60         # or point < self.min \
61         # or point > self.max :
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 - x1
69             ya = point - y1
70             xb = x2 - x1
71             yb = y2 - y1
72             if xa * yb < xb * ya:
73                 return False
74         return True
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.
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.replace('</coordinates>', '').replace('\n', '').replace('\r', '')
91     coordinates = [ xyz for xyz in coordinates.split(' ') if xyz ]
92     assert coordinates == coordinates[-1], \
93         'First and last coordinates of %s should be the same: %s, %s' % \
94         (filename, coordinates, coordinates[-1])
95     assert len(coordinates)>3, 'polygon should have 3 edges minimum'
97     area = Area()
98     for xyz in coordinates[0:-1]:
99         x, y, z = xyz.split(',')
101     assert area.check(), 'Polygon should be counter-clockwise and convex.'
102     return area
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"