python - 通过高程网格查找地面大致面积的算法？

Question

我有一个 10mX10m“像素”的 200X200 网格/阵列(即 40,000 个值)，其值表示该点土地的平均海拔。

网格上有广阔的连接区域，其中高程值为 0 米，因为它们代表实际的海洋。

问题:有没有快速求出土地大概面积的算法？据我所知，我可能会乘以 200^2*10^2 以获得面积的粗略近似值，但有些值差异很大。

我想我知道一种相当昂贵的方法，即通过对顶点位于高程的所有三角形求和。但是有更快/更简单的方法吗？

Answer 1

NumPy和 SciPy是解决这类问题的工具。这是一个 200×200 的合成景观，点位于 10 米的网格上，高度范围可达海拔 40 米:

>>> import numpy as np
>>> xaxis = yaxis = np.arange(0, 2000, 10)
>>> x, y = np.meshgrid(xaxis, yaxis)
>>> z = np.maximum(40 * np.sin(np.hypot(x, y) / 350), 0)

我们可以在 Matplotlib 中查看这个:

>>> import matplotlib.pyplot as plt
>>> import mpl_toolkits.mplot3d.axes3d as axes3d
>>> _, axes = plt.subplots(subplot_kw=dict(projection='3d'))
>>> axes.plot_surface(x, y, z, cmap=plt.get_cmap('winter'))
>>> plt.show()

现在，陆地上的点数(即高度大于 0)计算起来很简单，您可以将其乘以网格正方形的大小(在您的问题中为 100 平方米)来估算土地面积:

>>> (z > 0).sum() * 100
1396500

但是从这个问题中，我了解到您想要一个更准确的估计，一个考虑到土地坡度的估计。一种方法是制作覆盖土地的三角形网格，然后将三角形的面积相加。

首先，将坐标数组转换为点数组(point cloud):

>>> points = np.vstack((x, y, z)).reshape(3, -1).T
>>> points
array([[  0.000000e+00,   0.000000e+00,   0.000000e+00],
       [  1.000000e+01,   0.000000e+00,   1.142702e+00],
       [  2.000000e+01,   0.000000e+00,   2.284471e+00],
       ..., 
       [  1.970000e+03,   1.990000e+03,   3.957136e+01],
       [  1.980000e+03,   1.990000e+03,   3.944581e+01],
       [  1.990000e+03,   1.990000e+03,   3.930390e+01]])

其次，使用scipy.spatial.Delaunay在二维中进行三角剖分，得到表面网格:

>>> from scipy.spatial import Delaunay
>>> tri = Delaunay(points[:,:2])
>>> len(tri.simplices)
79202
>>> tri.simplices
array([[39698, 39899, 39898],
       [39899, 39698, 39699],
       [39899, 39700, 39900],
       ..., 
       [19236, 19235, 19035],
       [19437, 19236, 19237],
       [19436, 19236, 19437]], dtype=int32)

三角剖分中每个三角形的值是三角形中三个点的 points 数组中的索引。

第三，选择其中有一些土地的三角形:

>>> land = (points[tri.simplices][...,2] > 0).any(axis=1)
>>> triangles = tri.simplices[land]
>>> len(triangles)
27858

第四，计算这些三角形的面积:

>>> v0, v1, v2 = (points[triangles[:,i]] for i in range(3))
>>> areas = np.linalg.norm(np.cross(v1 - v0, v2 - v0), axis=1) / 2
>>> areas
array([ 50.325028,  50.324343,  50.32315 , ...,  50.308673,  50.313157, 50.313649])

最后，把它们加起来:

>>> areas.sum()
1397829.2847141961

这与最初的估计差别不大，这是意料之中的，因为斜坡很浅。