CODEWITHSUNDEEP
rgisr-sfcentroid
Unstack
Unstack

Reputation: 561

How do I get gCentroid to work at the pole?

I am struggling with gCentroid, because it doesn't seem -- to me -- to give the 'right' answer near a pole of the Earth.

For instance:

library(rgeos)

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

does not give me the North Pole, it gives:

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

How do I get gCentroid to work on the surface of the Earth?

Upvotes: 1

Views: 58

Answers (1)

Jindra Lacko
Jindra Lacko

Reputation: 8699

The GEOS library is limited to planar geometry operations; this can bring issues in edge cases / the poles being a notorious example.

For the centroid via GEOS to work as intended you need to transform your coordinates from WGS84 to a coordinate reference system appropriate to polar regions; for Arctic regions I suggest 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

Also note that your workflow - while not wrong in principle - is a bit dated, and the {rgeos} package is approaching its end of life.

It may be good time to give a strong consideration to {sf} package, which is newer, actively developed and can, via interface to s2 library from Google, handle spherical geometry operations.

For an example of {sf} based workflow consider this code; the result (centroid = North Pole) is equivalent to the sp / rgeos one.

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)

Upvotes: 1

Related Questions