寻找地面近似面积的算法，基于高程网格？

Question

我有一个200x200的网格/数组（即40,000个值），每个"像素"的大小为10m x 10m，其值表示该点的平均海拔高度。

在整个网格中，有广阔的连通区域，其高度值为0m，因为它们代表实际上的海洋。

问题：是否存在一种快速算法来获取陆地的近似面积？我知道可以将200^2 x 10^2相乘以获得大致的面积近似值，但其中一些值变化很大。

我认为我知道一种相当昂贵的方法，即求出所有顶点高程相同的三角形的总面积。但是否有更快/更简单的方法？

- bjd2385

---

@lmiguelvargasf 这不是一个课堂作业。我是数学专业的，而不是计算机科学专业的学生。我正在为一架无人机（UAV俱乐部）编写程序/脚本，并且我已经使用另一个我编写的脚本通过Google Maps API获取了我提到的数据。这是一个（Python）类的一部分，它返回有关无人机飞行地面的信息。 - bjd2385

我无法确定任务是排除海洋瓦片（因为海洋！=“陆地”，否则为什么要提到海洋），还是计算所有瓦片的“三维面积”，或者两者兼而有之。如果您想要3D区域，请记住，除非您对地面的陡峭程度有一个界限，否则任何近似都可能是任意糟糕的。 - j_random_hacker

那么您是想简单地计算出数组中土地瓦片（level > 0）的数量，并将其乘以10平方米吗？这将近似于在投影到地面平面后的面积（即忽略高度，仅用于确定是否为水域）。 - Byte Commander

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

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

首先，将坐标数组转换为点数组（点云）：

>>> 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

这与原始估计值并没有太大差异，因为斜率很缓。

- Gareth Rees

---

太完美了！非常感谢！对于最后的计算（即计算面积），你似乎使用了鞋带公式（也在Widder的高等微积分中），这太棒了；我找到了它，但我怀疑我不会像你一样高效地计算出来。 - bjd2385

Answer 2

首先是一些有用的测试附加内容：

# a function to create a random map as simulated input for testing:
def get_map(x_size, y_size, h_min, h_max):
    import random
    return [[random.randint(h_min, h_max) for x in range(x_size)] for y in range(y_size)]

# a function to nicely print the map for debug and visualization
def print_map(hmap):
    print(*hmap, sep="\n")

然后我们编写实际的土地面积计算器：

# calculate approximate land area where the height is greater than zero
# map is a list of lists, tile_size is in m², min_level is the sea level
def calc_land_area(hmap, tile_size=100, min_level=0):
    land_tiles = sum(len([tile for tile in row if tile>min_level]) for row in hmap)
    return tile_size * land_tiles

现在进行一项测试：

hmap = get_map(5, 5, 0, 3)
print_map(hmap)
print("land area:", calc_land_area(hmap), "m²")

那可能会导致例如这个随机示例输出：

[2, 0, 3, 0, 2]
[3, 0, 0, 0, 2]
[1, 0, 0, 1, 2]
[3, 3, 3, 2, 1]
[3, 1, 1, 3, 0]
land area: 1700 m²

你在 25 个地图瓦片上看到了 8 个海洋瓦片，因此有 800 平方米是海洋，1700 平方米是陆地。

点击此处在 ideone.com 上查看此代码的运行结果

- Byte Commander

---

不要使用len([tile for tile in row if tile>min_level])这种方式，它会创建一个列表然后立即丢弃。可以改为使用sum(tile > min_level for tile in row)的方式。 - Gareth Rees

@GarethRees 好的，看起来这是一个改进。谢谢！ - Byte Commander