```
library(sfnetworks)
library(sf)
library(tidygraph)
library(tidyverse)
library(osmdata)
library(slopes)
library(terra)
library(units)
```

# Analyzing real-world geospatial networks in R

Geospatial networks are graphs embedded in geographical space. Space is relevant in their analysis, since the graph topology alone does not contain all information. To facilitate such an integrated analysis, we created the R package sfnetworks, which combines the power of the sf package for spatial data science and the tidygraph package for standard graph analysis. At the UseR! Conference in Salzburg, we presented a simple example of how sfnetworks in combination with other R packages can be used to analyze a real-world spatial network. This notebook contains the code belonging to that talk.

## Analysis

The analysis aims to find out how many restaurants we can reach from the conference venue, if we would travel by bicycle. It consists of five steps. First, we obtain the centerlines of streets and the locations of restaurants in Salzburg from OpenStreetMap, using the osmdata package. Then, we estimate a cycling travel time for each street segment, taking into account the gradient which we compute with the slopes package. We also assign each street segment a cycling suitability level, describing how well it is designed for cycling. The street centerlines are then converted into a routable network structure with the sfnetworks package. Finally, we compute a simple cumulative opportunities metric to quantify the accessibility of restaurants from the conference venue.

### Prepare

We will use the following packages in our analysis. To be clear which function comes from which package, we will explicitly use namespaces in all code chunks.

### 1. Get OpenStreetMap data

In this step, we use the osmdata package to retrieve data from OpenStreetMap. This package basically provides an R interface to the Overpass API of OpenStreetMap. First, we set the conference venue and draw a circular buffer around it with a radius of 3km, which will serve as our area of interest.

If your queries are larger or more frequent, you should probably not overload the Overpass API. Instead, you can download an OpenStreetMap extract (e.g., from Geofabrik) and read the required data into R using the osmextract package.

```
= sf::st_sfc(sf::st_point(c(13.04205, 47.81415)), crs = 4326)
venue = sf::st_buffer(venue, units::set_units(3, "km")) area
```

Now we query all streets within our area of interest from OpenStreetMap. In OpenStreetMap, streets are those features that have a value for the *highway* key. In R, they will be represented by an sf data frame with linestring geometries. Since some streets may be drawn as a closed polygon (like a square or a roundabout) we convert polygon borders into linestrings, before extracting only the linestrings from the API response. Finally, we cut streets at the edge of our area of interest and make sure all of them are single linestrings, meaning that the different parts of multi-linestrings will be seen as different street segments (note that to do this properly we first need to cast everything to multi-linestrings).

```
= osmdata::opq(sf::st_bbox(area)) |>
streets ::add_osm_feature(key = "highway") |>
osmdata::osmdata_sf() |>
osmdata::osm_poly2line() |>
osmdata::pluck("osm_lines") |>
purrr::st_intersection(area) |>
sf::st_cast("MULTILINESTRING") |>
sf::st_cast("LINESTRING") sf
```

In a similar way, we also query all restaurant locations within our area of interest. They will serve later on as the destinations in our accessibility analysis. In OpenStreetMap, restaurants can be identified as those features having the value *restaurant* for the *amenity* key. They are represented as points instead of linestrings.

```
= osmdata::opq(sf::st_bbox(area)) |>
restaurants ::add_osm_feature(key = "amenity", value = "restaurant") |>
osmdata::osmdata_sf() |>
osmdata::pluck("osm_points") |>
purrr::st_intersection(area) sf
```

```
par(mar = c(0, 0, 0, 0))
plot(sf::st_geometry(streets))
plot(sf::st_geometry(restaurants), col = "orange", pch = 20, cex = 0.5, add = TRUE)
```