В настоящий момент я работаю над проектом с событиями точечных паттернов в линейной сети (автомобильные аварии) и читаю главу 17 spatstat
книги: «Пространственные точечные паттерны: методология и приложения с R».
Авторы книги объясняют, что они определили новый класс объектов под названием lpp
для анализа точечных паттернов в линейной сети. Каркас каждого lpp
объекта - это объект linnet
, и есть несколько функций для создания объекта linnet
. Для моего приложения соответствующими функциями являются linnet
и as.linnet
. Функция linnet
создает линейный сетевой объект из пространственного положения каждой вершины и информации о том, какие вершины соединены ребром, в то время как функция as.linnet
может применяться к psp
объекту, который преобразуется в linnet
объекты, определяя связность сети, используя заданный порог расстояния.
Причина, по которой я задаю этот вопрос, заключается в том, что я не знаю, как эффективно создать объект linnet, начиная с объекта sf
с геометрией LINESTRING
. Насколько мне известно, можно преобразовать объект sf
в объект sp
(т.е. объект SpatialLines
), затем я могу преобразовать объект sp
в объект psp
(с помощью функции as.psp
), а затем я могу преобразовать объект psp
в объект объект linnet с помощью функции as.psp.linnet
(которая определена в пакете maptools). Основная проблема с этим подходом (как заявили авторы пакета в своей книге) заключается в том, что предполагаемая сеть неверна каждый раз, когда в моих сетевых данных возникает эстакада или подземный переход, поскольку соответствующая сеть будет создавать искусственные пересечения в новой сети. Более того, как написали авторы в своей книге, код становится экспоненциально медленнее.
Следующий код представляет собой упрощенную версию того, что я делал до сих пор, но я думаю, что должен быть более простой и лучший способ создать объект linnet из объекта sf
. Я бы использовал функцию linnet
, но проблема в том, что я не знаю, как создать (разреженную) матрицу смежности для соответствующих вершин сети или матрицу связей между краями сети.
# packages
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: nlme
#> Loading required package: rpart
#>
#> spatstat 1.61-0 (nickname: 'Puppy zoomies')
#> For an introduction to spatstat, type 'beginner'
#>
#> Note: spatstat version 1.61-0 is out of date by more than 11 weeks; a newer version should be available.
library(maptools)
#> Loading required package: sp
#> Checking rgeos availability: TRUE
library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright
# download data
iow_polygon <- getbb("Isle of Wight, South East, England", format_out = "sf_polygon", featuretype = "state") %>%
st_transform(crs = 27700)
iow_highways <- st_read("https://download.geofabrik.de/europe/great-britain/england/isle-of-wight-latest.osm.pbf", layer = "lines") %>%
st_transform(crs = 27700)
#> Reading layer `lines' from data source `https://download.geofabrik.de/europe/great-britain/england/isle-of-wight-latest.osm.pbf' using driver `OSM'
#> Simple feature collection with 44800 features and 9 fields
#> geometry type: LINESTRING
#> dimension: XY
#> bbox: xmin: -5.716262 ymin: 43.35489 xmax: 1.92832 ymax: 51.16517
#> epsg (SRID): 4326
#> proj4string: +proj=longlat +datum=WGS84 +no_defs
# subset the data otherwise the code takes ages
iow_highways <- iow_highways[iow_polygon, ] %>%
subset(grepl(pattern = c("primary|secondary|tertiary"), x = highway))
# transform as sp
iow_highways_sp <- as(iow_highways %>% st_geometry(), "Spatial")
# transform as psp
iow_highways_psp <- as.psp(iow_highways_sp)
# transform as linnet
iow_highways_linnet <- as.linnet.psp(iow_highways_psp, sparse = TRUE)
Я могу извлечь координаты каждой вершины сети
stplanr::line2points(iow_highways)
#> Simple feature collection with 2814 features and 1 field
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 430780.7 ymin: 75702.05 xmax: 464851.7 ymax: 96103.72
#> epsg (SRID): 27700
#> proj4string: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
#> First 10 features:
#> id geometry
#> 1 1 POINT (464851.7 87789.73)
#> 2 1 POINT (464435.4 88250.85)
#> 3 2 POINT (464390.9 87412.27)
#> 4 2 POINT (464851.7 87789.73)
#> 5 3 POINT (462574.6 88987.62)
#> 6 3 POINT (462334.6 88709.92)
#> 7 4 POINT (464066.9 87576.84)
#> 8 4 POINT (464390.9 87412.27)
#> 9 5 POINT (464420 88227.79)
#> 10 5 POINT (464398.7 88225.33)
но тогда я не знаю, как построить матрицу смежности.
Создано 2 декабря 2019 г. пакетом REPEX (v0.3.0)