CODEWITHSUNDEEP
numpyshapely
asdflkjwqerqwr
asdflkjwqerqwr

Reputation: 1189

Removing numpy meshgrid points outside of a Shapely polygon

I have a 10 x 10 grid that I would like to remove points outside of a shapely Polygon:

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

gridX, gridY = np.mgrid[0.0:10.0, 0.0:10.0]
poly = Polygon([[1,1],[1,7],[7,7],[7,1]])

#plot original figure
fig = plt.figure()
ax = fig.add_subplot(111)
polyp = PolygonPatch(poly)
ax.add_patch(polyp)
ax.scatter(gridX,gridY)
plt.show()

Here is the resulting figure: original fig

And what I want the end result to look like: end

I know that I can reshape the array to a 100 x 2 array of grid points:

stacked = np.dstack([gridX,gridY])
reshaped = stacked.reshape(100,2)

I can see if the point lies within the polygon easily:

for i in reshaped:
    if Point(i).within(poly):
         print True

But I am having trouble taking this information and modifying the original grid

Upvotes: 4

Views: 2233

Answers (2)

hamflow
hamflow

Reputation: 321

with regards to @pseudocubic answer, I want to add an answer where you can include boundary condition too. the method adds, using shapely >> touches and instead of appending eligible points, excluding ineligible points, and at the end removing all empty points for both dimensions:

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

gridX, gridY = np.mgrid[0.0:10.0, 0.0:10.0]
poly = Polygon([[1,1],[1,7],[7,7],[7,1]])
stacked = np.dstack([gridX,gridY])
reshaped = stacked.reshape(100,2)
points = reshaped.tolist()
for i, point in enumerate(points):
        point_geom = Point([point])
        if not poly.contains(point_geom) and not poly.touches(point_geom): ## (x != x_min and x != x_max and y != y_min and y != y_max):  # second condition ensures within boundary
            points[i] = ([],[])  # Mark the point as None if it's outside the polygon

mod_points = [[coord for coord in point if coord] for point in points]  
mod_points = [point for point in mod_points if point is not None and point != []]

And for plotting,

#plot original points
fig = plt.figure()
ax = fig.add_subplot(111)
# Extract the x and y coordinates of polygon
poly_x, poly_y = poly.exterior.xy
ax.plot(poly_x, poly_y, c = "green")
#ploting modified points
ax.scatter(gridX,gridY)
mod_x, mod_y = zip(*mod_points)
ax.scatter(mod_x,mod_y , c = 'red')
plt.show()

points inside polygon

Upvotes: 0

pseudocubic
pseudocubic

Reputation: 1039

You're pretty close already; instead of printing True, you could just append the points to a list.

output = []
for i in reshaped:
    if Point(i).within(poly):
        output.append(i)

output = np.array(output)
x, y = output[:, 0], output[:, 1]

It seems that Point.within doesn't consider points that lie on the edge of the polygon to be "within" it though.

Upvotes: 3

Related Questions