OSM，rgeos，osmar，面积计算，不累加

Question

我正在尝试使用osmar从OSM获取多边形的大小，以下载数据。但是，理智检查告诉我这是不对的。

以下是我的意思的示例。

（1）伦敦海德公园附近的地理区域。提取标记为“公园”的所有方式和关系。

devtools::install_github('osmdatar/osmdata')
library(osmdata)
library(osmar)
library(sp)
library(sf)
library(rgeos)

osmO <- get_osm(center_bbox(-0.167919, 51.5072682, 2000, 2000))

ids_relations <- osmO$relations$tags[osmO$relations$tags$v=="park","id"]
ids_ways <- osmO$ways$tags[osmO$ways$tags$v=="park","id"]
ids_sub <- find_down(osmO, way(c(ids_relations, ids_ways)))
sp_sub_park <- as_sp(subset(osmO, ids = ids_sub), "polygons")

现在，我想知道每个“公园”的面积（图1）（中间最大的是海德公园）。

spplot(sp_sub_park, c("version"), colorkey = FALSE, col.regions=c('green'))

有两种方法：

1）在多边形本身中使用插槽“区域”。

a1 <- sapply(sp_sub_park@polygons, function(x) x@area)

2）用指定的投影计算面积。

bg_poly_t <- spTransform(sp_sub_park, CRS("+proj=longlat +datum=WGS84"))
a2 <- rgeos::gArea(bg_poly_t, byid=TRUE)

这两个给我相同的结果[图2]（请注意，两个最大的区域是海德公园，被一条道路分成两部分）

plot(a1*1000000, a2*1000000)

但是，大小不是我期望的。该面积以平方公里（绘制的平方米）为单位返回。据称，海德公园的两个部分总计约300平方米，相当于一个大公寓的大小，但不是公园（海德公园〜1.420.000平方米）。

有任何想法吗？

Answer 1

由于您的“地图”采用纬度/经度坐标，因此要计算区域，您必须将其转换为“公制”投影（如@Phil建议），或使用球面几何（例如，在包中实现geosphere）。

st_area幸运的是，该包的功能sf可以计算多边形的面积，geosphere以防对象位于地理坐标中。因此，您可以简单地执行以下操作：

sf_sub_park <- st_as_sf(sp_sub_park)
areas <- st_area(sf_sub_park)
sum(areas)

，给予：

2551269 米^2

，这与@Phil 结果非常接近。

华泰