Merge pull request #2 from ClementineCttn/main

Lesson 4 into Carpentries formatting

config.yaml

@@ -14,7 +14,7 @@ carpentry: 'incubator'
 title: 'Geospatial Data Carpentry for Urbanism'
 
 # Date the lesson was created (YYYY-MM-DD, this is empty by default)
-created:
+created: 
 
 # Comma-separated list of keywords for the lesson
 keywords: 'software, data, lesson, The Carpentries, geospatial'
@@ -33,7 +33,7 @@ source: 'https://github.com/carpentries-incubator/r-geospatial-urban'
 branch: 'main'
 
 # Who to contact if there are any issues
-contact: 'c.forgaci@tudelft.nl'
+contact: 'rbanism@tudelft.nl'
 
 # Navigation ------------------------------------------------
 #
@@ -58,6 +58,7 @@ contact: '[email protected]'
 # - another-learner.md
 
 # Order of episodes in your lesson
+
 episodes:
 - 01-intro-to-r.Rmd
 - 02-data-structures.Rmd
@@ -73,6 +74,8 @@ episodes:
 - 15-reproject-raster-data.Rmd
 - 16-raster-calculations.Rmd
 - 17-work-with-multi-band-rasters.Rmd
+- 18-import-and-visualise-osm-data.Rmd
+- 19-basic-gis-with-r-sf.Rmd
 
 # Information for Learners
 learners:

episodes/18-import-and-visualise-osm-data.Rmd

@@ -0,0 +1,241 @@
+---
+title: 'Import and Visualise OSM Data'
+teaching: 45
+exercises: 25
+---
+
+:::::::::::::::::::::::::::::::::::::: questions 
+
+- How to import and work with vector data from OpenStreetMap?
+
+::::::::::::::::::::::::::::::::::::::::::::::::
+
+::::::::::::::::::::::::::::::::::::: objectives
+
+After completing this episode, participants should be able to…
+
+- Import OSM vector data from the API
+- Select and manipulate OSM vector data
+- Visualise and map OSM Vector data
+- Use Leaflet for interactive mapping
+
+::::::::::::::::::::::::::::::::::::::::::::::::
+
+```{r setup, include=FALSE} 
+knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
+```
+
+## What is OpenStreetMap?
+
+OpenStreetMap (OSM) is a collaborative project which aims at mapping the world and sharing geospatial data in an open way. Anyone can contribute, by mapping geographical objects they encounter, by adding topical information on existing map objects (their name, function, capacity, etc.), or by mapping buildings and roads from satellite imagery (cf. [HOT: Humanitarian OpenStreetMap Team](https://www.hotosm.org/)).
+
+This information is then validated by other users and eventually added to the common "map" or information system. This ensures that the information is accessible, open, verified, accurate and up-to-date.
+
+The result looks like this:
+![](fig/OSM1.png)
+
+The geospatial data underlying this interface is made of geometrical objects (i.e. points, lines, polygons) and their associated tags (#building #height, #road #secondary #90kph, etc.).
+
+## How to extract geospatial data from OpenStreetMap?
+
+```{r message=FALSE}
+library(tidyverse)
+library(sf)
+assign("has_internet_via_proxy", TRUE, environment(curl::has_internet))
+
+```
+
+
+### Bounding box
+
+The first thing to do is to define the area within which you want to retrieve data, aka the *bounding box*. This can be defined easily using a place name and the package `nominatimlite` to access the free Nominatim API provided by OpenStreetMap. 
+
+We are going to look at *Brielle* together, but you can also work with the small cities of *Naarden*, *Geertruidenberg*, *Gorinchem*, *Enkhuizen* or *Dokkum*.
+
+
+
+
+
+We first geocode our spatial text search and extract the corresponding polygon (`geo_lite_sf`) and then extract its bounding box (`st_bbox`).
+
+```{r}
+#install.packages("nominatimlite")
+library(nominatimlite)
+
+nominatim_polygon <- geo_lite_sf(address = "Brielle", points_only = FALSE)
+bb <- st_bbox(nominatim_polygon)
+bb
+```
+
+::::::::::::::::::::::::::::::::::::: callout
+
+If you encounter an error linked to your internet proxy ("Error: Overpass query unavailable without internet R"), run this line of code. It might not be needed, but ensures that your machine knows it has internet.
+
+```{r}
+assign("has_internet_via_proxy", TRUE, environment(curl::has_internet))
+```
+
+::::::::::::::::::::::::::::::::::::::::::::::::
+
+
+### A word of caution
+
+There might be multiple responses from the API query, corresponding to different objects at the same location, or different objects at different locations.
+For example: Brielle (Netherlands) and Brielle (New Jersey)
+
+![Brielle, Netherlands](fig/Brielle_NL.jpeg){width=40%}
+
+![Brielle, New Jersey](fig/Brielle_NJ.jpeg "Brielle, New Jersey"){width=40%}
+
+By default, `geo_lite_sf()` from the `nominatimlite` package returns the first item. This means that regardless of the number of returned locations with the given name, the function will return a bounding box and it might be that we are not looking for the first item. We should therefore try to be as unambiguous as possible by adding a country code or district name.
+
+```{r}
+nominatim_polygon <- geo_lite_sf(address = "Brielle, NL", points_only = FALSE)
+bb <- st_bbox(nominatim_polygon)
+bb
+```
+
+
+
+:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: instructor
+
+If this does not work for some reason, the `nominatim_polygon` can be found in the data folder: "episodes/data/boundingboxBrielle.shp".
+
+::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+
+
+## Extracting features
+
+A [feature](https://wiki.openstreetmap.org/wiki/Map_features) in the OSM language is a category or tag of a geospatial object. Features are described by general keys (e.g. "building", "boundary", "landuse", "highway"), themselves decomposed into sub-categories (values) such as "farm", "hotel" or "house" for `buildings`, "motorway", "secondary" and "residential" for `highway`. This determines how they are represented on the map.
+
+
+### Searching documentation
+
+Let's say we want to download data from OpenStreetMap and we know there is a package for it named `osmdata`, but we don't know which function to use and what arguments are needed. Where should we start?
+
+Let's check the documentation [online](https://docs.ropensci.org/osmdata/):
+
+![The OSMdata Documentation page](fig/osmdata.png){width=80%}
+
+It appears that there is a function to extract features, using the Overpass API. This function is `opq()` (for OverPassQuery) which, in combination with `add_osm_feature()`, seems to do the job. However, it might not be crystal clear how to apply it to our case. Let's click on the function name in the documentation to find out more.
+
+![The Overpass Query Documentation page](fig/opq.png){width=80%}
+
+
+
+On this page we can read about the arguments needed for each function: a bounding box for `opq()` and some `key` and `value` for `add_osm_feature()`. Thanks to the examples provided, we can assume that these keys and values correspond to different levels of tags from the OSM classification. In our case, we will keep it at the first level of classification, with "buildings" as `key`, and no value. We also see from the examples that another function is needed when working with the `sf` package: `osmdata_sf()`. This ensures that the type of object is suited for `sf`. With these tips and examples, we can write our feature extraction function as follows:
+
+
+```{r}
+#install.packages("osmdata")
+library(osmdata)
+
+
+x <- opq(bbox = bb) %>%
+   add_osm_feature(key = 'building') %>%
+   osmdata_sf()
+```
+
+
+
+### Structure of objects
+
+What is this `x` object made of? It is a data frame of all the buildings contained in the bounding box, which gives us their OSM id, their geometry and a range of attributes, such as their name, building material, building date, etc. The completion level of this data frame depends on user contributions and open resources. Here, for instance, the national [BAG](https://data.overheid.nl/dataset/10491-bag) dataset was used and it is quite complete, but that is different in other countries.
+
+
+
+```{r}
+str(x$osm_polygons)
+```
+
+
+
+## Mapping 
+
+
+Let's map the building age of post-1900 Brielle buildings.
+
+
+### Projections
+
+First, we are going to select the polygons and reproject them with the Amersfoort/RD New projection, suited for maps centred on the Netherlands. This code for this projection is: 28992.
+
+```{r}
+buildings <- x$osm_polygons %>%
+  st_transform(.,crs=28992)
+```
+
+:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: instructor
+
+If this does not work for some reason, the `buildings` can be found in the data folder: "episodes/data/dataBrielle.shp".
+
+::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+
+
+### Visualisation
+
+Then we create a new variable using the threshold at 1900. Every date before 1900 will be recoded as `1900`, so that buildings older than 1900 will be represented with the same shade.
+
+Then we use the `ggplot()` function to visualise the buildings by age. The specific function to represent information as a map is `geom_sf()`. The rest works like other graphs and visualisation, with `aes()` for the aesthetics.
+
+```{r}
+start_date <- as.numeric(buildings$start_date)
+buildings$build_date <- if_else(start_date < 1900, 1900, start_date)
+
+ ggplot(data = buildings) +
+   geom_sf(aes(fill = build_date, colour=build_date))  +
+   scale_fill_viridis_c(option = "viridis")+
+   scale_colour_viridis_c(option = "viridis")
+```
+
+So this reveals the historical centre of Brielle (or the city you chose) and the various urban extensions through time.
+Anything odd? What? Around the centre? Why these limits / isolated points?
+
+
+
+::::::::::::::::::::::::::::::::::::: challenge 
+
+## Challenge: import an interactive basemap layer under the buildings with `Leaflet` (20min)
+
+- Check out the [leaflet package documentation](https://rstudio.github.io/leaflet/)
+- Plot a basemap in Leaflet and try different tiles in the [basemap documentation](https://rstudio.github.io/leaflet/basemaps.html)
+- Transform the buildings into WGS84 projection and add them to the basemap layer with the `addPolygons()` function.
+- Have the `fillColor` of these polygons represent the `build_date` variable. See the [choropleth documentation](https://rstudio.github.io/leaflet/choropleths.html) for use of colors. Tip: use the examples given in the documentation and replace the variable names where needed.
+
+:::::::::::::::::::::::: solution 
+
+## One solution
+
+```{r}
+#install.packages("leaflet")
+library(leaflet)
+
+
+buildings2 <- buildings %>%
+  st_transform(.,crs=4326)
+
+leaflet(buildings2) %>%
+# addTiles()
+  addProviderTiles(providers$CartoDB.Positron) %>%
+  addPolygons(color = "#444444", weight = 0.1, smoothFactor = 0.5,
+  opacity = 0.2, fillOpacity = 0.8,
+  fillColor = ~colorQuantile("YlGnBu", -build_date)(-build_date),
+  highlightOptions = highlightOptions(color = "white", weight = 2,
+    bringToFront = TRUE))
+```
+
+:::::::::::::::::::::::::::::::::
+
+::::::::::::::::::::::::::::::::::::: 
+
+::::::::::::::::::::::::::::::::::::: keypoints 
+
+- Use the `Nominatim` and `Overpass` APIs within R
+- Use the `osmdata` and `nominatimlite` packages to retrieve geospatial data
+- Select features and attributes among OSM tags
+- Use the `ggplot`, `sf` and `leaflet` packages to map data
+
+::::::::::::::::::::::::::::::::::::::::::::::::
+