OSM，rgeos，osmar，面积计算，不相加
OSM, rgeos, osmar, area calculation, does not add up
我正在尝试从 OSM 获取多边形的大小，使用 osmar 下载数据。然而，健全性检查告诉我这些是不对的。
下面是我的意思的一个例子。
(1) 伦敦海德公园周围的地理区域。提取标记为 'park'.
的所有方式和关系
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")
现在，我想知道每个'parks'[图1]的面积（中间那个大的是海德公园）。
spplot(sp_sub_park, c("version"), colorkey = FALSE, col.regions=c('green'))
有两种方法：
1) 在多边形本身中使用插槽 'area'。
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平方米）。
有什么想法吗？
我已将数据转换为使用米作为长度单位的英国国家网格（WGS84 可能使用度数，IDK？）。那么面积是250万平方米，这似乎更合理？
sp_sub_park <- spTransform(sp_sub_park, CRS("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs "))
gArea(sp_sub_park)
# [1] 2550387
这比您估计的约 140 万平方米要高，但是海德公园在您的地图上都是很大的区域（即肯辛顿花园是分开的吗？）。
由于您的 "map" 在 lat/lon 坐标中，要计算面积，您必须将其转换为 "metric" 投影（如@Phil 所建议），或使用球形几何（例如，在包 geosphere 中实现）。
幸运的是，如果对象在地理坐标中，sf 包的函数 st_area 使用 geosphere 计算多边形的面积。因此，你可以简单地做：
sf_sub_park <- st_as_sf(sp_sub_park)
areas <- st_area(sf_sub_park)
sum(areas)
，给予：
2551269 m^2
，这与@Phil 的结果非常接近。
HTH
我正在尝试从 OSM 获取多边形的大小，使用 osmar 下载数据。然而，健全性检查告诉我这些是不对的。
下面是我的意思的一个例子。
(1) 伦敦海德公园周围的地理区域。提取标记为 'park'.
的所有方式和关系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")
现在，我想知道每个'parks'[图1]的面积（中间那个大的是海德公园）。
spplot(sp_sub_park, c("version"), colorkey = FALSE, col.regions=c('green'))
有两种方法：
1) 在多边形本身中使用插槽 'area'。
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平方米）。
有什么想法吗？
我已将数据转换为使用米作为长度单位的英国国家网格（WGS84 可能使用度数，IDK？）。那么面积是250万平方米，这似乎更合理？
sp_sub_park <- spTransform(sp_sub_park, CRS("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs "))
gArea(sp_sub_park)
# [1] 2550387
这比您估计的约 140 万平方米要高，但是海德公园在您的地图上都是很大的区域（即肯辛顿花园是分开的吗？）。
由于您的 "map" 在 lat/lon 坐标中，要计算面积，您必须将其转换为 "metric" 投影（如@Phil 所建议），或使用球形几何（例如，在包 geosphere 中实现）。
幸运的是，如果对象在地理坐标中，sf 包的函数 st_area 使用 geosphere 计算多边形的面积。因此，你可以简单地做：
sf_sub_park <- st_as_sf(sp_sub_park)
areas <- st_area(sf_sub_park)
sum(areas)
，给予：
2551269 m^2
，这与@Phil 的结果非常接近。
HTH