如何从具有 LINESTRING 几何列的 sf 对象开始创建 linnet 对象？
How can I create a linnet object starting from an sf object with LINESTRING geometry column?
目前我正在处理线性网络（车祸）上的点模式事件的项目，我正在阅读 spatstat 书的第 17 章："Spatial Point Patterns: Methodology and Applications with R"。
本书的作者解释说，他们定义了一个名为 lpp 的新 class 对象，用于分析线性网络上的点模式。每个 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-02 创建
只是确认spatstat包没有提供处理其他格式的功能；我们的期望是 maptools 或其他包将提供格式转换代码；这对于 sf 对象格式尚不可用，大概是因为 sf 相对较新。
关键问题是具有 LINESTRING 几何形状的 sf 对象是否包含足够的信息来确定网络的连通性。如果是这样，那么我建议你制作一个 2 列矩阵，列出所有由边连接的顶点对，并调用 spatstat::linnet。如果不是，则可用数据不足...
最后请注意，GitHub 提供的 spatstat(1.61-0.061) 的当前开发版本比线性网络上的许多操作的当前版本 (1.61-0) 快得多.即将公开发布。
我不确定你为什么在去 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)
目前我正在处理线性网络（车祸）上的点模式事件的项目，我正在阅读 spatstat 书的第 17 章："Spatial Point Patterns: Methodology and Applications with R"。
本书的作者解释说，他们定义了一个名为 lpp 的新 class 对象，用于分析线性网络上的点模式。每个 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-02 创建只是确认spatstat包没有提供处理其他格式的功能；我们的期望是 maptools 或其他包将提供格式转换代码；这对于 sf 对象格式尚不可用，大概是因为 sf 相对较新。
关键问题是具有 LINESTRING 几何形状的 sf 对象是否包含足够的信息来确定网络的连通性。如果是这样，那么我建议你制作一个 2 列矩阵，列出所有由边连接的顶点对，并调用 spatstat::linnet。如果不是，则可用数据不足...
最后请注意，GitHub 提供的 spatstat(1.61-0.061) 的当前开发版本比线性网络上的许多操作的当前版本 (1.61-0) 快得多.即将公开发布。
我不确定你为什么在去 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)