1

I have a list of boundary coordinates for a polygon. I want to extract some interior points from this polygon based on some sort of grid spacing that I can choose.

For example, say I have this:

from shapely.geometry import MultiPoint
Polygon([(0,0),
        (1,0),
        (1,1),
        (0,1)])

enter image description here

This yields a simple square. However, we can imagine this square is actually a grid filled with an infinite number of points. For this example, let's say I want a list of points inside the square at 0.1 grid spacing. Is there a way to extract a list of these points?

For bonus points, how would we not only get the interior points, but also all of the boundary points at 0.1 grid spacing?

Edit: This question was mainly posed for non-rectangular shapes and has the purpose of finding points for shapes of counties/provinces/districts, etc. I just used a rectangle as a simple example. However, I accepted the proposed solution. Will give out upvotes to anyone who provides a more complex solution for non-rectangular objects

Eli Turasky
  • 981
  • 2
  • 11
  • 28
  • 1
    (Bonus points) In the last paragraph, they call it bounties here. ;p – swatchai Jan 22 '21 at 18:35
  • 1
    Does this answer your question? [Get all lattice points lying inside a Shapely polygon](https://stackoverflow.com/questions/44399749/get-all-lattice-points-lying-inside-a-shapely-polygon) – Georgy Jan 23 '21 at 15:57
  • 1
    do you need something that would generalize to an arbitrary shape? What would the solution for a hexagon look like? – Paul H Jan 24 '21 at 16:44
  • Yes. It would be nice to have a solution that works for non-rectangular shapes. – Eli Turasky Jan 25 '21 at 13:29

1 Answers1

1
  • Regarding to the extraction of the interior points, based on a grid shape, it can be done with simple math, like so :
from shapely.geometry import Polygon, MultiPoint
from math import ceil

# The setup
pol = Polygon([(0,0), (1,0), (1,1), (0,1)])
cell_size = 0.1

# Extract coordinates of the bounding box:
minx, miny, maxx, maxy = pol.bounds

# How many points in each dimension ?
rows = int(ceil((maxy - miny) / cell_size))
cols = int(ceil((maxx - minx) / cell_size))

# Actually generate the points:
pts = []
x, y = minx, miny
for countcols in range(cols):
    for countrows in range(rows):
        pts.append((x, y))
        y += cell_size
    x += cell_size
    y = miny

# Create the MultiPoint from this list
result_pts_interior = MultiPoint(pts)

Result:

Result1

But note that if the polygon does not have a rectangular shape, it will also generate points which will be outside the polygon (you should then test if they intersect the polygon with pol.intersects(Point(x, y)) before adding them to the list).


  • Regarding to the extraction of boundary points, it can be done using the interpolate method of shapely LineStrings, like so :
from shapely.geometry import Polygon, MultiPoint
import numpy as np

# The setup
pol = Polygon([(0,0), (1,0), (1,1), (0,1)])
distance_between_pts = 0.1

boundary = pol.boundary # Boundary of polygon as a linestring
boundary_length = boundary.length # Its length
# Build a list of points spaced by 0.1 along this linestring:
pts_boundary = [
    boundary.interpolate(n, False) for n
    in np.linspace(0, boundary_length, int(boundary_length / distance_between_pts) + 1)
]
# Create the MultiPoint from this list
result_pts_boundary = MultiPoint(pts_boundary)

Result:

Result2

mgc
  • 5,223
  • 1
  • 24
  • 37