rgeospatialr-sfcentroid

How to compute the weighted centroid of multiple st_points with the spatialEco package?

Input data

looks like:

Simple feature collection with 4 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -1.752458 ymin: 49.31726 xmax: 1.11638 ymax: 49.87318
Geodetic CRS:  Hu Tzu Shan 1950
# A tibble: 4 × 2
  weight             geometry
   <dbl>          <POINT [°]>
1     12 (-1.752458 49.53935)
2      8  (1.095099 49.41049)
3      3   (1.11638 49.31726)
4     15 (0.8276434 49.87318)

I am using the wt.centroid function from the spatialEco package:

wt.centroid(cities$geometry, cities$weight)

But I get this error:

Error in wt.centroid(cities$geometry, cities$weight) : 
  Projection must be in distance units, not lat/long

Questions

  1. What should I do to use distance units?
  2. Is there another function that can compute centroids without the transformation?

dput() for reproductibility

structure(list(weight = c(12, 8, 3, 15), geometry = structure(list(
    structure(c(-1.75245754632, 49.5393529841), class = c("XY", 
    "POINT", "sfg")), structure(c(1.09509946531, 49.4104887444
    ), class = c("XY", "POINT", "sfg")), structure(c(1.11637971624, 
    49.3172603768), class = c("XY", "POINT", "sfg")), structure(c(0.827643384885, 
    49.8731818649), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = -1.75245754632, 
ymin = 49.3172603768, xmax = 1.11637971624, ymax = 49.8731818649
), class = "bbox"), crs = structure(list(input = "EPSG:4236", 
    wkt = "GEOGCRS[\"Hu Tzu Shan 1950\",\n    DATUM[\"Hu Tzu Shan 1950\",\n        ELLIPSOID[\"International 1924\",6378388,297,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Geodesy.\"],\n        AREA[\"Taiwan, Republic of China - onshore - Taiwan Island, Penghu (Pescadores) Islands.\"],\n        BBOX[21.87,119.25,25.34,122.06]],\n    ID[\"EPSG\",4236]]"), class = "crs"), n_empty = 0L)), row.names = c(NA, 
-4L), sf_column = "geometry", agr = structure(c(weight = NA_integer_), levels = c("constant", 
"aggregate", "identity"), class = "factor"), class = c("sf", 
"tbl_df", "tbl", "data.frame"))

needed packages

library(spatialEco)
library(tidyverse)
library(sf)
Solution

  • There seems to be a problem in wt.centroid if you pass it a tibble such as cities rather than a data frame. It is to do with the differing behaviour of st_drop_geometry (used in the wt.centroid function to extract the column of weights), which returns a vector of weights for a normal df, but returns a single column tibble for a tibble. The latter fails is.numeric, but the former passes. [Update - I have raised this as an issue with the developer of spatialEco].

    So the trick is to convert cities to a df, then transform it to a projected crs, then calculate the centroid and convert back to EPSG:4326...

    cities <- cities %>% 
     as.data.frame() %>%    #convert to df
     st_as_sf() %>%         #reinstate as sf
     st_transform(2154)     #convert to projected crs

cities.centroid <- wt.centroid(cities, p = 'weight', spatial = TRUE) %>% 
     st_transform(4326)     #convert back to lat/long
 
cities.centroid

Simple feature collection with 1 feature and 1 field
Attribute-geometry relationship: 1 constant, 0 aggregate, 0 identity
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 0.09100739 ymin: 49.63301 xmax: 0.09100739 ymax: 49.63301
Geodetic CRS:  WGS 84
  ID                    geometry
1  1 POINT (0.09100739 49.63301)