10

I have seen this question asked but have not really been able to find a full response. I have a simple shapely polygon, called polygon. I would like to extract this polygon as a binary mask (ideally a numpy array). How would I go about doing this?

I have also managed to convert from shapely to geopandas as shown here so extracting a mask from geopandas would work as well, but I have not really been able to find a thread on this unfortunately.

EDIT: To be clear, if I am to instead use a coordinate grid, my grid contains x and y cartesian coordinates (unordered) corresponding to the points that make the contour of the shape. These are floats, so solutions that require int inputs won't quite work. Ideally I would like the starting point to be a shapely polygon instead of a set of points, but I can use an unordered set of points instead if this is preferable (or alternately somehow extract clockwise vertices from a shapely polygon)

I have tried Yusuke's method described here But the mask I get does not quite make sense.

Yusuke's method:

#%% create grid and plot
nx, ny = 100, 100
poly_verts = Plane1verts #this is a list of tuples containing cartesian coordinate pairs of the shape contour in x and y
# Create vertex coordinates for each grid cell...
# (<0,0> is at the top left of the grid in this system)
x, y = np.meshgrid(np.arange(nx), np.arange(ny))
x, y = x.flatten(), y.flatten()

points = np.vstack((x,y)).T

path = Path(poly_verts)
grid = path.contains_points(points)
grid = grid.reshape((ny,nx))

plt.imshow(grid)
plt.title('Grid plot')
plt.show()

the resulting plot of the mask is enter image description here

Which is not what I expected. Whereas plotting from geopandas as described below shows the correct shape.

#%% create shapely and plot for comparison
from shapely.geometry import Polygon
#convert the sets of points dict to a shapely object
polygon1_plane1=Polygon(Plane1vert_tuple)

p = gpd.GeoSeries(polygon1_plane1)
p.plot()
plt.show()

Resulting in the plot enter image description here

EDIT2: here is a copy of the coordinate grid I am using as a list of tuples

[(-8.982, -12.535), (-7.478, -12.535), (-5.975, -12.535), (-4.471, -12.535), (-4.471, -12.535), (-2.967, -11.031), (-1.463, -11.031), (-1.463, -11.031), (0.041, -9.527), (0.041, -9.527), (1.544, -8.023), (3.048, -8.023), (4.552, -8.023), (4.552, -8.023), (6.056, -6.52), (7.559, -6.52), (7.559, -6.52), (7.559, -5.016), (9.063, -3.512), (10.567, -3.512), (10.567, -3.512), (10.567, -2.008), (10.567, -0.505), (10.567, 0.999), (10.567, 2.503), (10.567, 4.007), (10.567, 4.007), (9.063, 5.51), (9.063, 5.51), (7.559, 7.014), (7.559, 7.014), (6.056, 8.518), (6.056, 8.518), (4.552, 10.022), (4.552, 11.526), (4.552, 11.526), (3.048, 11.526), (1.544, 11.526), (1.544, 11.526), (1.544, 10.022), (0.041, 8.518), (0.041, 8.518), (0.041, 7.014), (-1.463, 5.51), (-2.967, 5.51), (-4.471, 5.51), (-4.471, 5.51), (-5.975, 4.007), (-7.478, 4.007), (-8.982, 4.007), (-10.486, 4.007), (-11.99, 4.007), (-13.493, 4.007), (-13.493, 4.007), (-14.997, 2.503), (-14.997, 2.503), (-16.501, 0.999), (-18.005, 0.999), (-18.005, 0.999), (-18.005, -0.505), (-19.508, -2.008), (-19.508, -2.008), (-19.508, -3.512), (-19.508, -5.016), (-19.508, -5.016), (-18.005, -6.52), (-18.005, -8.023), (-18.005, -8.023), (-16.501, -9.527), (-16.501, -9.527), (-14.997, -9.527), (-13.493, -11.031), (-13.493, -11.031), (-11.99, -11.031), (-10.486, -12.535), (-10.486, -12.535)]
hex93
  • 319
  • 3
  • 15
  • Does this answer your question? [SciPy Create 2D Polygon Mask](https://stackoverflow.com/questions/3654289/scipy-create-2d-polygon-mask) – Georgy May 26 '21 at 15:56
  • I have seen that thread but unfortunately not as my coordinate grid is of floats instead of ints |(including coordinates in negative x and y directions). I have a set of cartesian coordinates, not pixel coordinates, and so it seems that the suggested solution is not working. My plot using Yusuke's method vs directly plotting with shapely does not agree. For some reason I am getting a very different shape. – hex93 May 27 '21 at 10:27
  • are you working with geopandas? – j-i-l May 27 '21 at 10:44
  • I am, though currently I am only doing this to plot my contour and verify it is behaving properly. – hex93 May 27 '21 at 10:56

2 Answers2

13

rasterio.features.rasterize sounds like exactly what you are looking for.

from shapely.geometry import Polygon
import rasterio.features
import matplotlib.pyplot as plt

poly = Polygon([(0, 50), (10, 10), (30, 0), (45, 45), (0, 50)])
img = rasterio.features.rasterize([poly], out_shape=(60, 50))
plt.imshow(img)

rasterized image

Martin Valgur
  • 5,793
  • 1
  • 33
  • 45
  • Would this work if the points on the polygon are not in any particular order? Mine are unordered unfortunately and I am having the same issue as in my original post, i.e. the gpd plot looks identical to using Yusuke's method. – hex93 May 27 '21 at 10:55
  • I have also tried feeding in the polygon directly to rasterize but I get the same strange looking shape from the first plot of my original post. – hex93 May 27 '21 at 11:02
  • You could try applying `.convex_hull` to the polygons / set of points? – Martin Valgur May 27 '21 at 11:10
  • Gave that a go, the shape changes a bit but I still have the same problem that it only takes up a small portion of the output at the top left and the shape does not match that on geopandas. Strange! (thanks for the suggestions so far by the way) – hex93 May 27 '21 at 11:13
  • You can try using the `transform` param of `rasterize()` to transform the polygon to the image coordinate range. The rasterized polygon appears in the top left because that's where the coordinates of your polygon vertices are located. – Martin Valgur May 27 '21 at 11:17
  • That might work, what would I feed to this param? The default from docs looks to be transform=Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), I'm not sure what I would change it to in my case. – hex93 May 27 '21 at 11:21
  • This will work for a single polygon, good job. For multi-polygon or a whole GeoPandas while you would need a wrapper or another approach. – Олег Місько Jan 20 '23 at 16:50
1

You can use OpenCV with cv2.fillPoly()

Example:

import numpy as np
import cv2
from shapely.geometry import box


mask = np.zeros([960, 1280])
# create an example bounding box polygon
x1, y1, x2, y2 = 480, 540, 780, 840
polygon = box(x1, y1, x2, y2)
points = [[x, y] for x, y in zip(*polygon.boundary.coords.xy)]

mask = cv2.fillPoly(mask, np.array([points]).astype(np.int32), color=1)
biendltb
  • 1,149
  • 1
  • 13
  • 20