时间戳

首页 > 解决方案 > 在 sf 对象中按距离聚合特征

r - 在 sf 对象中按距离聚合特征

问题描述

我有两个sf对象：多边形county（注意：这是一个多多边形，即许多县）和点monitor2。

county如下所示。汉字不能正常显示，但问题不大。

Simple feature collection with 6 features and 4 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 113.15 ymin: 20.58265 xmax: 124.5656 ymax: 40.10793
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
                      City                 District                 Province   Code                       geometry
1 <U+53F0><U+6E7E><U+7701> <U+53F0><U+6E7E><U+7701> <U+53F0><U+6E7E><U+7701> 710000 MULTIPOLYGON (((116.7346 20...
2 <U+5317><U+4EAC><U+5E02> <U+671D><U+9633><U+533A> <U+5317><U+4EAC><U+5E02> 110105 MULTIPOLYGON (((116.4834 40...
3 <U+4E0A><U+6D77><U+5E02> <U+666E><U+9640><U+533A> <U+4E0A><U+6D77><U+5E02> 310107 MULTIPOLYGON (((121.3562 31...
4 <U+4E0A><U+6D77><U+5E02> <U+5B9D><U+5C71><U+533A> <U+4E0A><U+6D77><U+5E02> 230506 MULTIPOLYGON (((121.4855 31...
5 <U+5E7F><U+5DDE><U+5E02> <U+767D><U+4E91><U+533A> <U+5E7F><U+4E1C><U+7701> 440111 MULTIPOLYGON (((113.4965 23...
6 <U+798F><U+5DDE><U+5E02> <U+9F13><U+697C><U+533A> <U+798F><U+5EFA><U+7701> 320106 MULTIPOLYGON (((119.2611 26...

monitor2如下所示。

Simple feature collection with 6 features and 5 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 116.17 ymin: 39.8673 xmax: 116.473 ymax: 40.2865
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
# A tibble: 6 x 6
  code  name     city  ref   value          geometry
  <chr> <chr>    <chr> <chr> <dbl>       <POINT [°]>
1 1001A 万寿西宫 北京  N      47.8 (116.366 39.8673)
2 1002A 定陵     北京  Y      45.9  (116.17 40.2865)
3 1003A 东四     北京  N      42.2 (116.434 39.9522)
4 1004A 天坛     北京  N      51.2 (116.434 39.8745)
5 1005A 农展馆   北京  N      46.9 (116.473 39.9716)
6 1006A 官园     北京  N      49.5 (116.361 39.9425)

第一个任务是将value特征加入monitor2到county. 我用st_is_within_distanceand做到了这一点st_join。请参阅下面的代码。我将距离设置为 50 公里。新多边形中的某些县可能具有来自 50 公里缓冲区内多个点的值。

new = st_join(county, monitor2,
            join = st_is_within_distance, dist = 50)

第二个任务来了。我需要从 50 公里缓冲区内的不同点到县质心的距离汇总值。我如何完成这项任务？

欢迎任何意见。

标签： rsf

解决方案

如果没有可重复的数据，很难确切地知道您想要什么，但这里尝试展示您如何做到这一点。

获取样本数据。我们在这里从纬度/经度重新投影到带有米的东西，这样我们就可以进行基于距离的空间操作。我们将使用样本数据中的 3 个县，并使用中间县作为我们要测量距离的主要多边形，并添加散布在三个县的点的随机样本。

library(sf)
nc <- st_read(system.file("shape/nc.shp", package="sf"))
nc <- st_transform(nc, 32119) # NC state plane projection in metres
county = st_cast(nc[2,],"POLYGON")
p1 = st_as_sf(st_sample(nc[1:3, ], 200)) # random points

# Visualize
plot(st_geometry(nc)[1:3])
plot(county, col = "grey80", add = TRUE)

在此处输入图像描述

我们只想关注距目标县一定距离内的点。让我们通过使用st_buffer.

plot(st_buffer(county, dist = 10000), col = NA, border = "red", lty = 3, add = TRUE)

在此处输入图像描述 我们可以对中心县 10000m 范围内的点进行子集化，这与与对象st_is_within_distance相交的效果相同。st_buffer

p1_10 <- p1[st_is_within_distance(county,p1,dist = 10000, sparse = FALSE),]

测量质心和该子集的每个元素之间的距离是直截了当的。然后，我们可以将距离测量值分配为子空间对象中的变量。

p1_10$distance_to_centroid <- as.vector(st_distance(st_centroid(county), p1_10))

这是完全绘制的样子

plot(st_geometry(nc)[1:3])
plot(county, col = "grey80", add = TRUE)
plot(p1, add = TRUE, pch = 19)
plot(st_buffer(county, dist = 10000), col = NA, border = "red", lty = 3, add = TRUE)
plot(st_centroid(county), col = "red", pch = 15, cex = 1, axes = TRUE, add = TRUE)
plot(p1_10["distance_to_centroid"], add = TRUE, pch = 19)

在此处输入图像描述 这是 p1_10 obj 在这里的样子：

> p1_10
Simple feature collection with 78 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 389967.6 ymin: 293489.4 xmax: 448197.1 ymax: 315140.7
CRS:            EPSG:32119
First 10 features:
                           x distance_to_centroid
1  POINT (437228.1 294079.7)           21703.5425
2  POINT (425029.8 305656.7)            5868.4917
3  POINT (425131.4 309137.8)            6665.0253
4  POINT (409851.2 294971.7)           14549.0585
5  POINT (393070.6 303879.7)           26207.5651
6  POINT (436666.3 296282.2)           20070.5879
7  POINT (442623.8 295976.3)           25549.5662
8  POINT (400517.2 307897.4)           18746.6918
9    POINT (418763.7 306728)             724.6165
10 POINT (405001.4 294845.7)           18125.0738

因此，您可以从这里使用您想要的任何方法按距离聚合您的特征。在dplyr中，它非常简单。例如，假设我想以 5 公里的间隔聚合。

library(dplyr)
p1_10 %>% 
  mutate(dist_group = ceiling(distance_to_centroid/5000)) %>% 
  group_by(dist_group) %>% 
  tally() %>% # stop here if you want the result to retain geography
  as_tibble() %>% 
  select(dist_group, n)

# A tibble: 7 x 2
  dist_group     n
       <dbl> <int>
1          1     7
2          2    15
3          3    22
4          4    13
5          5    11
6          6     9
7          7     1

推荐阅读