How to handle datasets that don't list elements in a manner that reflects directionality #145

andtheWings · 2021-03-20T22:42:33Z

andtheWings
Mar 20, 2021

If you execute the following code and compare the 'from' and 'to' attributes of the edges to the 'ONE_WAY_DIR' variable, it seems like the elements of each LINESTRING are not listed in such a way to reflect proper directionality.

If I'm correct, how would you go about fixing? Here's a proposed workflow, but I'm not sure how to execute in code (assuming you are in the Eastern and Northern hemispheres):

Check value of 'ONE_WAY_DIR'
- If 'NORTH', check latitude of 'from' and 'to'
  - If 'from' is less than 'to', leave as is
  - If 'from' is greater than 'to', apply st_reverse
- If 'SOUTH', check latitude of 'from' and 'to'
  - If 'from' is less than 'to', apply st_reverse
  - If 'from' is greater than 'to', leave as is
- If 'EAST', check longitude of 'from' and 'to'
  - If 'from' is less than 'to', leave as is
  - If 'from' is greater than 'to', apply st_reverse
- If 'WEST', check longitude of 'from' and 'to'
  - If 'from' is less than 'to', apply st_reverse
  - If 'from is greater than 'to', leave as is

library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.1.3, PROJ 7.1.0
library(sfnetworks)
library(tidyverse)

Indy_streets <- st_read("https://opendata.arcgis.com/datasets/fa9ec663cf25407a9b7645ff14334a7f_13.geojson") 
#> Reading layer `Street_Centerlines' from data source `https://opendata.arcgis.com/datasets/fa9ec663cf25407a9b7645ff14334a7f_13.geojson' using driver `GeoJSON'
#> Simple feature collection with 69475 features and 53 fields
#> geometry type:  MULTILINESTRING
#> dimension:      XY
#> bbox:           xmin: -86.35518 ymin: 39.62654 xmax: -85.93584 ymax: 39.94967
#> geographic CRS: WGS 84

one_way_Indiana_ave_streets <- Indy_streets %>%
    # Crop to Indiana Ave District
    st_crop(c(xmin = -86.189233, ymin = 39.770670, xmax = -86.160613, ymax = 39.789350)) %>%
    # Convert geometries from MULTILINESTRINGS to LINESTRINGS
    st_cast("LINESTRING") %>%
    # Only keep directionality variable
    select(ONE_WAY_DIR, geometry) %>%
    # Filter for one-way streets
    filter(ONE_WAY_DIR %in% c("EAST", "WEST", "NORTH", "SOUTH")) %>%
    # Convert to network object
    as_sfnetwork()
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
#> Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only
    
edges <- one_way_Indiana_ave_streets %>%
    activate("edges") %>%
    st_as_sf()

nodes <- one_way_Indiana_ave_streets %>%
    activate("nodes") %>%
    st_as_sf()

^{Created on 2021-03-20 by the reprex package (v1.0.0)}

agila5 · 2021-03-21T20:11:05Z

agila5
Mar 21, 2021
Collaborator

Hi @andtheWings and thanks for your comment. If the POINTS composing each LINESTRING are specified using the same direction as in ONE_WAY_DIR, then I think you can adopt the to_spatial_directed morpher. For example:

# packages
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(sfnetworks)
library(dplyr)
library(tidygraph)

# data
Indy_streets <- st_read(
  "https://opendata.arcgis.com/datasets/fa9ec663cf25407a9b7645ff14334a7f_13.geojson"
)
#> Reading layer `Street_Centerlines' from data source `https://opendata.arcgis.com/datasets/fa9ec663cf25407a9b7645ff14334a7f_13.geojson' using driver `GeoJSON'
#> Simple feature collection with 69475 features and 53 fields
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -86.35518 ymin: 39.62654 xmax: -85.93584 ymax: 39.94967
#> Geodetic CRS:  WGS 84

# Crop to Indiana Ave District
Indy_streets <- Indy_streets %>%
  st_crop(c(xmin = -86.189233, ymin = 39.770670, xmax = -86.160613, ymax = 39.789350))
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

# Convert geometries from MULTILINESTRINGS to LINESTRINGS
Indy_streets <- Indy_streets %>% 
  st_cast("MULTILINESTRING") %>% 
  st_cast("LINESTRING")
#> Warning in st_cast.sf(., "LINESTRING"): repeating attributes for all sub-
#> geometries for which they may not be constant

# Only keep directionality variable and filter for oneway streets
Indy_streets_oneway <- Indy_streets %>% 
  select(ONE_WAY_DIR) %>% 
  filter(ONE_WAY_DIR %in% c("EAST", "WEST", "NORTH", "SOUTH"))

# create an UNDIRECTED sfnetwork object
indy_sfn <- as_sfnetwork(Indy_streets_oneway, directed = FALSE)

# conver to a spatially directed network. 
indy_sfn %>% 
  convert(to_spatial_directed, .clean = TRUE)
#> # A sfnetwork with 148 nodes and 143 edges
#> #
#> # CRS:  WGS 84 
#> #
#> # A directed simple graph with 10 components with spatially explicit edges
#> #
#> # Node Data:     148 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: -86.18923 ymin: 39.77067 xmax: -86.16061 ymax: 39.78935
#>               geometry
#>            <POINT [°]>
#> 1 (-86.16398 39.77146)
#> 2 (-86.16694 39.77152)
#> 3 (-86.16854 39.77437)
#> 4  (-86.17004 39.7744)
#> 5 (-86.16639 39.78734)
#> 6 (-86.16639 39.78856)
#> # ... with 142 more rows
#> #
#> # Edge Data:     143 x 4
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: -86.18923 ymin: 39.77067 xmax: -86.16061 ymax: 39.78935
#>    from    to ONE_WAY_DIR                                               geometry
#>   <int> <int> <chr>                                             <LINESTRING [°]>
#> 1     1     2 EAST        (-86.16398 39.77146, -86.16452 39.77147, -86.16672 39~
#> 2     3     4 WEST                       (-86.16854 39.77437, -86.17004 39.7744)
#> 3     5     6 NORTH       (-86.16639 39.78734, -86.16637 39.78775, -86.16639 39~
#> # ... with 140 more rows

^{Created on 2021-03-21 by the reprex package (v1.0.0)}

Session info

sessioninfo::session_info()
#> - Session info ---------------------------------------------------------------
#>  setting  value                       
#>  version  R version 4.0.4 (2021-02-15)
#>  os       Windows 10 x64              
#>  system   x86_64, mingw32             
#>  ui       RTerm                       
#>  language (EN)                        
#>  collate  Italian_Italy.1252          
#>  ctype    Italian_Italy.1252          
#>  tz       Europe/Berlin               
#>  date     2021-03-21                  
#> 
#> - Packages -------------------------------------------------------------------
#>  package         * version    date       lib
#>  abind             1.4-5      2016-07-21 [1]
#>  assertthat        0.2.1      2019-03-21 [1]
#>  backports         1.2.1      2020-12-09 [1]
#>  class             7.3-18     2021-01-24 [2]
#>  classInt          0.4-3      2020-04-07 [1]
#>  cli               2.3.1      2021-02-23 [1]
#>  colorspace        2.0-0      2020-11-11 [1]
#>  crayon            1.4.1      2021-02-08 [1]
#>  DBI               1.1.1      2021-01-15 [1]
#>  deldir            0.2-10     2021-02-16 [1]
#>  digest            0.6.27     2020-10-24 [1]
#>  dplyr           * 1.0.5      2021-03-05 [1]
#>  e1071             1.7-5      2021-03-15 [1]
#>  ellipsis          0.3.1      2020-05-15 [1]
#>  evaluate          0.14       2019-05-28 [1]
#>  fansi             0.4.2      2021-01-15 [1]
#>  fs                1.5.0      2020-07-31 [1]
#>  generics          0.1.0      2020-10-31 [1]
#>  ggplot2           3.3.3      2020-12-30 [1]
#>  glue              1.4.2      2020-08-27 [1]
#>  goftest           1.2-2      2019-12-02 [1]
#>  gtable            0.3.0      2019-03-25 [1]
#>  highr             0.8        2019-03-20 [1]
#>  htmltools         0.5.1.1    2021-01-22 [1]
#>  igraph            1.2.6      2020-10-06 [1]
#>  KernSmooth        2.23-18    2020-10-29 [2]
#>  knitr             1.31       2021-01-27 [1]
#>  lattice           0.20-41    2020-04-02 [2]
#>  lifecycle         1.0.0      2021-02-15 [1]
#>  lwgeom            0.2-5      2020-06-12 [1]
#>  magrittr          2.0.1      2020-11-17 [1]
#>  Matrix            1.3-2      2021-01-06 [2]
#>  mgcv              1.8-33     2020-08-27 [2]
#>  munsell           0.5.0      2018-06-12 [1]
#>  nlme              3.1-152    2021-02-04 [2]
#>  pillar            1.5.1      2021-03-05 [1]
#>  pkgconfig         2.0.3      2019-09-22 [1]
#>  polyclip          1.10-0.001 2021-03-16 [1]
#>  proxy             0.4-25     2021-03-05 [1]
#>  purrr             0.3.4      2020-04-17 [1]
#>  R6                2.5.0      2020-10-28 [1]
#>  Rcpp              1.0.6      2021-01-15 [1]
#>  reprex            1.0.0      2021-01-27 [1]
#>  rlang             0.4.10     2020-12-30 [1]
#>  rmarkdown         2.7        2021-02-19 [1]
#>  rpart             4.1-15     2019-04-12 [2]
#>  scales            1.1.1      2020-05-11 [1]
#>  sessioninfo       1.1.1      2018-11-05 [1]
#>  sf              * 0.9-8      2021-03-17 [1]
#>  sfheaders         0.4.0      2020-12-01 [1]
#>  sfnetworks      * 0.5.0.9000 2021-03-21 [1]
#>  spatstat          2.0-1      2021-03-13 [1]
#>  spatstat.core     2.0-0      2021-03-16 [1]
#>  spatstat.data     2.1-0      2021-03-16 [1]
#>  spatstat.geom     2.0-0      2021-03-18 [1]
#>  spatstat.linnet   2.0-0      2021-03-18 [1]
#>  spatstat.sparse   2.0-0      2021-03-16 [1]
#>  spatstat.utils    2.1-0      2021-03-16 [1]
#>  stringi           1.5.3      2020-09-09 [1]
#>  stringr           1.4.0      2019-02-10 [1]
#>  styler            1.3.2      2020-02-23 [1]
#>  tensor            1.5        2012-05-05 [1]
#>  tibble            3.1.0      2021-02-25 [1]
#>  tidygraph       * 1.2.0      2020-05-12 [1]
#>  tidyr             1.1.3      2021-03-03 [1]
#>  tidyselect        1.1.0      2020-05-11 [1]
#>  units             0.7-0      2021-02-25 [1]
#>  utf8              1.2.1      2021-03-12 [1]
#>  vctrs             0.3.6      2020-12-17 [1]
#>  withr             2.4.1      2021-01-26 [1]
#>  xfun              0.22       2021-03-11 [1]
#>  yaml              2.2.1      2020-02-01 [1]
#>  source                                    
#>  CRAN (R 4.0.3)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.3)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.3)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  Github (baddstats/polyclip@55623e8)       
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  local                                     
#>  CRAN (R 4.0.4)                            
#>  Github (spatstat/spatstat.core@22b0c87)   
#>  Github (spatstat/spatstat.data@1611456)   
#>  Github (spatstat/spatstat.geom@ecdc5ca)   
#>  Github (baddstats/spatstat.linnet@4d457b9)
#>  CRAN (R 4.0.4)                            
#>  Github (spatstat/spatstat.utils@840166e)  
#>  CRAN (R 4.0.3)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.3)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.4)                            
#>  CRAN (R 4.0.3)                            
#> 
#> [1] C:/Users/Utente/Documents/R/win-library/4.0
#> [2] C:/Program Files/R/R-4.0.4/library

A few comments:

I'm not 100% sure that I understand your question, so please check the output of the reprex and check if the nodes are specified following the direction specified in the ONE_WAY_DIR column.
If you need to convert a mixture of LINESTRING and MULTILINESTRING to LINESTRING (as it occurs in this example), then I think you should 1) convert everything to MULTILINESTRING and 2) convert everything to LINESTRING. If you skip the first step, then you may lose some geometries. For example:

# packages
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

# fake data
mixture_lines <- st_sfc(
  st_multilinestring(
    list(rbind(c(0, 0), c(0, 1)), rbind(c(0, 1), c(1, 1))), 
  ), 
  st_linestring(rbind(c(1, 1), c(2, 2)))
)

# plot
par(mar = rep(0, 4))
plot(mixture_lines)

# skip first step
plot(st_cast(mixture_lines, "LINESTRING"))
#> Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

# better approach
plot(mixture_lines %>% st_cast("MULTILINESTRING") %>% st_cast("LINESTRING"))

^{Created on 2021-03-21 by the reprex package (v1.0.0)}

2 replies

andtheWings Mar 26, 2021
Author

Thank you for the detailed response! Unfortunately, the POINTS composing each LINESTRING are not specified using the same direction as in ONE_WAY_DIR. Edge 1 is a good example where its ONE_WAY_DIR value is EAST, but if you plot the from and to nodes, they're actually directed West.

Also thanks for the tip on converting all elements to MULTILINESTRING before converting to LINESTRING!

agila5 Mar 28, 2021
Collaborator

Hi @andtheWings, a small update. I just found that you can estimate the angle between two points specified on the sphere using geosphere::bearing. For example:

# packages
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(tidygraph)
#> 
#> Attaching package: 'tidygraph'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(sfnetworks)
library(geosphere)

# data
Indy_streets <- st_read(
  "https://opendata.arcgis.com/datasets/fa9ec663cf25407a9b7645ff14334a7f_13.geojson"
)
#> Reading layer `Street_Centerlines' from data source `https://opendata.arcgis.com/datasets/fa9ec663cf25407a9b7645ff14334a7f_13.geojson' using driver `GeoJSON'
#> Simple feature collection with 69477 features and 53 fields
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -86.35518 ymin: 39.62654 xmax: -85.93584 ymax: 39.94967
#> Geodetic CRS:  WGS 84

# Crop to Indiana Ave District
Indy_streets <- Indy_streets %>%
  st_crop(c(xmin = -86.189233, ymin = 39.770670, xmax = -86.160613, ymax = 39.789350))
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

# Convert geometries from MULTILINESTRINGS to LINESTRINGS
Indy_streets <- Indy_streets %>% 
  st_cast("MULTILINESTRING") %>% 
  st_cast("LINESTRING")
#> Warning in st_cast.sf(., "LINESTRING"): repeating attributes for all sub-
#> geometries for which they may not be constant

# Only keep directionality variable and filter for oneway streets
Indy_streets_oneway <- Indy_streets %>% 
  select(ONE_WAY_DIR) %>% 
  filter(ONE_WAY_DIR %in% c("EAST", "WEST", "NORTH", "SOUTH"))

Now we can estimate the angle between each pair of startpoint and endpoint using geosphere::bearing. Check ?geosphere::bearing for more details.

(my_bearing <- bearing(
  p1 = lwgeom::st_startpoint(Indy_streets_oneway) %>% st_coordinates(), 
  p2 = lwgeom::st_endpoint(Indy_streets_oneway) %>% st_coordinates()
))
#>   [1]  -88.44454205  -88.18603238   -0.01888859    1.71529260    1.71561712
#>   [6]    1.71529371  -90.32558456  -88.32793795    1.71539175   -0.28888808
#>  [11]   -0.28753259   -1.23811504    1.71019725    1.29459307  -50.18765764
#>  [16]  -93.57405364  -87.94928986  -89.18559072  -90.40482168  -90.95525975
#>  [21]  -92.10191415  124.31591932  -90.95536007   -0.28973679    1.71545389
#>  [26]  -88.18523670   30.69540708  -90.96250690   -0.28741107  -75.02173326
#>  [31]  -84.72831774  -84.91007920    2.97263040 -177.02975429  -89.04190690
#>  [36] -119.19412137  -44.56864390  -91.43828841  -90.95565169  -90.36973310
#>  [41]    1.71531282    9.96945495    1.71563451 -104.63975530  -88.40667852
#>  [46]  -88.18935323  -87.64458163    1.89018094  -87.85911313  -94.40913487
#>  [51]  -94.67527139  -89.18686799  -84.79893610  -88.18691115   99.55331499
#>  [56]    1.71556057  -88.62972236    1.71514135  -90.85228877  -87.70032972
#>  [61]    1.71528901    1.71555674   13.09037022 -173.20392891 -171.80191553
#>  [66]  -96.48259660  -93.57263414  -90.33965323 -172.98666787   -0.28963888
#>  [71]  -92.23531129  -90.90180955  -88.18615338  -80.45829614  -18.10265804
#>  [76]  -88.32368012   92.99391792  -87.85994480    1.53175498   -0.11004324
#>  [81]   93.23571949   89.43505701    1.71654467   96.25401821   37.23735095
#>  [86]  -88.12388879  -70.62192311  -89.18771493  -76.38893408  -87.98133970
#>  [91]  -68.38775021  -85.28671741   -0.55172032   97.80832472   -9.03834386
#>  [96]   -0.46362358    1.01720782    1.71615382  171.04380061  -87.08525235
#> [101]  -88.32263005  -88.32309468  -88.18495248    1.71505639   -3.37721384
#> [106]  -90.45552506  -87.86270097  -88.18730996  -17.04192519 -112.48297227
#> [111]  -90.45703277  112.16162023 -178.77320796  -90.63602189  -90.43095547
#> [116]    1.71561671   65.50625444   30.53907603  -88.19139861 -111.94071865
#> [121]  -54.05097305  -85.26655838  -75.49092939  -16.64923893  -90.58324090
#> [126]  -88.17936505  -90.64088980  -19.17020788 -177.34117404  -88.38350146
#> [131]  -89.15866199  -90.34132525   34.52599073  -90.68877054  -77.60193592
#> [136]   -0.09955898   -0.47917108  -88.18329150  -90.10913100  -90.16806387
#> [141]  -81.55845541 -112.48411274  -90.16607011

How would you deal with segments that are specified in a non-standard direction? For example:

par(mfrow = c(2, 2), mar = rep(0, 4))
plot(Indy_streets_oneway[15, ] %>% st_geometry())
plot(Indy_streets_oneway[37, ] %>% st_geometry())
plot(Indy_streets_oneway[117, ] %>% st_geometry())
plot(Indy_streets_oneway[121, ] %>% st_geometry())

^{Created on 2021-03-28 by the reprex package (v1.0.0)}

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

How to handle datasets that don't list elements in a manner that reflects directionality #145

Uh oh!

{{title}}

Uh oh!

Replies: 1 comment 2 replies

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Select a reply

Uh oh!

How to handle datasets that don't list elements in a manner that reflects directionality #145

Uh oh!

andtheWings Mar 20, 2021

Replies: 1 comment · 2 replies

Uh oh!

agila5 Mar 21, 2021 Collaborator

Uh oh!

andtheWings Mar 26, 2021 Author

Uh oh!

agila5 Mar 28, 2021 Collaborator

andtheWings
Mar 20, 2021

Replies: 1 comment 2 replies

agila5
Mar 21, 2021
Collaborator

andtheWings Mar 26, 2021
Author

agila5 Mar 28, 2021
Collaborator