如何使用 sf 包从 3D 数据中提取建筑面积

How to extract floor area from 3D data using the sf package

我有以下来自瑞士联邦地形局 swisstopo 的建筑体积 3D 数据:

https://shop.swisstopo.admin.ch/en/products/landscape/build3D

要下载数据,您需要点击阅读更多,向下滚动并点击示例数据。

这只是一个示例。请注意,我使用 GDB 是因为我的真实数据集采用这种格式。 然后我把数据读入R:

library(sf)

data <- st_read("swissbuildings3dlv03/GDB/swissBUILDINGS3d_10.gdb")

在下一步中,我删除三维并计算面积:

data %>% 
 st_zm() %>% 
 st_area()

然而,与这些建筑的官方数据相比,面积总是比预期的大两倍。我想,我将尺寸从 3D 缩小到 2D 的方式有问题。任何帮助或提示将不胜感激!

如果您仔细观察第一个要素的几何形状:

> st_as_text(st_geometry(map[1,]))
GEOMETRYCOLLECTION Z ( TIN Z .....
[...]
MULTIPOLYGON Z (((601608.1 197988.1 525.57, 601604 197974.8 525.57, 601593.7 197977.9 525.57, 601597.8 197991.3 525.57, 601608.1 197988.1 525.57)), ((601608.1 197988.1 517.9725, 601597.8 197991.3 517.9725, 601593.7 197977.9 517.9725, 601604 197974.8 517.9725, 601608.1 197988.1 517.9725))))"

MULTIPOLYGON 位是相关的。它由两个 POLYGONS 组成:

> st_as_text(st_cast(st_collection_extract(st_geometry(map)[1]),"POLYGON"))
[1] "POLYGON Z ((601608.1 197988.1 525.57,
                 601604 197974.8 525.57,
                 601593.7 197977.9 525.57,
                 601597.8 197991.3 525.57,
                 601608.1 197988.1 525.57))"          
[2] "POLYGON Z ((601608.1 197988.1 517.9725,
                 601597.8 197991.3 517.9725,
                 601593.7 197977.9 517.9725,
                 601604 197974.8 517.9725,
                 601608.1 197988.1 517.9725))"

这看起来像两个相同的多边形在不同的 Z 高度。通过 st_zm 过滤将产生两个相同的 2D 多边形。因此,您将获得双倍的面积。

如果你保证每个特征都是这样表示的,那么除以二得到足迹面积是有效的。如果可能存在其他复杂性(瑞士有金字塔吗?或者形状像 Toblerone 的建筑物?)顶部多边形与底部多边形不同,那么您需要提取最低的多边形(如果您想要足迹在地面上)或所有多边形的联合,如果你想要从上面看到的区域覆盖。

下面是第一个要素 st_union 通过合并顶部和底部多边形将面积减半的示例:

> st_area(st_union(st_collection_extract(st_geometry(map)[1])))
150.5087 [m^2]

与原版相比:

> st_area(st_collection_extract(st_geometry(map)[1]))
301.0174 [m^2]