Python确定数组中网格上多边形内的平均值并提取最大值

Question

我有 3 个 NumPy 数组，其中包含 UTM-X(256) 和 UTM-Y(256) 坐标，以及 UTM 中天气雷达 256x256 (km) 的累积降雨量 (65536)。

我在网格边界内还有一个多边形，它是 UTM 中的集水边界。

我需要确定流域多边形（RADAR 数据的裁剪子集）的平均降雨量、最大值和最大值的位置。我已经确定了整个雷达网格的平均值。

所以问题是：如何对由 Polygon 确定的 NumPy 数组的子集执行分析？我原以为这将是一个非常常见的操作，但没有找到任何 Python 脚本来执行此操作。

以下是数据集的示意图：

Answer 1

这是一种可能的方法的概述。

首先找到界定流域边界的多边形。假设您知道您的全套点的 UTM 坐标中的哪一个形成了该集水区边界，假设是这样的，

catchment = (UTM_X, UTM_Y) 点元组的 np.array

您可以使用 scipy.spatial.ConvexHull

找到该点集的边界

boundary=scipy.spatial.ConvexHull(catchment)

接下来，对于您的降雨数据数组，您必须测试坐标是落在凸包边界之内还是之外。

This previous SO question 有一些很好的答案来解释如何进行坐标测试。

最后，您将收集那些通过边界内测试的降雨数据点，并使用适当的 NumPy/SciPy 统计函数执行您想要执行的任何统计计算。

Answer 2

假设边界是作为多边形顶点列表给出的，您可以让 matplotlib 在数据坐标上为您生成一个掩码，然后使用该掩码仅对轮廓内的值求和。

换句话说，当您有一系列坐标来定义标记感兴趣区域的多边形边界时，让 matplotlib 生成一个布尔掩码，指示该多边形内的所有坐标。然后可以使用此掩码仅提取轮廓内有限的降雨数据集。

以下简单示例向您展示了这是如何完成的：

import numpy as np
from matplotlib.patches import PathPatch
from matplotlib.path import Path
import matplotlib.pyplot as plt

# generate some fake data
xmin, xmax, ymin, ymax = -10, 30, -4, 20 
y,x = np.mgrid[ymin:ymax+1,xmin:xmax+1]
z = (x-(xmin+xmax)/2)**2 + (y-(ymin + ymax)/2)**2
extent = [xmin-.5, xmax+.5, ymin-.5, ymax+.5]
xr, yr = [np.random.random_integers(lo, hi, 3) for lo, hi
         in ((xmin, xmax), (ymin, ymax))] # create a contour

coordlist = np.vstack((xr, yr)).T  # create an Nx2 array of coordinates
coord_map = np.vstack((x.flatten(), y.flatten())).T # create an Mx2 array listing all the coordinates in field

polypath = Path(coordlist)
mask = polypath.contains_points(coord_map).reshape(x.shape) # have mpl figure out which coords are within the contour

f, ax = plt.subplots(1,1)
ax.imshow(z, extent=extent, interpolation='none', origin='lower', cmap='hot')
ax.imshow(mask, interpolation='none', extent=extent, origin='lower', alpha=.5, cmap='gray')
patch = PathPatch(polypath, facecolor='g', alpha=.5)
ax.add_patch(patch)
plt.show(block=False)
print(z[mask].sum())  # prints out the total accumulated

在此示例中，x 和 y 代表您的 UTM-X 和 UTM-Y 数据范围。 z 代表天气降雨量数据，但在这种情况下是一个矩阵，与平均降雨量的单列视图不同（很容易重新映射到网格上）。

在最后一行，我总结了轮廓内z 的所有值。如果您想要平均值，只需将 sum 替换为 mean。