首页 > 解决方案 > 如何让 gCentroid 在极点工作?

问题描述

我在gCentroid上苦苦挣扎,因为在我看来,在地球极点附近给出“正确”的答案似乎并没有。

例如:

library(rgeos)

gCentroid(SpatialPoints(coords=data.frame(longitude=c(-135,-45,45,135),latitute=c(80,80,80,80)),proj4string = CRS('EPSG:4326')))

没有给我北极,它给了我:

> SpatialPoints: 
>   x  y 
> 1 0 80 
> Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84 +no_defs

如何让gCentroid在地球表面工作?

标签: rgissfcentroid

解决方案


GEOS 库仅限于平面几何操作;这可能会在极端情况下带来问题/极点是一个臭名昭著的例子。

为了使通过 GEOS 的质心按预期工作,您需要将坐标从 WGS84 转换为适合极地的坐标参考系统;对于北极地区,我建议EPSG:3995

library(sp)
library(dplyr)
library(rgeos)

points_sp <- SpatialPoints(coords=data.frame(longitude=c(-135,-45,45,135),latitute=c(80,80,80,80)),proj4string = CRS('EPSG:4326'))

points_updated <- points_sp %>% 
  spTransform(CRS("EPSG:3995")) # a projected CRS apropriate for Arctic regions
  
centroid <- gCentroid(points_updated) %>% 
  spTransform(CRS("EPSG:4326")) # back to safety of WGS84!

centroid # looks better now...
# SpatialPoints:
#  x  y
# 1 0 90
# Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84 +no_defs 

另请注意,您的工作流程 - 虽然原则上没有错误 - 有点过时,并且{rgeos}包即将结束。

现在可能是认真考虑{sf}package 的好时机,它是较新的、积极开发的,并且可以通过与 Google 的 s2 库的接口来处理球面几何操作。

对于{sf}基于工作流的示例,请考虑此代码;结果(质心 = 北极)相当于 sp / rgeos 之一。

library(sf)

points_sf <- points_sp %>% # declared earlier
  st_as_sf()


centroid_sf <- points_sf %>% 
  st_union() %>%  # unite featrues / from 4 points >> 1 multipoint
  st_centroid()

centroid_sf # the North Pole in a slightly different (sf vs sp) format
# Geometry set for 1 feature 
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 0 ymin: 90 xmax: 0 ymax: 90
# Geodetic CRS:  WGS 84 (with axis order normalized for visualization)
# POINT (0 90)

推荐阅读