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:
A base class for creating sample points located in a given region of interest,
from shapely.geometry import Polygon
def __init__(self, polygon = ''):
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.
self.polygon = polygon
self.polygon = Polygon()
self.samples = list()
self.sample_count = 0
self.prepared = False
def add_polygon(self, polygon):
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
Print all sample points using their WKT representation.
for sample_pt in self.samples:
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.prepared = True
Create a stub for the actual sampling procedure.