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

【问题标题】：How can I create a linnet object starting from an sf object with LINESTRING geometry column?如何从具有 LINESTRING 几何列的 sf 对象开始创建 linnet 对象？
【发布时间】：2019-12-02 18:30:28
【问题描述】：

目前，我正在研究一个线性网络上的点模式事件（车祸）项目，我正在阅读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函数），然后我可以使用 as.psp.linnet 函数（在 maptools 包中定义）将 psp 对象转换为 linnet 对象。这种方法的主要问题（正如该软件包的作者在他们的书中所说）是每次在我的网络数据中出现立交桥或地下通道时，推断的网络都是错误的，因为相应的 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 (v0.3.0) 于 2019 年 12 月 2 日创建}

【问题讨论】：

标签： r sf spatstat

【解决方案1】：

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

iow_highways_linnet <- as.linnet.SpatialLines(iow_highways_sp)

这会将SpatialLines 直接转换为linnet 并融合线路共享一个顶点。我认为地下通道不会与除非两条线在交点处都有顶点。请参见下面的示例：

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)

【讨论】：

对不起，我没有看到你的答案。我从未使用过as.linnet.SpatialLines 函数，因为我不知道它存在。无论如何，我刚刚在 IOW 示例上测试了您的方法，并且运行顺利。谢谢！

【解决方案2】：

只是确认spatstat包不提供处理其他格式的功能；我们的期望是maptools或其他包会提供格式转换代码；这还不适用于sf 对象格式，大概是因为sf 相对较新。

关键问题是具有 LINESTRING 几何形状的 sf 对象是否包含足够的信息来确定网络的连通性。如果是这样，那么我建议您制作一个 2 列矩阵，列出所有由边连接的顶点对，然后调用 spatstat::linnet。如果没有，那么可用的数据还不够……

最后请注意，GitHub 提供的spatstat(1.61-0.061) 的当前开发版本比线性网络上的许多操作的当前版本 (1.61-0) 快得多。很快就会公开发布。

【讨论】：