本文旨在解决在Python 2D NumPy数组中通过多边形区域填充数值时常见的坐标系混淆问题。我们将深入探讨“真实世界”坐标与数组索引坐标之间的差异，并提供一种统一且高效的方法，通过在同一坐标系统（通常是数组索引）中定义多边形顶点并检查点，从而正确地修改数组特定区域的值。文章将通过示例代码演示如何避免常见错误，并讨论性能优化及进阶库的使用。
理解问题：2D数组中的坐标系统与多边形定义
在Python中使用NumPy数组处理二维数据时，我们通常会遇到两种主要的坐标系统，这两种系统在概念上有所不同，但在实际操作中很容易混淆：
- “真实世界”坐标系统 (Real-world Coordinates): 这通常与物理尺寸、空间位置或连续几何空间相关联，例如一个从 (0, 0) 延伸到 (Lx, Ly) 的平面。在代码中，这类坐标常通过 np.linspace 生成的等间隔点和 np.meshgrid 创建的网格 X, Y 来表示。
- 数组索引坐标系统 (Array Index Coordinates): 这直接对应于NumPy数组的离散行和列索引。例如，一个 (Ny+1, Nx+1) 大小的数组，其索引范围从 (0, 0) 到 (Ny, Nx)。
当目标是在2D NumPy数组中定义一个几何形状（如不规则多边形）并基于此形状修改数组元素时，一个常见的错误是混合使用这两种坐标系统。例如，多边形的顶点可能在一个坐标系中定义（例如，使用数组索引），而用于判断点是否在多边形内部的待检查点却来自另一个坐标系（例如，使用“真实世界”坐标）。这种坐标系统的不一致性会导致多边形无法正确地在数组中被识别和填充。
考虑以下示例，我们希望在一个初始化为零的 s 2D数组中，将一个由数组索引定义的矩形区域内的值设置为 100：
立即学习“Python免费学习笔记（深入）”；
import numpy as np
import matplotlib.path as mpath
import pandas as pd
# 变量声明：定义“真实世界”空间
Lx = 1.7
Ly = 1.7
Nx = 170
Ny = 170
x = np.linspace(0, Lx, Nx + 1)
y = np.linspace(0, Ly, Ny + 1)
X, Y = np.meshgrid(x, y) # 生成“真实世界”坐标网格
# s 数组初始化
s = np.zeros((Ny + 1, Nx + 1))
# 定义多边形顶点：这里使用数组索引 (列, 行)
polygon_vertices = [(10, 10), (50, 10), (50, 70), (10, 70)]
polygon_path = mpath.Path(polygon_vertices)
# 尝试修改 s 数组：存在坐标系混淆
for i in range(Ny + 1):
for j in range(Nx + 1):
# 问题所在：多边形顶点是基于数组索引定义的，
# 但此处却使用了“真实世界”坐标 (X[i, j], Y[i, j]) 进行点判断
if polygon_path.contains_point((X[i, j], Y[i, j])):
s[i, j] = 100
# df = pd.DataFrame(s)
# df.to_excel('s_incorrect.xlsx', index=False) # 导出查看结果，会发现 s 数组仍为全零在上述代码中，polygon_vertices 列表中的元组 (10, 10)、(50, 10) 等显然是意图表示数组的列和行索引。然而，在循环内部，contains_point 方法却接收了 (X[i, j], Y[i, j])，这些值是根据 Lx, Ly 和 np.meshgrid 生成的“真实世界”坐标。由于 matplotlib.path.Path 实例 polygon_path 是基于数组索引定义的，当它尝试判断一个“真实世界”坐标点是否包含在其内部时，会因坐标系不匹配而始终返回 False，导致 s 数组无法被正确修改。
解决方案：统一坐标系统进行多边形填充
解决此问题的核心在于保持坐标系统的一致性。如果多边形的顶点是基于数组索引定义的，那么在判断数组中的每个点是否位于多边形内部时，也应该使用其对应的数组索引。这意味着我们可以直接将数组的 (列索引, 行索引) 作为点传递给 contains_point 方法。
以下是修正后的代码示例：
import numpy as np
import matplotlib.path as mpath
import pandas as pd
# 变量声明
# 当仅使用数组索引进行操作时，Lx, Ly, x, y, X, Y 等“真实世界”坐标相关的变量不再需要
Nx = 170
Ny = 170
# s 数组初始化
s = np.zeros((Ny + 1, Nx + 1))
# 定义多边形顶点：保持使用数组索引 (列索引, 行索引)
polygon_vertices = [(10, 10), (50, 10), (50, 70), (10, 70)]
polygon_path = mpath.Path(polygon_vertices)
# 修改 s 数组：使用统一的数组索引坐标系
for i in range(Ny + 1): # 遍历行索引
for j in range(Nx + 1): # 遍历列索引
# 修正：直接使用数组索引 (j, i) 进行点判断
# 注意：mpath.Path 默认接受 (x, y) 形式，对应于数组的 (列索引, 行索引)
if polygon_path.contains_point((j, i)):
s[i, j] = 100 # 将指定区域值设置为 100
# df = pd.DataFrame(s)
# df.to_excel('s_correct.xlsx', index=False) # 导出查看正确结果，会发现多边形区域已被填充通过将 if polygon_path.contains_point((X[i, j], Y[i, j])) 修改为 if polygon_path.contains_point((j, i))，我们确保了多边形定义和点检查都发生在相同的数组索引坐标系中。这样，多边形区域内的 s 数组元素将被正确地设置为 100。
代码解析与优化
- 坐标系统一致性： matplotlib.path.Path 接受的顶点列表 [(x1, y1), (x2, y2), ...] 默认是笛卡尔坐标系中的 (x, y)。在NumPy数组的上下文里，这通常对应于 (列索引, 行索引)。因此，当遍历 s[i, j] 时，j 是列索引（对应x轴），i 是行索引（对应y轴）。所以 polygon_path.contains_point((j, i)) 是正确的用法。
- 冗余变量移除： 由于我们完全在数组索引坐标系中操作，原先用于定义“真实世界”坐标的 Lx, Ly, x, y, X, Y 等变量变得不再必要，可以移除以简化代码，提高可读性。
- 循环效率： 尽管上述解决方案功能正确，但对于大型 Ny 和 Nx 值，嵌套的 for 循环和 contains_point 的逐点调用可能会非常慢。matplotlib.path.Path.contains_points 方法可以一次性检查多个点，但即便如此，生成所有点的坐标并传递给它仍然涉及大量计算。
性能考量与进阶方法
对于需要更高性能的场景，尤其是在处理大型2D数组时，可以考虑使用专门的图像处理库，它们通常提供了更优化的算法来处理几何形状的绘制和填充。
-
使用 skimage.draw.polygon (Scikit-image):skimage.draw 模块提供了多种高效的函数用于在图像（NumPy数组）上绘制形状。skimage.draw.polygon 函数可以直接返回多边形内部所有像素的行和列索引，效率远高于逐点检查。
from skimage.draw import polygon import numpy as np import pandas as pd Ny, Nx = 171, 171 # 数组尺寸，与 (Ny+1, Nx+1) 对应 s_optimized = np.zeros((Ny, Nx)) # 定义多边形顶点 (行索引列表, 列索引列表) # 注意：skimage.draw.polygon 期望 (行索引列表, 列索引列表)， # 而我们之前使用的 (x, y) 形式是 (列索引, 行索引)，所以这里需要转换 r_coords = np.array([10, 10,
