在 R 中找到周长与周围邻居的比率

Question

有一个形状文件

我将 R 中的数据读取为名为“SD”的 sf 对象，并包含“GEOID”和“geometry”列以及其他列。我在名为“POP”的数据中创建了一个新列

SD$POP <- sample(2000:8000, 35)

首先，我想测量每个“GEOID”的面积和周长。

其次，我需要找出每个 GEOID 的周长与其周围的 GEOID 共享的比例。例如：假设 GEOID 09 被 GEOID 06、08、10、13 包围。我想找出 GEOID 09 的周长与 06、08、10 和 13 分别共享的百分比。然后创建新列“POP_PERCENT”和“ADJ_POP”这将是

ADJ-POP = SUM OF POP for GEOID 06, 08, 10, and 13' .

POP_PERCENT =  Percent share with 06 + Percent share with 08 + Percent share with 10 +
Percent share with 13

对于冗长的问题，我深表歉意。如果有任何帮助，我将不胜感激。

Answer 1

我认为通常您需要的是 st_intersection 加上一些数据操作。数据集没有 GEOID 列，我不完全清楚 POP_PERCENT 的预期输出是什么（相邻人口按周长重叠百分比加权？）。所以这不是一个完全完整的答案。

这是使用 sf 和其他一些 tidyverse 函数的开始。这将计算每个多边形的面积、周长和交点。我从 rownames 创建了一个 GEOID 列。希望这足以让您朝着正确的方向前进。

先算出每个的面积和周长

library(tidyverse)
library(sf)

SD_ap <- SD %>%
  st_transform(5070) %>% #convert from lat/long to projected coord system
  mutate(area = st_area(.),
         length_perim = lwgeom::st_perimeter(.))

#-----
  GEOID  POP                       geometry              area         perim
1     1 6582 MULTIPOLYGON (((860088.4 12...   357945281 [m^2]  139932.2 [m]
2     2 3199 MULTIPOLYGON (((843395.6 11...   137754181 [m^2]  100879.4 [m]

然后计算重叠周长。将多边形转换为线，然后运行 st_intersection。这 returns 141 条线段不是 35 个多边形，而是以原始重叠多边形 (origins)

区分

SD_overlap <- SD_ap %>%
  st_cast('MULTILINESTRING') %>%
  st_intersection() %>%
  mutate(length_intersection = st_length(.)) %>%
  filter(length_intersection > units::as_units(0,'m'))

#----
     GEOID  POP              area  length_perim n.overlaps origins length_intersection
1        1 6559   357945281 [m^2]  139932.2 [m]          2    1, 2        6761.611 [m]
2        2 4781   137754181 [m^2]  100879.4 [m]          2    2, 3       24201.037 [m]
2.2      2 4781   137754181 [m^2]  100879.4 [m]          2    2, 4       31628.198 [m]
3        3 3424   295676500 [m^2]  212930.4 [m]          2    3, 4       22967.725 [m]

绘制前 4 GEOID 秒以检查是否按预期运行

ggplot(data = filter(SD, GEOID %in% 1:4),
       aes(fill = as.factor(GEOID),
           label = GEOID)) +
  geom_sf(alpha = .2) +
  geom_sf(data = filter(SD_overlap, map_lgl(origins, ~all(.x %in% 1:4))), 
          aes(color = map_chr(origins, ~paste(.x, collapse = ','))),
          fill = NA,
          lwd = 2,
          show.legend = 'line') +
  geom_sf_label(fill = "white",  # override the fill from aes()
                fun.geometry = sf::st_centroid) +
  guides(fill = 'none')+
  labs(color='Overlap')

要回答哪些种群相邻的问题，您可以尝试这样的方法。

sum_adj_pops <- function(GEOID_i){
  
  ID_adj <- SD_overlap %>%
    st_drop_geometry() %>%
    filter(map_lgl(origins, ~GEOID_i %in% .x)) %>%
    unnest_longer(origins) %>%
    filter(origins != GEOID_i) %>%
    pull(origins)
  
  SD %>%
    st_drop_geometry() %>%
    filter(GEOID %in% ID_adj) %>%
    summarize(POP_adj = sum(POP)) %>%
    pull(POP_adj)
}

SD %>%
  mutate(POP_adj = map_dbl(GEOID, sum_adj_pops))

#----
   GEOID  POP                       geometry POP_adj
1      1 6559 MULTIPOLYGON (((-86.64623 3...   24024
2      2 4781 MULTIPOLYGON (((-86.861 33....   27132
3      3 3424 MULTIPOLYGON (((-86.97402 3...   21907
4      4 4230 MULTIPOLYGON (((-86.74407 3...   22259
5      5 7054 MULTIPOLYGON (((-87.78275 3...   22021

数据

url_shp <- 'https://github.com/ahadzaman1002/example_data/raw/main/fe_2007_01_sldu00.zip'

tmp_zip <- tempfile(fileext = '.zip')
tmp_dir <- tempdir()
download.file(url = url_shp, tmp_zip)
unzip(tmp_zip, exdir = tmp_dir)

path_shp <- list.files(tmp_dir, '.shp$', full.names = TRUE)

set.seed(415) #for reproducibility with sample()

SD <- st_read(path_shp) %>%
  mutate(POP = sample(2000:8000, 35)) %>%
  rowid_to_column('GEOID')

在 R 中找到周长与周围邻居的比率

Finding ratio of perimeter with surrounding neighbors in R

gis

r

spatial

sf