I needed to generate a set of random points within a polygon, including convex and concave. The need arouse in a geospatial domain where polygons are rather small (on a geo-scale) and wouldn’t span more than say 10 miles, though the benefit of employing more complex algorithms to deal with spheroid properties are negligible. Plane geometry provided enough to meet this requirement. Point-in-Polygon tests are rather simple and are used to test whether a point exists in a polygon. The test is performed using a Ray casting algorithm which test the intersections of a ray across the x-axis starting from the point in question.
Another concept is the Minimum Bounding Rectangle (Bounding Box), which is the minimal rectangle needed to enclose a geographical object (i.e. polygon).
So, one can generate random points within a polygon by…
- Generating a bounding box
- Generating a point within the bounding box. This is a simple algorithm.
- Using Point-in-Polygon to test whether this point exists within the polygon.
Because of the random sampling nature and false positives from step 2, which must be tested in step 3, the above must be performed in a loop until the Point-in-Polygon test passes.
This works quite well for generating test data, as there are no tight bounds on the performance characteristics of random generation. One could also use the above algorithm in production as long as the ration of polygon to bounding box is rather large, which is usually the case for convex polygons. The ratio might be too small convex polygons, though causing a more than acceptable number of false positives in step #2.
To utilize the geo-utils to generate random points within a polygon, you would do the following:
from vtown import geo from vtown.geo.polygon import Polygon polygon = Polygon( geo.LatLon(42.39321,-82.92114), geo.LatLon(42.39194,-82.91669), geo.LatLon(42.39147,-82.91796), geo.LatLon(42.39090,-82.91974), geo.LatLon(42.39321,-82.92114)) point = polygon.random_point()
The above polygon is generated using lat/lon coordinates, but you can generate them using simple x/y coordinates with geo.Point(x,y)
Here are some code snippets from the implementation. I only pasted the relevant parts. For boilerplate and relevant data structures, see the geo-utils package.
class BoundingBox(object): def __init__(self, *points): """docstring for __init__""" xmin = ymin = float('inf') xmax = ymax = float('-inf') for p in points: if p.x < xmin: xmin = p.x if p.y < ymin: ymin = p.y if p.x > xmax: xmax = p.x if p.y > ymax: ymax = p.y self.interval_x = Interval(xmin, xmax) self.interval_y = Interval(ymin, ymax) def random_point(self): x = self.interval_x.random_point() y = self.interval_y.random_point() return Point(x, y) class Polygon: ## __init__ omitted here... def contains(self, point): seg_counter = private.SegmentCounter(point) for i in range(1, len(self.points)): line = Line(*self.points[i-1:i+1]) if seg_counter.process_segment(line): return True return seg_counter.crossings % 2 == 1 def random_point(self): bb = BoundingBox(*self.points) while True: print("GENERATING RANDOM POINT...") p = bb.random_point() if self.contains(p): return p class SegmentCounter(object): def __init__(self, point): self.point = point self.crossings = 0 def process_segment(self, line): p, p1, p2 = self.point, line.point1, line.point2 if p1.x < p.x and p2.x < p.x: return False if (p.x == p2.x and p.y == p2.y): return True if p1.y == p.y and p2.y == p.y: minx = p1.x maxx = p2.x if minx > maxx: minx = p2.x maxx = p1.x if p.x >= minx and p.x <= maxx: return True return False if ((p1.y > p.y) and (p2.y <= p.y)) \ or ((p2.y > p.y) and (p1.y <= p.y)): x1 = p1.x - p.x y1 = p1.y - p.y x2 = p2.x - p.x y2 = p2.y - p.y det = numpy.linalg.det([[x1, y1], [x2, y2]]) if det == 0.0: return True if y2 < y1: det = -det if det > 0.0: self.crossings += 1