How can I create a linnet object starting from an sf object with LINESTRING geometry column?

Question

At the moment I'm working on a project with point pattern events on a linear network (car crashes) and I'm reading chapter 17 of spatstat book: "Spatial Point Patterns: Methodology and Applications with R".

The authors of the book explain that they defined a new class of objects called lpp for analyzing point patterns on a linear network. The skeleton of each lpp object is a linnet object and there are several functions to create a linnet object. For my application the relevant functions are linnet and as.linnet. The function linnet creates a linear network object from the spatial location of each vertex and information about which vertices are joined by an edge, while the as.linnet function can be applied to a psp object which is transformed into linnet objects inferring the connectivity of the network using a specified distance threshold.

The reason why I'm asking this question is that I don't know how to efficiently create a linnet object starting from a sf object with a LINESTRING geometry. As far as I know, it's possible to transform the sf object into an sp object (i.e. a SpatialLines object), then I can transform the sp object into a psp object (using as.psp function) and then I can transform the psp object into a linnet object using the as.psp.linnet function (which is defined in the maptools package). The main problem with this approach (as the authors of the package said in their book) is that the inferred network is wrong every time an overpass or an underpass occurs in my network data since the corresponding linnet will create artificial intersections in the nework. Moreover, as the authors said in their book, the code gets exponentially slower.

The following code is a simplified version of what I did so far but I think that there must be an easier and better way to create a linnet object from an sf object. I would use the linnet function but the problem is that I don't know how to create a (sparse) adjacency matrix for the corresponding vertices of the network or a matrix of links between the edges of the network.

# 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)

I can extract the coordinates of each vertex of the network

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)

but then I don't know how to build the adjacency matrix.

^{Created on 2019-12-02 by the reprex package (v0.3.0)}

Ege Rubak Ege Rubak · Accepted Answer · 2019-12-04T14:31:50

I’m not sure why you go through the psp format on the way to linnet. Try to replace the last two lines of your first code chunk by:

iow_highways_linnet <- as.linnet.SpatialLines(iow_highways_sp)

This converts SpatialLines directly to linnet and fuses lines that share a vertex. I don’t think an underpass will be fused to an overpass unless both lines have a vertex at the intersection point. See example below:

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)

How can I create a linnet object starting from an sf object with LINESTRING geometry column?

2 Answers