creating_gridded_accessibility_raster.R

# this script provides 1km gridded data on how accessible the grid square is by the proportion of 100m grid cells within it are accessible by footpath.

#rationale:
#OSM data is really rich and detailed but sometimes you need gridded summaries in raster format. This example shows to generate a raster with the proportional coverage within 1km grid squares containing footpaths. And a second example with the save proportional coverage but for car parks as an indication of how much parking there is in a 1km grid cell.
#These are metrics/indices and act as coarser approximations

#load packages
library(osmextract)
library(magrittr) # for %>%

#for spatial operations
library(sf)
library(terra)

#set the extent of the area to get data from (use a geofabrik region name)
#download_area <- "Great Britain"
download_area <- "North Yorkshire" # use a smaller area for testing

#EXAMPLE 1: AREAS ACCESSIBLE BY FOOT
#Access lines
vt_opts_1 = c(
  "-select", "osm_id, highway, designation, footway, sidewalk",
  "-where", "highway IN ('footway', 'path', 'residential','unclassified','tertiary','sidewalk') OR designation IN ('public_footpath','byway_open_to_all_traffic','restricted_byway','public_bridleway','access_land') OR footway = 'sidewalk' OR sidewalk IN ('both','left','right')"
)

#download the data for the area and create a file called access_lines.gpkg which contains all the layers
oe_get(
  download_area,
  layer = "lines", #have to specify here if it's lines or polygons
  provider = "geofabrik",
  match_by = "name",
  max_string_dist = 1,
  level = NULL,
  force_download = F,
  vectortranslate_options = vt_opts_1,
  extra_tags = c("designation","footway","sidewalk"),
  force_vectortranslate = T,
  quiet = FALSE
) %>% st_write(file.path("access_lines.gpkg"),delete_layer = T) #save file to local directory

#EXAMPLE2: CAR PARKS
#Access polygon
vt_opts_2 = c(
  "-select", "osm_id, amenity",
  "-where", "amenity = 'parking'"
)

#download the data for the area and create a file called access_lines.gpkg which contains all the layers
oe_get(
  download_area,
  layer = "multipolygons", #have to specify here if it's lines or polygons
  provider = "geofabrik",
  match_by = "name",
  max_string_dist = 1,
  level = NULL,
  force_download = F,
  vectortranslate_options = vt_opts_2,
  extra_tags = c("amenity"),
  force_vectortranslate = T,
  quiet = FALSE
) %>% st_write(file.path("parking_areas.gpkg"),delete_layer = T) #save file to local directory

#load in the layer back in
#choose whether to do footpaths or parking here by commenting out the line you don't want
access_layer <- st_read(file.path("access_lines.gpkg")) #FOOTPATHS
access_layer <- st_read(file.path("parking_areas.gpkg")) # CAR PARKs

#transform to OSGB and convert to spatvect for using with terra
access_layer <- st_transform(access_layer,crs = 27700)
access_layer <- vect(access_layer)

#define the gridded area (raster)
x <- rast(resolution = 100, # in metres 1000 = 1km grid
     nlyrs=1, 
     xmin=432000, 
     xmax=452000, 
     ymin=471000, 
     ymax=491000,
     vals =1) #set all values to 1 to start

#set crs of grid
crs(x) <- "epsg:27700"

#view the objects make sure things are looking good
print(x)
print(access_layer)

#crop the spatvect to the extent of the raster
access_layer <- crop(access_layer,x)

# for each 1km set to 0 if no access features present
grid <- rasterize(access_layer,x,
                 background = 0,
                 touches = T) #set the values without target feature to 0

# aggregate to 1km grid squares
grid <- aggregate(grid,fact = 10,fun= "mean")

#visualise
plot(grid)
lines(access_layer) # can take a while to add the lines if plotting for a large area

#export raster
writeRaster(grid,"output_raster.grd")