convert lesson 4 into Rmd format with working R chunks

episodes/basic-gis-operations-with-sf.Rmd

@@ -0,0 +1,226 @@
+---
+title: 'Basic GIS operations with `sf`'
+teaching: 45
+exercises: 25
+---
+
+:::::::::::::::::::::::::::::::::::::: questions 
+
+- How to perform basic GIS operations with the `sf` package?
+
+::::::::::::::::::::::::::::::::::::::::::::::::
+
+::::::::::::::::::::::::::::::::::::: objectives
+
+- Perform unions, joins and intersection operations
+- Compute the area of spatial polygons
+- Create buffers and centroids 
+- Map the results
+
+::::::::::::::::::::::::::::::::::::::::::::::::
+
+## Why `sf` for GIS?
+
+`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 unite polygons, create buffers and centroids, etc. cf. the `sf` [cheatsheet](https://raw.githubusercontent.com/rstudio/cheatsheets/main/sf.pdf).
+
+
+### Conservation in Brielle, NL
+
+Let's focus on old buildings 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 100m around pre-1800 buildings.
+
+Let's select them and see where they are.
+
+
+```{r}
+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")
+ 
+```
+As conservationists, we want to create a zone around historical buildings where building regulation will have special restrictions to preserve historical buildings. 
+
+## Buffers
+
+Let's say this zone should be 100 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}
+ distance <- 100 # 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}
+ single_old_buffer <- st_union(buffer_old_buildings) %>%
+   st_cast(to = "POLYGON") %>%
+  st_as_sf() 
+
+single_old_buffer<- single_old_buffer %>%
+  mutate("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}
+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}
+ centroids_buffers <- 
+  st_intersection(centroids_old,single_old_buffer) %>%
+   mutate(n=1)
+
+ centroid_by_buffer <- centroids_buffers %>%
+   group_by(ID) %>%
+   summarise(n = sum(n))
+ single_buffer <- single_old_buffer %>%
+   mutate(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}
+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") 
+
+```
+
+::::::::::::::::::::::::::::::::::::: challenge 
+
+## Challenge: Conservation rules have changed. 
+
+The historical threshold now applies to all pre-war buildings, but the distance to these building is reduced to 10m. Can you map the number of all buildings per 10m fused buffer?
+
+
+:::::::::::::::::::::::: solution 
+
+```{r}
+ 
+ old <- 1939 
+ distance <- 10
+ #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 %>%
+   mutate("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) %>%
+   mutate(n=1)
+ centroid_by_buffer <- centroids_buffers %>% 
+   group_by(ID) %>%
+   summarise(
+   n = sum(n)
+   )
+ single_buffer <- single_old_buffer %>% 
+   mutate(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}
+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") 
+ 
+```
+
+
+:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: instructor
+
+
+::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+
+
+::::::::::::::::::::::::::::::::::::: keypoints 
+
+- Use the `st_*` functions from `sf` for basic GIS operations
+- Perform unions, joins and intersection operations
+- Compute the area of spatial polygons
+
+::::::::::::::::::::::::::::::::::::::::::::::::
+