You need to sign up for the openrouteservice to perform the following exercises. This allows you to create an API key that allows the use of the ORS. See: https://openrouteservice.org/dev/#/signup

The free API comes with some restrictions with respect to the number of operations you are allowed to execute per day and per minute (see https://openrouteservice.org/plans/). But this will not harm us here.

1 Helper functions

A helper function that eases the communication with the API a bit. It handles the processing of the data to ensure that a properly formated data.frame is returned.

getOsomeStat <- function(uri, valueFieldName="value", ...)
{
  results <- postForm(uri, binary=FALSE, .encoding = "utf-8", ...)
  
  #browser()
  resultList <- RJSONIO::fromJSON(results, simplify = TRUE)
  resultsDf <- data.frame(do.call("rbind", (resultList$result)))
  # make sure the right data types are used
  # for users we have  fromTimestamp  and  toTimestamp fields, not timestampe
  if(length(grep(x=names(resultsDf), pattern = "timestamp"))> 0)
  {
    resultsDf$timestamp <- parse_datetime( as.character(resultsDf$timestamp))  
  }
  if(length(grep(x=names(resultsDf), pattern = "fromTimestamp"))> 0)
  {
    resultsDf$fromTimestamp <- parse_datetime( as.character(resultsDf$fromTimestamp))  
  }
  if(length(grep(x=names(resultsDf), pattern = "toTimestamp"))> 0)
  {
    resultsDf$toTimestamp <- parse_datetime( as.character(resultsDf$toTimestamp))  
  }
  # rename value field
  resultsDf$value <- as.numeric(as.character(resultsDf$value))
  idxValueField <- which(names(resultsDf)=="value")
  names(resultsDf)[idxValueField] <- valueFieldName
  
  
  return(resultsDf)
}
adminLimaCallao <- st_read(dsn="data_day5/lima.gpkg", layer="limaCallao")
## Reading layer `limaCallao' from data source `/home/slautenb/Documents/lehre/HD/ws_2020_21/heikaLab/R4UrbanDataLab_2020/data_day5/lima.gpkg' using driver `GPKG'
## Simple feature collection with 2 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -77.25156 ymin: -12.51993 xmax: -76.62082 ymax: -11.57244
## geographic CRS: WGS 84
limaBbox <- st_bbox(adminLimaCallao)
limaBboxStr <- paste(limaBbox, collapse = ", ")
response <- RCurl::postForm("https://api.ohsome.org/v1/elements/centroid", bboxes = limaBboxStr, filter = "amenity=hospital", time = "2020-05-01",  binary=FALSE, properties="tags")

content <- geojsonio::as.json(response)
## Registered S3 method overwritten by 'geojson':
##   method        from     
##   print.geojson geojsonsf
hospitalsLima <- geojsonio::geojson_sf(content)
nrow(hospitalsLima)
## [1] 130
st_crs(hospitalsLima) # needs to be WGS84 - which it is (no surprise)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

2 Calculate isochrones

theCoords <- data.frame(st_coordinates(hospitalsLima))
names(theCoords) <-c("lon", "lat")
head(theCoords)
## # A tibble: 6 x 2
##     lon   lat
##   <dbl> <dbl>
## 1 -77.0 -11.9
## 2 -77.0 -12.2
## 3 -77.1 -12.0
## 4 -77.1 -12.0
## 5 -77.1 -12.1
## 6 -76.9 -12.2

We need a loop since 5 is the maximum number of locations ors_isochrones accept. Would need to stop at 500 since this is the imit per day. We also need to pause every 20 requests snce this is the limit per minute.

I run the code and saved the data to avoid performing the same operation over and over again. If you look into the RMakrdown code you will see that I used the eval=FALSE and include=FALSE chunk options to hide this a bit from the resulting HTML document.

isocCar <- NULL
nHospitals <- nrow(hospitalsLima)
# nrow(theCoords)
for(i in 1:nHospitals)
{
  if(i%%20==0)
  {
    cat(paste(round(i/nrow(theCoords)*100), "% processed.\n"))
    Sys.sleep(70)
  }
  
  testRes <- try(
    {
      res <- ors_isochrones(locations= theCoords[i,], range = c(300, 600, 1200, 2400, 3600), output = "sf", profile="driving-car", range_type="time",  attributes=c("area", "reachfactor", "total_pop"))  
    })
  if(inherits(testRes, "try-error"))
    next
  res$osm_id <- hospitalsLima$X.osmId[i]
  if(is.null(isocCar))
  { 
    isocCar <- res
  } else {
    isocCar <- rbind(isocCar, res)
  }
  
  
}

dim(isocCar)
head(st_drop_geometry(isocCar))

2.1 Export

save(isocCar, file= "data_day5/isocCar_hospitalsLima.Rdata")
st_write(isocCar, dsn="data_day5/isochronesLima.gpkg", layer="iisocCar_hospitalsLima", delete_layer=TRUE)

3 Analysis of the isochrones

Since the isocrones overlap each other it is useful to separate them and plot them in the right order. Another alternative would be to clip them.

isocCar300 <- isocCar %>% filter(value == 300)
isocCar600 <- isocCar %>% filter(value == 600)
isocCar1200 <- isocCar %>% filter(value == 1200)
isocCar2400 <- isocCar %>% filter(value == 2400)
isocCar3600 <- isocCar %>% filter(value == 3600)
theColors <- RColorBrewer::brewer.pal(5, "YlOrBr")
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(isocCar3600) + tm_fill(col=theColors[5]) + 
  tm_shape(isocCar2400) + tm_fill(col=theColors[4]) + 
  tm_shape(isocCar1200) + tm_fill(col=theColors[3]) + 
  tm_shape(isocCar600) + tm_fill(col=theColors[2]) + 
  tm_shape(isocCar300) + tm_fill(col=theColors[1]) + 
  tm_shape(hospitalsLima) + tm_dots()+
  tm_scale_bar()

3.1 Population data

The ORS retunred already the population data per isochrone area based on WorldPop data. However, we need to be clear about that the population cannot be added across isochrones related to different hospital since this would involve heavy double counting: overlapping isochrones share the population in the overlapping areas. What we could do however is to look at the population of the isochrones of each hospital.

isocCar %>% ggplot(mapping=aes(x=osm_id, y=total_pop)) + 
  geom_point(size=.5) + facet_grid(~value) + 
  xlab("Hospital") + ylab("Population") +
  labs(title= "Population in isochrones of hospitals in Lima", subtitle = "Isochrones are in minutes driving distance by car")