Subdivision morpher does nothing on my network

[Jean-Romain]

I'm learning sfnetwork and can't reproduce with my data some examples I previously reproduced with the vignettes data or other data of mine. I have a network that seems ok to me but for which the edges must be subdivided. However to_spatial_subdivision() does nothing and returns 48 edges like the original input. I'm expecting something lke 60 or 70 edges. Yet, the linestrings are intersecting at rounded and equal coordinates so start and end points coordinates are matching with other linestring nodes. If there is something wrong in my input ? You can get the data here. Thanks.

library(sf)
library(sfnetworks)
library(tidygraph)
col = sample(rainbow(100))

net = st_read("network.gpkg")
net = st_geometry(net)
net = as_sfnetwork(net)
subdivision = convert(net, to_spatial_subdivision)
subdivision

par(mfrow = c(1,2), mar = c(0,0,0,0))
plot(st_geometry(net, "edges"), col = col, lwd =3)
plot(st_geometry(net, "nodes"), pch = 20, cex = 1, add = TRUE)
plot(st_geometry(subdivision, "edges"), col = col, lwd = 3)
plot(st_geometry(subdivision, "nodes"), pch = 20, cex = 1, add = TRUE)

[agila5]

Hi @Jean-Romain. I think you may have some problems with the input data.

1. The LINESTRINGS stored in network.gpkg have a Z dimension which explains why to_spatial_subdivision does not split any segment (since the common points probably have a different value for the Z covariate);

# packages
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE

seg = st_read("C:/Users/Andrea Gilardi/Downloads/network.gpkg", quiet = TRUE)
seg
#> Simple feature collection with 48 features and 1 field
#> Geometry type: LINESTRING
#> Dimension:     XYZ
#> Bounding box:  xmin: 269011.3 ymin: 5154036 xmax: 277985.3 ymax: 5161983
#> z_range:       zmin: 1 zmax: 177
#> Projected CRS: NAD83(CSRS) / MTM zone 9
#> First 10 features:
#>    level                           geom
#> 1      1 LINESTRING Z (269011.3 5160...
#> 2      1 LINESTRING Z (275391.5 5154...
#> 3      1 LINESTRING Z (276392.2 5154...
#> 4      1 LINESTRING Z (277985.3 5158...
#> 5      2 LINESTRING Z (270575 516073...
#> 6      2 LINESTRING Z (270983 516040...
#> 7      2 LINESTRING Z (271255 516030...
#> 8      2 LINESTRING Z (272537 516004...
#> 9      2 LINESTRING Z (272827 515979...
#> 10     2 LINESTRING Z (273639 515923...

^{Created on 2022-06-22 by the reprex package (v2.0.1)}

Do you know what the Z coordinate represents?

2. If we ignore the Z dimension, then your input LINESTRINGS have several repeated internal points (which I think it's weird):

# packages
options(width = 100)
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE

# data
seg = st_read("C:/Users/Andrea Gilardi/Downloads/network.gpkg", quiet = TRUE)
seg = st_geometry(seg)

# I focus on an even simpler example
seg = seg[1]
st_cast(seg, "POINT")[32]
#> Geometry set for 1 feature 
#> Geometry type: POINT
#> Dimension:     XYZ
#> Bounding box:  xmin: 269181 ymin: 5160641 xmax: 269181 ymax: 5160641
#> z_range:       zmin: 2 zmax: 2
#> Projected CRS: NAD83(CSRS) / MTM zone 9
#> POINT Z (269181 5160641 2)
st_cast(seg, "POINT")[33]
#> Geometry set for 1 feature 
#> Geometry type: POINT
#> Dimension:     XYZ
#> Bounding box:  xmin: 269181 ymin: 5160641 xmax: 269181 ymax: 5160641
#> z_range:       zmin: 3 zmax: 3
#> Projected CRS: NAD83(CSRS) / MTM zone 9
#> POINT Z (269181 5160641 3)

# or even
sprintf("%.15f", st_coordinates(st_cast(seg, "POINT")[32]))
#> [1] "269181.000000000000000"  "5160641.000000000000000" "2.000000000000000"
sprintf("%.15f", st_coordinates(st_cast(seg, "POINT")[33]))
#> [1] "269181.000000000000000"  "5160641.000000000000000" "3.000000000000000"

^{Created on 2022-06-22 by the reprex package (v2.0.1)}

What's going on?

[Jean-Romain]

Yes it has Z but Z can be removed if it is an issue. Z is expected to be the elevation of the roads (but here it contains dummy data). I removed Z this way

net = st_read("network.gpkg")
net = st_zm(net)

But in this case to_spatial_subdivision subdivides on every single node so I get 3000 edges

[Jean-Romain]

Ok I see the duplicated points. This may be related to previous coordinate rounding step I performed to solve my issue with subdivision. I'll regenerate a clean network to see if I still have duplicated points.

[Jean-Romain]

I fixed my code that generates the lines to avoid duplicated node every 100 m and the subdivision morpher now works like charm. Thank you for identifying the issue 👍