1 Στόχοι της άσκησης

2 Links

3 Βιβλιοθήκες

Εκτελούμε το παρακάτω μονο μια φορά. Μετα τοποθετουμε comments # και δεντο ξανα-εκτελουμε:

#library(devtools)
#devtools::install_github("dimitrisk/goal", force=T)

library(remotes)
remotes::install_github("dimitrisk/goal", force=T)
#remotes::install_version("rgeos", version = "0.6-4")
#remotes::install_version("simplevis", version = "7.1.0")
#devtools::install_github("dimitrisk/goal", force=T)

Εγκατάσταση ολων των βιβλιωθηκών.

install.packages(c(
  "sf", "sfnetworks", "tidygraph", "igraph", "dplyr",
  "tibble", "ggplot2", "units", "tmap", "osmdata", "nabor", "sfheaders"
))

Το παρακάτω μέρος, φορτώνει τις βιβλιοθηκες:

library(sf)
library(sfnetworks)
library(tidygraph)
library(igraph)
library(dplyr)
library(tibble)
library(ggplot2)
library(units)
library(tmap)
library(osmdata)
#library(link2GI)
library(nabor)
library(goal)

4 Δεδομένα δικτύου

Εδώ ορίζουμε το BoundingBox (ΒΒ) της δικής μας περιοχής μελέτης! Κατεβάζει οδικό δίκτυο εντός του BB της περιοχής μελέτης.

q=c(26.545029,39.088569,26.570177,39.116810) # Μυτιλήνη
net2 = goal::osm.getRoads(q, withBB=TRUE, outcrs=4326)


Using bbox

poly = goal::osm.bb_2_pol(q, outcrs =  4326) # bbox σε spatial polygon
net3 = goal::osm.ClipSFnetwork_with_poly(net2, poly) # clip network by spatial polygon

plot(net3,col="grey", main="Clipped sfnetwork of Mytilini")
plot(poly,add=T)

# Create a random point
#gps = sfheaders::sf_point(data.frame(y = 26.55257, x = 39.10575, ID=-1))  %>% st_set_crs(4326)

# nearest edge (road) to the point. The network must have edges activated.
#near_edge = st_nearest_feature(gps, net3 %>% st_as_sf())

#near_edge
#st_as_sf(net3)[near_edge,]

#p3 = ggplot() +
  #geom_sf(data = st_as_sf(net3), color = 'black') +
  #geom_sf(data = gps, color = 'red') +
  #geom_sf(data = st_as_sf(net3)[near_edge,], color = 'orange')
#p3
#plot(net3)
#net3

Non-Isolated Nodes

net3 %>%  sfnetworks::activate("nodes") %>% dplyr::filter(!tidygraph::node_is_isolated())

# A sfnetwork with 1292 nodes and 744 edges
#
# CRS:  EPSG:4326 
#
# A directed acyclic simple graph with 549 components with spatially explicit edges
#
# Node data: 1,292 × 1 (active)
             geometry
          <POINT [°]>
1  (26.55257 39.1052)
2 (26.54515 39.10575)
3   (26.561 39.10523)
4 (26.55839 39.10629)
5 (26.55782 39.09519)
6 (26.54682 39.09666)
# ℹ 1,286 more rows
#
# Edge data: 744 × 6
   from    to osm_id   name             highway                         geometry
  <int> <int> <chr>    <chr>            <chr>                   <LINESTRING [°]>
1     1     2 39428970 Ζωοδόχου Πηγής   primary (26.55257 39.1052, 26.55211 39.…
2     3   630 40847716 8ης Νοεμβρίου    primary (26.561 39.10523, 26.56096 39.1…
3   631   632 40847725 Αργύρη Εφταλιώτη primary (26.56032 39.10481, 26.561 39.1…
# ℹ 741 more rows

Isolated Nodes

net3 %>% sfnetworks::activate("nodes") %>% dplyr::filter(tidygraph::node_is_isolated()) %>% st_as_sf() %>% nrow()

[1] 0

Loop Edges

net3 %>% sfnetworks::activate("edges") %>% dplyr::filter(tidygraph::edge_is_loop()) %>% st_as_sf() %>% nrow()

[1] 0

Multiple Edges (has any parallel siblings)

net3 %>% sfnetworks::activate("edges") %>% dplyr::filter(tidygraph::edge_is_multiple()) %>% st_as_sf() %>% nrow()

[1] 0

Simplify network.

net3b = net3 %>% activate("edges") %>%
  filter(!edge_is_multiple()) %>%
  filter(!edge_is_loop())

5 Ενσωμάτωση των σημείων

Κατασκευή τοπολογικά ορθού δικτύου (net4) από τα κατεβασμένα δεδομένα της περιοχής μελέτης. Subdivide network Smooth pseudo nodes

net4 = tidygraph::convert(net3b, to_spatial_subdivision)
net4 = tidygraph::convert(net4, to_spatial_smooth)
net4 = tidygraph::convert(net4, to_spatial_simple)

Πόσοι κόμβοι στο δίκτυο μας (nodes)?

n_nodes = net4 %>% activate(nodes) %>% st_as_sf()%>% nrow()

Σε κάθε δρόμο (edges) πρόσθεσε μια νέα στήλη “length” με το μήκος του κάθε δρόμου.

net4 = net4 %>% activate(edges) %>%
  mutate(length = edge_length())%>% st_set_crs(4326)

Σε κάθε κόμβο (nodes), πρόσθεσε μια νέα στήλη “ID” αλλά και μια νέα στήλη “ishouse”=1.

net4 = net4 %>% activate(nodes) %>%
  mutate(ID = 1:n_nodes)%>%
  mutate(ishouse = 1)%>% st_set_crs(4326)

Αυτά είναι τα επιλεγμένα σημεία καταφυγής μας. Τα βρήκαμε από το GoogleΜaps και μπορεί να είναι πάρκα, εκκλησίες, σχολεία, γήπεδα.

gps1 = sfheaders::sf_point(data.frame(y = 26.55295, x = 39.10692))  %>% st_set_crs(4326)
gps2 = sfheaders::sf_point(data.frame(y = 26.55200, x = 39.10100 ))  %>% st_set_crs(4326)
gps3 = sfheaders::sf_point(data.frame(y = 26.560898277367066, x = 39.090741799541824 ))  %>% st_set_crs(4326)
gps4 = sfheaders::sf_point(data.frame(y = 26.553421104184043, x = 39.11297446115607))  %>% st_set_crs(4326)


mycol = c("blue","green",'red',"purple")
evac_points = rbind(gps1, gps2, gps3,gps4) #%>% rowid_to_column() # Νέα στήλη 'rowid'

Ενσωμάτωση των χώρων καταφυγής, στο δίκτυο μας ως 4 νέους κόμβους (nodes):

blended = st_network_blend(net4, evac_points  )

Άλλαξε τις τιμές των στηλών: ‘isevac’, ‘ishouse’ ανάλογα:

blended = blended %>% activate(nodes) %>%
  mutate(isevac = ifelse(is.na(ishouse), 1, 0) ) %>%
  mutate(ishouse = ifelse(is.na(ishouse), 0, 1) )

Οι 10 τελευταίοι κόμβοι (nodes) του δικτύου μας μέχρι στιγμής. Φαίνονται οι 4 νέοι κόμβοι (nodes) isevac = 1

blended %>%  activate("nodes") %>% st_as_sf() %>% as.data.frame() %>% tail(10)

                      geometry .tidygraph_node_index   ID ishouse isevac
1010  POINT (26.5514 39.11576)                  1010 1010       1      0
1011 POINT (26.56026 39.10863)                  1011 1011       1      0
1012 POINT (26.55188 39.10789)                  1012 1012       1      0
1013  POINT (26.5592 39.10647)                  1013 1013       1      0
1014 POINT (26.55063 39.11499)                  1014 1014       1      0
1015 POINT (26.55493 39.10828)                  1015 1015       1      0
1016 POINT (26.56105 39.09064)                    NA   NA       0      1
1017 POINT (26.55217 39.10101)                    NA   NA       0      1
1018 POINT (26.55342 39.11296)                    NA   NA       0      1
1019 POINT (26.55322 39.10683)                    NA   NA       0      1

Αφαίρεση των disconnected

table(components(blended)$membership)


   1    2    3    4    5    6    7    8 
1003    2    2    2    2    2    2    4

blended = blended %>% activate('nodes')%>%
  filter(components(blended)$membership == 1)

6 Ανάθεση

Το δίκτυο μας μέχρι στιγμής:

ggplot() +
  geom_sf(data = st_as_sf(net4%>%  activate("edges")), color = 'grey') +
  geom_sf(data = st_as_sf(net4%>%  activate("nodes")), color = 'grey')+
  geom_sf(data = evac_points, color = mycol, cex=3, pch=17)

Τα ‘rowids’ των χώρων καταφυγής isevac=1 :

rowids_evac = blended %>%  activate("nodes") %>% as.data.frame() %>%  rowid_to_column() %>% filter(isevac == 1) %>% pull(rowid)
tail(rowids_evac)

[1] 1000 1001 1002 1003

Τα ‘rowids’ των σπιτιών ishouse=1 :

rowids_houses = blended %>%  activate("nodes") %>% as.data.frame() %>%  rowid_to_column() %>% filter(ishouse == 1) %>% pull(rowid)
tail(rowids_houses)

[1] 994 995 996 997 998 999

Δημιουργία ‘Spatial object’ των σπιτιών, χώρων καταφυγής. Θα μας χρειαστούν πιο μετά για εξαγωγή σε Shapefiles.

evac_sf = blended %>%  activate("nodes") %>% st_as_sf() %>% filter(isevac==1)
houses_sf = blended %>%  activate("nodes") %>% st_as_sf() %>% filter(ishouse==1)   #%>% rowid_to_column()

Πόσοι δρόμοι (edges)?

n_edges =  blended %>% activate(edges)  %>% st_as_sf()%>%  nrow()

Δημιουργία στήλης σε κάθε δρόμο (edge) με τίτλο στήλης IDedge

blended = blended %>% activate(edges) %>% mutate(IDedge = 1:n_edges)

Υπολογισμός του Distance Matrix:

dm = st_network_cost(blended, from =rowids_houses , to =rowids_evac ,  direction="all")
head(dm)

Units: [m]
         [,1]      [,2]      [,3]     [,4]
[1,] 2000.417  498.8066  968.3447 192.5096
[2,] 2418.773  917.1625 1386.7005 610.8654
[3,] 2539.549 1037.9381 1507.4761 731.6410
[4,] 2594.360 1092.7495 1562.2875 786.4524
[5,] 2640.730 1139.1193 1608.6574 832.8223
[6,] 2620.220 1118.6091 1645.8889 870.0538

Για κάθε σπίτι, εύρεση ποιου από τους 4 χώρους καταφυγής είναι πιο κοντά. Για κάθε σπίτι, κατασκευή δυο νέων στηλών:

closest_index: Ποιος χώρος καταφυγής είναι πιο κοντά?
closest_index_dist: Σε τι απόσταση είναι ο πιο κοντινός χώρος καταφυγής?

houses_sf$closest_index = apply(dm, 1, function(x) which(x == min(x))[1])
houses_sf$closest_index_dist = apply(dm, 1, function(x)  min(x)[1])

Οπτικοποίηση τελικής ανάθεσης στους 4 χώρους καταφυγής:

plot(blended, col="grey")
plot(st_geometry(houses_sf), cex=1.5, col=mycol[houses_sf$closest_index], pch=21, add=T)
plot(st_geometry(evac_sf) , cex=2, pch=17, add=T, col=mycol)
plot(poly, add=T)

6.1 Στατιστικά ανάθεσης

Πόσοι κόμβοι ανατέθηκαν σε κάθε χώρο καταφυγής?

table(houses_sf$closest_index)


  1   2   3   4 
152 294 199 354

Απόσταση μεταξύ κόμβων και χώρων καταφυγής (Ελάχιστη, Μέγιστη, Μέσος Όρος):

houses_sf %>% as.data.frame()%>% group_by(closest_index) %>%
  summarise( min_dist=min(closest_index_dist),
             max_dist=max(closest_index_dist),
             mean_dist=mean(closest_index_dist) )

# A tibble: 4 × 4
  closest_index min_dist max_dist mean_dist
          <int>    <dbl>    <dbl>     <dbl>
1             1     22.1    1805.      673.
2             2     34.5    1363.      555.
3             3     12.6    1004.      398.
4             4     51.0    1183.      522.

7 Περιορισμός απόστασης

Ορισμός μιας μεγίστης απόστασης περπατήματος.

apostasi = 700

dm2 = dm
dm2 = units::drop_units(dm2) #  Αφαίρεση μονάδων μέτρησης από το Distance Matrix
dm2[dm2>=apostasi] = NA


houses_sf$closest_index2 = apply(dm2, 1, function(x) which(x == min(x, na.rm=T))[1])
houses_sf$closest_index_dist2 = apply(dm2, 1, function(x)  min(x, na.rm=T)[1])

Οπτικοποίηση τελικής ανάθεσης στους 4 χώρους καταφυγής:

plot(blended, col="grey", main = sprintf("Με περιορισμό απόστασης %sμ", apostasi))
plot(st_geometry(houses_sf), cex=1.5, col=mycol[houses_sf$closest_index2],  pch=21, add=T)
plot(st_geometry(evac_sf) , cex=2, pch=17, add=T, col=mycol)
plot(poly, add=T)
plot(st_geometry(houses_sf),cex=0.5, add=T, col="grey", pch=20)

Πόσοι κόμβοι δεν εξυπηρετούνται?

table(houses_sf$closest_index2,  useNA="ifany")


   1    2    3    4 <NA> 
  75  207  182  264  271

8 Περιορισμός πλήθους

Κάθε χώρος καταφυγής μπορεί να εξυπηρετήσει μόνο μέχρι 55 άτομα.

dm3 = dm

df = as.data.frame(dm3)
df$whichMin = apply(dm3, 1, which.min)
df$minDistance = apply(dm3, 1, FUN=min, na.rm=T)

Ιεράρχηση ανά γκρούπ

library(dplyr)
df3 = df %>%
  group_by(whichMin) %>%
  mutate(my_ranks = order(order(minDistance, decreasing=F)))


df3

# A tibble: 999 × 7
# Groups:   whichMin [4]
      V1    V2    V3    V4 whichMin minDistance my_ranks
     [m]   [m]   [m]   [m]    <int>       <dbl>    <int>
 1 2000.  499.  968.  193.        4        193.       14
 2 2419.  917. 1387.  611.        4        611.      238
 3 2540. 1038. 1507.  732.        4        732.      277
 4 2594. 1093. 1562.  786.        4        786.      297
 5 2641. 1139. 1609.  833.        4        833.      313
 6 2620. 1119. 1646.  870.        4        870.      323
 7 2652. 1151. 1682.  906.        4        906.      334
 8 2564. 1422. 1385. 1000.        4       1000.      349
 9 2445. 1303. 1266.  881.        4        881.      325
10 2532. 1389. 1382.  967.        4        967.      344
# ℹ 989 more rows

df3$whichMin2 = NA
df3[df3$whichMin==1 & df3$my_ranks %in% c(1:55),]$whichMin2 = 1
df3[df3$whichMin==2 & df3$my_ranks %in% c(1:55),]$whichMin2 = 2
df3[df3$whichMin==3 & df3$my_ranks %in% c(1:55),]$whichMin2 = 3
df3[df3$whichMin==4 & df3$my_ranks %in% c(1:55),]$whichMin2 = 4

houses_sf$closest_index3 = df3$whichMin2

plot(blended, col="grey", main = sprintf("Με περιορισμό πλήθους %s ατόμων", 55))
plot(st_geometry(houses_sf), cex=1.5, col=mycol[houses_sf$closest_index3],  pch=21, add=T)
plot(st_geometry(evac_sf) , cex=2, pch=17, add=T, col=mycol)
plot(poly, add=T)
plot(st_geometry(houses_sf),cex=0.5, add=T, col="grey", pch=20)

9 Ερωτήματα

Μέσα στο .zip αυτού του εργαστηρίου, υπάρχει ένα έτοιμο RScript.r με το Bounding-Box της πόλης των Χανίων και τις συντεταγμένες για τους 4 χώρους καταφυγής. Χρησιμοποιήστε το και επαναλάβετε όλη τη παραπάνω διαδικασία και απαντήστε στα ακόλουθα ερωτήματα:

Πόσους κόμβους (σημεία) και ακμές (γραμμές) έχει το δίκτυο των Χανίων?
Ποια είναι τα 4 ‘rowids’ των 4 χώρων καταφυγής των Χανίων?

Αφού αναθέσετε όλους τους κόμβους στους 4 χώρους καταφυγής με βάση την απόσταση:

Οπτικοποίηση
Πόσοι κόμβοι ανατέθηκαν σε κάθε χώρο καταφυγής?
Ποιες οι τιμές για: “Ελάχιστη”, “Μέγιστη”, “Μέσος Όρος” απόστασης μεταξύ κόμβων και χώρων καταφυγής?

Αφού αναθέσετε κόμβους στους 4 χώρους καταφυγής με περιορισμό τη μέγιστη απόσταση 800μ:

Οπτικοποίηση
Πόσοι κόμβοι ανατέθηκαν σε κάθε χώρο καταφυγής?
Ποιες οι τιμές για: “Ελάχιστη”, “Μέγιστη”, “Μέσος Όρος” απόστασης μεταξύ κόμβων και χώρων καταφυγής?

A3: Spatial Network

Δρ. Δημήτρης Καβρουδάκης