In this edition of our series dedicated to polygon sampling techniques we will look into the process of creating regularly gridded sample points with non-uniform intervals.

Previously published parts of this series include:

• Introduction
• Part 1: Centroids and representative points

Regular grid sampling

Using a single point to represent a whole polygon geometry may satisfy only the most basic sampling demands. Another – and certainly more applicable – way to arrange sample points is to create virtual grids of such points that stretch out over the whole area of the polygon to be processed. To do so we need a regular interval between sample points that may be defined identical (i.e. uniform) in both x- and y-direction. Here we will go the more generic route and implement separately adjustable (i.e. non-uniform) intervals for x and y.

The value range of sample coordinates is of course laid out by the extent of the source polygon. Using Shapely the extent of a polygon is called forth by the property bounds:

>>> from shapely.geometry import Polygon
>>> polygon = Polygon([(0,0), (6,0), (0,6)])
>>> print polygon.bounds
(0.0, 0.0, 6.0, 6.0)

Given two intervals in x and y it is now easy to create points at regular gridded positions laid out over the extent of the source polygon. Additionally it should be assured that created points are actually within the polygon to be sampled.

>>> from shapely.geometry import Point
>>> bounds = polygon.bounds
>>> ll = bounds[:2] # lower left coordinate pair of polygon's extent
>>> ur = bounds[2:] # upper right                  ~
>>> x_interval = 1.5
>>> y_interval = 2.0
>>> for x in floatrange(ll[0], ur[0], x_interval):
...     for y in floatrange(ll[1], ur[1], y_interval):
...             point = Point(x, y)
...             if point.within(polygon):
...                     print point
...
POINT (1.5 2)
POINT (1.5 4)
POINT (3 2)

Now real-world spatial data is different yet again as it rarely comes with extents fitting to full meters. It is still possible to have regularly gridded sample points that have *nicely* looking coordinates by creating extent ceilings and floors used as base for the sampling process. We can actually use our defined intervals to create these values:

>>> from shapely.wkt import loads
>>> polygon = loads(open(r"real_world_data.wkt").read())
>>> polygon.bounds
(366392.3774440361, 5630693.4900143575, 373404.5164361303, 5641855.842006282)
>>> ll = polygon.bounds[:2]
>>> ur = polygon.bounds[2:]
>>> x_interval = 100
>>> y_interval = 200
>>> low_x = int(ll[0]) / x_interval * x_interval
>>> upp_x = int(ur[0]) / x_interval * x_interval + x_interval
>>> low_y = int(ll[1]) / y_interval * y_interval
>>> upp_y = int(ur[1]) / y_interval * y_interval + y_interval
>>> print low_x, upp_x, low_y, upp_y
366300 373500 5630600 5642000

Putting it all together

To extend our previously introduced PolygonPointSampler, we can combine all our findings in a new sub class RegularGridSampler. This one will also make use of the possibility of creating a separate constructor as there is the need to define the sampling intervals.

class RegularGridSampler(PolygonPointSampler):
    def __init__(self, polygon = '', x_interval = 100, y_interval = 100):
        super(self.__class__, self).__init__(polygon)
        self.x_interval = x_interval
        self.y_interval = y_interval
    
    def perform_sampling(self):
        u"""
        Perform sampling by substituting the polygon with a regular grid of
        sample points within it. The distance between the sample points is
        given by x_interval and y_interval.
        """
        if not self.prepared:
            self.prepare_sampling()
        ll = self.polygon.bounds[:2]
        ur = self.polygon.bounds[2:]
        low_x = int(ll[0]) / self.x_interval * self.x_interval
        upp_x = int(ur[0]) / self.x_interval * self.x_interval + self.x_interval
        low_y = int(ll[1]) / self.y_interval * self.y_interval
        upp_y = int(ur[1]) / self.y_interval * self.y_interval + self.y_interval
        
        for x in floatrange(low_x, upp_x, self.x_interval):
            for y in floatrange(low_y, upp_y, self.y_interval):
                p = Point(x, y)
                if p.within(self.polygon):
                    self.samples.append(p)

Using our real-world example and a (uniform) sampling interval of 1,000 meters we arrive at the following result:

We may also use non-uniform sampling intervals of 1,000 meters in x- and 500 meters in y-direction:

A multi-part polygon may also be used to apply the sampler, for example with a uniform sampling interval of 250 m:

WTH is a floatrange?

In the meantime you most likely have already realized that we’re using a non-standard range function to iterate over a range of float values. Based on a StackExchange suggestion I have defined an according routine in a utility function:

def floatrange(start, stop, step):
    while start < stop:
        yield start
        start += step

Following the declaration of my PolygonPointSampler’s base class in the previous post, I will now start to implement multiple classes each representing a different method to derive sample points for a given polygon.

Centroid sampling

Let’s start with the simplest polygon-to-point sampling method available: the creation of a centroid point. The centroid represents the geometric center of a polygon. If we’re looking at a triangle it emerges as the intersection of the triangle’s median lines, hence it is sometimes also called the median point. Shapely allows to derive a polygon’s centroid by querying the corresponding attribute:

>>> from shapely.geometry import Polygon
>>> polygon = Polygon([(0,0), (6,0), (0,6)])
>>> centroid = polygon.centroid
>>> print centroid
POINT (2 2)

Applying this method to some real-word data will often lead to the phenomenon visible below: The centroid comes to lie outside of the original polygon’s interior.

Representative point sampling

To deal with this problem, Shapely equips geometry objects with a method called representative_point() that – as the documentation reads – “returns a cheaply computed point that is guaranteed to be within the geometric object”. Judging from the Shapely source code this method goes back to the PointOnSurface function provided by GEOS. I haven’t been able to find out how exactly it is guaranteed that the sample point is an interior point, but from the looks of it, it is most likely that the algorithm described in a fitting thread at gis.stackexchange.com has been applied to it. In any way for our Python example the code would look like the following:

>>> representative_point = polygon.representative_point()
>>> print representative_point
POINT (1.5 3)

Please note that querying for the representative point of a polygon actually calls a Python function – as indicated by the brackets following the according declaration – while the centroid is a property of the Shapely polygon object. Using real-world data we arrive at the situation displayed below:

It is also possible to apply the sampler on multi-part polygons:

Putting it all together

Based on the knowledge laid out above, it is now possible to furnish the previously created PolygonPointSampler base object with extensions that create sample points by either using the centroid or representative point method. I have called them CentroidSampler and RepresentativePointSampler, respectively:

class CentroidSampler(PolygonPointSampler):
    def perform_sampling(self):
        u"""
        Perform sampling by reprensenting each source polygon with its centroid.
        """
        if not self.prepared:
            self.prepare_sampling()
        for src in self.src:
            self.samples.append(src.centroid)

class RepresentativePointSampler(PolygonPointSampler):
    def perform_sampling(self):
        u"""
        Perform sampling by representing each source polygon with a
        representative point whose coordinates are guaranteed to be within
        the polygon's geometry.
        """
        if not self.prepared:
            self.prepare_sampling()
        for src in self.src:
            self.samples.append(src.representative_point())