Creating sample points for a given region of interest is a common task in geospatial analysis. It is therefore logically consistent that there is already a number of ways available to create some including the Create Random Points tool of the ArcToolbox in ArcGIS or the corresponding components of the fTools plugin for QGIS. As I’m trying to follow an explicit under-the-hood-philosophy, I will – starting with this very first “real” entry at Portolan – implement my very own sampling methodology for anyone willing to follow.

As it is my (current) programming language of choice I will use Python to accomplish my task. Adding to the fact that I have plenty of working experience with it, Python has the advantage of being very well positioned in the realm of geospatial processing. This is courtesy of a wide range of libraries dealing with corresponding tasks including two that I will use extensively, namely OGR (as part of GDAL) and Shapely. While OGR serves as a well suited toolbox for all things vector – including export and import to external file and/or database formats, basic dataset creation and editing as well as more sophisticated procedures such as generalization – I have found Shapely (basically representing a Pythonic interface to libgeos that is also used by OGR itself) to be a more direct link to the topological operations that are bread and butter for any kind of geographical information system.

As Python explicitly encourages the object-oriented programming paradigm, I will follow that and implement my very own PolygonPointSampler in compliance to that paradigm. Mind you, I’m not an explicitly stated computer scientist, but a cartographer by trait that somehow turned into a mostly self-taught specialist for geoinformatics. My theoretical ramblings with regard to programming may be a little off from time to time, however all of the things presented here do actually work in practice – which is most important for me. And maybe for the reader as well.

Corresponding to these prerequisites a base class for a PolygonPointSampler could be implemented as set out in the following listing:

u"""
A base class for creating sample points located in a given region of interest,
i.e. polygon.
"""

from shapely.geometry import Polygon

class PolygonPointSampler(object):

    def __init__(self, polygon = ''):
        u"""
        Initialize a new PolygonPointSampler object using the specified polygon
        object (as allocated by Shapely). If no polygon is given a new empty
        one is created and set as the base polygon.
        """
        if polygon:
            self.polygon = polygon
        else:
            self.polygon = Polygon()
        self.samples = list()
        self.sample_count = 0
        self.prepared = False

    def add_polygon(self, polygon):
        u"""
        Add another polygon entity to the base polygon by geometrically unifying
        it with the current one.
        """
        self.polygon = self.polygon.union(polygon)
        self.prepared = False

    def print_samples(self):
        u"""
        Print all sample points using their WKT representation.
        """
        for sample_pt in self.samples:
            print sample_pt

    def prepare_sampling(self):
        u"""
        Prepare the actual sampling procedure by splitting up the specified base
        polygon (that may consist of multiple simple polygons) and appending its
        compartments to a dedicated list.
        """
        self.src = list()
        if hasattr(self.polygon, 'geoms'):
            for py in self.polygon:
                self.src.append(py)
        else:
            self.src.append(self.polygon)
        self.prepared = True

    def perform_sampling(self):
        u"""
        Create a stub for the actual sampling procedure.
        """
        raise NotImplementedError