data-raw/test-rnet_merge.md

Test out network merging functionality

# Get zip in working directory
list.files("~/Downloads/", pattern = "zip")
f = list.files("~/Downloads/", pattern = "OneDrive_1_19-09-2023.zip", full.names = TRUE)
f
# Unzip
unzip(f, exdir = "./data-raw/csna_data")

Read-in data:

library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.3     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.3     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
list.files("./data-raw/csna_data/", pattern = "shp")
character(0)
csna_data = sf::read_sf("csna_data/Thurrock Roads_polyline.shp")
csna_data = sf::st_transform(csna_data, "EPSG:4326")

# Get OSM data for area
csna_convex_hull = sf::st_convex_hull(sf::st_union(csna_data))
csna_centroid = sf::st_centroid(csna_convex_hull)
plot(csna_convex_hull)

osm_data = osmextract::oe_get_network(csna_centroid, "driving", extra_tags = c("highway", "name", "maxspeed"), boundary = csna_convex_hull, boundary_type = "clipsrc")
The input place was matched with Essex. 
The chosen file was already detected in the download directory. Skip downloading.
Starting with the vectortranslate operations on the input file!

0...10...20...30...40...50...60...70...80...90...100 - done.

Warning in CPL_gdalvectortranslate(source, destination, options, oo, doo, :
GDAL Message 1: A geometry of type MULTILINESTRING is inserted into layer lines
of geometry type LINESTRING, which is not normally allowed by the GeoPackage
specification, but the driver will however do it. To create a conformant
GeoPackage, if using ogr2ogr, the -nlt option can be used to override the layer
geometry type. This warning will no longer be emitted for this combination of
layer and feature geometry type.

Finished the vectortranslate operations on the input file!

Reading layer `lines' from data source 
  `/home/robin/data/osm/geofabrik_essex-latest.gpkg' using driver `GPKG'
Simple feature collection with 8993 features and 12 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 0.216299 ymin: 51.45155 xmax: 0.5495789 ymax: 51.57737
Geodetic CRS:  WGS 84
plot(osm_data$geometry)

# Length of networks
sum(sf::st_length(csna_data)) |>
  units::set_units("km")
927.9664 [km]
sum(sf::st_length(osm_data)) |>
  units::set_units("km")
1180.52 [km]

We’ll take a subset in a 500 m buffer in the centre of the study area.

buffer = sf::st_buffer(csna_centroid, 500)
osm_subset = sf::st_intersection(osm_data, buffer)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
csna_subset = sf::st_intersection(csna_data, buffer)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries

Let’s plot the two networks:

plot(csna_subset$geometry, col = "red")
plot(osm_subset$geometry, col = "blue", add = TRUE)

summary(csna_subset$CSNA_Level)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.000   2.000   2.000   2.517   3.000   6.000
funs = list(CSNA_level = max)
osm_joined = stplanr::rnet_merge(osm_subset, csna_subset, funs = funs)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries

Warning in st_cast.sf(sf::st_cast(x, "MULTILINESTRING"), "LINESTRING"):
repeating attributes for all sub-geometries for which they may not be constant

Joining with `by = join_by(osm_id)`
Joining with `by = join_by(osm_id)`
waldo::compare(osm_joined$geometry, csna_subset$geometry)
`old` is length 107
`new` is length 118

`attr(old, 'bbox')`: "((0.3597854,51.50738),(0.3718281,51.51591))"
`attr(new, 'bbox')`: "((0.3597861,51.50738),(0.3717286,51.51591))"

`attr(old, 'crs')$input`: "WGS 84"   
`attr(new, 'crs')$input`: "EPSG:4326"

`dim(old[[1]])`: 6 2
`dim(new[[1]])`: 4 2

`old[[1]]`: "LINESTRING (0.3672797 51.50..."
`new[[1]]`: "LINESTRING (0.3682288 51.51..."

`dim(old[[2]])`: 23 2
`dim(new[[2]])`:  2 2

`old[[2]]`: "LINESTRING (0.360145 51.510..."
`new[[2]]`: "LINESTRING (0.3633168 51.51..."

`old[[3]]`: "LINESTRING (0.370691 51.509..."
`new[[3]]`: "LINESTRING (0.3654734 51.51..."

`dim(old[[4]])`:  9 2
`dim(new[[4]])`: 12 2

`old[[4]]`: "LINESTRING (0.3661647 51.51..."
`new[[4]]`: "LINESTRING (0.361379 51.512..."

And 209 more differences ...
summary(osm_joined$geometry %in% osm_subset$geometry)
   Mode    TRUE 
logical     107
osm_joined
Simple feature collection with 107 features and 14 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 0.3597854 ymin: 51.50738 xmax: 0.3718281 ymax: 51.51591
Geodetic CRS:  WGS 84
First 10 features:
    osm_id        name      highway waterway aerialway barrier man_made access
1  4384852 Rowley Road  residential     <NA>      <NA>    <NA>     <NA>   <NA>
2  4384853   High Road    secondary     <NA>      <NA>    <NA>     <NA>   <NA>
3  4384855 School Lane  residential     <NA>      <NA>    <NA>     <NA>   <NA>
4  4384857  Pound Lane  residential     <NA>      <NA>    <NA>     <NA>   <NA>
5  4568525    Fen Lane unclassified     <NA>      <NA>    <NA>     <NA>   <NA>
6  4568593   Mill Lane unclassified     <NA>      <NA>    <NA>     <NA>   <NA>
7  4568594   Mill Lane unclassified     <NA>      <NA>    <NA>     <NA>   <NA>
8  4568595   Mill Lane      service     <NA>      <NA>    <NA>     <NA>   <NA>
9  4568596   The Green  residential     <NA>      <NA>    <NA>     <NA>   <NA>
10 4568597        <NA>      service     <NA>      <NA>    <NA>     <NA>   <NA>
   service maxspeed z_order
1     <NA>     <NA>       3
2     <NA>     <NA>       6
3     <NA>     <NA>       3
4     <NA>     <NA>       3
5     <NA>     <NA>       3
6     <NA>     <NA>       3
7     <NA>     <NA>       3
8     <NA>     <NA>       0
9     <NA>     <NA>       3
10    <NA>     <NA>       0
                                                               other_tags
1                                                                    <NA>
2                                                           "ref"=>"B188"
3                                 "sidewalk"=>"both","surface"=>"asphalt"
4                                                                    <NA>
5                                                                    <NA>
6                                                         "oneway"=>"yes"
7                                                                    <NA>
8  "designation"=>"public_bridleway","horse"=>"yes","motorcar"=>"private"
9                                                                    <NA>
10                                                                   <NA>
   CSNA_level                       geometry  length_x
1           3 LINESTRING (0.3672797 51.50... 119.95760
2           3 LINESTRING (0.360145 51.510... 509.84284
3          NA LINESTRING (0.370691 51.509... 241.89090
4           2 LINESTRING (0.3661647 51.51... 229.95311
5           3 LINESTRING (0.3606277 51.51...  69.71302
6          NA LINESTRING (0.3623346 51.51... 355.78291
7          NA LINESTRING (0.3636589 51.50... 125.19843
8          NA LINESTRING (0.3642193 51.50...  19.90700
9           2 LINESTRING (0.3665946 51.51... 194.09989
10         NA LINESTRING (0.368502 51.511...  67.64835

Let’s plot the result:osm_subset

osm_joined |>
  select(CSNA_level) |>
  plot()

Plot with tmap

library(tmap)
The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
     (status 2 uses the sf package in place of rgdal)

Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')
brks = 1:5
m1 = tm_shape(osm_joined) +
  tm_lines(col = "CSNA_level", palette = "viridis", breaks = brks)
m2 = tm_shape(csna_subset) +
  tm_lines(col = "CSNA_Level", palette = "viridis", breaks = brks)
tmap_arrange(m2, m1)
Warning: Values have found that are higher than the highest break

Warning: Values have found that are higher than the highest break

We can try to improve the fit by changing the arguments that feed into the rnet_merge function.

osm_joined2 = stplanr::rnet_merge(osm_subset, csna_subset, funs = funs, dist = 15)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries

Warning in st_cast.sf(sf::st_cast(x, "MULTILINESTRING"), "LINESTRING"):
repeating attributes for all sub-geometries for which they may not be constant

Joining with `by = join_by(osm_id)`
Joining with `by = join_by(osm_id)`
m1 = tm_shape(osm_joined2) +
  tm_lines(col = "CSNA_level", palette = "viridis", breaks = brks)
m2 = tm_shape(csna_subset) +
  tm_lines(col = "CSNA_Level", palette = "viridis", breaks = brks)
tmap_arrange(m2, m1)
Warning: Values have found that are higher than the highest break

Warning: Values have found that are higher than the highest break

Warning: Values have found that are higher than the highest break

Warning: Values have found that are higher than the highest break

We’ll set a max segment length.

args(stplanr::rnet_join)
function (rnet_x, rnet_y, dist = 5, length_y = TRUE, key_column = 1, 
    subset_x = TRUE, dist_subset = NULL, segment_length = 0, 
    endCapStyle = "FLAT", contains = TRUE, max_angle_diff = NULL, 
    ...) 
NULL
osm_joined3 = stplanr::rnet_merge(osm_subset, csna_subset, funs = funs, dist = 15, segment_length = 10)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries

Warning in st_cast.sf(sf::st_cast(x, "MULTILINESTRING"), "LINESTRING"):
repeating attributes for all sub-geometries for which they may not be constant

Joining with `by = join_by(osm_id)`
Linking to GEOS 3.11.1, GDAL 3.6.4, PROJ 9.1.1; sf_use_s2() is TRUE
Joining with `by = join_by(osm_id)`
m1 = tm_shape(osm_joined3) +
  tm_lines(col = "CSNA_level", palette = "viridis", breaks = brks)
m2 = tm_shape(csna_subset) +
  tm_lines(col = "CSNA_Level", palette = "viridis", breaks = brks)
tmap_arrange(m2, m1)
Warning: Values have found that are higher than the highest break

Warning: Values have found that are higher than the highest break

Warning: Values have found that are higher than the highest break

Warning: Values have found that are higher than the highest break

Finally we’ll set a max angle difference.

osm_joined4 = stplanr::rnet_merge(osm_subset, csna_subset, funs = funs, dist = 7, segment_length = 10, max_angle_diff = 30)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries

Warning in st_cast.sf(sf::st_cast(x, "MULTILINESTRING"), "LINESTRING"):
repeating attributes for all sub-geometries for which they may not be constant

Joining with `by = join_by(osm_id)`

Warning: st_centroid assumes attributes are constant over geometries

Joining with `by = join_by(osm_id)`
m1 = tm_shape(osm_joined4) +
  tm_lines(col = "CSNA_level", palette = "viridis", breaks = brks)
m2 = tm_shape(csna_subset) +
  tm_lines(col = "CSNA_Level", palette = "viridis", breaks = brks)
tmap_arrange(m2, m1)
Warning: Values have found that are higher than the highest break

Warning: Values have found that are higher than the highest break

Warning: Values have found that are higher than the highest break

Warning: Values have found that are higher than the highest break

Save the output as follows.

sf::st_write(osm_joined4, "osm_joined4.geojson")
Writing layer `osm_joined4' to data source 
  `osm_joined4.geojson' using driver `GeoJSON'
Writing 107 features with 14 fields and geometry type Line String.


ropensci/stplanr documentation built on Aug. 26, 2024, 9:48 p.m.