英文:
How to compute the weighted centroid of multiple st_points with the spatialEco package?
问题
以下是翻译好的内容:
### 输入数据
看起来像这样:
```R
具有4个要素和1个字段的简单要素集
几何类型:点
维度:XY
边界框:xmin:-1.752458 ymin:49.31726 xmax:1.11638 ymax:49.87318
地理坐标系:胡子山1950
# 一个小表格: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)
我正在使用spatialEco
包中的wt.centroid
函数:
wt.centroid(cities$geometry, cities$weight)
但是我遇到了这个错误:
Error in wt.centroid(cities$geometry, cities$weight) :
投影必须使用距离单位,而不是经纬度单位
问题
- 为了使用距离单位,我应该怎么做?
- 是否有另一个函数可以在不进行转换的情况下计算质心?
<details>
<summary>英文:</summary>
### 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)
</details>
# 答案1
**得分**: 1
在将 `cities` 转换为数据框之后,将其转换为投影坐标系,然后计算质心并转换回 EPSG:4326...
```R
cities <- cities %>%
as.data.frame() %>%
st_as_sf() %>%
st_transform(2154) # 转换为投影坐标系
cities.centroid <- wt.centroid(cities, p = 'weight', spatial = TRUE) %>%
st_transform(4326) # 转换回经纬度坐标系
cities.centroid
简单特征集,包含1个特征和1个字段
属性-几何关系:1个常数,0个聚合,0个身份
几何类型:点
维度:XY
边界框:xmin: 0.09100739 ymin: 49.63301 xmax: 0.09100739 ymax: 49.63301
大地测量坐标系:WGS 84
ID geometry
1 1 POINT (0.09100739 49.63301)
<details>
<summary>英文:</summary>
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)
</details>
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论