如何让 gCentroid 在极点工作?
How do I get gCentroid to work at the pole?
我正在与 gCentroid 作斗争,因为在我看来它似乎无法在地球的一极附近给出 'right' 答案。
例如:
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 在地球表面工作?
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}
包的好时机,它更新、积极开发并且可以通过 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)
我正在与 gCentroid 作斗争,因为在我看来它似乎无法在地球的一极附近给出 'right' 答案。
例如:
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 在地球表面工作?
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}
包的好时机,它更新、积极开发并且可以通过 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)