построчная работа с sf - PullRequest
Новая
       25

построчная работа с sf

0 голосов
/ 11 ноября 2019

Я пытаюсь получить список соседних округов в другом столбце. Я использую sf и границы округа от «USAboundaries». Я написал следующий код, но я чувствую, что должен быть вычислительно более эффективный способ. Любая идея?

Заранее спасибо. 

library(USAboundaries)
library(sf)

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$neighbor[i] <- county$geoid[unlist(st_is_within_distance(county[i,], county, 70000, sparse = TRUE ))]

}

1 Ответ

2 голосов
/ 11 ноября 2019

Более эффективным подходом может быть:

  1. для буферизации слоя county по соседнему выбранному вами порогу, например 70000 м,
  2. , использование st_intersects на«Буферизованный» слой по отношению к исходному слою, чтобы выяснить, какие графства пересекаются в каждом «буферизованном» округе,
  3. и вставить список идентификаторов, чтобы они поместились в столбец 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

Чтобы убедиться, что результат имеет смысл, следующий код также вычисляет соседа count: 

county = st_transform(county, 2163)
x = st_buffer(county, 70000)
int = st_intersects(x, county)
county$n = sapply(int, length)

А вот график соседей: 

plot(county[, "n"])

...