r - 如何从具有 LINESTRING 几何列的 sf 对象开始创建 linnet 对象？

Question

目前，我正在开发一个涉及线性网络上的点模式事件(车祸)的项目，并且我正在阅读 spatstat 一书的第 17 章:“空间点模式:方法论和R 的应用程序”。

本书的作者解释说，他们定义了一类名为 lpp 的新对象，用于分析线性网络上的点模式。每个lpp对象的骨架是一个linnet对象，并且有几个函数可以创建linnet对象。对于我的应用程序，相关函数是 linnet 和 as.linnet。函数 linnet 根据每个顶点的空间位置以及有关哪些顶点由边连接的信息创建线性网络对象，而 as.linnet 函数可应用于一个 psp 对象，该对象转换为 linnet 对象，使用指定的距离阈值推断网络的连接性。

我问这个问题的原因是我不知道如何从具有 LINESTRING 几何形状的 sf 对象开始有效地创建 linnet 对象。据我所知，可以将 sf 对象转换为 sp 对象(即 SpatialLines 对象)，然后我可以将sp 对象转换为 psp 对象(使用 as.psp 函数)，然后我可以将 psp 对象转换为使用 as.psp.linnet 函数(在 maptools 包中定义)的 linnet 对象。这种方法的主要问题(正如该包的作者在他们的书中所说)是，每次网络数据中出现立交桥或地下 channel 时，推断的网络都是错误的，因为相应的 linnet 会在网络中创建人工交叉点。此外，正如作者在书中所说，代码速度呈指数级下降。

以下代码是我迄今为止所做的简化版本，但我认为必须有一种更简单、更好的方法从 sf 对象创建 linnet 对象。我会使用 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)  

但是我不知道如何构建邻接矩阵。

^{由reprex package于2019年12月2日创建(v0.3.0)}

Answer 1

我不知道你为什么要浏览 psp格式在路上linnet 。尝试将第一个代码块的最后两行替换为:

iow_highways_linnet <- as.linnet.SpatialLines(iow_highways_sp)

这会转换SpatialLines直接联系linnet和 fuse 线共享一个顶点。我不认为地下 channel 会与立交桥，除非两条线在交点处都有顶点。请参阅下面的示例:

l1 <- sf::st_linestring(matrix(c(-1,1,-1,1,1,1,-1,-1), ncol = 2))
l2 <- sf::st_linestring(matrix(c(-1,-1,1,1,2,1,-1,-2), ncol = 2))
l_sf <- sf::st_sf(id = 1:2, geom = sf::st_sfc(l1,l2))
l_sp <- sf::as_Spatial(l_sf)
l <- maptools::as.linnet.SpatialLines(l_sp)
plot(l)