routing and isochrones in R

code/classroom/QgisinR.Rmd

@@ -45,7 +45,7 @@ TRIPSgeo_freg = st_read("geo/TRIPSgeo_freg.gpkg", quiet = TRUE)
 
 Represent Transport Zones with Total, and with Car %.
 
-```{r}
+```{r carper}
 # create Car_per variable
 TRIPSgeo_mun = TRIPSgeo_mun |> mutate(Car_per = Car / Total * 100)
 TRIPSgeo_freg = TRIPSgeo_freg |> mutate(Car_per = Car / Total * 100)
@@ -59,7 +59,7 @@ mapview(TRIPSgeo_freg, zcol = "Car_per", col.regions = rev(hcl.colors(9, "-Infer
 
 ## Geometric centroids
 
-```{r}
+```{r geocentroid}
 CENTROIDSgeo = st_centroid(TRIPSgeo_mun)
 # mapview(CENTROIDSgeo)
 ```
@@ -77,7 +77,7 @@ It is not that easy to estimate weighted centroids with R. See [here](https://wz
 We will make a bridge connection to QGIS to use its native function of mean coordinates.  
 
 
-```{r}
+```{r weighted}
 library(qgisprocess)
 # qgis_search_algorithms("mean") # search the exact function name
 # qgis_get_argument_specs("native:meancoordinates") |> select(name, description) # see the required inputs
@@ -178,11 +178,14 @@ mapview(SURVEYeuclidean)
 SURVEY$distance = st_length(SURVEYeuclidean) # compute distance and add directly in the first layer
 
 summary(SURVEY$distance) # in meters
+
+# to remove the units - can be useful
+SURVEY$distance = units::drop_units(SURVEY$distance)
 ```
 
 The same function can be used to find the closest GIRA station to each survey home location. And also check where are the ones that are far away from GIRA.
 
-```{r}
+```{r girahub}
 GIRA = st_read("geo/GIRA2023.geojson", quiet = TRUE) # we can also read geojson with this function!
 
 nearest = st_nearest_feature(SURVEY, GIRA) # creates an index of the closest GIRA station id
@@ -195,9 +198,220 @@ mapview(SURVEY, zcol = "distanceGIRA") +
 
 ## Routing distance
 
+We use the [openrouteservice-r](https://giscience.github.io/openrouteservice-r/) package. For that you need to [create an account](https://openrouteservice.org/dev/#/signup) and get a Token / api key.
+
+> **Note on how to store credentials in R**  
 Using an API key from [OpenRouteService](https://openrouteservice.org/dev/#/signup), you should store it at your computer and **never show it directly on code**.  
-For that `usethis::edit_r_environ()` and paste your token as `ORS_API_KEY="xxxxxxxxxxxxxxxxxxxxxx"` (replace with your token). Save the .Renviron file, and press `Ctrl+Shift+F10` to restart R so it can take effect.
+For that `usethis::edit_r_environ()` and paste your token as `ORS_API_KEY="xxxxxxxxxxxxxxxxxxxxxx"` (replace with your token).    
+Save the .Renviron file, and press `Ctrl+Shift+F10` to restart R so it can take effect.
+
+See the [documentation](https://giscience.github.io/openrouteservice-r/articles/openrouteservice.html) for more details.
+
+### Distances 1 point to many points
+
+Estimate the time and distance by `foot-waking` and `driving-car`, `fastest` mode, from survey locations (under 2 km^[for speed-up purposes - api request limit up to 40 / minute]) to IST.
+
+```{r routingdataprep}
+# devtools::install_github("GIScience/openrouteservice-r")
+library(openrouteservice)
+# ors_api_key(Sys.getenv("ORS_API_KEY")) # one time setup
+
+# get coordinates variable
+SURVEY$coordinates = st_coordinates(SURVEY)
+IST$coordinates = st_coordinates(IST)
+
+# Filter only the locations up to 2km euclidean
+SURVEYsample = SURVEY |> filter(distance <= 2000)
+# nrow(SURVEYsample) # 95
+```
+
+
+Although it is the same algorithm, here it works differently from QGIS.
+
+There are many ways of doing this. If we want to know only time and distance, and **not the route** itself, we can use the `ors_matrix()`. See example [here](https://web.tecnico.ulisboa.pt/~rosamfelix/gis/trips/timedistancematrix.html#Distance_and_time_matrix).  
+If we need the route, we should use the function `ors_directions()`. This one is not that easy to set-up because the function is prepared to retrieve only one result per request :( So we do a loop. Don't worry, it is not that
+
+```{r orsloop1, eval=FALSE, include=TRUE}
+ROUTES_foot = data.frame() # initial empty data frame
+
+# loop - the origin (i) is the survey location, and the IST is always the same destination
+for (i in 1:nrow(SURVEYsample)) {
+  ROUTES1 = ors_directions(
+    data.frame(
+      lon = c(SURVEYsample$coordinates[i, 1], IST$coordinates[1, 1]),
+      lat = c(SURVEYsample$coordinates[i, 2], IST$coordinates[1, 2])
+    ),
+    profile = "foot-walking", # or driving-car cycling-regular cycling-electric
+    preference = "fastest", # or shortest
+    output = "sf"
+  )
+  ROUTES1$distance = ROUTES1$summary[[1]]$distance # extract these values from summary
+  ROUTES1$duration = ROUTES1$summary[[1]]$duration
+  
+  ROUTES_foot = rbind(ROUTES_foot, ROUTES1) # to keep adding in the same df
+}
+
+ROUTES_foot = ROUTES_foot |>
+  select(distance, duration, geometry) |> # discard unnecessary variables
+  mutate(ID = SURVEYsample$ID) # cbind with syrvey ID
+```
+
+Repeat the same for `car-driving`.
+
+```{r orsloop2, eval=FALSE, include=TRUE}
+ROUTES_car = data.frame() # initial empty data frame
+
+# loop - the origin (i) is the survey location, and the IST is always the same destination
+for (i in 1:nrow(SURVEYsample)) {
+  ROUTES1 = ors_directions(
+    data.frame(
+      lon = c(SURVEYsample$coordinates[i, 1], IST$coordinates[1, 1]),
+      lat = c(SURVEYsample$coordinates[i, 2], IST$coordinates[1, 2])
+    ),
+    profile = "driving-car", # or cycling-regular cycling-electric foot-walking
+    preference = "fastest", # or shortest
+    output = "sf"
+  )
+  ROUTES1$distance = ROUTES1$summary[[1]]$distance # extract these values from summary
+  ROUTES1$duration = ROUTES1$summary[[1]]$duration
+  
+  ROUTES_car = rbind(ROUTES_car, ROUTES1) # to keep adding in the same df
+}
+
+ROUTES_car = ROUTES_car |>
+  select(distance, duration, geometry) |> # discard unnecessary variables
+  mutate(ID = SURVEYsample$ID) # cbind with syrvey ID
+```
+
+```{r importexport1, include=FALSE}
+# st_write(ROUTES_foot, "original/routes_foot.geojson")
+# st_write(ROUTES_car, "original/routes_car.geojson")
+
+ROUTES_foot = st_read("original/routes_foot.geojson", quiet = TRUE)
+ROUTES_car = st_read("original/routes_car.geojson", quiet = TRUE)
+```
+
+## Compare distances
+
+We can compare the euclidean and routing distances that we estimated for the survey locations under 2 km.
+
+```{r distancessummary}
+summary(SURVEYsample$distance) # Euclidean
+summary(ROUTES_foot$distance) # Walk
+summary(ROUTES_car$distance) # Car
+```
+
+## Vizualise routes
+
+Visualize with transparency of 30%
+
+```{r maproutes}
+mapview(ROUTES_foot, alpha = 0.3)
+mapview(ROUTES_car, alpha = 0.3, color = "red")
+```
+
+We can also use the `overline()` [function from stplanr package](https://docs.ropensci.org/stplanr/reference/overline.html) to break up the routes when they overline, and add them up.
+
+```{r overline, message=FALSE, warning=FALSE}
+library(stplanr)
+
+# we create a value that we can later sum, it also could be the number of trips represented by this route. in this case is only one respondent per route
+ROUTES_foot$trips = 1 
+
+ROUTES_foot_overline = overline(
+  ROUTES_foot,
+  attrib = "trips",
+  fun = sum
+)
+
+mapview(ROUTES_foot_overline, zcol = "trips", lwd = 3)
+```
+
+*   How many people are entering IST by the stairs near *Bar de Civil*?  
+*   And by the North gate?
+*   And from Alameda stairs?
 
-<!-- Use [openrouteservice-r](https://giscience.github.io/openrouteservice-r/) package devtools::install_github("GIScience/openrouteservice-r") -->
+# Buffers vs. Isochones and Service Areas
+
+## Buffer
+
+Represent a buffer of 500 m and 2000 m from IST^[Here I selected only the first variable because now we also have the coordinates information (unecessary for this procedure)].
+
+```{r bufferIST}
+
+# BUFFERist500 = st_buffer(IST, dist = 500) # non  projected - results may be weird
+BUFFERist500 = geo_buffer(IST[1], dist = 500) # from stplnar, to make sure it is in meters.
+BUFFERist2000 = geo_buffer(IST[1], dist = 2000)
+
+mapview(BUFFERist500) + mapview(BUFFERist2000, alpha.regions = 0.5)
+```
+
+## Isochrone
+
+### Isochrone from 1 point - distance
+
+We use again the `openrouteservice` r package.
+
+```{r isoch1}
+ISOCist = ors_isochrones(
+  IST$coordinates,
+  profile = "foot-walking",
+  range_type = "distance", # or time
+  range = c(500, 1000, 2000),
+  output = "sf"
+)
+
+ISOCist = arrange(ISOCist, -value) # to make the larger polygons on top of the table so the are displayed behind.
+
+mapview(ISOCist, zcol = "value", alpha.regions = 0.5)
+```
+
+As you can see, the distance buffer of 500m is larger than the isochrone of 500m.
+Actually we can measure their area of reach.
+
+```{r arearatio}
+ISOCist$area = st_area(ISOCist)
+BUFFERist500$area = st_area(BUFFERist500)
+BUFFERist2000$area = st_area(BUFFERist2000)
+
+ratio1 = BUFFERist500$area / ISOCist$area[ISOCist$value == 500] # 1.71
+ratio2 = BUFFERist2000$area / ISOCist$area[ISOCist$value == 2000] # 1.22
+```
+
+The euclidean buffer of 500m is `r round(ratio1, 2)` times larger than its isochrone, and the buffer of 2000m is `r round(ratio2, 2)` times larger than its isochrone.
+
+
+### Isochrone from more than 1 point - time
+
+For this purpose we will use the high schools dataset.
+
+```{r getschools}
+# import schools
+SCHOOLS = st_read("geo/SCHOOLS_basicsec.gpkg", quiet = TRUE)
+
+SCHOOLS$coordinates = st_coordinates(SCHOOLS) # create coordinate variable
+
+SCHOOLShigh = SCHOOLS |> filter(Nivel == "Secundario") # filter the high schools
+
+```
+
+And proceed with the time isochrones.
+
+```{r isoch2, eval=FALSE, include=FALSE}
+coor = data.frame(lon = SCHOOLS$coordinates[,1], lat = SCHOOLS$coordinates[,2])
+coor = coor[1:5,] # error of api, get a loop to overpass this!
+
+ISOCist = ors_isochrones(
+  coor,
+  profile = "foot-walking",
+  range_type = "time", # or distance
+  range = c(5, 20),
+  output = "sf"
+)
+
+ISOCist = arrange(ISOCist, -value) # to make the larger polygons on top of the table so the are displayed behind.
+
+mapview(ISOCist, zcol = "value", alpha.regions = 0.5)
+```
 
 *Work In Progress*