更有效的方法可能是:
- 通过您选择的相邻阈值缓冲
county
图层,例如 70000 m,
- 在原始图层的“缓冲”层上使用
st_intersects
,找出每个“缓冲”县相交的县,
- 并粘贴 ID 列表,使其适合
character
列
这是一个可重现的示例:
library(USAboundaries)
library(sf)
library(dplyr)
library(stringr)
county <- us_counties() %>%
filter(str_detect(statefp,"02|72|78|15|66")==FALSE) %>%
mutate(geoid = as.integer(geoid))
# for(i in 1:nrow(county)){
# county$neighbor1[i] <- county$geoid[unlist(st_is_within_distance(county[i,], county, 70000, sparse = TRUE ))]
#
# }
county = st_transform(county, 2163) # Transform to EPSG:2163
x = st_buffer(county, 70000) # Calculate buffered county layer
int = st_intersects(x, county) # Find intersections
int = lapply(int, function(x) county$geoid[x]) # Translate intersection indices to IDs (i.e., geoid)
int = sapply(int, paste, collapse = "|") # Paste into single string
county$neighbors = int # Bind to the county layer
结果打印在下面。新neighbors
列最后包含邻居geoid
s:
head(county)
## Simple feature collection with 6 features and 13 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 95466.62 ymin: -1693172 xmax: 1475761 ymax: 24168.64
## epsg (SRID): 2163
## proj4string: +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
## statefp countyfp countyns affgeoid geoid name lsad
## 1 39 131 01074078 0500000US39131 39131 Pike 06
## 2 46 003 01266983 0500000US46003 46003 Aurora 06
## 3 55 035 01581077 0500000US55035 55035 Eau Claire 06
## 4 48 259 01383915 0500000US48259 48259 Kendall 06
## 5 40 015 01101795 0500000US40015 40015 Caddo 06
## 6 19 093 00465235 0500000US19093 19093 Ida 06
## aland awater state_name state_abbr jurisdiction_type
## 1 1140324458 9567612 Ohio OH state
## 2 1834813753 11201379 South Dakota SD state
## 3 1652211310 18848512 Wisconsin WI state
## 4 1715747531 1496797 Texas TX state
## 5 3310745124 30820525 Oklahoma OK state
## 6 1117599859 1406461 Iowa IA state
## geometry
## 1 MULTIPOLYGON (((1424415 -50...
## 2 MULTIPOLYGON (((95466.62 -1...
## 3 MULTIPOLYGON (((656691.6 17...
## 4 MULTIPOLYGON (((104718.2 -1...
## 5 MULTIPOLYGON (((124977.8 -1...
## 6 MULTIPOLYGON (((348649.9 -2...
## neighbors
## 1 39131|39047|39027|39097|21043|39129|39045|54053|39023|39049|54011|39165|21135|21069|39025|39057|39079|39087|39145|21019|39141|21205|39073|21161|39015|21023|54099|39001|39009|39163|21089|39127|39071|39053|39105
## 2 46003|46043|46067|46077|46015|31015|46061|46087|46073|31103|31089|46135|46009|46035|46005|46111|46059|46069|46053|46085|46097|31107|46023|46065|46017|46123
## 3 55035|55107|55141|55017|55005|27169|27109|55119|55081|55121|55033|55011|55053|55019|55093|55099|55057|55063|55095|55091|27157|55073|27049|55109
## 4 48259|48163|48171|48385|48463|48029|48055|48453|48265|48187|48325|48493|48013|48019|48319|48209|48267|48299|48053|48031|48091
## 5 40015|40149|40011|40109|40049|40141|40039|40073|40043|40031|40067|40009|40083|40051|40087|40017|40027|40093|40033|40075|40137|40065|40129|40055|40019
## 6 19093|19147|19165|19009|19193|19141|19077|19027|19041|31043|19133|19035|19021|19167|19151|31173|19085|19149|19047|19025|31177|31021|46127|19073|19161
为了确保结果有意义,以下代码还计算了邻居计数:
county = st_transform(county, 2163)
x = st_buffer(county, 70000)
int = st_intersects(x, county)
county$n = sapply(int, length)
这是一个邻居计数图:
plot(county[, "n"])