14  Buffers vs. Isochones

We will use the stplnar package to create buffers in meters, and the openrouteservice package to create isochrones.

See the ors setup instructions for more details.

14.1 Buffer

Represent a buffer of 500 m and 2000 m from IST1.

IST = st_sfc(st_point(c(-9.1397404, 38.7370168)), crs = 4326) # create a point
IST$coordinates = st_coordinates(IST)

# 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)

14.2 Isochrone

Isochrone from 1 point - distance

We use again the openrouteservice r package.

ISOCist = ors_isochrones(
  IST$coordinates,
  profile = "foot-walking",
  range_type = "distance", # or time
  range = c(500, 1000, 2000),
  output = "sf",
  api_key = Sys.getenv("ORS_USHIFT")
)

ISOCist = arrange(ISOCist, -value) |>  # to make the larger polygons on top of the table so the are displayed behind.
          select(-center) # unnecessary variable
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.

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 1.7 times larger than its isochrone, and the buffer of 2000m is 1.23 times larger than its isochrone.

Isochrone from more than 1 point - time

For this purpose we will use the high schools dataset.

# 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
SCHOOLShigh = SCHOOLShigh |> mutate(id = 1:nrow(SCHOOLShigh)) # provide and id

# list of XY coordinates for ORS
coord_schools = data.frame(lon = SCHOOLShigh$coordinates[, 1],
                           lat = SCHOOLShigh$coordinates[, 2])
                           # id = SCHOOLShigh$id)

And proceed with the time isochrones, for a range of 20 min, with 5 min intervals.

ISOCschools = ors_isochrones(
  coord_schools,
  profile = "foot-walking",
  range_type = "time", # or distance
  range = 20*60, # 20 minutes in seconds
  interval = 5*60, # to have intervals of 5 minutes
  attributes = "area", #you can directly get area, population, and so on.
  output = "sf",
  api_key = Sys.getenv("ORS_USHIFT")
)

ISOCschools = arrange(ISOCschools, -value) |> # larger polygons on top of the table so the are displayed behind.
  select(-center)

And now merge this information with the schools’ names.

ISOCschools = ISOCschools |>
  mutate(id = group_index + 1, # because group_index starts at 0
         value = value/60) |> 
    left_join(SCHOOLShigh |>
                st_drop_geometry() |>
                select(id, INF_NOME, Alunos))
mapview(ISOCschools, zcol = "value", alpha.regions = 0.5)

And estimate areas for the 20min isochrones

summary(ISOCschools$area[ISOCschools$value==20])/1000000 # in km²
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  3.295   4.303   5.363   5.098   5.804   6.511

  1. Here I selected only the first variable because now we also have the coordinates information (unnecessary for this procedure)↩︎