From 3152adc5bc4f80b1e15c8e26614b2beda819c1dc Mon Sep 17 00:00:00 2001 From: Clementine Cottineau Date: Thu, 6 Feb 2025 11:46:10 +0100 Subject: [PATCH 1/4] change bbox function from nominatimelite to getbb in replicability example --- episodes/18-import-and-visualise-osm-data.Rmd | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/episodes/18-import-and-visualise-osm-data.Rmd b/episodes/18-import-and-visualise-osm-data.Rmd index 01dbb0a6..68812de1 100644 --- a/episodes/18-import-and-visualise-osm-data.Rmd +++ b/episodes/18-import-and-visualise-osm-data.Rmd @@ -199,8 +199,7 @@ We have produced a proof a concept on Brielle, but can we factorise our work to We might replace the name in the first line and run everything again. Or we can create a function. ```{r reproducibility} extract_buildings <- function(cityname, year=1900){ - nominatim_polygon <- geo_lite_sf(address = cityname, points_only = FALSE) - bb <- st_bbox(nominatim_polygon) + bb <- getbb(cityname) x <- opq(bbox = bb) %>% add_osm_feature(key = 'building') %>% From c941c4e31ab1481a4598d37a776e596c88478c2a Mon Sep 17 00:00:00 2001 From: Clementine Cottineau Date: Thu, 6 Feb 2025 11:51:33 +0100 Subject: [PATCH 2/4] add replicability example in slides --- instructors/4-gis-slides.html | 414 ++++++++++++++++++++-------------- instructors/4-gis-slides.qmd | 45 ++++ 2 files changed, 289 insertions(+), 170 deletions(-) diff --git a/instructors/4-gis-slides.html b/instructors/4-gis-slides.html index 94486636..6acdd9fb 100644 --- a/instructors/4-gis-slides.html +++ b/instructors/4-gis-slides.html @@ -443,7 +443,7 @@

Learning objectives

-

What is OpenStreetMap? 🗺

+

What is OpenStreetMap?

@@ -461,7 +461,7 @@

What is OpenStreetMap? 🗺

-

What is OpenStreetMap? 🗺

+

What is OpenStreetMap?

Anyone can contribute, by:

  • adding information on existing map objects

  • @@ -470,7 +470,7 @@

    What is OpenStreetMap? 🗺

-

OSM 🗾

+

OSM

The OSM information system relies on :

  • geometrical objects (i.e. points, lines, polygons)

  • @@ -486,7 +486,7 @@

    How to extract geospatial data from Open Street Map with R?

    You can also try with Naarden, Geertruidenberg, Gorinchem, Enkhuizen or Dokkum if you prefer.

-

The Bounding Box ⬜

+

The Bounding Box

We first geocode our spatial text search and extract the corresponding bounding box with the getbb() function

library(tidyverse)
@@ -503,7 +503,7 @@ 

The Bounding Box ⬜

-

The Bounding Box 🔲

+

The Bounding Box

A Problem with download? Try:

assign("has_internet_via_proxy", TRUE, environment(curl::has_internet))
@@ -715,12 +715,56 @@ 

Mapping urbanisation in Brielle

+
+

Replicability

+

We have produced a proof a concept on Brielle, but can we factorise our work to be replicable with other small fortified cities? You can use any of the following cities: Naarden, Geertruidenberg, Gorinchem, Enkhuizen or Dokkum.

+

We might replace the name in the first line and run everything again… or we can create a function.

+
+
+

Replicability

+
+
extract_buildings <- function(cityname, year=1900){
+  bb <- getbb(cityname)
+  
+  x <- opq(bbox = bb) %>%
+     add_osm_feature(key = 'building') %>%
+     osmdata_sf()
+     
+  buildings <- x$osm_polygons %>%
+    st_transform(.,crs=28992)
+    
+  start_date <- as.numeric(buildings$start_date)
+  
+  buildings$build_date <- if_else(start_date < year, year, 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") +
+     ggtitle(paste0("Old buildings in ",cityname)) +
+     coord_sf(datum = st_crs(28992))
+}
+
+
+
+

Test on Brielle

+
+
extract_buildings("Brielle, NL")
+ +
+
+
+

Test on Naarden

+
+
extract_buildings("Naarden, NL")
+ +
+
-

⏰ Challenge

+

Challenge

Import an interactive basemap layer under the buildings with Leaflet, using the package documentation.

-
+
20:00
@@ -733,28 +777,24 @@

⏰ Challenge

-
-

⏰ To submit your work in Zoom:

- -
-

⏰ One solution

+

One solution

-
library(leaflet)
-
-buildings2 <- buildings %>%
-  st_transform(.,crs=4326)
-
-leaflet(buildings2) %>%
-  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))
+
library(leaflet)
+
+buildings2 <- buildings %>%
+  st_transform(.,crs=4326)
+
+leaflet(buildings2) %>%
+  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))
-
- +
+
@@ -776,6 +816,30 @@

Part 2. How to perform basic GIS operations with the sf package

Objectives:

+

By the end of this session, you should be able to:
+## Objectives: By the end of this session, you should be able to:

+
    +
  • Perform geoprocessing operations such as union, join, and intersection
    +
  • +
  • Create buffers and centroids
    +
  • +
  • Compute the area of spatial polygons
    +
  • +
  • Calculate density within spatial units
    +
  • +
  • Map the results
    +
  • +
  • Perform geoprocessing operations such as union, join, and intersection
    +
  • +
  • Create buffers and centroids
    +
  • +
  • Compute the area of spatial polygons
    +
  • +
  • Calculate density within spatial units
    +
  • +
  • Map the results
    +## Objectives:
  • +

By the end of this session, you should be able to:

  • Perform unions, joins and intersection operations

  • @@ -785,30 +849,30 @@

    Objectives:

-

the ‘sf’ package 📤

+

the ‘sf’ package

-

the ‘sf’ cheatsheet 🗺

+

the ‘sf’ cheatsheet

-

Conservation in Brielle, NL 🏫

+

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.

If you were using pen&paper or a traditional GIS software, what would your workflow be?

-

Conservation in Brielle, NL 🏫

+

Conservation in Brielle, NL

-
bb <- osmdata::getbb("Brielle, NL")
-x <- opq(bbox = bb) %>%
-    add_osm_feature(key = 'building') %>%
-    osmdata_sf()
-
-buildings <- x$osm_polygons %>%
-  st_transform(.,crs=28992)
-
-summary(buildings$start_date)
+
bb <- osmdata::getbb("Brielle, NL")
+x <- opq(bbox = bb) %>%
+    add_osm_feature(key = 'building') %>%
+    osmdata_sf()
+
+buildings <- x$osm_polygons %>%
+  st_transform(.,crs=28992)
+
+summary(buildings$start_date)
   Length     Class      Mode 
     10702 character character 
@@ -816,26 +880,32 @@

Conservation in Brielle, NL 🏫

-

Conservation in Brielle, NL 🏭

+

Conservation in Brielle, NL

-
old <- 1800 # year prior to which you consider a building old
-
-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") +
-   coord_sf(datum = st_crs(28992)) 
+
old <- 1800 # year prior to which you consider a building old
+
+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") +
+   coord_sf(datum = st_crs(28992)) 

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 100 meters. In GIS terms, we want to create a buffer around polygons. The corresponding function sf is st_buffer, with 2 arguments: : the radius of the buffer.

+ +

Buffer example

+
+

Buffer

+

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: - “x”: the polygons around which to create buffers - “dist”: the radius of the buffer.

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 100 meters. In GIS terms, we want to create a buffer around polygons. The corresponding function sf is st_buffer, with 2 arguments:

    @@ -843,14 +913,14 @@

    Buffer

  • “dist”: the radius of the buffer.
-
+

Buffer

-
 distance <- 100 # in meters 
- 
-# First, we check that the "old_buildings" 
-# layer projection is measured in meters:
-st_crs(old_buildings)
+
 distance <- 100 # in meters 
+ 
+# First, we check that the "old_buildings" 
+# layer projection is measured in meters:
+st_crs(old_buildings)
Coordinate Reference System:
   User input: EPSG:28992 
@@ -896,27 +966,27 @@ 

Buffer

-
+

Buffer

-
#then we use `st_buffer()`
-buffer_old_buildings <- 
-  st_buffer(x = old_buildings, dist = distance)
- 
-ggplot(data = buffer_old_buildings) + 
-  geom_sf() +   
-  coord_sf(datum = st_crs(28992))
+
#then we use `st_buffer()`
+buffer_old_buildings <- 
+  st_buffer(x = old_buildings, dist = distance)
+ 
+ggplot(data = buffer_old_buildings) + 
+  geom_sf() +   
+  coord_sf(datum = st_crs(28992))

Union

Now, we have 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().

+

Instead of overlapping buffers, we want to create continuous conservation zones by merging adjacent areas into a single polygon. This process is called a union, and the corresponding function in sf is st_union(). Union dissolves boundaries between overlapping geometries.

-
single_old_buffer <- st_union(buffer_old_buildings) %>%
-   st_cast(to = "POLYGON") %>%
-   st_as_sf() 
+
single_old_buffer <- st_union(buffer_old_buildings) %>%
+   st_cast(to = "POLYGON") %>%
+   st_as_sf() 
@@ -924,9 +994,9 @@

Union

We also used st_cast() to explicit the type of the resulting object (POLYGON instead of the default MULTIPOLYGON) and st_as_sf() to transform the polygon into an sf object. With this function, we ensure that we end up with an sf object, which was not the case after we forced the union of old buildings into a POLYGON format.

Then we create unique IDs to identify the new polygons.

-
single_old_buffer<- single_old_buffer %>%
-  mutate("ID"=as.factor(1:nrow(single_old_buffer))) %>%
-  st_transform(.,crs=28992) 
+
single_old_buffer<- single_old_buffer %>%
+  mutate("ID"=as.factor(1:nrow(single_old_buffer))) %>%
+  st_transform(.,crs=28992) 
@@ -937,99 +1007,103 @@

Centroids

Centroids

-
sf::sf_use_s2(FALSE)  # s2 works with geographic projections, so to calculate centroids in projected CRS units (meters), we need to disable it
-
-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) +
-  coord_sf(datum = st_crs(28992))
+
sf::sf_use_s2(FALSE)  # s2 works with geographic projections, so to calculate centroids in projected CRS units (meters), we need to disable it
+
+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) +
+  coord_sf(datum = st_crs(28992))

Intersection

-

Now, 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().

+

Now, 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 dissolved buffer polygon contains.

+

This operation means intersecting the polygon layer with the point layer. The corresponding function is st_intersection(). The corresponding function is st_intersection().

Intersection

st_intersection here adds the attributes of the intersected polygon buffers to the data table of the centroids. This means we will now know about each centroid, the ID of its intersected polygon-buffer, and a variable called “n” which is population with 1 for everyone. This means that all centroids will have the same weight when aggregated.

+
+
<img src="[IMAGE_URL_1](https://gisgeography.com/wp-content/uploads/2020/10/Intersect-Table-2048x1170.png)" alt="Point to Polygon" width="45%">
+<img src="https://gisgeography.com/wp-content/uploads/2020/10/Intersect-Tool.png" alt="Polygon to Polygon" width="45%">
+
-
centroids_buffers <- 
-  st_intersection(centroids_old,single_old_buffer) %>%
-  mutate(n = 1)
+
centroids_buffers <- 
+  st_intersection(centroids_old,single_old_buffer) %>%
+  mutate(n = 1)

Intersection

We aggregate them by ID number (group_by(ID)) and sum the variable n to know how many centroids are contained in each polygon-buffer.

-
centroid_by_buffer <- centroids_buffers %>%
-   group_by(ID) %>%
-   summarise(n = sum(n))
- 
-single_buffer <- single_old_buffer %>%
-   mutate(n_buildings = centroid_by_buffer$n)
+
centroid_by_buffer <- centroids_buffers %>%
+   group_by(ID) %>%
+   summarise(n = sum(n))
+ 
+single_buffer <- single_old_buffer %>%
+   mutate(n_buildings = centroid_by_buffer$n)

Mapping the number of buildings

-
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")
+
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")

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

+
+

Calculating area and density of spatial features

Let’s compute the density of old buildings per buffer zone (= number of buildings / area of buffer).

-
single_buffer$area <- sf::st_area(single_buffer)  %>% 
-  units::set_units(., km^2)
-
-single_buffer$old_buildings_per_km2 <- as.numeric(
-  single_buffer$n_buildings / single_buffer$area
-  )
+
single_buffer$area <- sf::st_area(single_buffer)  %>% 
+  units::set_units(., km^2)
+
+single_buffer$old_buildings_per_km2 <- as.numeric(
+  single_buffer$n_buildings / single_buffer$area
+  )

Final Output

Let’s map this layer over the initial map of individual buildings and save the result.

-
p <- 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") 
-
-ggsave(filename = "fig/ConservationBrielle.png", 
-       plot = p)
+
p <- 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") 
+
+ggsave(filename = "fig/ConservationBrielle.png", 
+       plot = p)

Final Output

-
p 
+
p 
-

⏰Challenge: Conservation rules have changed!

+

Challenge: Conservation rules have changed!

The historical threshold now applies to all pre-WWII buildings, but the distance to these building is reduced to 10m.

Can you map the density of old buildings per 10m fused buffer?

-
+
10:00
@@ -1039,65 +1113,65 @@

⏰Challenge: Conservation rules have changed!

One solution

-
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)
-  
-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)
-
-p <- 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") 
+
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)
+  
+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)
+
+p <- 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") 

One solution

-
p
+
p

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.

+

We have seen how to create spatial buffers and centroids, how to merge overlapping geometries using union, how to intersect vector data, and how to retrieve the area of polygons. Additionally, we have explored how to calculate density intersect vector data and how retrieve the area of polygons.

In short:

  • Use the st_* functions from sf for basic GIS operations

  • @@ -1106,7 +1180,7 @@

    Summary and keypoints

-

🦜 A few words of caution

+

A few words of caution

We have taught you to think, type, try, test, read the documentation. This is not only the old fashion way, but the foundation.

When you encounter a bug, a challenge, a question that we have not covered, you could always make use of search engines and AI… but!

Be careful to keep it as a help, tool and support whose quality you can still assess.

@@ -1154,7 +1228,7 @@

Thank you

This workshop is offered on a yearly basis. You are welcome to join the team for future events as helpers and/then as instructors. This strengthens the community… and it can bring you GS Credits.

-

🐰The end…

+

The end…