Создание объекта inla.graph с правильными индексными метками из объекта sf linestring

Question

Создание объекта inla.graph с правильными индексными метками из объекта sf linestring

У меня есть шейп-файл, содержащий linestringsкоторый описывает связь между городами в Бразилии. Я хотел бы преобразовать эти соединения в объект соседства с кодом города, установленным в качестве имени строки, чтобы сделать его совместимым с моим фреймом данных:

      > head(regic_link)
Simple feature collection with 6 features and 6 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: -45.7931 ymin: -21.1472 xmax: -41.3903 ymax: -15.9032
proj4string:   +proj=longlat +ellps=GRS80 +no_defs 
     id origin_code            nome_ori dest_code        nome_dest   dist_km                       geometry
1 14016     3123304      Dores do Turvo   3165701  Senador Firmino 11.732462 LINESTRING (-43.1884 -20.97...
2 14117     3124708   Estrela do Indaiá   3166600 Serra da Saudade  9.080594 LINESTRING (-45.7878 -19.52...
3 15205     3138658              Lontra   3135357         Japonvar 11.104687 LINESTRING (-44.3029 -15.90...
4 17147     3163300  São José do Divino   3144904      Nova Módica 13.066889 LINESTRING (-41.3903 -18.48...
5 17151     3163409 São José do Goiabal   3121803         Dionísio 12.078370 LINESTRING (-42.7077 -19.92...
6 12463     3102100       Alto Rio Doce   3121506 Desterro do Melo 17.982702 LINESTRING (-43.412 -21.025...

Таким образом, имя строки будет установлено как, а сосед - в списке, и наоборот (в конечном итоге это будет изменено на индекс, который я создаю, но это упрощает проверку). По сути, мне нужен linestring эквивалент для следующего кода, используемого для многоугольников:

      nb.orig <- poly2nb(as_Spatial(shp), row.names = shp$index)
names(nb.orig) <- attr(nb.orig, "region.id")

nb2INLA("output/nb_orig.graph",  nb.orig)

( shp - шейп-файл, состоящий из многоугольников, а index переменная присутствует как в шейп-файле, так и во фрейме данных).

До сих пор я использовал sfnetworks а также igraph пакеты для создания объекта окрестности, но не смогли присоединить значения индекса к именам строк:

      regic_net <- as_sfnetwork(regic_link, directed = F) 

net_adj <- as_adjacency_matrix(regic_net, names = T)
nb_net <- mat2listw(net_adj)$neighbours

Функция as_sfnetwork по умолчанию присваивает номер каждому из узлов в сети (в граничных данных это переменные 'from' и 'to', которые нельзя изменить с помощью функции mutate):

      > regic_net
# A sfnetwork with 827 nodes and 3786 edges
#
# CRS:  NA 
#
# An undirected simple graph with 1 component with spatially explicit edges
#
# Edge Data:     3,786 x 7 (active)
# Geometry type: LINESTRING
# Dimension:     XY
# Bounding box:  xmin: -50.6938 ymin: -22.855 xmax: -39.9496 ymax: -14.2696
   from    to origin_code nome_ori           dest_code nome_dest                                                                                    geometry
  <int> <int>       <dbl> <chr>                  <dbl> <chr>                                                                                <LINESTRING [°]>
1     1     2     3123304 Dores do Turvo       3165701 Senador Firmino    (-43.1884 -20.975, -43.18106 -20.97017, -43.17373 -20.96534, -43.16639 -20.9605, …
2     3     4     3163409 São José do Goiab…   3121803 Dionísio           (-42.7077 -19.9255, -42.71344 -19.91851, -42.71917 -19.91152, -42.72491 -19.90453…
3     5     6     3167707 Sobrália             3162609 São João do Orien… (-42.0972 -19.2363, -42.1016 -19.2437, -42.106 -19.2511, -42.11039 -19.2585, -42.…
4     5     7     3167707 Sobrália             3168408 Tarumirim          (-42.0972 -19.2363, -42.08899 -19.24038, -42.08079 -19.24447, -42.07258 -19.24855…
5     8     9     3115474 Catuti               3141009 Mato Verde         (-42.9607 -15.3612, -42.95243 -15.36428, -42.94415 -15.36735, -42.93588 -15.37043…
6    10    11     3130606 Inconfidentes        3108305 Borda da Mata      (-46.328 -22.3174, -46.31896 -22.31505, -46.30993 -22.3127, -46.30089 -22.31034, …
# … with 3,780 more rows
#
# Node Data:     827 x 1
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -50.6938 ymin: -22.855 xmax: -39.9496 ymax: -14.2696
             geometry
          <POINT [°]>
1  (-43.1884 -20.975)
2  (-43.1004 -20.917)
3 (-42.7077 -19.9255)
# … with 824 more rows

Эти значения затем используются как имена строк на следующем шаге в nb объект, созданный с использованием as_adjacency_matrix а также mat2listw. Кто-нибудь знает, как прикрепить мой собственный индекс или даже извлечь индекс, который был создан, чтобы я мог согласовать его между этим и фреймом данных? Я попытался создать таблицу преобразования, взяв / origin_code а также to/ dest_code но они не совпадают, from value является самым низким из двух индексов и поэтому не всегда является источником моих данных.

Также, чтобы усложнить задачу, на 57 узлов больше, чем ожидалось, и я не могу понять, почему это произошло. Я проверил дубликаты, они имеют одинаковые координаты и код / имя в данных, но обрабатываются отдельно с другим номером индекса.

Таким образом, я хотел бы преобразовать объект линии в объект соседства, который можно использовать в модели INLA, где города, соединенные линиями, считаются соседями. Мне нужно прикрепить индекс к этим городам, чтобы объект окрестности установил индекс в качестве имени строки, и его можно было объединить с данными в модели. Надеюсь, это имеет смысл, дайте мне знать, если нет!

2

r igraph r-sf sfnetwork r-inla

Источник

24 авг '21 в 12:46

1 ответ

Другие вопросы по тегам r igraph r-sf sfnetwork r-inla

user10280411 25 авг '21 в 16:12 2021-08-25 16:12 · Answer 1 · 2021-08-25 16:12

Я представлю решение на упрощенном примере всего с 5 LINESTRING(s). Если следующий код соответствует желаемому решению, его будет легко перевести на реальные случаи.

Сначала определите некоторые функции и загрузите соответствующие пакеты:

      '%!in%' <- function(x,y) !('%in%'(x,y))
Sys.setenv(`_R_S3_METHOD_REGISTRATION_NOTE_OVERWRITES_` = "false")
suppressPackageStartupMessages({
  library(sf)
  library(lwgeom)
  library(tidygraph)
  library(sfnetworks)
})

Загрузите данные и выберите только соответствующие столбцы + только пять строк

      regic_link <- st_read(
  dsn = "city_link/REGIC2018_Ligacoes_entre_Cidades.shp",
  query = "SELECT cod_ori, cod_dest FROM REGIC2018_Ligacoes_entre_Cidades LIMIT 5"
)
#> Reading query `SELECT cod_ori, cod_dest FROM REGIC2018_Ligacoes_entre_Cidades LIMIT 5' from data source `C:\Users\Utente\Desktop\city_link\REGIC2018_Ligacoes_entre_Cidades.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 5 features and 2 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -63.0338 ymin: -13.4945 xmax: -60.1348 ymax: -9.912
#> Geodetic CRS:  SIRGAS 2000

Отфильтровать повторяющиеся соединения

      regic_link <- regic_link %>%
  filter(paste0(cod_ori, cod_dest) %!in% paste0(cod_dest, cod_ori))

Создайте объект sfnetwork

      regic_sfn <- as_sfnetwork(regic_link, directed = FALSE)

Извлеките узлы (т.е. уникальные граничные точки для этих линий)

      regic_nodes <- st_as_sf(regic_sfn, "nodes")

Каждый узел должен соответствовать граничной точке LINESTRINGS, определенной в объекте. regic_link. Как вы заметили, мы не можем просто использовать столбцы из / в, поскольку они не сохраняют порядок отправления и назначения. Следовательно, мы должны вручную сопоставить узлы с начальными и конечными точками (предполагая, что каждый раз они имеют одинаковые координаты, тогда они также имеют одинаковый идентификатор). Более того, каждый узел может соответствовать нескольким (идентичным) граничным точкам, но нам просто нужно взять первое совпадение (т.е. [1] обозначения ниже):

      start <- lapply(
  X = st_equals(regic_nodes, st_startpoint(regic_link)),
  FUN = function(x) regic_link$cod_ori[x[1]]
)

Проверьте вывод:

      start
#> [[1]]
#> [1] "1100015"
#> 
#> [[2]]
#> [1] NA
#> 
#> [[3]]
#> [1] NA
#> 
#> [[4]]
#> [1] NA
#> 
#> [[5]]
#> [1] "1100023"
#> 
#> [[6]]
#> [1] NA
#> 
#> [[7]]
#> [1] "1100031"
#> 
#> [[8]]
#> [1] NA

Это означает, что первый узел соответствует ID 1100015седьмой узел соответствует ID 1100031и так далее. С другой стороны, другие узлы не соответствуют какой-либо начальной точке. Повторите для конечных точек:

      end <- lapply(
  X = st_equals(regic_nodes, st_endpoint(regic_link)),
  FUN = function(i) regic_link$cod_dest[i[1]]
)

Проверьте вывод, который должен быть противоположен start и с теми же интерпретациями:

      end
#> [[1]]
#> [1] NA
#> 
#> [[2]]
#> [1] "1100049"
#> 
#> [[3]]
#> [1] "1100189"
#> 
#> [[4]]
#> [1] "1100296"
#> 
#> [[5]]
#> [1] NA
#> 
#> [[6]]
#> [1] "1100114"
#> 
#> [[7]]
#> [1] NA
#> 
#> [[8]]
#> [1] "1100304"

Объедините два объекта. dplyr::coalescence используется для получения идентификаторов, отличных от NA

      idxs <- mapply(dplyr::coalesce, start, end)

Добавьте новые идентификаторы в таблицу узлов:

      regic_sfn <- regic_sfn %N>%
  mutate(name = idxs)

Оцените матрицу смежности между городами

      regic_adj <- igraph::as_adjacency_matrix(regic_sfn, names = TRUE)
regic_adj
#> 8 x 8 sparse Matrix of class "dgCMatrix"
#>         1100015 1100049 1100189 1100296 1100023 1100114 1100031 1100304
#> 1100015       .       1       1       1       .       .       .       .
#> 1100049       1       .       .       .       .       .       .       .
#> 1100189       1       .       .       .       .       .       .       .
#> 1100296       1       .       .       .       .       .       .       .
#> 1100023       .       .       .       .       .       1       .       .
#> 1100114       .       .       .       .       1       .       .       .
#> 1100031       .       .       .       .       .       .       .       1
#> 1100304       .       .       .       .       .       .       1       .

Если мы проверим первую строку, мы заметим, что город соответствует правильным столбцам:

      regic_link %>% filter(cod_ori == 1100023)
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -63.0338 ymin: -10.4323 xmax: -62.4721 ymax: -9.912
#> Geodetic CRS:  SIRGAS 2000
#>   cod_ori cod_dest                       geometry
#> 1 1100023  1100114 LINESTRING (-63.0338 -9.912...

С другой стороны

      regic_link %>% filter(cod_ori == 1100049)
#> Simple feature collection with 0 features and 2 fields
#> Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
#> Geodetic CRS:  SIRGAS 2000
#> [1] cod_ori  cod_dest geometry
#> <0 rows> (or 0-length row.names)

поскольку это соответствует одному из cod_dest.

      regic_link %>% filter(cod_dest == 1100049)
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -62.0041 ymin: -11.9342 xmax: -61.4512 ymax: -11.4356
#> Geodetic CRS:  SIRGAS 2000
#>   cod_ori cod_dest                       geometry
#> 1 1100015  1100049 LINESTRING (-62.0041 -11.93...

Наконец, мы можем преобразовать матрицу adj в inla.graph объект (но я думаю, что это не обязательно для функций INLA):

      INLA::inla.read.graph(regic_adj)
#> $n
#> [1] 8
#> 
#> $nnbs
#> [1] 3 1 1 1 1 1 1 1
#> 
#> $nbs
#> $nbs[[1]]
#> [1] 2 3 4
#> 
#> $nbs[[2]]
#> [1] 1
#> 
#> $nbs[[3]]
#> [1] 1
#> 
#> $nbs[[4]]
#> [1] 1
#> 
#> $nbs[[5]]
#> [1] 6
#> 
#> $nbs[[6]]
#> [1] 5
#> 
#> $nbs[[7]]
#> [1] 8
#> 
#> $nbs[[8]]
#> [1] 7
#> 
#> 
#> $graph.file
#> [1] NA
#> 
#> $cc
#> $cc$id
#> [1] 1 1 1 1 2 2 3 3
#> 
#> $cc$n
#> [1] 3
#> 
#> $cc$nodes
#> $cc$nodes[[1]]
#> [1] 1 2 3 4
#> 
#> $cc$nodes[[2]]
#> [1] 5 6
#> 
#> $cc$nodes[[3]]
#> [1] 7 8
#> 
#> 
#> 
#> attr(,"class")
#> [1] "inla.graph"

^{Создано 2021-08-25 пакетом REPEX (v2.0.0)}

Надеюсь, это поможет. Вы можете проверить полный пример удаления LIMIT 5 из аргумента запроса в st_read. Могут быть более простые подходы, но сейчас я не могу придумать ни одного.