Accessibility Map Example - Burkina Faso
s5joschi
July 18, 2019
Introduction
This tutorial showcases the creation of Accessibility Maps also referred to as Travel Time Maps in R using the AccessibilityMaps package. We will create maps that exhibit travel times from any given location in Burkina Faso to cities with more than 1 Million inhabitants inside the country. We will create two distinct maps, one for the dry-season where unpaved roads function better for car transport and one for the rainy season where travel speed on unpaved roads is dramatically reduced.
The AccessibilityMaps package provides mostly wrapper functions to existing packages and FOSS GIS-libraries in and outside of R. Most geodata processing is done with GDAL which is more efficient and less depended on RAM in contrast to the raster package in R. This allows you to process large raster data-sets with several million raster-cells to create high resolution maps for larger research areas. The Travel Time Maps are calculated using the r.cost function from GRASS (for a documentation see here).
Installation
In order to use the package you need to have GRASS GIS and GDAL installed. Both are automatically installed if you install Quantum GIS which is also a good software solution for visualizing the results and creating good locking, publishable maps. If you prefer to install GRASS and GDAL as standalone versions please refer to https://gdal.org/ and https://grass.osgeo.org/download/ to obtain the latest stable versions. It is highly recommended to use the stable versions instead of the latest, but unstable releases. Please make sure that you have GRASS and GDAL installed before you proceed.
The easiest way to install the AccessibilityMaps package is using the devtools package in R which allows you to download and install everything within R. The installation command should look like this:
library("devtools")
install_github("s5joschi/accessibility")
# once installed just activate the package with
library(AccessibilityMaps)The AccessibilityMaps package will install automatically a set of R-packages in the background, if they are not already installed on your system Those packages are raster, rgdal, gdalUtils, rgrass7 and sp.
We will proceed by setting up the basic working environment for this tutorial and by downloading the input data. The input-data consist of a vector layer with roads in Burkina Faso from the Open Streetmap project. This data was downloaded from the website of Geofabrik. Furthermore we will use a prototype of land-cover data for Africa in 2016 which is provided in 10 meters resolution from the ESA copernicus mission (found here). You do not need to download the data from these websites. Just use the script below to download a copy from this repository.
# create and set working directory
dir.create("accessibility_example")
setwd("accessibility_example/")
# download data and unzip it
download.file(
"https://github.com/s5joschi/accessibility/raw/master/example/input.zip",
destfile = "input.zip"
)
unzip("input.zip")
# create an output folder for the results
dir.create("output")Processing Scheme
Travel Time Maps are created by first producing a friction surface which shows the amount of travel time that is necessary to cross a raster cell in either horizontal or vertical direction. This friction surface depends on the underlying land-cover which is included usually from vector data containing information on the infrastructure and raster data containing information on the landcover. The friction surface can be corrected for the effect of slope on travel velocity. The corrected friction surface is than aggregated in a second step to produce the Travel Time Map. The whole processing logic can be seen in the following image.
The Friction Map
A friction map is a raster layer that shows how much time is required to cross each raster cell in vertical or horizontal space depending on the underlying surface. Cells that contain roads or navigable rivers on flat surfaces have a very low friction i.e. they are quick to be crossed. In contrast a forest or bush-lands in mountains areas have an extremely high friction. The friction map together with a set of travel destinies are the two main ingredients to create travel time maps.
Let’s First load the input-data to create the friction map into R:
# landcover
r_landcover <-
raster(
"input/landuse/esa-copernicus-prototype_0005dergrees.tif"
)
# digital elevation model (DEM)
r_dem<-
raster(
"input/dem/bf_elevation.tif"
)
# administrational boundaries
spodf_admin<-
readOGR("input/admin/",
"gis_osm_places_a_free_1"
)
# OSM roads
spodf_roads<-
readOGR(
"input/roads/",
"gis_osm_roads_free_1_mainroads"
)
# OSM places
spodf_sources<-
readOGR(
"input/admin/",
"gis_osm_places_free_1"
)
# convert popuplation data from OSM to numeric
spodf_sources@data$population<-as.numeric(as.character(spodf_sources@data$population))Let’s plot some of the input-data to create a map with land-cover classes, the road network and nonadministrative boundaries of Burkina Faso:
plot(
r_landcover
)
plot(
spodf_admin,
border="red",
add=T)
plot(
spodf_roads,
col="black",
add=T
)
We will start now to process the layers that will be later used to create the friction map. To prepare our road layer for processing we first need to create a column in the attribute table that contains travel speeds in km/h based on the surface type of the road. If it is a road with the category primary or primary_link (for details see the OSM documentation) we assume it is paved and assign it an average travel speed of 60 km/h. Otherwise we assume that the roads are not paved and assign it a travel speed of 40 km/h in the dry season- and only 10km/h in the rainy season. We use our land-cover map as a base-layer to define the projection system, extent and resolution of the layers to be created. We use the function acc_vec2fric to process the vector input-data.
# 3.1 roads
# create columns with travelspeeds for dry and wetseason
spodf_roads$s_dry <-
ifelse(spodf_roads$fclass %in% c("primary", "primary_link"),
60,
40)
spodf_roads$s_wet <-
ifelse(spodf_roads$fclass %in% c("primary", "primary_link"),
60,
10)
# create friction input layers from roads
r_roads_dry <-
acc_vec2fric(my_input = spodf_roads,
my_baselayer = r_landcover,
my_speedfield = "s_dry")
r_roads_wet <-
acc_vec2fric(my_input = spodf_roads,
my_baselayer = r_landcover,
my_speedfield = "s_wet")
# plot the rasterized roads
plot(r_roads_dry,
col=c("orange","black"))
The next step now is to process the land-cover layer. In order to process the raster input data we use the function acc_ras2fric. Our original land-cover data has 11 classes from 0 to 10 (for the original classes see the documentation on the link above). Our assumptions on travel speeds on the different land surfaces are listed in the code below. We provide the land-cover input classes as a vector with the sequence 0,1,2,3,4,5,6,7,8,9,10 and another vector containing the corresponding travel speeds in the same order.
# reclassification values.
# no data will be reclassified to 0 km/h,
# tree cover and shrublands and similar to 3 km/h
# open landscapes to 10 km/h
# urban areas to 30 km/h
# water to 20 km/h (boat transportation)
r_landcover_reclass <-
acc_ras2fric(
my_input = r_landcover,
my_baselayer = r_landcover,
my_reclass_inputvalues = 0:10,
my_reclass_outputvalues = c(0, 3, 3, 10, 10, 3, 10, 10, 30, 2, 20)
)
# plot the output raster
plot(r_landcover_reclass)
Now we will correct the input layers for the effect that slope exerts on them. Our main assumption is that steep slopes reduce travel speed. In order to do so we first need to convert a Digital Elevation Model (DEM) into a map that contains information about the slope measure in radians. We use the function acc_radians to create this map assuring that it has the same spatial characteristics as all other layers using the land-cover map as a base-layer. Once we have a radians map we can use the function acc_slopecorr to correct the input layers.
# Create radiansmap
r_radians <-
acc_radians(
my_input = r_dem,
my_baselayer = r_landcover,
resampling_method = "near"
)
# correct for slope effects
r_roads_dry_corr <-
acc_slopecorr(my_input = r_roads_dry,
my_radians = r_radians)
r_roads_wet_corr <-
acc_slopecorr(my_input = r_roads_wet,
my_radians = r_radians)
r_landcover_reclass_corr <-
acc_slopecorr(my_input = r_landcover_reclass,
my_radians = r_radians)The final step to create the friction map is to use all pre-processed input layers with the function acc_friction. We create two distinct friction maps: one for the rainy and one for the dry season. We us an output resolution of 500 meters to facilitate the creation of the travel times maps, which demands quite some computational resources. The processing time will ultimately not only depend on your hardware specification but more importantly on the extent of your research area and the spatial resolution.
r_friction_dry <-
acc_friction(list(r_roads_dry_corr, r_landcover_reclass_corr),
my_outputresolution = 500)
r_friction_wet <-
acc_friction(list(r_roads_wet_corr, r_landcover_reclass_corr),
my_outputresolution = 500)
# plot friction map of dry_season
plot(r_friction_dry)
The friction map contains values that define how much time is needed (in seconds) to cross each grid-cell in horizontal or vertical space. Clearly the roads can be identified in white color to be the cells with less friction. The friction map is the main input data for the travel time map together with the travel destinies (sources) which are provided as a vector input layer.
The Travel Time Map
In the last step we will create travel time maps that show the time necessary for any place in Burkina Faso to reach the next city with more than 1 Mio. inhabitants. Once you have the friction maps to your liking you basically can now create travel time maps for any type of source with the code below.
Note, that we use 3 MB of RAM with the max_ram parameter. You can adjust this parameter. It is recommended that you do not use much more than around 50-60 % of currently available RAM to not slow down your computer excessively. Also note that you will have to define here the path to the installed GRASS binaries. This path varies depending on your local Operation System and the version of GRASS installed. Helpful information about the location of GRASS on your system can be found here: https://grasswiki.osgeo.org/wiki/R_statistics/rgrass7#GRASS_within_R. Another easy way to find the binaries is to use the function findGRASS from the the package link2GI. Consider installing it if you are unable to find it manually. If you encounter an error regarding iconv.dll or other configuration errors for GRASS on Windows, this CRAN page might help.
# subset the cities for such that have more than 1 Mio inhabitants.
spodf_sources<-spodf_sources[spodf_sources@data$population >100000, ]
r_accessibility_dry <-
acc_accessibility(
my_friction = r_friction_dry,
my_sources = spodf_sources,
knightsmove = T,
grassbin = "/usr/lib/grass72", # might change depending on OS. See function documentation.
max_ram = 3000
)
r_accessibility_wet <-
acc_accessibility(
my_friction = r_friction_wet,
my_sources = spodf_sources,
knightsmove = T,
grassbin = "/usr/lib/grass72",
max_ram = 3000
)
# convert output (seconds) to minutes
r_accessibility_dry <-
r_accessibility_dry / 60
r_accessibility_wet <-
r_accessibility_wet / 60
# plot both maps
plot(r_accessibility_dry, breaks=0:1200,col = topo.colors(1200),legend=F)
legend("topright", legend = c("< 1hr","5 hrs","10 hrs","15 hrs","20 hrs"), fill = topo.colors(5))
plot(r_accessibility_wet, breaks=0:1200,col = topo.colors(1200),legend=F)
legend("topright", legend = c("< 1hr","5 hrs","10 hrs","15 hrs","20 hrs"), fill = topo.colors(5))
# export rasters
writeRaster(
r_accessibility_dry,
filename = "output/r_accessibility_dry.tif",
datatype = "INT2U",
overwrite = T
)
writeRaster(
r_accessibility_wet,
filename = "output/r_accessibility_wet.tif",
datatype = "INT2U",
overwrite = T
)We can clearly detect differences in wet- and dry-season accessibility for several regions in Burkina Faso which is mainly the result of different road qualities. These differences might have an important impact on agricultural value chains or socio-economic welfare for people living in remote rural areas.
Final Note
Please note that it may make also sense to export the layers created during processing to inspect the results in a proper GIS software. If you use OSM data please note, that the quality varies greatly from country to country. Depending on your research objective you might want to always complement it with data from governmental agencies. Also it makes sense to include additional infrastructure layers such as waterways or railroads which we did not do in this example for the sake of simplicity. We hope that you can make use of this package and if you encounter any errors make sure to file a bug report here or ask a question on Stackoverflow.