Несоответствия геометрии после вычисления пересечения

Question

Я новичок в геопространственной стороне вещей, и я пытаюсь разделить шейп-файл на части сетки, а затем вычислить пересечение по сетке. Тем не менее, я сталкиваюсь с проблемами согласованности с выходными данными. Вот минимальный пример:

x <- st_read("https://data.sfgov.org/api/geospatial/rarb-5ahf?method=export&format=GeoJSON")
grids <- grids <- sf::st_make_grid(x, n = 40)
int <- st_intersection(x, grids)

Я (как и ожидалось) получаю следующее сообщение:

, хотя координаты долгота / широта, st_intersection предполагает, что они плоские

Затем выполнение st_area дает следующий вывод:

> sum(st_area(x))
600592318 m^2
> sum(st_area(int))
600594187 m^2

Полученная площадь пересечения. В зависимости от того, где я нахожусь в мире, он иногда теряет область.

Сценарий реального мира, с которым я использую это, включает шейп-файл и растр, содержащий весь мир, спроецированный в longlat. Мне также нужен мой вывод в longlat.

Answer 1

Комментарий @Spacedman представляется правильным. Вы можете проверить, уплотнив полигоны, используя функцию densify в пакете smoothr:

library(sf)
library(smoothr)

x <- st_read("https://data.sfgov.org/api/geospatial/rarb-5ahf?method=export&format=GeoJSON")
grids <- grids <- sf::st_make_grid(x, n = 40)
int <- st_intersection(x, grids)

#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant
#> throughout all geometries

sum(st_area(x))
#> 600592318 [m^2]
sum(st_area(int))
#> 600594187 [m^2]

grids_dense <- smoothr::densify(grids, n = 60)
x_dense     <- smoothr::densify(x, n = 60)
int_dense   <- st_intersection(x_dense, grids_dense)

#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant
#> throughout all geometries

sum(st_area(x_dense))
#> 600594404 [m^2]
sum(st_area(int_dense))
#> 600594404 [m^2]

Answer 2

Спасибо всем, кто ответил. Решение densify сработало, но с моей стороны было слишком дорого в вычислительном отношении. Я попробовал другие библиотеки и нашел решение с rgeos.

library(rgeos)
library(sf)

x <- sf::st_read("https://data.sfgov.org/api/geospatial/rarb-5ahf?method=export&format=GeoJSON")
grids <- grids <- sf::st_make_grid(x, n = 40)
int <- sf::st_intersection(x, grids)

> sum(st_area(x))
600592318 m^2
> sum(st_area(int))
600594187 m^2

> rgeos::gArea(sf::as_Spatial(x))
[1] 0.06140787
> rgeos::gArea(sf::as_Spatial(int))
[1] 0.06140787

. Это работает для меня, так как меня больше всего волнует пропорция площади.