Fastest way to produce a grid of points that fall within a polygon or shape?

Question

I am using shapely in python and trying to generate evenly spaced points in a grid that fall within a shape in the fastest O(n) time. The shape may be any closed polygon, not just a square or circle. My current approach is:

1. Find min/max y & x to build a rectangle.
2. Build a grid of points given a spacing parameter (resolution)
3. Verify one-by-one if the points fall within the shape.

Is there a faster way to do this?

# determine maximum edges
polygon = shape(geojson['features'][i]['geometry'])
latmin, lonmin, latmax, lonmax = polygon.bounds

# construct a rectangular mesh
points = []
for lat in np.arange(latmin, latmax, resolution):
    for lon in np.arange(lonmin, lonmax, resolution):
        points.append(Point((round(lat,4), round(lon,4))))

# validate if each point falls inside shape
valid_points.extend([i for i in points if polygon.contains(i)])

Frankfort

Are the shapes always convex? If they are you can probably just check the outer points. – Displode

Does that polygon.contains() accept numpy arrays? If it does, you can use numpy.meshdrid to create two matrix for x and y coordinates, witch would be faster than iterate – Mascot

np.meshgrid IS faster, however shape.contains does NOT take numpy arrays it would seem. you have sped up my software. – Frankfort

Related questions: Get all lattice points lying inside a Shapely polygon, Extract interior points of polygon given the exterior coordinates. – Poulenc

what library the shape function belongs to ? – Accipiter

Answer 1

I saw that you answered your question (and seems to be happy with using intersection) but also note that shapely (and the underlying geos library) have prepared geometries for more efficient batch operations on some predicates (contains, contains_properly, covers, and intersects). See Prepared geometry operations.

Adapted from the code in your question, it could be used like so:

from shapely.prepared import prep

# determine maximum edges
polygon = shape(geojson['features'][i]['geometry'])
latmin, lonmin, latmax, lonmax = polygon.bounds

# create prepared polygon
prep_polygon = prep(polygon)

# construct a rectangular mesh
points = []
for lat in np.arange(latmin, latmax, resolution):
    for lon in np.arange(lonmin, lonmax, resolution):
        points.append(Point((round(lat,4), round(lon,4))))

# validate if each point falls inside shape using
# the prepared polygon
valid_points.extend(filter(prep_polygon.contains, points))

Caspar

Answer 2

Oh why hell yes. Use the intersection method of shapely.

polygon = shape(geojson['features'][i]['geometry'])
lonmin, latmin, lonmax, latmax = polygon.bounds

# construct rectangle of points
x, y = np.round(np.meshgrid(np.arange(lonmin, lonmax, resolution), np.arange(latmin, latmax, resolution)),4)
points = MultiPoint(list(zip(x.flatten(),y.flatten())))

# validate each point falls inside shapes
valid_points.extend(list(points.intersection(polygon)))

Frankfort

Answer 3

The best i can think is do this:

X,Y = np.meshgrid(np.arange(latmin, latmax, resolution),
              np.arange(lonmin, lonmax, resolution))

#create a iterable with the (x,y) coordinates
points = zip(X.flatten(),Y.flatten())

valid_points.extend([i for i in points if polygon.contains(i)])

Mascot

Answer 4

if you want to generate n points in a shapely.geometry.Polygon, there is a simple iterative function to do it. Manage tol (tolerance) argument to speed up the points generation.

import numpy as np
from shapely.geometry import Point, Polygon


def gen_n_point_in_polygon(self, n_point, polygon, tol = 0.1):
    """
    -----------
    Description
    -----------
    Generate n regular spaced points within a shapely Polygon geometry
    -----------
    Parameters
    -----------
    - n_point (int) : number of points required
    - polygon (shapely.geometry.polygon.Polygon) : Polygon geometry
    - tol (float) : spacing tolerance (Default is 0.1)
    -----------
    Returns
    -----------
    - points (list) : generated point geometries
    -----------
    Examples
    -----------
    >>> geom_pts = gen_n_point_in_polygon(200, polygon)
    >>> points_gs = gpd.GeoSeries(geom_pts)
    >>> points_gs.plot()
    """
    # Get the bounds of the polygon
    minx, miny, maxx, maxy = polygon.bounds    
    # ---- Initialize spacing and point counter
    spacing = polygon.area / n_point
    point_counter = 0
    # Start while loop to find the better spacing according to tolérance increment
    while point_counter <= n_point:
        # --- Generate grid point coordinates
        x = np.arange(np.floor(minx), int(np.ceil(maxx)), spacing)
        y = np.arange(np.floor(miny), int(np.ceil(maxy)), spacing)
        xx, yy = np.meshgrid(x,y)
        # ----
        pts = [Point(X,Y) for X,Y in zip(xx.ravel(),yy.ravel())]
        # ---- Keep only points in polygons
        points = [pt for pt in pts if pt.within(polygon)]
        # ---- Verify number of point generated
        point_counter = len(points)
        spacing -= tol
    # ---- Return
    return points

Tadio