如何使用spatialEco包计算多个st_points的加权质心?

huangapple go评论70阅读模式
英文:

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) : 
  投影必须使用距离单位,而不是经纬度单位

问题

  1. 为了使用距离单位,我应该怎么做?
  2. 是否有另一个函数可以在不进行转换的情况下计算质心?

<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[&quot;Hu Tzu Shan 1950&quot;,\n DATUM[&quot;Hu Tzu Shan 1950&quot;,\n ELLIPSOID[&quot;International 1924&quot;,6378388,297,\n LENGTHUNIT[&quot;metre&quot;,1]]],\n PRIMEM[&quot;Greenwich&quot;,0,\n ANGLEUNIT[&quot;degree&quot;,0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[&quot;geodetic latitude (Lat)&quot;,north,\n ORDER[1],\n ANGLEUNIT[&quot;degree&quot;,0.0174532925199433]],\n AXIS[&quot;geodetic longitude (Lon)&quot;,east,\n ORDER[2],\n ANGLEUNIT[&quot;degree&quot;,0.0174532925199433]],\n USAGE[\n SCOPE[&quot;Geodesy.&quot;],\n AREA[&quot;Taiwan, Republic of China - onshore - Taiwan Island, Penghu (Pescadores) Islands.&quot;],\n BBOX[21.87,119.25,25.34,122.06]],\n ID[&quot;EPSG&quot;,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 &lt;- cities %&gt;% 
         as.data.frame() %&gt;%    #convert to df
         st_as_sf() %&gt;%         #reinstate as sf
         st_transform(2154)     #convert to projected crs
    
    cities.centroid &lt;- wt.centroid(cities, p = &#39;weight&#39;, spatial = TRUE) %&gt;% 
         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>



huangapple
  • 本文由 发表于 2023年5月25日 23:51:26
  • 转载请务必保留本文链接:https://go.coder-hub.com/76334169.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定