Определить отдельные пространственные линии на карте - PullRequest
Новая
       34

Определить отдельные пространственные линии на карте

0 голосов
/ 15 апреля 2020

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

Пример: 

dat <- structure(list(geometry = structure(list(structure(c(169.023627307075, 
  169.02315299228, -45.3068517761089, -45.3081870363656), .Dim = c(2L, 
    2L), class = c("XY", "LINESTRING", "sfg")), structure(c(169.01854883529, 
      169.018930977689, 169.02315299228, -45.3083004691879, -45.3083134946477, 
      -45.3081870363656), .Dim = 3:2, class = c("XY", "LINESTRING", 
        "sfg")), structure(c(169.02315299228, 169.0330663306, -45.3081870363656, 
          -45.3144702778175), .Dim = c(2L, 2L), class = c("XY", "LINESTRING", 
            "sfg")), structure(c(169.015997396195, 169.022945130719, -45.3119974282578, 
              -45.3168289670259), .Dim = c(2L, 2L), class = c("XY", "LINESTRING", 
                "sfg")), structure(c(169.022945130719, 169.032555385154, -45.3168289670259, 
                  -45.3163448854193), .Dim = c(2L, 2L), class = c("XY", "LINESTRING", 
                    "sfg")), structure(c(169.01868555271, 169.022945130719, -45.3174947235968, 
                      -45.3168289670259), .Dim = c(2L, 2L), class = c("XY", "LINESTRING", 
                        "sfg"))), n_empty = 0L, crs = structure(list(input = "4326", 
                          wkt = "GEOGCS[\"GCS_WGS_1984\",\n    DATUM[\"WGS_1984\",\n        SPHEROID[\"WGS_84\",6378137,298.257223563]],\n    PRIMEM[\"Greenwich\",0],\n    UNIT[\"Degree\",0.017453292519943295],\n    AUTHORITY[\"EPSG\",\"4326\"]]"), class = "crs"), class = c("sfc_LINESTRING", 
                            "sfc"), precision = 0, bbox = structure(c(xmin = 169.015997396195, 
                              ymin = -45.3174947235968, xmax = 169.0330663306, ymax = -45.3068517761089
                            ), class = "bbox"))), row.names = c(NA, 6L), class = c("sf", 
                              "data.frame"), sf_column = "geometry", agr = structure(integer(0), .Label = c("constant", 
                                "aggregate", "identity"), class = "factor"))

library(ggplot2)
library(sf)

ggplot() + 
  geom_sf(data = dat, size = 1) + 
  coord_sf()

enter image description here

In в этом тривиальном примере есть две «реки», но у меня нет способа идентифицировать их в кадре данных. Мне нужен такой столбец: 

dat$ID <- c("A", "A", "A", "B", "B", "B")

Единственный способ, с помощью которого я могу думать, - это сравнение значений координат концов линий, но это не очень эффективно. Любая помощь высоко ценится.

1 Ответ

1 голос
/ 15 апреля 2020

Расчет идентификаторов кластеров на основе пересечения можно выполнить с помощью igraph следующим образом: 

library(sf)
library(igraph)

dat <- structure(...)

# Detect clusters
m = st_intersects(dat, sparse = FALSE)
g = graph_from_adjacency_matrix(m)
dat$id = clusters(g)$membership
dat

## Simple feature collection with 6 features and 1 field
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 169.016 ymin: -45.31749 xmax: 169.0331 ymax: -45.30685
## geographic CRS: WGS 84
##                         geometry id
## 1 LINESTRING (169.0236 -45.30...  1
## 2 LINESTRING (169.0185 -45.30...  1
## 3 LINESTRING (169.0232 -45.30...  1
## 4 LINESTRING (169.016 -45.312...  2
## 5 LINESTRING (169.0229 -45.31...  2
## 6 LINESTRING (169.0187 -45.31...  2

# Plot
plot(st_geometry(dat))
text(st_coordinates(st_centroid(dat)), as.character(dat$id))

enter image description here

...