bin/area.py

   1 #!/usr/bin/env python
   2 # -*- encoding: utf-8
   3
   4 __all__ = [ 'Area', 'load_area_from_kml_polygon' ]
   5
   6 import sys
   7
   8 class Area:
   9     """
  10     That class defines an area (on the Earth)
  11     It provides testing whether a point is inside or not
  12     """
  13     def __init__(self, points=[]):
  14         self.points = []
  15         for p in points:
  16             self.addpoint(p)
  17     def addpoint(self, point):
  18         self.points.append(point)
  19         if len(self.points)==1:
  20             self.min = point
  21             self.max = point
  22             return
  23         if point[0] < self.min[0]:
  24             self.min = (point[0], self.min[1])
  25         elif point[0] > self.max[0]:
  26             self.max = (point[0], self.max[1])
  27         if point[1] < self.min[1]:
  28             self.min = (self.min[0], point[1])
  29         elif point[1] > self.max[1]:
  30             self.max = (self.max[0], point[1])
  31
  32     def check(self):
  33         for point in self.points:
  34             if not self.contains(point):
  35                 return False
  36         return True
  37
  38     def contains(self, point):
  39         if not self.points:
  40             return False
  41         # first test the bounding box
  42         #if point[0] < self.min[0] \
  43         # or point[0] > self.max[0] \
  44         # or point[1] < self.min[1] \
  45         # or point[1] > self.max[1] :
  46         #    return False
  47         for i in range(len(self.points)):
  48             p1 = self.points[i]
  49             x1, y1 = p1
  50             p2 = self.points[(i+1)%len(self.points)]
  51             x2, y2 = p2
  52             xa = point[0] - x1
  53             ya = point[1] - y1
  54             xb = x2 - x1
  55             yb = y2 - y1
  56             if xa * yb < xb * ya:
  57                 return False
  58         return True
  59
  60 def load_area_from_kml_polygon(filename):
  61     file = open(filename)
  62     coordinates_lines = [ line for line in file.readlines() if '</coordinates>' in line ]
  63     if len(coordinates_lines) != 1:
  64         print >> sys.stderr, 'There should be exactly one line with coordinates in', filename
  65         sys.exit(1)
  66     coordinates = coordinates_lines[0].replace('</coordinates>', '').replace('\n', '').replace('\r', '')
  67     coordinates = [ xyz for xyz in coordinates.split(' ') if xyz ]
  68     if coordinates[0] != coordinates[-1]:
  69         print >> sys.stderr, 'First and last coordinates of', filename, 'should be the same'
  70         print >> sys.stderr, coordinates[0]
  71         print >> sys.stderr, coordinates[-1]
  72         sys.exit(1)
  73     assert len(coordinates)>3, 'polygon should have 3 edges minimum'
  74
  75     area = Area()
  76     for xyz in coordinates[0:-1]:
  77         x,y,z = xyz.split(',')
  78         area.addpoint((float(y),float(x)))
  79     assert area.check(), 'Polygon should be counter-clockwise and convex.'
  80     return area
  81
  82
  83 #if __name__ == '__main__':
  84 # counter clock-wise : Positive
  85 #pelagos = Area([
  86 #    (42.91, 12.5),
  87 #    (45.3612930132714, 10.01843703552244),
  88 #    (43.6,5.5),
  89 #    (40.57,8.6)
  90 #    ])
  91 #for p in [
  92 #    (42,9),
  93 #    (41,5),
  94 #    (40,12),
  95 #    (45,13),
  96 #    (45,7),
  97 #    ]:
  98 #    print "testing", p
  99 #    if pelagos.contains(p):
 100 #        print "inside"
 101 #    else:
 102 #        print"outside"
 103
 104