定义2D对象并将其区域用作布尔值

Question

我定义了两个空间维度（x和z），并且我能够手动“绘制”一个对象以将其用作解决方程的布尔值。我将其定义如下：

A = np.zeros((nz,nx))
object = np.ones_like(A)
object[ int(5/dz):int(10/dz) , int(5/dx):int(10/dz) ] = 2
object = object == 2

通过这样做，我想在z维度中定义一个5x10的正方形，在x维度中定义一个5x10的正方形，然后应用将其理解为一个区域的algorythim。但是，当要绘制复杂区域时，最终很难通过小正方形和矩形来完成。

因此，我想通过单击鼠标来自动生成区域，并且希望能够将该区域用作布尔值。

我能够使用：绘制多边形：

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import Polygon

fig, ax = plt.subplots()

object = np.array(plt.ginput(n=-100,mouse_stop=2))
p = Polygon(object, alpha=0.5)
plt.gca().add_artist(p)


plt.draw()
plt.show()

但是这会输出顶点的z和x坐标，我试图将其用作boleean，但我无法编写它，以便python将其理解为这些点定义的区域。

这个问题容易解决吗？

Answer 1

[如果只想计算普通多边形的面积，则可以使用Shapely python包，例如：

import numpy as np
import matplotlib.pyplot as plt
from shapely.ops import Polygon
from matplotlib.patches import Polygon as PltPolygon

# Get the coordinate input
canvas_size = np.array([1, 1])
canvas_lim = np.array([[0, canvas_size[0]], [0, canvas_size[1]]])
fig, ax = plt.subplots()
plt.xlim(canvas_lim[0])
plt.ylim(canvas_lim[1])
ax.set_aspect("equal")
coordinates = np.array(plt.ginput(n=-100, mouse_stop=2))

# Use shapely.ops.Polygon to calculate the area
poly = Polygon(coordinates)

area = poly.area
print("The area is {} units^2".format(area))

# Draw the polygon
p = PltPolygon(coordinates, alpha=0.5)
ax.add_artist(p)
plt.show()

如果确实需要遮罩，这是使用numpy和matplotlib.path对其进行光栅化的一种方法。有关详细信息，请参见代码中的注释：

import numpy as np
import matplotlib.path as mpltPath
import matplotlib.pyplot as plt

# Define the limits of our polygon
canvas_desired_size = np.array([110, 100])
# The pixel size with which we calculate (number of points to consider)
# The higher this number, the more we have to calculate, but the
# closer the approximation will be
pixel_size = 0.1

# Cacluate the actual size of the canvas
num_pxiels = np.ceil(canvas_desired_size / pixel_size).astype(int)
canvas_actual_size = num_pxiels * pixel_size

# Let's create a grid where each pixel's value is it's position in our 2d image
x_coords = np.linspace(
    start=0,
    stop=canvas_actual_size[0],
    endpoint=False,
    num=canvas_desired_size[0] / pixel_size,
)
y_coords = np.linspace(
    start=0,
    stop=canvas_actual_size[1],
    endpoint=False,
    num=canvas_desired_size[1] / pixel_size,
)
# Since it makes more sense to check if the middle of the pixel is in the
# polygion, we shift everything with half pixel size
pixel_offset = pixel_size / 2
x_centers = x_coords + pixel_offset
y_centers = y_coords + pixel_offset

xx, yy = np.meshgrid(x_centers, y_centers, indexing="ij")

# Flatten our xx and yy matrixes to an N * 2 array, which contains
# every point in our grid
pixel_centers = np.array(
    list(zip(xx.flatten(), yy.flatten())), dtype=np.dtype("float64")
)

# Now prompt for the imput shape
canvas_lim = np.array([[0, canvas_actual_size[0]], [0, canvas_actual_size[1]]])
fig, ax = plt.subplots()
plt.xlim(canvas_lim[0])
plt.ylim(canvas_lim[1])
ax.set_aspect("equal")
shape_points = np.array(plt.ginput(n=-100, mouse_stop=2))

# Create a Path object
shape = mpltPath.Path(shape_points)
# Use Path.contains_points to calculate if each point is
# within our shape
shape_contains = shape.contains_points(pixel_centers)

# Reshape the result to be a matrix again
mask = np.reshape(shape_contains, num_pxiels)

# Calculate area
print(
    "The shape area is roughly {} units^2".format(
        np.sum(shape_contains) * pixel_size ** 2
    )
)

# Show the rasterized shape to confirm it looks correct
plt.imshow(np.transpose(mask), aspect="equal", origin="lower")
plt.xlim([0, num_pxiels[0]])
plt.ylim([0, num_pxiels[1]])
plt.show()

或者，更简单的解决方案是using your plot as an image并将其阈值化以获得布尔掩码。应该在Google上提供许多示例说明。