Как извлечь площадь пола из 3D-данных с помощью пакета sf

Question

У меня есть следующие трехмерные данные из швейцарского Федерального управления топографии Швейцарии по объемам зданий:

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.Любая помощь или подсказка будет высоко оценен!

Answer 1

Если вы внимательно посмотрите на геометрию первого объекта:

> 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 имеет значение.Он состоит из двух многоугольников:

> 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-полигонам.Следовательно, вы получите удвоенную площадь.

Если вы гарантируете, что каждая функция представлена ​​таким образом, то деление на два для получения площади зоны действия является действительным.Если может быть какая-то другая сложность (есть ли пирамиды в Швейцарии? Или, может быть, здания в форме Тоблерона?), Где верхний многоугольник не совпадает с нижним, тогда вам нужно будет извлечь самый нижний многоугольник (если вы хотите иметь место)на земле) или объединение всех многоугольников, если вы хотите, чтобы покрытие области было видно сверху.

Вот пример 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]