Session4b.Rmd

---
source: Rmd
title: "Basic GIS operations with the `sf` package"
author: "Clémentine Cottineau"
teaching: 40
exercises: 20
---


```{r, include=FALSE}
source("bin/chunk-options.R")
knitr_fig_path("04b-")

```

# Question: "How can I use the `sf` package to handle spatial data?"

## Objectives:
- Show how to create spatial buffers and centroids
- Demonstrate how to intersect vector data
- Present the function to retrieve the area of polygons




```{r setup, include=FALSE}
knitr::opts_chunk$set(
	echo = TRUE,
	message = FALSE,
	warning = FALSE
)
```

Packages needed 
```{r lib}
library(dplyr)
library(sf)
library(tibble)
library(ggplot2)
library(remotes)
remotes::install_github('ropensci/osmdata')
library(osmdata)
library(osmextract)
```

---

# Introduction

## What is `sf`?

`sf` is a package which supports simple features (sf), ["a standardized way to
encode spatial vector data."](https://cran.r-project.org/web/packages/sf/sf.pdf).
It contains a large set of functions to achieve all the operations on vector spatial data for which you might use traditional GIS software: change the coordinate system, join layers, intersect or unit polygons, create buffers and centroids, etc. cf. the `sf` [cheatsheet](https://raw.githubusercontent.com/rstudio/cheatsheets/main/sf.pdf).

We are going to go through some of these basics with the case study of Rotterdam buildings.



## Conservation in Rotterdam

Let's focus on really old building and imagine we're in charge of their conservation. We want to know how much of the city would be affected by a non-construction zone of 500m around pre-1800 buildings.

Let's select them and see where they are.

```{r old-buildings}
assign("has_internet_via_proxy", TRUE, environment(curl::has_internet))

bb <- getbb('Rotterdam', format_out = 'sf_polygon')
x <- opq(bbox = bb) %>% 
   add_osm_feature(key = 'building') %>%
    osmdata_sf ()
buildings <- x$osm_polygons %>% st_transform(.,crs=28992)

old <- 1800 # year prior to which you consider a building old

summary(buildings$start_date)
buildings$start_date <- as.numeric(as.character(buildings$start_date))

old_buildings <- buildings |>
  filter(start_date <= old)

 ggplot(data = old_buildings) + geom_sf(colour="red")
```

# Basic GIS operations

As conservationists, we want to create a zone around historical buildings where building regulation will have special restrictions to preserve historical buildings. 

## Buffer

Let's say this zone should be 500 meters. In GIS terms, we want to create a buffer around polygons. The corresponding function `sf` is `st_buffer`, with 2 arguments: the polygons around which to create buffers, and the radius of the buffer.

```{r buffers}
 distance <- 500 # in meters 
 
buffer_old_buildings <- 
  st_buffer(x = old_buildings, dist = distance)
 
 ggplot(data = buffer_old_buildings) + geom_sf()
  
```

## Union

Now, we have a lot of overlapping buffers. We would rather create a unique conservation zone rather than overlapping ones in that case. So we have to fuse the overlapping buffers into one polygon. This operation is called union and the corresponding function is `st_union`.

```{r union}
 single_old_buffer <- st_union(buffer_old_buildings) |> 
   st_cast(to = "POLYGON") |> 
  st_as_sf() 

single_old_buffer<- single_old_buffer |>
  add_column("ID"=as.factor(1:nrow(single_old_buffer)))  |> 
  st_transform(.,crs=28992) 
  
```

We also use `st_cast()` to explicit the type of the resulting object (*POLYGON* instead of the default *MULTIPOLYGON*), `st_as_sf()` to transform the polygon into an `sf` object.


## Centroids
For the sake of visualisation speed, we would like to represent buildings by a single point (for instance: their geometric centre) rather than their actual footprint. This operation means defining their centroid and the corresponding function is `st_centroid`.

```{r centroids}
sf::sf_use_s2(FALSE)
centroids_old <- st_centroid(old_buildings) |> 
  st_transform(.,crs=28992)  

ggplot() + 
    geom_sf(data = single_old_buffer, aes(fill=ID)) +
    geom_sf(data = centroids_old)
```

# Intersection
Now what we would like to distinguish conservation areas based on the number of historic buildings they contain. In GIS terms, we would like to know how many centroids each fused buffer polygon contains. This operation means intersecting the layer of polygons with the layer of points and the corresponding function is `st_intersection`.


```{r intersection}
 centroids_buffers <- st_intersection(centroids_old,single_old_buffer) |> 
   add_column(n=1)
 centroid_by_buffer <- centroids_buffers |> 
   group_by(ID) |>
   summarise(n = sum(n))
 single_buffer <- single_old_buffer |> 
   add_column(n_buildings = centroid_by_buffer$n)

  ggplot() + 
   geom_sf(data = single_buffer, aes(fill=n_buildings)) +
   scale_fill_viridis_c(alpha = 0.8,
                        begin = 0.6,
                        end = 1,
                        direction = -1,
                        option = "B")
```

## Final output:

Let's map this layer over the initial map of individual buildings.

```{r visu-1800}

ggplot() + 
   geom_sf(data = buildings) +
   geom_sf(data = single_buffer, aes(fill=n_buildings), colour = NA) +
   scale_fill_viridis_c(alpha = 0.6,
                        begin = 0.6,
                        end = 1,
                        direction = -1,
                        option = "B") 

```
 
### Exercise
 
## Challenge: Conservation rules have changed. The historical threshold now applies to all pre-war buildings, but the distance to these building is reduced to 100m. Can you map the number of all buildings per 100m fused buffer?

---

## Solution
 
```{r prewar-replication}
 
 old <- 1939 
 distance <- 100
 #select
 old_buildings <- buildings |>
   filter(start_date <= old)
 #buffer
 buffer_old_buildings <- st_buffer(old_buildings, dist = distance)
  #union
 single_old_buffer <- st_union(buffer_old_buildings) |> 
   st_cast(to = "POLYGON") |> 
   st_as_sf()  
 
 single_old_buffer <- single_old_buffer |>
   add_column("ID"=1:nrow(single_old_buffer))  |> 
   st_transform(single_old_buffer,crs=4326) 
 #centroids
 centroids_old <- st_centroid(old_buildings) |> st_transform(.,crs=4326)  
  #intersection
  centroids_buffers <- st_intersection(centroids_old,single_old_buffer) |> 
   add_column(n=1)
 centroid_by_buffer <- centroids_buffers |> 
   group_by(ID) |>
   summarise(
   n = sum(n)
   )
 single_buffer <- single_old_buffer |> 
   add_column(n_buildings = centroid_by_buffer$n)
  ggplot() + 
   geom_sf(data = buildings) +
   geom_sf(data = single_buffer, aes(fill=n_buildings), colour = NA) +
   scale_fill_viridis_c(alpha = 0.6,
                        begin = 0.6,
                        end = 1,
                        direction = -1,
                        option = "B") 
   
```


*Problem: there are many pre-war buildings and the buffers are large so the number of old buildings is not very meaningful. Let's compute the density of old buildings per buffer zone.*

## Area

```{r area}
single_buffer$area <- st_area(single_buffer, )  %>% units::set_units(., km^2)

single_buffer$old_buildings_per_km2 <- as.numeric(single_buffer$n_buildings / single_buffer$area)

 ggplot() + 
   geom_sf(data = buildings) +
   geom_sf(data = single_buffer, aes(fill=old_buildings_per_km2), colour = NA) +
   scale_fill_viridis_c(alpha = 0.6,
                        begin = 0.6,
                        end = 1,
                        direction = -1,
                        option = "B") 
 
```


---


# Summary and Keypoints:
We have seen how to create spatial buffers and centroids, how to intersect vector data and how retrieve the area of polygons.

In short:

- Use the `sf` package to treat geospatial data
- Use the `st_*` functions for basic GIS operations
- Use the `ggplot` package to map the results