From 28a0b646a5cf2e5da4ef873ae963160c5efcfcd5 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 6 Nov 2021 23:03:25 +0000 Subject: [PATCH 001/104] Split-out pkg loading and reading for #669 --- 04-spatial-operations.Rmd | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index b28a364b8..509aad0a4 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -9,10 +9,16 @@ library(sf) library(terra) library(dplyr) library(spData) +``` + +- You also need to read in a couple of datasets as follows for Section \@ref(spatial-ras) + +```{r 04-spatial-operations-1-1} elev = rast(system.file("raster/elev.tif", package = "spData")) grain = rast(system.file("raster/grain.tif", package = "spData")) ``` + ## Introduction Spatial operations are a vital part of geocomputation\index{geocomputation}. From 6c5d3c3c3b0a6958f2fce3615e6c7a3298396f46 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 8 Nov 2021 11:55:44 +0000 Subject: [PATCH 002/104] Refactor index.Rmd text about bookdown, update links --- index.Rmd | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/index.Rmd b/index.Rmd index 1b75df6a4..47dcffba2 100644 --- a/index.Rmd +++ b/index.Rmd @@ -34,12 +34,11 @@ This is the online home of *Geocomputation with R*, a book on geographic data an The geocompr book cover **Note**: The first edition of the book has been published by CRC Press in the [R Series](https://www.routledge.com/Chapman--HallCRC-The-R-Series/book-series/CRCTHERSER). -You can buy the book from [CRC Press](https://www.routledge.com/9781138304512), or [Amazon](https://www.amazon.com/Geocomputation-R-Robin-Lovelace-dp-0367670577/dp/0367670577/), and see the archived first edition on the open source book platform [bookdown.org](https://bookdown.org/robinlovelace/geocompr/). +You can buy the book from [CRC Press](https://www.routledge.com/9781138304512), or [Amazon](https://www.amazon.com/Geocomputation-R-Robin-Lovelace-dp-0367670577/dp/0367670577/), and see the archived **First Edition** hosted on [bookdown.org](https://bookdown.org/robinlovelace/geocompr/). -Inspired by [**bookdown**](https://github.com/rstudio/bookdown) and the Free and Open Source Software for Geospatial ([FOSS4G](https://foss4g.org/)) movement, this book is open source. -This ensures its contents are reproducible and publicly accessible for people worldwide. - -The online version of the book is hosted at [geocompr.robinlovelace.net](https://geocompr.robinlovelace.net) and kept up-to-date by [GitHub Actions](https://github.com/Robinlovelace/geocompr/actions), which provides information on its 'build status' as follows: +Inspired by [**bookdown**](https://bookdown.org/) and the Free and Open Source Software for Geospatial ([FOSS4G](https://foss4g.org/)) movement, the code and prose underlying this book is open source, ensuring it's reproducible, accessible and modifiable (e.g. in case you find an inevitable typo) for the benefit of people worldwide. +The online version of the book is hosted at [geocompr.robinlovelace.net](https://geocompr.robinlovelace.net) and kept up-to-date by [GitHub Actions](https://github.com/Robinlovelace/geocompr/actions). +Its current 'build status' as follows: [![Actions](https://github.com/Robinlovelace/geocompr/workflows/Render/badge.svg)](https://github.com/Robinlovelace/geocompr/actions) ``` From 4f98482659f1330a4c8ae594483f0f289852bd81 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 9 Nov 2021 08:25:45 +0000 Subject: [PATCH 003/104] Update ex location --- 03-attribute-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/03-attribute-operations.Rmd b/03-attribute-operations.Rmd index b696c0155..d129f9235 100644 --- a/03-attribute-operations.Rmd +++ b/03-attribute-operations.Rmd @@ -764,6 +764,6 @@ This, however, may lead to a restricted usability of packages depending on the d ## Exercises ```{r, echo=FALSE, results='asis'} -res = knitr::knit_child('_05-ex.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) +res = knitr::knit_child('_03-ex.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) cat(res, sep = '\n') ``` From d1f809763515ef54df48589349acc7bf2562efe4 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 9 Nov 2021 08:26:16 +0000 Subject: [PATCH 004/104] Update ex location --- 03-attribute-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/03-attribute-operations.Rmd b/03-attribute-operations.Rmd index b696c0155..d129f9235 100644 --- a/03-attribute-operations.Rmd +++ b/03-attribute-operations.Rmd @@ -764,6 +764,6 @@ This, however, may lead to a restricted usability of packages depending on the d ## Exercises ```{r, echo=FALSE, results='asis'} -res = knitr::knit_child('_05-ex.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) +res = knitr::knit_child('_03-ex.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) cat(res, sep = '\n') ``` From 400daec0d4eeee380dcf208afdedb6410a248624 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 11 Nov 2021 18:04:09 +0000 Subject: [PATCH 005/104] Define spatial ops upfront --- 04-spatial-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 509aad0a4..739a23c49 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -21,7 +21,7 @@ grain = rast(system.file("raster/grain.tif", package = "spData")) ## Introduction -Spatial operations are a vital part of geocomputation\index{geocomputation}. +Spatial operations, including spatial joins between vector datasets and local and focal operations on raster datasets, are a vital part of geocomputation\index{geocomputation}. This chapter shows how spatial objects can be modified in a multitude of ways based on their location and shape. The content builds on the previous chapter because many spatial operations have a non-spatial (attribute) equivalent. This is especially true for *vector* operations: Section \@ref(vector-attribute-manipulation) on vector attribute manipulation provides the basis for understanding its spatial counterpart, namely spatial subsetting (covered in Section \@ref(spatial-subsetting)). From 0e4672c73dd5aa21cac987615613ca25ec150eef Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 11 Nov 2021 18:41:39 +0000 Subject: [PATCH 006/104] Update prose in Section 4.1 --- 04-spatial-operations.Rmd | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 739a23c49..40ef57f1b 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -23,27 +23,25 @@ grain = rast(system.file("raster/grain.tif", package = "spData")) Spatial operations, including spatial joins between vector datasets and local and focal operations on raster datasets, are a vital part of geocomputation\index{geocomputation}. This chapter shows how spatial objects can be modified in a multitude of ways based on their location and shape. -The content builds on the previous chapter because many spatial operations have a non-spatial (attribute) equivalent. +Many spatial operations have a non-spatial (attribute) equivalent, so concepts such as subsetting and joining datasets demonstrated in the previous chapter are applicable here. This is especially true for *vector* operations: Section \@ref(vector-attribute-manipulation) on vector attribute manipulation provides the basis for understanding its spatial counterpart, namely spatial subsetting (covered in Section \@ref(spatial-subsetting)). Spatial joining (Section \@ref(spatial-joining)) and aggregation (Section \@ref(spatial-aggr)) also have non-spatial counterparts, covered in the previous chapter. Spatial operations differ from non-spatial operations in some ways, however. To illustrate the point, imagine you are researching road safety. Spatial joins can be used to find road speed limits related with administrative zones, even when no zone ID is provided. -But this raises the question: should the road completely fall inside a zone for its values to be joined? -Or is simply crossing or being within a certain distance sufficient? -When posing such questions, it becomes apparent that spatial operations differ substantially from attribute operations on data frames: +But this raises the question: should the road completely fall inside a zone for its values to be joined, or is simply crossing or being within a certain distance sufficient? +When posing such questions, it becomes apparent that spatial operations differ substantially from attribute operations on data frames and that spatial operations are often more complex: the *type* of spatial relationship between objects must be considered. -These are covered in Section \@ref(topological-relations), on topological relations. - +Various types of spatial operations are covered in Section \@ref(topological-relations), on topological relations between vector features. \index{spatial operations} Another unique aspect of spatial objects is distance. -All spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship. These are covered in Section \@ref(distance-relations). +All spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship, as described in the context of vector data in Section \@ref(distance-relations). -Spatial operations also apply to raster objects. -Spatial subsetting of raster objects is covered in Section \@ref(spatial-raster-subsetting); merging several raster 'tiles' into a single object is covered in Section \@ref(merging-rasters). -For many applications, the most important spatial operation on raster objects is *map algebra*, as we will see in Sections \@ref(map-algebra) to \@ref(global-operations-and-distances). -Map algebra is also the prerequisite for distance calculations on rasters, a technique which is covered in Section \@ref(global-operations-and-distances). +Spatial operations on raster objects include subsetting --- covered in Section \@ref(spatial-raster-subsetting) --- and merging several raster 'tiles' into a single object, as demonstrated in Section \@ref(merging-rasters). +*Map algebra* cover a range of operations that modify raster cell values, with or without reference to surrounding cell values, is is vital for many applications. +The concept of map algebra is introduced in Section \@ref(map-algebra); local, focal and zonal map algebra operations are covered in sections \@ref(local-operations), \@ref(focal-operations), and \@ref(zonal-operations), respectively. Global map algebra operations, which generate summary statistics representing an entire raster dataset, and distance calculations on rasters, are discussed in Section \@ref(global-operations-and-distances). +In the final section before the exercises (\@ref(merging-rasets)) the process of merging two raster datasets is discussed and demonstrated with reference to a reproducible example. ```{block2 04-spatial-operations-2, type='rmdnote'} It is important to note that spatial operations that use two spatial objects rely on both objects having the same coordinate reference system, a topic that was introduced in Section \@ref(crs-intro) and which will be covered in more depth in Chapter \@ref(reproj-geo-data). @@ -607,7 +605,7 @@ The next subsection explores these and related operations in more detail. ### Map algebra \index{map algebra} -Map algebra makes raster processing really fast. +Map algebra operations tend to be fast. This is because raster datasets only implicitly store coordinates. To derive the coordinate of a specific cell, we have to calculate it using its matrix position and the raster resolution and origin. For the processing, however, the geographic position of a cell is barely relevant as long as we make sure that the cell position is still the same after the processing. From 17990b1b8bcc3e47b0bf8e07d3aff3839e4a93f4 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 08:21:02 +0000 Subject: [PATCH 007/104] Define map algebra --- 04-spatial-operations.Rmd | 3 +++ 1 file changed, 3 insertions(+) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 40ef57f1b..a5d41ae27 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -605,6 +605,9 @@ The next subsection explores these and related operations in more detail. ### Map algebra \index{map algebra} +The term 'map algebra' was coined in the late 1970s to describe a "set of conventions, capabilities, and techniques" for the analysis of geographic data [@tomlin_1994_map]. + + Map algebra operations tend to be fast. This is because raster datasets only implicitly store coordinates. To derive the coordinate of a specific cell, we have to calculate it using its matrix position and the raster resolution and origin. From ac9bfa3e75df6049bb5b97bd7d3ea7b87c4f4696 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 08:59:03 +0000 Subject: [PATCH 008/104] Update map algebra section --- 04-spatial-operations.Rmd | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index a5d41ae27..ab433c076 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -39,7 +39,7 @@ Another unique aspect of spatial objects is distance. All spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship, as described in the context of vector data in Section \@ref(distance-relations). Spatial operations on raster objects include subsetting --- covered in Section \@ref(spatial-raster-subsetting) --- and merging several raster 'tiles' into a single object, as demonstrated in Section \@ref(merging-rasters). -*Map algebra* cover a range of operations that modify raster cell values, with or without reference to surrounding cell values, is is vital for many applications. +*Map algebra* covers a range of operations that modify raster cell values, with or without reference to surrounding cell values, is is vital for many applications. The concept of map algebra is introduced in Section \@ref(map-algebra); local, focal and zonal map algebra operations are covered in sections \@ref(local-operations), \@ref(focal-operations), and \@ref(zonal-operations), respectively. Global map algebra operations, which generate summary statistics representing an entire raster dataset, and distance calculations on rasters, are discussed in Section \@ref(global-operations-and-distances). In the final section before the exercises (\@ref(merging-rasets)) the process of merging two raster datasets is discussed and demonstrated with reference to a reproducible example. @@ -605,20 +605,21 @@ The next subsection explores these and related operations in more detail. ### Map algebra \index{map algebra} -The term 'map algebra' was coined in the late 1970s to describe a "set of conventions, capabilities, and techniques" for the analysis of geographic data [@tomlin_1994_map]. +The term 'map algebra' was coined in the late 1970s to describe a "set of conventions, capabilities, and techniques" for the analysis of geographic raster *and* (although less prominently) vector data [@tomlin_1994_map]. +Although the concept never became widely adopted, the term usefully encapsulates and helps classify the range operations that can be undertaken on raster datasets. +In this context we define map algebra more narrowly, as operations that modify or summarise raster cell values, with reference to surrounding cells, zones, or statistical functions that apply to every cell. - -Map algebra operations tend to be fast. -This is because raster datasets only implicitly store coordinates. -To derive the coordinate of a specific cell, we have to calculate it using its matrix position and the raster resolution and origin. +Map algebra operations tend to be fast, because raster datasets only implicitly store coordinates (hence the oversimplifying phrase "raster is faster but vector is corrector"). +The location of cells in raster datasets can be calculated it using its matrix position and the resolution and origin of the dataset (stored in the header). For the processing, however, the geographic position of a cell is barely relevant as long as we make sure that the cell position is still the same after the processing. Additionally, if two or more raster datasets share the same extent, projection and resolution, one could treat them as matrices for the processing. -This is exactly what map algebra is doing in R. -First, the **terra** package checks the headers of the rasters on which to perform any algebraic operation, and only if they are correspondent to each other, the processing goes on. -And secondly, map algebra retains the so-called one-to-one locational correspondence. -This is where it substantially differs from matrix algebra which changes positions when for example multiplying or dividing matrices. -Map algebra (or cartographic modeling) divides raster operations into four subclasses [@tomlin_geographic_1990], with each working on one or several grids simultaneously: +This is the way that map algebra works with the **terra** package. +First, the headers of the raster datasets are queried and (in cases where map algebra operations work on more than one dataset) checked to ensure the datasets are compatible. +Second, map algebra retains the so-called one-to-one locational correspondence, meaning that cells cannot move. +This differs from matrix algebra, in which values change position, for example when multiplying or dividing matrices. + +Map algebra (or cartographic modeling with raster data) divides raster operations into four subclasses [@tomlin_geographic_1990], with each working on one or several grids simultaneously: 1. *Local* or per-cell operations 2. *Focal* or neighborhood operations. From aadc73e27eacfd02467655bce9be7aa8b3964b33 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 20:28:19 +0000 Subject: [PATCH 009/104] Shorten and simplify description of spatial vs attribute operations --- 04-spatial-operations.Rmd | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index ab433c076..af0c2fca4 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -27,14 +27,10 @@ Many spatial operations have a non-spatial (attribute) equivalent, so concepts s This is especially true for *vector* operations: Section \@ref(vector-attribute-manipulation) on vector attribute manipulation provides the basis for understanding its spatial counterpart, namely spatial subsetting (covered in Section \@ref(spatial-subsetting)). Spatial joining (Section \@ref(spatial-joining)) and aggregation (Section \@ref(spatial-aggr)) also have non-spatial counterparts, covered in the previous chapter. -Spatial operations differ from non-spatial operations in some ways, however. -To illustrate the point, imagine you are researching road safety. -Spatial joins can be used to find road speed limits related with administrative zones, even when no zone ID is provided. -But this raises the question: should the road completely fall inside a zone for its values to be joined, or is simply crossing or being within a certain distance sufficient? -When posing such questions, it becomes apparent that spatial operations differ substantially from attribute operations on data frames and that spatial operations are often more complex: -the *type* of spatial relationship between objects must be considered. -Various types of spatial operations are covered in Section \@ref(topological-relations), on topological relations between vector features. -\index{spatial operations} +Spatial operations differ from non-spatial operations in in a number of ways, however: +Spatial joins, for example, can be done in a number of ways --- including matching entities that intersect with or are within a certain distance of the target dataset --- while the attribution joins discussed in Section \@ref(vector-attribute-joining) in the previous chapter can only be done in one way (except when using fuzzy joins, as described in the documentation of the [**fuzzyjoin**](https://cran.r-project.org/package=fuzzyjoin) package). +The *type* of spatial relationship between objects must be considered when undertaking spatial operations, as described in Section \@ref(topological-relations), on topological relations between vector features. +\index{spatial operations}. Another unique aspect of spatial objects is distance. All spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship, as described in the context of vector data in Section \@ref(distance-relations). From 4ae179b5f47ecd0493f75e3707dbb7ab04687d7f Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 21:03:28 +0000 Subject: [PATCH 010/104] Minor tweak to c4 --- 04-spatial-operations.Rmd | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index af0c2fca4..8dced00db 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -31,8 +31,7 @@ Spatial operations differ from non-spatial operations in in a number of ways, ho Spatial joins, for example, can be done in a number of ways --- including matching entities that intersect with or are within a certain distance of the target dataset --- while the attribution joins discussed in Section \@ref(vector-attribute-joining) in the previous chapter can only be done in one way (except when using fuzzy joins, as described in the documentation of the [**fuzzyjoin**](https://cran.r-project.org/package=fuzzyjoin) package). The *type* of spatial relationship between objects must be considered when undertaking spatial operations, as described in Section \@ref(topological-relations), on topological relations between vector features. \index{spatial operations}. -Another unique aspect of spatial objects is distance. -All spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship, as described in the context of vector data in Section \@ref(distance-relations). +Another unique aspect of spatial objects is distance: all spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship, as described in the context of vector data in Section \@ref(distance-relations). Spatial operations on raster objects include subsetting --- covered in Section \@ref(spatial-raster-subsetting) --- and merging several raster 'tiles' into a single object, as demonstrated in Section \@ref(merging-rasters). *Map algebra* covers a range of operations that modify raster cell values, with or without reference to surrounding cell values, is is vital for many applications. From 55f20af61525dc25083b6056976b73be19d56af3 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 21:05:09 +0000 Subject: [PATCH 011/104] It uses the terra not raster package --- 04-spatial-operations.Rmd | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 8dced00db..fe956b937 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -44,7 +44,8 @@ It is important to note that spatial operations that use two spatial objects rel ## Spatial operations on vector data {#spatial-vec} -This section provides an overview of spatial operations on vector geographic data represented as simple features in the **sf** package before Section \@ref(spatial-ras), which presents spatial methods using the **raster** package. +This section provides an overview of spatial operations on vector geographic data represented as simple features in the **sf** package. +Section \@ref(spatial-ras) presents spatial operations on raster datasets using functions from the **terra** package. ### Spatial subsetting From 8eb32b6305b916b54b9654f64450cea9191517b0 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 21:22:26 +0000 Subject: [PATCH 012/104] Introduce terra pkg in context of raster in c2 --- 02-spatial-data.Rmd | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/02-spatial-data.Rmd b/02-spatial-data.Rmd index 7e8c08615..51160cda4 100644 --- a/02-spatial-data.Rmd +++ b/02-spatial-data.Rmd @@ -723,11 +723,23 @@ rast_table = tibble::tribble( ### An introduction to terra -The **terra** package supports raster objects in R. -It provides an extensive set of functions to create, read, export, manipulate and process raster datasets. +The **terra** package supports raster objects in R. +Like its predecessor **raster** (created by the same developer, Robert Hijmans), it provides an extensive set of functions to create, read, export, manipulate and process raster datasets. +**terra**'s functionality is largely the same as the more mature **raster** package, but there are some differences: **terra** functions are usually more computationally efficient than **raster** equivalents. + +On the other hand, the **raster** class system is popular and used by many other packages; as with **sf** and **sp**, the good news is you can seamlessly translate between the two types of object to ensure backwards compatibility with older scripts and packages, for example, with the functions `as.raster(my_rast)` . + +```{r, echo=FALSE, eval=FALSE} +# test raster/terra conversions +my_raster = as.raster(my_rast) +library(raster) +my_rast2 = rast(my_raster) # how to convert back to terra? fails for me (RL 2021-11) +identical(my_rast, my_rast2) +``` + Aside from general raster data manipulation, **terra** provides many low-level functions that can form the basis to develop more advanced raster functionality. \index{terra (package)|see {terra}} -**terra** also lets you work on large raster datasets that are too large to fit into the main memory. +**terra** also lets you work on large raster datasets that are too large to fit into the main memory. In this case, **terra** provides the possibility to divide the raster into smaller chunks, and processes these iteratively instead of loading the whole raster file into RAM. For the illustration of **terra** concepts, we will use datasets from the **spDataLarge**. @@ -738,6 +750,7 @@ First, let's create a `SpatRaster` object named `my_rast`: ```{r 02-spatial-data-37, message=FALSE} raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge") my_rast = rast(raster_filepath) +class(my_rast) ``` Typing the name of the raster into the console, will print out the raster header (dimensions, resolution, extent, CRS) and some additional information (class, data source, summary of the raster values): From 61c04f47259ad211b9474eec6af3ebc38f55a436 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 21:25:48 +0000 Subject: [PATCH 013/104] Cross reference c1 --- 02-spatial-data.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/02-spatial-data.Rmd b/02-spatial-data.Rmd index 51160cda4..258ac9614 100644 --- a/02-spatial-data.Rmd +++ b/02-spatial-data.Rmd @@ -727,7 +727,7 @@ The **terra** package supports raster objects in R. Like its predecessor **raster** (created by the same developer, Robert Hijmans), it provides an extensive set of functions to create, read, export, manipulate and process raster datasets. **terra**'s functionality is largely the same as the more mature **raster** package, but there are some differences: **terra** functions are usually more computationally efficient than **raster** equivalents. -On the other hand, the **raster** class system is popular and used by many other packages; as with **sf** and **sp**, the good news is you can seamlessly translate between the two types of object to ensure backwards compatibility with older scripts and packages, for example, with the functions `as.raster(my_rast)` . +On the other hand, the **raster** class system is popular and used by many other packages; as with **sf** and **sp**, the good news is you can seamlessly translate between the two types of object to ensure backwards compatibility with older scripts and packages, for example, with the functions `as.raster(my_rast)` (see the previous chapter for more on the evolution of R packages for working with geographic data). ```{r, echo=FALSE, eval=FALSE} # test raster/terra conversions From 942bba3d14633988ef0f14cc7ac954cee5951e98 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 00:48:59 +0000 Subject: [PATCH 014/104] Reword 4.2.1, make it clearer (hopefully) --- 04-spatial-operations.Rmd | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index fe956b937..acc5821e3 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -45,21 +45,22 @@ It is important to note that spatial operations that use two spatial objects rel ## Spatial operations on vector data {#spatial-vec} This section provides an overview of spatial operations on vector geographic data represented as simple features in the **sf** package. -Section \@ref(spatial-ras) presents spatial operations on raster datasets using functions from the **terra** package. +Section \@ref(spatial-ras) presents spatial operations on raster datasets using classes and functions from the **terra** package. ### Spatial subsetting -Spatial subsetting is the process of selecting features of a spatial object based on whether or not they in some way *relate* in space to another object. -It is analogous to *attribute subsetting* (covered in Section \@ref(vector-attribute-subsetting)) and can be done with the base R square bracket (`[`) operator or with the `filter()` function from the **tidyverse**\index{tidyverse (package)}. +Spatial subsetting is the process of taking a spatial object and returning a new object containing only features that *relate* in space to another object. +Analogous to *attribute subsetting* (covered in Section \@ref(vector-attribute-subsetting)), subsets of `sf` data frames can be created with square bracket (`[`) operator using the syntax `x[y, , op = st_intersects]`, where `x` is an `sf` object from which a subset of rows will be returned, `y` is the 'subsetting object' and `, op = st_intersects` is an optional argument that specifies the topological relation (also known as the binary predicate) used to do the subsetting. +The default topological relation used when an `op` argument is not provided is `st_intersects()`: the command `x[y, ]` is identical to `x[y, , op = st_intersects]` shown above but not `x[y, , op = st_disjoint]` (the meaning of these and other topological relations is described in the next section). +The `filter()` function from the **tidyverse**\index{tidyverse (package)} can also be used but this approach more verbose, as we will see in the examples below. \index{vector!subsetting} \index{spatial!subsetting} -An example of spatial subsetting is provided by the `nz` and `nz_height` datasets in **spData**. -These contain projected data on the 16 main regions and 101 highest points in New Zealand, respectively (Figure \@ref(fig:nz-subset)). -The following code chunk first creates an object representing Canterbury, then uses spatial subsetting to return all high points in the region: +To demonstrate spatial subsetting, we will use the `nz` and `nz_height` datasets in the **spData** package, which contain geographic data on the 16 main regions and 101 highest points in New Zealand, respectively (Figure \@ref(fig:nz-subset)), in a projected coordinate system. +The following code chunk creates an object representing Canterbury, then uses spatial subsetting to return all high points in the region: ```{r 04-spatial-operations-3} canterbury = nz %>% filter(Name == "Canterbury") @@ -80,8 +81,9 @@ p_hpnz2 = tm_shape(nz) + tm_polygons(col = "white") + tmap_arrange(p_hpnz1, p_hpnz2, ncol = 2) ``` -Like attribute subsetting `x[y, ]` subsets features of a *target* `x` using the contents of a *source* object `y`. -Instead of `y` being of class `logical` or `integer` --- a vector of `TRUE` and `FALSE` values or whole numbers --- for spatial subsetting it is another spatial (`sf`) object. +Like attribute subsetting, the command `x[y, ]` subsets features of a *target* `x` using the contents of a *source* object `y`. +Instead of `y` being a vector of class `logical` or `integer`, however, for spatial subsetting both `x` and `y` must be geographic objects. +Specifically, objects used for spatial subsetting in this way must have the class `sf` or `sfc`. Various *topological relations* can be used for spatial subsetting. These determine the type of spatial relationship that features in the target object must have with the subsetting object to be selected, including *touches*, *crosses* or *within* (see Section \@ref(topological-relations)). From 613fe320392a523cfe2bc414f795672080e6c228 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 01:14:47 +0000 Subject: [PATCH 015/104] Complete review of Section 4.2.1 --- 04-spatial-operations.Rmd | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index acc5821e3..552df3e73 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -81,15 +81,14 @@ p_hpnz2 = tm_shape(nz) + tm_polygons(col = "white") + tmap_arrange(p_hpnz1, p_hpnz2, ncol = 2) ``` -Like attribute subsetting, the command `x[y, ]` subsets features of a *target* `x` using the contents of a *source* object `y`. +Like attribute subsetting, the command `x[y, ]` (equivalent to `nz_height[canterbury, ]`) subsets features of a *target* `x` using the contents of a *source* object `y`. Instead of `y` being a vector of class `logical` or `integer`, however, for spatial subsetting both `x` and `y` must be geographic objects. -Specifically, objects used for spatial subsetting in this way must have the class `sf` or `sfc`. +Specifically, objects used for spatial subsetting in this way must have the class `sf` or `sfc`: both `nz` and `nz_height` are geographic vector data frames and have the class `sf`, and the result of the operation returns another `sf` object representing the features in the target `nz_height` object that intersect with (in this case high points that are located within) the `canterbury` region. -Various *topological relations* can be used for spatial subsetting. -These determine the type of spatial relationship that features in the target object must have with the subsetting object to be selected, including *touches*, *crosses* or *within* (see Section \@ref(topological-relations)). -*Intersects* is the default spatial subsetting operator, a default that returns `TRUE` for many types of spatial relations, including *touches*, *crosses* and *is within*. -These alternative spatial operators can be specified with the `op =` argument, a third argument that can be passed to the `[` operator for `sf` objects. -This is demonstrated in the following command which returns the opposite of `st_intersects()`, points that do not intersect with Canterbury (see Section \@ref(topological-relations)): +Various *topological relations* can be used for spatial subsetting which determine the type of spatial relationship that features in the target object must have with the subsetting object to be selected. +These include *touches*, *crosses* or *within*, as we will see shortly in Section \@ref(topological-relations). +The default setting `st_intersects` is a 'catch all' topological relation that will return features in the target that *touch*, *cross* or are *within* the source 'subsetting' object. +As indicated above, alternative spatial operators can be specified with the `op =` argument, as demonstrated in the following command which returns the opposite of `st_intersects()`, points that do not intersect with Canterbury (see Section \@ref(topological-relations)): ```{r 04-spatial-operations-4, eval=FALSE} nz_height[canterbury, , op = st_disjoint] @@ -101,16 +100,17 @@ One can use this to change the subsetting operation in many ways. `nz_height[canterbury, 2, op = st_disjoint]`, for example, returns the same rows but only includes the second attribute column (see `` sf:::`[.sf` `` and the `?sf` for details). ``` -For many applications, this is all you'll need to know about spatial subsetting for vector data. -In this case, you can safely skip to Section \@ref(topological-relations). - +For many applications, this is all you'll need to know about spatial subsetting for vector data: it just works. +If you are impatient to learn about more topological relations, beyond `st_intersects()` and `st_disjoint()`, skip to the next section (\@ref(topological-relations)). If you're interested in the details, including other ways of subsetting, read on. -Another way of doing spatial subsetting uses objects returned by *topological operators*. -This is demonstrated in the first command below: + +Another way of doing spatial subsetting uses objects returned by topological operators. +These objects can be useful in their own right, for example when exploring the graph network of relationships between contiguous regions, but they can also be used for subsetting, as demonstrated in the code chunk below: ```{r 04-spatial-operations-6} sel_sgbp = st_intersects(x = nz_height, y = canterbury) class(sel_sgbp) +sel_sgbp sel_logical = lengths(sel_sgbp) > 0 canterbury_height2 = nz_height[sel_logical, ] ``` @@ -125,7 +125,7 @@ Note: another way to return a logical output is by setting `sparse = FALSE` (mea Note: the solution involving `sgbp` objects is more generalisable though, as it works for many-to-many operations and has lower memory requirements. ``` -It should be noted that a logical can also be used with `filter()` as follows (`sparse = FALSE` is explained in Section \@ref(topological-relations)): +It should be noted that a logical can also be used with `filter()` as follows (`sparse = FALSE` is explained in Section \@ref(topological-relations)): ```{r 04-spatial-operations-7b} canterbury_height3 = nz_height %>% From eba5e53b90c2f0a67e9c22d7d27e11d75224449a Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 07:57:33 +0000 Subject: [PATCH 016/104] Improve segway to topological relations section --- 04-spatial-operations.Rmd | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 552df3e73..324e89cc1 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -133,7 +133,8 @@ canterbury_height3 = nz_height %>% ``` At this point, there are three versions of `canterbury_height`, one created with spatial subsetting directly and the other two via intermediary selection objects. -To explore these objects and spatial subsetting in more detail, see the supplementary vignettes on `subsetting` and [`tidyverse-pitfalls`](https://geocompr.github.io/geocompkg/articles/). +To explore spatial subsetting in more detail, see the supplementary vignettes on `subsetting` and [`tidyverse-pitfalls`](https://geocompr.github.io/geocompkg/articles/) on the [geocompkg website](https://geocompr.github.io/geocompkg/articles/). +The next section explores different types of spatial relation, also known as binary predicates, that can be used to identify whether or not two features are spatially related or not. ### Topological relations From 95e5a8653cb63e530a614ba167ff17ed8b4bf75f Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 08:39:19 +0000 Subject: [PATCH 017/104] Update with reference for #677 --- 04-spatial-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 324e89cc1..774ec9f24 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -138,7 +138,7 @@ The next section explores different types of spatial relation, also known as bin ### Topological relations -Topological relations describe the spatial relationships between objects. +Topological relations describe the spatial relationships between objects [@egenhofer_1990_mathematical]. To understand them, it helps to have some simple test data to work with. Figure \@ref(fig:relation-objects) contains a polygon (`a`), a line (`l`) and some points (`p`), which are created in the code below. \index{topological relations} From 5ab7ab6ea17598cde40f2fbc237acb81a39cb0f8 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 13:34:57 +0000 Subject: [PATCH 018/104] Add new refs and update bib --- geocompr.bib | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/geocompr.bib b/geocompr.bib index c5c9b98b6..e7a9302d4 100644 --- a/geocompr.bib +++ b/geocompr.bib @@ -529,6 +529,14 @@ @article{eddelbuettel_extending_2018 keywords = {\#nosource} } +@inproceedings{egenhofer_mathematical_1990, + title = {A Mathematical Framework for the Definition of Topological Relations}, + booktitle = {Proc. the Fourth International Symposium on Spatial Data Handing}, + author = {Egenhofer, Max and Herring, John}, + date = {1990}, + pages = {803--813} +} + @article{essd-11-647-2019, title = {{{ICGEM}} – 15 Years of Successful Collection and Distribution of Global Gravitational Models, Associated Services, and Future Plans}, author = {Ince, E. S. and Barthelmes, F. and Reißland, S. and Elger, K. and Förste, C. and Flechtner, F. and Schuh, H.}, @@ -676,7 +684,7 @@ @article{hengl_random_2018 } @article{hesselbarth_opensource_2021, - title = {Open-Source Tools in r for Landscape Ecology}, + title = {Open-Source Tools in {{R}} for Landscape Ecology}, author = {Hesselbarth, Maximillian H.K. and Nowosad, Jakub and Signer, Johannes and Graham, Laura J.}, date = {2021-06}, volume = {6}, @@ -1434,7 +1442,8 @@ @article{pebesma_simple_2018 @book{pebesma_spatial_2022, title = {Spatial {{Data Science}} with Applications in {{R}}}, author = {Pebesma, Edzer and Bivand, Roger}, - date = {2022} + date = {2022}, + url = {https://r-spatial.org/book} } @book{perpinan_rastervis_2016, @@ -1700,6 +1709,19 @@ @book{tomlin_geographic_1990 keywords = {\#nosource,Cartography,Data processing,Geographic information systems} } +@article{tomlin_map_1994, + title = {Map Algebra: One Perspective}, + shorttitle = {Map Algebra}, + author = {Tomlin, C. Dana}, + date = {1994}, + journaltitle = {Landscape and Urban Planning}, + volume = {30}, + number = {1-2}, + pages = {3--12}, + publisher = {{Elsevier}}, + doi = {10/dm2qm2} +} + @book{usgs_geological_2016, title = {U.{{S}}. {{Geological Survey}} ({{USGS}}) {{Earth Resources Observation}} and {{Science}} ({{EROS}}) {{Center}}}, author = {{USGS}}, From 98abe142dce110e6c86abc328820768b2eed5cdb Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 14:37:36 +0000 Subject: [PATCH 019/104] Add new references on algebraic topology --- geocompr.bib | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/geocompr.bib b/geocompr.bib index e7a9302d4..dd9162031 100644 --- a/geocompr.bib +++ b/geocompr.bib @@ -472,6 +472,20 @@ @article{coppock_history_1991 keywords = {\#nosource,⛔ No DOI found,History of GIS} } +@book{dieck_algebraic_2008, + title = {Algebraic Topology}, + author = {tom Dieck, Tammo}, + date = {2008}, + series = {{{EMS}} Textbooks in Mathematics}, + publisher = {{European Mathematical Society}}, + location = {{Zürich}}, + url = {https://www.maths.ed.ac.uk/~v1ranick/papers/diecktop.pdf}, + isbn = {978-3-03719-048-7}, + pagetotal = {567}, + keywords = {Algebraic topology,Homology theory,Homotopy theory}, + annotation = {OCLC: ocn261176011} +} + @book{diggle_modelbased_2007, title = {Model-Based Geostatistics}, author = {Diggle, Peter and Ribeiro, Paulo Justiniano}, @@ -1594,6 +1608,17 @@ @book{sherman_desktop_2008 keywords = {\#nosource} } +@book{spanier_algebraic_1995, + title = {Algebraic Topology}, + author = {Spanier, Edwin Henry}, + date = {1995}, + edition = {1. corr. Springer ed}, + publisher = {{Springer}}, + location = {{New York Berlin Barcelona Budapest}}, + isbn = {978-0-387-94426-5 978-3-540-90646-9 978-0-387-90646-1}, + pagetotal = {528} +} + @book{talbert_ancient_2014, title = {Ancient {{Perspectives}}: Maps and {{Their Place}} in {{Mesopotamia}}, {{Egypt}}, {{Greece}}, and {{Rome}}}, shorttitle = {Ancient {{Perspectives}}}, From eb51533dbb0ab20ac5ca5e9b2cb3c005f0bfcf53 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 14:40:02 +0000 Subject: [PATCH 020/104] Update text with refs for #677 --- 04-spatial-operations.Rmd | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 774ec9f24..0f0a2f311 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -138,7 +138,14 @@ The next section explores different types of spatial relation, also known as bin ### Topological relations -Topological relations describe the spatial relationships between objects [@egenhofer_1990_mathematical]. +Topological relations describe the spatial relationships between objects. +Specifically, "binary topological relationships" are logical statements (in that the answer can only be `TRUE` or `FALSE`) about the way that discrete regions of space defined by ordered sets of points (typically forming points, lines and polygons) relate to each other in two or more dimensions [@egenhofer_1990_mathematical]. +That may sound rather abstract and, indeed, the definition and classification of topological relations is based on earlier mathematical foundations first published as a book in 1966 [@spanier_algebraic_1995].^[ +See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. +] +However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for different types of topological relationships, as illustrated in Figure \@ref(fig:relation-objects). +In `sf`, functions testing for different types of topological relations are as 'binary predicates', as described in the vignette *Manipulating Simple Feature Geometries*, which can be viewed with the command [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), and in the help page [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) [@pebesma_simple_2018]. + To understand them, it helps to have some simple test data to work with. Figure \@ref(fig:relation-objects) contains a polygon (`a`), a line (`l`) and some points (`p`), which are created in the code below. \index{topological relations} From de750fc6f34b3856fbd8deed68bf2c1032e857fd Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 21:02:33 +0000 Subject: [PATCH 021/104] Test ggplot2 code for #677 --- 04-spatial-operations.Rmd | 49 ++++++++++++++++++++++++++++++++++----- 1 file changed, 43 insertions(+), 6 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 0f0a2f311..9ba950a24 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -146,16 +146,53 @@ See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for different types of topological relationships, as illustrated in Figure \@ref(fig:relation-objects). In `sf`, functions testing for different types of topological relations are as 'binary predicates', as described in the vignette *Manipulating Simple Feature Geometries*, which can be viewed with the command [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), and in the help page [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) [@pebesma_simple_2018]. -To understand them, it helps to have some simple test data to work with. -Figure \@ref(fig:relation-objects) contains a polygon (`a`), a line (`l`) and some points (`p`), which are created in the code below. +Figure \@ref(fig:relation-objects) + \index{topological relations} -```{r 04-spatial-operations-8} +```{r 04-spatial-operations-8, echo=FALSE} +library(ggplot2) # create a polygon -a_poly = st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1)))) -a = st_sfc(a_poly) +p1 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) +p2 = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) +p3 = st_sfc(st_polygon(list(rbind(c(0.7, 0.3), c(0.7, 0.95), c(0.9, 0.95), c(0.7, 0.3))))) +p4 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(0.9, 0.5), c(0.6, 0.1))))) +p5 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) + +# convert to sf objects +psf1 = st_sf(data.frame(Object = "Polygon1"), geometry = p1) +psf2 = st_sf(data.frame(Object = "Polygon2"), geometry = p2) +psf3 = st_sf(data.frame(Object = "Polygon3"), geometry = p3) +psf4 = st_sf(data.frame(Object = "Polygon4"), geometry = p4) +psf5 = st_sf(data.frame(Object = "Polygon5"), geometry = p5) + +# plot them +ps1 = rbind(psf1, psf2) +ps2 = rbind(psf1, psf3) +# sf::st_intersects(psf1, psf2, sparse = FALSE) +# sf::st_disjoint(psf1, psf2, sparse = FALSE) +# sf::st_overlaps(psf1, psf2, sparse = FALSE) +# sf::st_relate(psf1, psf2, sparse = FALSE) +# sf::st_relate(psf1, psf3) +# tm_shape(ps1) + tm_polygons(col = "Object") + tm_layout(legend.position = c("left", "top")) + +theme_set(new = theme_void()) +g1 = ggplot(ps1) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) +# g1 + annotate("text", x = 0.3, y = 0.9, label = "st_intersects(Polygon1, Polygon2)") +g1 + annotate("text", x = 0.1, y = 0.95, label = "intersects TRUE\ndisjoint FALSE\ntouches TRUE\n", hjust = "left", vjust = "top") +# Try annotating only which type of relations apply +# g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches", hjust = "left", vjust = "top") +g1an = g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches\nDE-9IM: FF2F11212", hjust = "left", vjust = "top") + +g2 = ggplot(ps2) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) +g2an = g2 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects,\ntouches, overlaps\nDE-9IM: 212101212", hjust = "left", vjust = "top") + +library(patchwork) +g1an + g2an + + # create a line -l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 2)) +l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 3)) l = st_sfc(l_line) # create points p_matrix = matrix(c(0.5, 1, -1, 0, 0, 1, 0.5, 1), ncol = 2) From 2011877f93918b6a5eab4bd4cbeb41ad960fdb58 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 22:15:34 +0000 Subject: [PATCH 022/104] Add placeholder section on DE-9IM #677 --- 04-spatial-operations.Rmd | 18 ++++++++++++++++-- geocompr.bib | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 48 insertions(+), 2 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 9ba950a24..c1ba199be 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -139,11 +139,13 @@ The next section explores different types of spatial relation, also known as bin ### Topological relations Topological relations describe the spatial relationships between objects. -Specifically, "binary topological relationships" are logical statements (in that the answer can only be `TRUE` or `FALSE`) about the way that discrete regions of space defined by ordered sets of points (typically forming points, lines and polygons) relate to each other in two or more dimensions [@egenhofer_1990_mathematical]. +"Binary topological relationships", to give them their full name, are logical statements (in that the answer can only be `TRUE` or `FALSE`) about the spatial relationships between two objects defined by ordered sets of points (typically forming points, lines and polygons) in two or more dimensions [@egenhofer_1990_mathematical]. That may sound rather abstract and, indeed, the definition and classification of topological relations is based on earlier mathematical foundations first published as a book in 1966 [@spanier_algebraic_1995].^[ See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. ] -However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for different types of topological relationships, as illustrated in Figure \@ref(fig:relation-objects). +However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships defined in English (and other human) languages as illustrated in Figure \@ref(fig:relations). + + In `sf`, functions testing for different types of topological relations are as 'binary predicates', as described in the vignette *Manipulating Simple Feature Geometries*, which can be viewed with the command [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), and in the help page [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) [@pebesma_simple_2018]. Figure \@ref(fig:relation-objects) @@ -310,6 +312,18 @@ plot(l, add = TRUE) plot(p, add = TRUE) ``` +### DE-9IM strings + +Underlying the binary predicates demonstrated in the previous section is the Dimensionally Extended 9-Intersection Model (DE-9IM). +The model was originally labelled "DE + 9IM" by its inventors, referring to the "dimension of the intersections of' boundaries, interiors, and exteriors of two features" [@clementini_1995_comparison] + +The way the model works can be demonstrated with reference to two intersecting polygons, as illustrated in Figure \@ref(fig:de-9im). + +```{r} + +``` + + ### Spatial joining Joining two non-spatial datasets relies on a shared 'key' variable, as described in Section \@ref(vector-attribute-joining). diff --git a/geocompr.bib b/geocompr.bib index dd9162031..1f6a7e1d9 100644 --- a/geocompr.bib +++ b/geocompr.bib @@ -429,6 +429,22 @@ @incollection{cheshire_spatial_2015 keywords = {\#nosource} } +@article{clementini_comparison_1995a, + title = {A Comparison of Methods for Representing Topological Relationships}, + author = {Clementini, Eliseo and Di Felice, Paolino}, + date = {1995-05-01}, + journaltitle = {Information Sciences - Applications}, + shortjournal = {Information Sciences - Applications}, + volume = {3}, + number = {3}, + pages = {149--178}, + issn = {1069-0115}, + doi = {10.1016/1069-0115(94)00033-X}, + url = {https://www.sciencedirect.com/science/article/pii/106901159400033X}, + urldate = {2021-11-13}, + langid = {english} +} + @article{conrad_system_2015, title = {System for {{Automated Geoscientific Analyses}} ({{SAGA}}) v. 2.1.4}, author = {Conrad, O. and Bechtel, B. and Bock, M. and Dietrich, H. and Fischer, E. and Gerlitz, L. and Wehberg, J. and Wichmann, V. and Böhner, J.}, @@ -1600,6 +1616,22 @@ @article{schratz_performance_nodate keywords = {\#nosource,⛔ No DOI found,Computer Science - Machine Learning,Statistics - Machine Learning,Statistics - Methodology} } +@article{shen_classification_2018, + title = {Classification of Topological Relations between Spatial Objects in Two-Dimensional Space within the Dimensionally Extended 9-Intersection Model}, + author = {Shen, Jingwei and Chen, Min and Liu, Xintao}, + date = {2018}, + journaltitle = {Transactions in GIS}, + volume = {22}, + number = {2}, + pages = {514--541}, + issn = {1467-9671}, + doi = {10.1111/tgis.12328}, + url = {https://onlinelibrary.wiley.com/doi/abs/10.1111/tgis.12328}, + urldate = {2021-11-13}, + langid = {english}, + annotation = {\_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/tgis.12328} +} + @book{sherman_desktop_2008, title = {Desktop {{GIS}}: Mapping the {{Planet}} with {{Open Source Tools}}}, author = {Sherman, Gary}, From cf85720f9178ce81fb2d0a952b6206a911ed6376 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sun, 14 Nov 2021 11:09:12 +0000 Subject: [PATCH 023/104] Tweak comment and replace code with link to terra issue in c2 --- 02-spatial-data.Rmd | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/02-spatial-data.Rmd b/02-spatial-data.Rmd index 258ac9614..e3fd940f5 100644 --- a/02-spatial-data.Rmd +++ b/02-spatial-data.Rmd @@ -730,14 +730,11 @@ Like its predecessor **raster** (created by the same developer, Robert Hijmans), On the other hand, the **raster** class system is popular and used by many other packages; as with **sf** and **sp**, the good news is you can seamlessly translate between the two types of object to ensure backwards compatibility with older scripts and packages, for example, with the functions `as.raster(my_rast)` (see the previous chapter for more on the evolution of R packages for working with geographic data). ```{r, echo=FALSE, eval=FALSE} -# test raster/terra conversions -my_raster = as.raster(my_rast) -library(raster) -my_rast2 = rast(my_raster) # how to convert back to terra? fails for me (RL 2021-11) -identical(my_rast, my_rast2) +# # test raster/terra conversions +# See https://github.com/rspatial/terra/issues/399 ``` -Aside from general raster data manipulation, **terra** provides many low-level functions that can form the basis to develop more advanced raster functionality. +In addition to functions for raster data manipulation, **terra** provides many low-level functions that can form a foundation for developing new tools for working with raster datasets. \index{terra (package)|see {terra}} **terra** also lets you work on large raster datasets that are too large to fit into the main memory. In this case, **terra** provides the possibility to divide the raster into smaller chunks, and processes these iteratively instead of loading the whole raster file into RAM. From 0bb2df39837b2fbf8d26b066f8b5fa4f5fafe6a6 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sun, 14 Nov 2021 11:09:36 +0000 Subject: [PATCH 024/104] Test code for #677 --- 04-spatial-operations.Rmd | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index c1ba199be..aaa7369db 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -320,7 +320,35 @@ The model was originally labelled "DE + 9IM" by its inventors, referring to the The way the model works can be demonstrated with reference to two intersecting polygons, as illustrated in Figure \@ref(fig:de-9im). ```{r} +b = st_sfc(st_point(c(0, 1)), st_point(c(1, 1))) # create 2 points +b = st_buffer(b, dist = 1) # convert points to circles +bsf = sf::st_sf(data.frame(Object = c("a", "b")), geometry = b) +b9 = replicate(bsf, n = 9, simplify = FALSE) +b9sf = do.call(rbind, b9) +domains = c("Interior", "Boundary", "Exterior") +b9sf$domain_a = rep(rep(domains, 3), each = 2) +b9sf$domain_b = rep(rep(domains, each = 3), each = 2) +library(ggplot2) +ggplot(b9sf) + + geom_sf() + + facet_grid(domain_a ~ domain_b) + +plot(b9sf) +tmap_arrange( + tm_shape(b) + tm_polygons(alpha = 0.5) + tm_layout(title = "Interior-Interior"), + tm_shape(b) + tm_polygons(alpha = 0.5) + tm_layout(title = "Interior-Boundary"), + tm_shape(b) + tm_polygons(alpha = 0.5), + tm_shape(b) + tm_polygons(alpha = 0.5), + tm_shape(b) + tm_polygons(alpha = 0.5), + tm_shape(b) + tm_polygons(alpha = 0.5), + tm_shape(b) + tm_polygons(alpha = 0.5), + tm_shape(b) + tm_polygons(alpha = 0.5), + tm_shape(b) + tm_polygons(alpha = 0.5), + nrow = 3 +) +plot(b) +text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y")) # add text ``` From b1c2c4e8c736657c95d2163c702b0e6d87d45af1 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 20:02:35 +0000 Subject: [PATCH 025/104] Add start fun for #677 --- code/de_9im.R | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 code/de_9im.R diff --git a/code/de_9im.R b/code/de_9im.R new file mode 100644 index 000000000..7881dd6f3 --- /dev/null +++ b/code/de_9im.R @@ -0,0 +1,18 @@ +#' This function visualises sf objects and returns info on the +#' types of spatial relationship there is between them +#' +#' Context: [robinlovelace/geocompr#677](https://github.com/Robinlovelace/geocompr/issues/677) +#' +#' @examples +#' library(sf) +#' x = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) +#' y = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) +#' de_9im +de_9im = function(x, y, object_names = c("x", "y")) { + if(is(x, "sfc") && is(y, "sfc")) { + x = st_sf(data.frame(Object = object_names[1]), geometry = x) + y = st_sf(data.frame(Object = object_names[2]), geometry = y) + } + browser() + ps1 = rbind(psf1, psf2) +} \ No newline at end of file From af95b09af66d4d57d2be1b499d304cb42c498fab Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 20:16:45 +0000 Subject: [PATCH 026/104] MVP for de_9im fun --- code/de_9im.R | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/code/de_9im.R b/code/de_9im.R index 7881dd6f3..a1b491e6c 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -7,12 +7,17 @@ #' library(sf) #' x = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) #' y = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) -#' de_9im -de_9im = function(x, y, object_names = c("x", "y")) { +#' de_9im(x, y) +de_9im = function(x, y, object_names = c("x", "y"), funs = list("st_intersects", "st_disjoint"), sparse = FALSE) { + requireNamespace("sf", quietly = TRUE) if(is(x, "sfc") && is(y, "sfc")) { x = st_sf(data.frame(Object = object_names[1]), geometry = x) y = st_sf(data.frame(Object = object_names[2]), geometry = y) } - browser() - ps1 = rbind(psf1, psf2) -} \ No newline at end of file + xy = rbind(x, y) + funs_matched = lapply(funs, match.fun) + res = lapply(seq(length(funs)), function(i) { + funs_matched[[i]](x, y, sparse = sparse) + }) + unlist(res) +} From 552474a75bca852f6a8e20ddb77a40314871beb7 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 20:18:59 +0000 Subject: [PATCH 027/104] Move test code to function script, fix build --- 04-spatial-operations.Rmd | 33 --------------------------------- code/de_9im.R | 12 ++++++++++++ 2 files changed, 12 insertions(+), 33 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index aaa7369db..2a49736cf 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -153,7 +153,6 @@ Figure \@ref(fig:relation-objects) \index{topological relations} ```{r 04-spatial-operations-8, echo=FALSE} -library(ggplot2) # create a polygon p1 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) p2 = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) @@ -161,38 +160,6 @@ p3 = st_sfc(st_polygon(list(rbind(c(0.7, 0.3), c(0.7, 0.95), c(0.9, 0.95), c(0.7 p4 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(0.9, 0.5), c(0.6, 0.1))))) p5 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) -# convert to sf objects -psf1 = st_sf(data.frame(Object = "Polygon1"), geometry = p1) -psf2 = st_sf(data.frame(Object = "Polygon2"), geometry = p2) -psf3 = st_sf(data.frame(Object = "Polygon3"), geometry = p3) -psf4 = st_sf(data.frame(Object = "Polygon4"), geometry = p4) -psf5 = st_sf(data.frame(Object = "Polygon5"), geometry = p5) - -# plot them -ps1 = rbind(psf1, psf2) -ps2 = rbind(psf1, psf3) -# sf::st_intersects(psf1, psf2, sparse = FALSE) -# sf::st_disjoint(psf1, psf2, sparse = FALSE) -# sf::st_overlaps(psf1, psf2, sparse = FALSE) -# sf::st_relate(psf1, psf2, sparse = FALSE) -# sf::st_relate(psf1, psf3) -# tm_shape(ps1) + tm_polygons(col = "Object") + tm_layout(legend.position = c("left", "top")) - -theme_set(new = theme_void()) -g1 = ggplot(ps1) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) -# g1 + annotate("text", x = 0.3, y = 0.9, label = "st_intersects(Polygon1, Polygon2)") -g1 + annotate("text", x = 0.1, y = 0.95, label = "intersects TRUE\ndisjoint FALSE\ntouches TRUE\n", hjust = "left", vjust = "top") -# Try annotating only which type of relations apply -# g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches", hjust = "left", vjust = "top") -g1an = g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches\nDE-9IM: FF2F11212", hjust = "left", vjust = "top") - -g2 = ggplot(ps2) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) -g2an = g2 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects,\ntouches, overlaps\nDE-9IM: 212101212", hjust = "left", vjust = "top") - -library(patchwork) -g1an + g2an - - # create a line l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 3)) l = st_sfc(l_line) diff --git a/code/de_9im.R b/code/de_9im.R index a1b491e6c..9c15ae0bb 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -21,3 +21,15 @@ de_9im = function(x, y, object_names = c("x", "y"), funs = list("st_intersects", }) unlist(res) } + +# # Test code to functionalize: +# theme_set(new = theme_void()) +# g1 = ggplot(ps1) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) +# # g1 + annotate("text", x = 0.3, y = 0.9, label = "st_intersects(Polygon1, Polygon2)") +# g1 + annotate("text", x = 0.1, y = 0.95, label = "intersects TRUE\ndisjoint FALSE\ntouches TRUE\n", hjust = "left", vjust = "top") +# # Try annotating only which type of relations apply +# # g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches", hjust = "left", vjust = "top") +# g1an = g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches\nDE-9IM: FF2F11212", hjust = "left", vjust = "top") +# +# g2 = ggplot(ps2) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) +# g2an = g2 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects,\ntouches, overlaps\nDE-9IM: 212101212", hjust = "left", vjust = "top") From c9310395a8c506134debef7ee0376f52eba4cb3e Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 20:28:13 +0000 Subject: [PATCH 028/104] Working de_9im function --- code/de_9im.R | 28 +++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/code/de_9im.R b/code/de_9im.R index 9c15ae0bb..d05c425ee 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -1,16 +1,22 @@ -#' This function visualises sf objects and returns info on the +#' This function visualises sf objects and returns info on the #' types of spatial relationship there is between them -#' +#' #' Context: [robinlovelace/geocompr#677](https://github.com/Robinlovelace/geocompr/issues/677) -#' -#' @examples +#' +#' @examples #' library(sf) #' x = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) #' y = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) #' de_9im(x, y) -de_9im = function(x, y, object_names = c("x", "y"), funs = list("st_intersects", "st_disjoint"), sparse = FALSE) { +de_9im = function(x, + y, + object_names = c("x", "y"), + funs = list("st_intersects", "st_disjoint"), + sparse = FALSE, + output = "character" + ) { requireNamespace("sf", quietly = TRUE) - if(is(x, "sfc") && is(y, "sfc")) { + if (is(x, "sfc") && is(y, "sfc")) { x = st_sf(data.frame(Object = object_names[1]), geometry = x) y = st_sf(data.frame(Object = object_names[2]), geometry = y) } @@ -19,17 +25,21 @@ de_9im = function(x, y, object_names = c("x", "y"), funs = list("st_intersects", res = lapply(seq(length(funs)), function(i) { funs_matched[[i]](x, y, sparse = sparse) }) - unlist(res) + res = unlist(res) + if(output == "character") { + res = unlist(funs)[res] + } + res } # # Test code to functionalize: # theme_set(new = theme_void()) # g1 = ggplot(ps1) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) # # g1 + annotate("text", x = 0.3, y = 0.9, label = "st_intersects(Polygon1, Polygon2)") -# g1 + annotate("text", x = 0.1, y = 0.95, label = "intersects TRUE\ndisjoint FALSE\ntouches TRUE\n", hjust = "left", vjust = "top") +# g1 + annotate("text", x = 0.1, y = 0.95, label = "intersects TRUE\ndisjoint FALSE\ntouches TRUE\n", hjust = "left", vjust = "top") # # Try annotating only which type of relations apply # # g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches", hjust = "left", vjust = "top") # g1an = g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches\nDE-9IM: FF2F11212", hjust = "left", vjust = "top") -# +# # g2 = ggplot(ps2) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) # g2an = g2 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects,\ntouches, overlaps\nDE-9IM: 212101212", hjust = "left", vjust = "top") From fcf6bf026e79c5333973aaeb3461e388de0bcfcf Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 21:20:13 +0000 Subject: [PATCH 029/104] Update code, fix build mkII --- 04-spatial-operations.Rmd | 13 ++++++++++++- code/de_9im.R | 16 +++++++++++++++- 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 2a49736cf..6e8b43ee9 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -153,12 +153,16 @@ Figure \@ref(fig:relation-objects) \index{topological relations} ```{r 04-spatial-operations-8, echo=FALSE} -# create a polygon +source("code/de_9im.R") p1 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) p2 = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) p3 = st_sfc(st_polygon(list(rbind(c(0.7, 0.3), c(0.7, 0.95), c(0.9, 0.95), c(0.7, 0.3))))) p4 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(0.9, 0.5), c(0.6, 0.1))))) p5 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) +de_9im(p1, p2) +de_9im(p1, p3) +de_9im(p1, p4) +de_9im(p1, p5) # create a line l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 3)) @@ -169,6 +173,13 @@ p_multi = st_multipoint(x = p_matrix) p = st_cast(st_sfc(p_multi), "POINT") ``` +```{r} +a_poly = st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1)))) +a = st_sfc(a_poly) +l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 2)) +``` + + ```{r relation-objects, echo=FALSE, fig.cap="Points (p 1 to 4), line and polygon objects arranged to illustrate topological relations.", fig.asp=1, out.width="50%", fig.scap="Demonstration of topological relations."} par(pty = "s") plot(a, border = "red", col = "gray", axes = TRUE) diff --git a/code/de_9im.R b/code/de_9im.R index d05c425ee..f7760a39c 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -11,7 +11,21 @@ de_9im = function(x, y, object_names = c("x", "y"), - funs = list("st_intersects", "st_disjoint"), + funs = list( + "st_intersects", + "st_disjoint", + "st_touches", + "st_crosses", + "st_within", + "st_contains", + "st_contains_properly", + "st_overlaps", + "st_equals", + "st_covers", + "st_covered_by" + # , + # "st_equals_exact" # requuires par argument + ), sparse = FALSE, output = "character" ) { From d771041b499c1c50cb6c0b427c03911a0b404e48 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 21:51:52 +0000 Subject: [PATCH 030/104] Progress on #677 --- 04-spatial-operations.Rmd | 1 + code/de_9im.R | 25 +++++++++++++++++++++---- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 6e8b43ee9..7b79cc4dd 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -162,6 +162,7 @@ p5 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) de_9im(p1, p2) de_9im(p1, p3) de_9im(p1, p4) +de_9im(p4, p1) de_9im(p1, p5) # create a line diff --git a/code/de_9im.R b/code/de_9im.R index f7760a39c..c1c9a9d92 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -8,9 +8,14 @@ #' x = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) #' y = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) #' de_9im(x, y) +#' p3 = st_sfc(st_polygon(list(rbind(c(0.7, 0.3), c(0.7, 0.95), c(0.9, 0.95), c(0.7, 0.3))))) +#' p4 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(0.9, 0.5), c(0.6, 0.1))))) +#' p5 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) +#' de_9im(x, p3) de_9im = function(x, y, object_names = c("x", "y"), + plot = TRUE, funs = list( "st_intersects", "st_disjoint", @@ -27,7 +32,8 @@ de_9im = function(x, # "st_equals_exact" # requuires par argument ), sparse = FALSE, - output = "character" + output = "character", + collapse = "\n" ) { requireNamespace("sf", quietly = TRUE) if (is(x, "sfc") && is(y, "sfc")) { @@ -43,9 +49,22 @@ de_9im = function(x, if(output == "character") { res = unlist(funs)[res] } + if(plot) { + res_text = paste(res, collapse = collapse) + message("Object x has the following spatial relations to y: ", res_text) + res = de_9im_plot(xy, label = res_text) + } res } +de_9im_plot = function(xy, label = "test", alpha = 0.5, show.legend = FALSE, x = 0.1, y = 0.95) { + require("ggplot2", quietly = TRUE) + # browser() + ggplot(xy) + geom_sf(aes(fill = Object), alpha = alpha, show.legend = show.legend) + + annotate("text", x = 0.1, y = 0.95, label = label, hjust = "left", vjust = "top") + + ggplot2::theme_void() +} + # # Test code to functionalize: # theme_set(new = theme_void()) # g1 = ggplot(ps1) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) @@ -53,7 +72,5 @@ de_9im = function(x, # g1 + annotate("text", x = 0.1, y = 0.95, label = "intersects TRUE\ndisjoint FALSE\ntouches TRUE\n", hjust = "left", vjust = "top") # # Try annotating only which type of relations apply # # g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches", hjust = "left", vjust = "top") -# g1an = g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches\nDE-9IM: FF2F11212", hjust = "left", vjust = "top") +# g1an = g1 + # -# g2 = ggplot(ps2) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) -# g2an = g2 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects,\ntouches, overlaps\nDE-9IM: 212101212", hjust = "left", vjust = "top") From 5a7002acb726915f8df498259fcf65509f64cda7 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 21:59:00 +0000 Subject: [PATCH 031/104] Require sf --- code/de_9im.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/code/de_9im.R b/code/de_9im.R index c1c9a9d92..e3b553ce5 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -35,7 +35,7 @@ de_9im = function(x, output = "character", collapse = "\n" ) { - requireNamespace("sf", quietly = TRUE) + require("sf", quietly = TRUE) if (is(x, "sfc") && is(y, "sfc")) { x = st_sf(data.frame(Object = object_names[1]), geometry = x) y = st_sf(data.frame(Object = object_names[2]), geometry = y) From 76499a39aed6c4382cc3dfed1acedb62cdf42810 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 22:18:24 +0000 Subject: [PATCH 032/104] New fig 4.2 for #677 --- 04-spatial-operations.Rmd | 12 +++++++++--- code/de_9im.R | 5 +++-- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 7b79cc4dd..0fcd56e58 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -152,19 +152,25 @@ Figure \@ref(fig:relation-objects) \index{topological relations} -```{r 04-spatial-operations-8, echo=FALSE} -source("code/de_9im.R") +```{r 04-spatial-operations-8, echo=FALSE, fig.cap="Topological relations between vector geometries", fig.show='hold', out.width="33%", message=FALSE, fig.width=3} +source("https://github.com/Robinlovelace/geocompr/raw/c4-v2-updates-rl/code/de_9im.R") +library(sf) p1 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) p2 = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) p3 = st_sfc(st_polygon(list(rbind(c(0.7, 0.3), c(0.7, 0.95), c(0.9, 0.95), c(0.7, 0.3))))) p4 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(0.9, 0.5), c(0.6, 0.1))))) -p5 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) +p5 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(1, 0.5), c(0.6, 0.1))))) +p6 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) de_9im(p1, p2) de_9im(p1, p3) de_9im(p1, p4) de_9im(p4, p1) de_9im(p1, p5) +de_9im(p1, p6) +``` + +```{r} # create a line l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 3)) l = st_sfc(l_line) diff --git a/code/de_9im.R b/code/de_9im.R index e3b553ce5..f9a61f1b7 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -27,7 +27,8 @@ de_9im = function(x, "st_overlaps", "st_equals", "st_covers", - "st_covered_by" + "st_covered_by", + ... # , # "st_equals_exact" # requuires par argument ), @@ -35,7 +36,7 @@ de_9im = function(x, output = "character", collapse = "\n" ) { - require("sf", quietly = TRUE) + require("sf") if (is(x, "sfc") && is(y, "sfc")) { x = st_sf(data.frame(Object = object_names[1]), geometry = x) y = st_sf(data.frame(Object = object_names[2]), geometry = y) From 9286ab98117dde0d90c7142413767ceb5cfde344 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 23:04:52 +0000 Subject: [PATCH 033/104] Fix polygons, improve vis for #677 --- 04-spatial-operations.Rmd | 21 ++++++++++++--------- code/de_9im.R | 15 ++++++++++----- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 0fcd56e58..fbfc8068f 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -152,19 +152,22 @@ Figure \@ref(fig:relation-objects) \index{topological relations} -```{r 04-spatial-operations-8, echo=FALSE, fig.cap="Topological relations between vector geometries", fig.show='hold', out.width="33%", message=FALSE, fig.width=3} -source("https://github.com/Robinlovelace/geocompr/raw/c4-v2-updates-rl/code/de_9im.R") +```{r 04-spatial-operations-8, echo=FALSE, fig.cap="Topological relations between vector geometries", fig.show='hold', out.width="33%", message=FALSE, fig.width=3, fig.height=3} +# source("https://github.com/Robinlovelace/geocompr/raw/c4-v2-updates-rl/code/de_9im.R") +source("code/de_9im.R") library(sf) -p1 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) -p2 = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) -p3 = st_sfc(st_polygon(list(rbind(c(0.7, 0.3), c(0.7, 0.95), c(0.9, 0.95), c(0.7, 0.3))))) -p4 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(0.9, 0.5), c(0.6, 0.1))))) -p5 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(1, 0.5), c(0.6, 0.1))))) -p6 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) +p1 = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(1, 0.5), c(0, 0))))) +p2 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 0.5), c(0, 0))))) +p3 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 0.7), c(0, 0))))) +p4 = st_sfc(st_polygon(list(rbind(c(0.7, 0.8), c(0.7, 0.5), c(0.9, 0.5), c(0.7, 0.8))))) +p5 = st_sfc(st_polygon(list(rbind(c(0.6, 0.7), c(0.7, 0.5), c(1, 0.5), c(0.6, 0.7))))) +p6 = st_sfc(st_polygon(list(rbind(c(0.1, 0), c(1, 0), c(1, 0.3), c(0.1, 0))))) +p7 = st_sfc(st_polygon(list(rbind(c(0.05, 0.4), c(0.05, 0.97), c(0.6, 0.97), c(0.5, 0.4), c(0.05, 0.4))))) + de_9im(p1, p2) de_9im(p1, p3) de_9im(p1, p4) -de_9im(p4, p1) +de_9im(p7, p1) de_9im(p1, p5) de_9im(p1, p6) ``` diff --git a/code/de_9im.R b/code/de_9im.R index f9a61f1b7..d27f52aad 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -27,14 +27,14 @@ de_9im = function(x, "st_overlaps", "st_equals", "st_covers", - "st_covered_by", - ... + "st_covered_by" # , # "st_equals_exact" # requuires par argument ), + include_relate = TRUE, sparse = FALSE, output = "character", - collapse = "\n" + collapse = " ✓\n" ) { require("sf") if (is(x, "sfc") && is(y, "sfc")) { @@ -50,6 +50,11 @@ de_9im = function(x, if(output == "character") { res = unlist(funs)[res] } + if(include_relate) { + relation = sf::st_relate(x, y) + relate_text = paste0(" \nDE-9IM string: \n", relation) + res = c(res, relate_text) + } if(plot) { res_text = paste(res, collapse = collapse) message("Object x has the following spatial relations to y: ", res_text) @@ -58,12 +63,12 @@ de_9im = function(x, res } -de_9im_plot = function(xy, label = "test", alpha = 0.5, show.legend = FALSE, x = 0.1, y = 0.95) { +de_9im_plot = function(xy, label = "test", alpha = 0.5, show.legend = FALSE, x = 0.1, y = 0.95, theme = ggplot2::theme_void()) { require("ggplot2", quietly = TRUE) # browser() ggplot(xy) + geom_sf(aes(fill = Object), alpha = alpha, show.legend = show.legend) + annotate("text", x = 0.1, y = 0.95, label = label, hjust = "left", vjust = "top") + - ggplot2::theme_void() + theme } # # Test code to functionalize: From 22cbb41328d5cc1d82542e8e9e2d779f983a3d54 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 23:58:38 +0000 Subject: [PATCH 034/104] Finish updating s. 4.2.2, close #677 --- 04-spatial-operations.Rmd | 96 +++++++++++++++++++++++---------------- 1 file changed, 57 insertions(+), 39 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index fbfc8068f..09761b643 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -143,16 +143,11 @@ Topological relations describe the spatial relationships between objects. That may sound rather abstract and, indeed, the definition and classification of topological relations is based on earlier mathematical foundations first published as a book in 1966 [@spanier_algebraic_1995].^[ See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. ] -However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships defined in English (and other human) languages as illustrated in Figure \@ref(fig:relations). - - -In `sf`, functions testing for different types of topological relations are as 'binary predicates', as described in the vignette *Manipulating Simple Feature Geometries*, which can be viewed with the command [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), and in the help page [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) [@pebesma_simple_2018]. - -Figure \@ref(fig:relation-objects) - +However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships defined in English (and other human languages) as illustrated in Figure \@ref(fig:relations), which shows a variety of geometry pairs and their associated relations. +The third and fourth pairs in Figure \@ref(fig:relations) (from left to right and then down) demonstrate that, for some relations, order is important: while the relations *equals*, *intersects*, *crosses*, *touches* and *overlaps* are symmetrical, meaning that if `function(x, y)` is true, `function(y, x)` will also by true, relations in which the order of the geomtries are important such as *contains* and *within* are not. \index{topological relations} -```{r 04-spatial-operations-8, echo=FALSE, fig.cap="Topological relations between vector geometries", fig.show='hold', out.width="33%", message=FALSE, fig.width=3, fig.height=3} +```{r relations, echo=FALSE, fig.cap="Topological relations between vector geometries", fig.show='hold', out.width="33%", message=FALSE, fig.width=3, fig.height=3} # source("https://github.com/Robinlovelace/geocompr/raw/c4-v2-updates-rl/code/de_9im.R") source("code/de_9im.R") library(sf) @@ -170,56 +165,79 @@ de_9im(p1, p4) de_9im(p7, p1) de_9im(p1, p5) de_9im(p1, p6) +# todo: add 3 more with line/point relations? ``` +In `sf`, functions testing for different types of topological relations are called 'binary predicates', as described in the vignette *Manipulating Simple Feature Geometries*, which can be viewed with the command [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), and in the help page [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) [@pebesma_simple_2018]. +To see how topological relations work in practice, let's create a simple reproducible example, building on the relations illustrated in Figure \@ref(fig:relations). -```{r} -# create a line -l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 3)) -l = st_sfc(l_line) -# create points -p_matrix = matrix(c(0.5, 1, -1, 0, 0, 1, 0.5, 1), ncol = 2) -p_multi = st_multipoint(x = p_matrix) -p = st_cast(st_sfc(p_multi), "POINT") +```{r, echo=FALSE} +# alternative way of getting reader to data, less clear... +# polygon_matrix = matrix(c( +# 0, 0.0, +# 0, 1.0, +# 1, 1.0, +# 1, 0.5, +# 0, 0.0 +# ), +# ncol = 2, +# byrow = TRUE +# ) ``` + ```{r} -a_poly = st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1)))) -a = st_sfc(a_poly) -l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 2)) +polygon_df = read.csv(text = "x, y +0, 0.0 +0, 1.0 +1, 1.0 +1, 0.5 +0, 0.0") +polygon_matrix = as.matrix(polygon_df) +polygon = st_polygon(list(polygon_matrix)) +polygon_sfc = st_sfc(polygon) +line_df = read.csv(text = "x,y +0.1, 0 +1, 0.1") +line = st_linestring(x = as.matrix(line_df)) +line_sfc = st_sfc(line) +# create points +point_df = read.csv(text = "x,y +0.1,0 +0.7,0.2 +0.4,0.8") +point_sf = st_as_sf(point_df, coords = c("x", "y")) ``` - ```{r relation-objects, echo=FALSE, fig.cap="Points (p 1 to 4), line and polygon objects arranged to illustrate topological relations.", fig.asp=1, out.width="50%", fig.scap="Demonstration of topological relations."} par(pty = "s") -plot(a, border = "red", col = "gray", axes = TRUE) -plot(l, add = TRUE) -plot(p, add = TRUE, lab = 1:4) -text(p_matrix[, 1] + 0.04, p_matrix[, 2] - 0.06, 1:4, cex = 1.3) +plot(polygon_sfc, border = "red", col = "gray", axes = TRUE) +plot(line_sfc, add = TRUE) +plot(point_sf, add = TRUE, lab = 1:4) +text(point_df[, 1] + 0.02, point_df[, 2] + 0.03, 1:3, cex = 1.3) ``` -A simple query is: which of the points in `p` intersect in some way with polygon `a`? +A simple query is: which of the points in `point_sf` intersect in some way with polygon `polygon_sfc`? The question can be answered by inspection (points 1 and 2 are over or touch the triangle). It can also be answered by using a *spatial predicate* such as *do the objects intersect*? This is implemented in **sf** as follows: ```{r 04-spatial-operations-9, eval=FALSE} -st_intersects(p, a) +st_intersects(point_sf, polygon_sfc) #> Sparse geometry binary ..., where the predicate was `intersects' -#> 1: 1 -#> 2: 1 -#> 3: (empty) -#> 4: (empty) +#> 1: (empty) +#> 2: (empty) +#> 3: 1 ``` The contents of the result should be as you expected: -the function returns a positive (`1`) result for the first two points, and a negative result (represented by an empty vector) for the last two. +the function returns a positive (`1`) for the third point, and a negative result (represented by an empty vector) for the first two which are outside the polygon's border. What may be unexpected is that the result comes in the form of a list of vectors. This *sparse matrix* output only registers a relation if one exists, reducing the memory requirements of topological operations on multi-feature objects. As we saw in the previous section, a *dense matrix* consisting of `TRUE` or `FALSE` values for each combination of features can also be returned when `sparse = FALSE`: ```{r 04-spatial-operations-10} -st_intersects(p, a, sparse = FALSE) +st_intersects(point_sf, polygon_sfc, sparse = FALSE) ``` The output is a matrix in which each row represents a feature in the target object and each column represents a feature in the selecting object. @@ -231,21 +249,21 @@ Note that `st_intersects()` returns `TRUE` for the second feature in the object The opposite of `st_intersects()` is `st_disjoint()`, which returns only objects that do not spatially relate in any way to the selecting object (note `[, 1]` converts the result into a vector): ```{r 04-spatial-operations-11} -st_disjoint(p, a, sparse = FALSE)[, 1] +st_disjoint(point_sf, polygon_sfc, sparse = FALSE)[, 1] ``` `st_within()` returns `TRUE` only for objects that are completely within the selecting object. This applies only to the first object, which is inside the triangular polygon, as illustrated below: ```{r 04-spatial-operations-12} -st_within(p, a, sparse = FALSE)[, 1] +st_within(point_sf, polygon_sfc, sparse = FALSE)[, 1] ``` Note that although the first point is *within* the triangle, it does not *touch* any part of its border. For this reason `st_touches()` only returns `TRUE` for the second point: ```{r 04-spatial-operations-13} -st_touches(p, a, sparse = FALSE)[, 1] +st_touches(point_sf, polygon_sfc, sparse = FALSE)[, 1] ``` What about features that do not touch, but *almost touch* the selection object? @@ -255,7 +273,7 @@ Note that although point 4 is one unit of distance from the nearest node of `a` This is illustrated in the code chunk below, the second line of which converts the lengthy list output into a `logical` object: ```{r 04-spatial-operations-14} -sel = st_is_within_distance(p, a, dist = 0.9) # can only return a sparse matrix +sel = st_is_within_distance(point_sf, polygon_sfc, dist = 0.9) # can only return a sparse matrix lengths(sel) > 0 ``` @@ -269,9 +287,9 @@ You can learn more at https://www.r-spatial.org/r/2017/06/22/spatial-index.html. ```{r 04-spatial-operations-16, eval=FALSE, echo=FALSE} # other tests -st_overlaps(p, a, sparse = FALSE) -st_covers(p, a, sparse = FALSE) -st_covered_by(p, a, sparse = FALSE) +st_overlaps(point_sf, polygon_sfc, sparse = FALSE) +st_covers(point_sf, polygon_sfc, sparse = FALSE) +st_covered_by(point_sf, polygon_sfc, sparse = FALSE) ``` ```{r 04-spatial-operations-17, eval=FALSE, echo=FALSE} From 55e40d449926a1edc5347b8761965f68d49cfc68 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 00:39:26 +0000 Subject: [PATCH 035/104] Add chunk to get around failing terra code --- 04-spatial-operations.Rmd | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 09761b643..26e68c01a 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -611,6 +611,11 @@ For the reader's convenience, these datasets can be also found in the **spData** ### Spatial subsetting {#spatial-raster-subsetting} +```{r} +knitr::opts_chunk$set(eval = FALSE) # terra code failing locally +``` + + The previous chapter (Section \@ref(manipulating-raster-objects)) demonstrated how to retrieve values associated with specific cell IDs or row and column combinations. Raster objects can also be extracted by location (coordinates) and other spatial objects. To use coordinates for subsetting, one can 'translate' the coordinates into a cell ID with the **terra** function `cellFromXY()`. From df4b7cc0f96f92103d2cf8a56454f0c58a9b1c7e Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 00:49:02 +0000 Subject: [PATCH 036/104] Fix references in c4 --- 04-spatial-operations.Rmd | 15 +++++---------- geocompr.bib | 39 ++++++++++++++++----------------------- 2 files changed, 21 insertions(+), 33 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 26e68c01a..3aeab48fe 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -139,7 +139,7 @@ The next section explores different types of spatial relation, also known as bin ### Topological relations Topological relations describe the spatial relationships between objects. -"Binary topological relationships", to give them their full name, are logical statements (in that the answer can only be `TRUE` or `FALSE`) about the spatial relationships between two objects defined by ordered sets of points (typically forming points, lines and polygons) in two or more dimensions [@egenhofer_1990_mathematical]. +"Binary topological relationships", to give them their full name, are logical statements (in that the answer can only be `TRUE` or `FALSE`) about the spatial relationships between two objects defined by ordered sets of points (typically forming points, lines and polygons) in two or more dimensions [@egenhofer_mathematical_1990]. That may sound rather abstract and, indeed, the definition and classification of topological relations is based on earlier mathematical foundations first published as a book in 1966 [@spanier_algebraic_1995].^[ See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. ] @@ -321,7 +321,7 @@ plot(p, add = TRUE) ### DE-9IM strings Underlying the binary predicates demonstrated in the previous section is the Dimensionally Extended 9-Intersection Model (DE-9IM). -The model was originally labelled "DE + 9IM" by its inventors, referring to the "dimension of the intersections of' boundaries, interiors, and exteriors of two features" [@clementini_1995_comparison] +The model was originally labelled "DE + 9IM" by its inventors, referring to the "dimension of the intersections of' boundaries, interiors, and exteriors of two features" [@clementini_comparison_1995] The way the model works can be demonstrated with reference to two intersecting polygons, as illustrated in Figure \@ref(fig:de-9im). @@ -611,11 +611,6 @@ For the reader's convenience, these datasets can be also found in the **spData** ### Spatial subsetting {#spatial-raster-subsetting} -```{r} -knitr::opts_chunk$set(eval = FALSE) # terra code failing locally -``` - - The previous chapter (Section \@ref(manipulating-raster-objects)) demonstrated how to retrieve values associated with specific cell IDs or row and column combinations. Raster objects can also be extracted by location (coordinates) and other spatial objects. To use coordinates for subsetting, one can 'translate' the coordinates into a cell ID with the **terra** function `cellFromXY()`. @@ -636,7 +631,7 @@ terra::extract(elev, matrix(c(0.1, 0.1), ncol = 2)) Raster objects can also be subset with another raster object, as demonstrated in the code chunk below: -```{r 04-spatial-operations-35} +```{r 04-spatial-operations-35, eval=FALSE} clip = rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45, resolution = 0.3, vals = rep(1, 9)) elev[clip] @@ -701,7 +696,7 @@ The next subsection explores these and related operations in more detail. ### Map algebra \index{map algebra} -The term 'map algebra' was coined in the late 1970s to describe a "set of conventions, capabilities, and techniques" for the analysis of geographic raster *and* (although less prominently) vector data [@tomlin_1994_map]. +The term 'map algebra' was coined in the late 1970s to describe a "set of conventions, capabilities, and techniques" for the analysis of geographic raster *and* (although less prominently) vector data [@tomlin_map_1994]. Although the concept never became widely adopted, the term usefully encapsulates and helps classify the range operations that can be undertaken on raster datasets. In this context we define map algebra more narrowly, as operations that modify or summarise raster cell values, with reference to surrounding cells, zones, or statistical functions that apply to every cell. @@ -872,7 +867,7 @@ This is in contrast to focal operations which return a raster object (see previo For example, to find the mean elevation for each grain size class (Figure \@ref(fig:cont-raster)), we use the `zonal()` function. -```{r 04-spatial-operations-43} +```{r 04-spatial-operations-43, eval=FALSE} z = zonal(elev, grain, fun = "mean") z ``` diff --git a/geocompr.bib b/geocompr.bib index 1f6a7e1d9..be2214d87 100644 --- a/geocompr.bib +++ b/geocompr.bib @@ -156,29 +156,18 @@ @inproceedings{bivand_open_2000 keywords = {\#nosource,⛔ No DOI found} } -@article{bivand_progress_2020a, +@article{bivand_progress_2020, + ids = {bivand_progress_2020a,bivand_progress_2021}, title = {Progress in the {{R}} Ecosystem for Representing and Handling Spatial Data}, author = {Bivand, Roger S.}, - date = {2020-10}, - journaltitle = {Journal of Geographical Systems}, - publisher = {{Springer Science and Business Media LLC}}, - doi = {10/ghnwg3}, - url = {https://doi.org/10.1007/s10109-020-00336-0} -} - -@article{bivand_progress_2021, - title = {Progress in the {{R}} Ecosystem for Representing and Handling Spatial Data}, - author = {Bivand, Roger S.}, - date = {2021-10-01}, + date = {2020-10-16}, journaltitle = {Journal of Geographical Systems}, shortjournal = {J Geogr Syst}, - volume = {23}, - number = {4}, - pages = {515--546}, + publisher = {{Springer Science and Business Media LLC}}, issn = {1435-5949}, doi = {10/ghnwg3}, url = {https://doi.org/10.1007/s10109-020-00336-0}, - urldate = {2021-11-04}, + urldate = {2020-11-08}, langid = {english} } @@ -429,7 +418,7 @@ @incollection{cheshire_spatial_2015 keywords = {\#nosource} } -@article{clementini_comparison_1995a, +@article{clementini_comparison_1995, title = {A Comparison of Methods for Representing Topological Relationships}, author = {Clementini, Eliseo and Di Felice, Paolino}, date = {1995-05-01}, @@ -1049,6 +1038,7 @@ @article{loecher_rgooglemaps_2015 } @article{loidl_spatial_2016, + ids = {loidl_spatial_2016a}, title = {Spatial Patterns and Temporal Dynamics of Urban Bicycle Crashes—{{A}} Case Study from {{Salzburg}} ({{Austria}})}, author = {Loidl, Martin and Traun, Christoph and Wallentin, Gudrun}, date = {2016-04-01}, @@ -1106,16 +1096,19 @@ @book{lovelace_geocomputation_2019 isbn = {1-138-30451-4} } -@article{lovelace_open_2021a, +@article{lovelace_open_2021, + ids = {lovelace_open_2021a}, title = {Open Source Tools for Geographic Analysis in Transport Planning}, author = {Lovelace, Robin}, - date = {2021-01}, - volume = {23}, - number = {4}, - pages = {547--578}, + date = {2021-01-16}, + journaltitle = {Journal of Geographical Systems}, + shortjournal = {J Geogr Syst}, publisher = {{Springer Science and Business Media LLC}}, + issn = {1435-5949}, doi = {10/ghtnrp}, - url = {https://doi.org/10.1007/s10109-020-00342-2} + url = {https://doi.org/10.1007/s10109-020-00342-2}, + urldate = {2021-01-17}, + langid = {english} } @article{lovelace_propensity_2017, From cc533650f72f5851c6c571c9ee2da9542a64c711 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 00:51:57 +0000 Subject: [PATCH 037/104] Comment out experimental new de-9im section for now... --- 04-spatial-operations.Rmd | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 3aeab48fe..0d11bcc03 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -318,14 +318,15 @@ plot(l, add = TRUE) plot(p, add = TRUE) ``` -### DE-9IM strings + + -Underlying the binary predicates demonstrated in the previous section is the Dimensionally Extended 9-Intersection Model (DE-9IM). -The model was originally labelled "DE + 9IM" by its inventors, referring to the "dimension of the intersections of' boundaries, interiors, and exteriors of two features" [@clementini_comparison_1995] + + -The way the model works can be demonstrated with reference to two intersecting polygons, as illustrated in Figure \@ref(fig:de-9im). + -```{r} +```{r, echo=FALSE, eval=FALSE} b = st_sfc(st_point(c(0, 1)), st_point(c(1, 1))) # create 2 points b = st_buffer(b, dist = 1) # convert points to circles bsf = sf::st_sf(data.frame(Object = c("a", "b")), geometry = b) @@ -357,7 +358,6 @@ plot(b) text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y")) # add text ``` - ### Spatial joining Joining two non-spatial datasets relies on a shared 'key' variable, as described in Section \@ref(vector-attribute-joining). From c4a62c3ef3370eefc2e3e8a418e7a1147fd934ce Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 00:55:03 +0000 Subject: [PATCH 038/104] Split paragraph in s. 4.2.2 in 2 --- 04-spatial-operations.Rmd | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 0d11bcc03..031278611 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -143,7 +143,8 @@ Topological relations describe the spatial relationships between objects. That may sound rather abstract and, indeed, the definition and classification of topological relations is based on earlier mathematical foundations first published as a book in 1966 [@spanier_algebraic_1995].^[ See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. ] -However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships defined in English (and other human languages) as illustrated in Figure \@ref(fig:relations), which shows a variety of geometry pairs and their associated relations. + +Despite their mathematical origins, topological relations can be understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships defined in English (and other human languages) as illustrated in Figure \@ref(fig:relations), which shows a variety of geometry pairs and their associated relations. The third and fourth pairs in Figure \@ref(fig:relations) (from left to right and then down) demonstrate that, for some relations, order is important: while the relations *equals*, *intersects*, *crosses*, *touches* and *overlaps* are symmetrical, meaning that if `function(x, y)` is true, `function(y, x)` will also by true, relations in which the order of the geomtries are important such as *contains* and *within* are not. \index{topological relations} From 0302b0712a7dcf2fa1ee7b1a946c3c6c1cfb3cf9 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 00:59:57 +0000 Subject: [PATCH 039/104] Update text in 4.2.2, improve figure 4.2 caption --- 04-spatial-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 031278611..d949c5523 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -148,7 +148,7 @@ Despite their mathematical origins, topological relations can be understood intu The third and fourth pairs in Figure \@ref(fig:relations) (from left to right and then down) demonstrate that, for some relations, order is important: while the relations *equals*, *intersects*, *crosses*, *touches* and *overlaps* are symmetrical, meaning that if `function(x, y)` is true, `function(y, x)` will also by true, relations in which the order of the geomtries are important such as *contains* and *within* are not. \index{topological relations} -```{r relations, echo=FALSE, fig.cap="Topological relations between vector geometries", fig.show='hold', out.width="33%", message=FALSE, fig.width=3, fig.height=3} +```{r relations, echo=FALSE, fig.cap="Topological relations between vector geometries, inspired by Figures 1 and 2 in Egenhofer and Herring (1990). The relations for which the function(x, y) is true are printed for each geometry pair, with x represented in pink and y represented in blue.", fig.show='hold', out.width="33%", message=FALSE, fig.width=3, fig.height=3} # source("https://github.com/Robinlovelace/geocompr/raw/c4-v2-updates-rl/code/de_9im.R") source("code/de_9im.R") library(sf) From b314c83b9f3993aa3c3a139486c4fd24fee2ec4b Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 01:08:53 +0000 Subject: [PATCH 040/104] Break-up geometry creation code chunk --- 04-spatial-operations.Rmd | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index d949c5523..a4dfce107 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -170,7 +170,8 @@ de_9im(p1, p6) ``` In `sf`, functions testing for different types of topological relations are called 'binary predicates', as described in the vignette *Manipulating Simple Feature Geometries*, which can be viewed with the command [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), and in the help page [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) [@pebesma_simple_2018]. -To see how topological relations work in practice, let's create a simple reproducible example, building on the relations illustrated in Figure \@ref(fig:relations). +To see how topological relations work in practice, let's create a simple reproducible example, building on the relations illustrated in Figure \@ref(fig:relations) and consolidating knowledge of how vector geometries are represented from a previous chapter (Section \@ref(geometry)). +Note that to create tabular data representing coordinates (x and y) of the polygon vertices, we use the base R function `read.csv()` and then convert the result into a matrix, a `POLYGON` and finally an `sfc` object: ```{r, echo=FALSE} # alternative way of getting reader to data, less clear... @@ -197,6 +198,12 @@ polygon_df = read.csv(text = "x, y polygon_matrix = as.matrix(polygon_df) polygon = st_polygon(list(polygon_matrix)) polygon_sfc = st_sfc(polygon) +``` + +We will create additional geometries to demonstrate spatial relations with the following commands. +Note the use of the function `st_as_sf()` and the argument `coords` to efficiently convert from a data frame containing columns representing coordinates to an `sf` object containing points: + +```{r} line_df = read.csv(text = "x,y 0.1, 0 1, 0.1") From 87eccc45b15cf2f505532ed52bf005a69c4ee155 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 01:20:44 +0000 Subject: [PATCH 041/104] Finish review of Section 4.2.2 --- 04-spatial-operations.Rmd | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index a4dfce107..b1f16d1c6 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -200,7 +200,7 @@ polygon = st_polygon(list(polygon_matrix)) polygon_sfc = st_sfc(polygon) ``` -We will create additional geometries to demonstrate spatial relations with the following commands. +We will create additional geometries to demonstrate spatial relations with the following commands which, when plotted on top of the polygon created above, relate in space to one another, as shown in Figure \@ref(fig:relation-objects). Note the use of the function `st_as_sf()` and the argument `coords` to efficiently convert from a data frame containing columns representing coordinates to an `sf` object containing points: ```{r} @@ -217,7 +217,7 @@ point_df = read.csv(text = "x,y point_sf = st_as_sf(point_df, coords = c("x", "y")) ``` -```{r relation-objects, echo=FALSE, fig.cap="Points (p 1 to 4), line and polygon objects arranged to illustrate topological relations.", fig.asp=1, out.width="50%", fig.scap="Demonstration of topological relations."} +```{r relation-objects, echo=FALSE, fig.cap="Points (`point_df` 1 to 3), line and polygon objects arranged to illustrate topological relations.", fig.asp=1, out.width="50%", fig.scap="Demonstration of topological relations."} par(pty = "s") plot(polygon_sfc, border = "red", col = "gray", axes = TRUE) plot(line_sfc, add = TRUE) @@ -249,11 +249,11 @@ st_intersects(point_sf, polygon_sfc, sparse = FALSE) ``` The output is a matrix in which each row represents a feature in the target object and each column represents a feature in the selecting object. -In this case, only the first two features in `p` intersect with `a` and there is only one feature in `a` so the result has only one column. +In this case, only the third feature in `point_sf` intersects with `polygon_sfc`. +There is only one feature in `polygon_sfc` so the result has only one column. The result can be used for subsetting as we saw in Section \@ref(spatial-subsetting). -Note that `st_intersects()` returns `TRUE` for the second feature in the object `p` even though it just touches the polygon `a`: *intersects* is a 'catch-all' topological operation which identifies many types of spatial relation. - +Note: `st_intersects()` returns `TRUE` even in cases where the features just touch: *intersects* is a 'catch-all' topological operation which identifies many types of spatial relation, as illustrated in Figure \@ref(fig:relate). The opposite of `st_intersects()` is `st_disjoint()`, which returns only objects that do not spatially relate in any way to the selecting object (note `[, 1]` converts the result into a vector): ```{r 04-spatial-operations-11} @@ -261,14 +261,14 @@ st_disjoint(point_sf, polygon_sfc, sparse = FALSE)[, 1] ``` `st_within()` returns `TRUE` only for objects that are completely within the selecting object. -This applies only to the first object, which is inside the triangular polygon, as illustrated below: +This also only applies to the third point, which is inside the polygon, as illustrated below: ```{r 04-spatial-operations-12} st_within(point_sf, polygon_sfc, sparse = FALSE)[, 1] ``` Note that although the first point is *within* the triangle, it does not *touch* any part of its border. -For this reason `st_touches()` only returns `TRUE` for the second point: +For this reason `st_touches()` only returns `FALSE` for all points: ```{r 04-spatial-operations-13} st_touches(point_sf, polygon_sfc, sparse = FALSE)[, 1] @@ -277,11 +277,12 @@ st_touches(point_sf, polygon_sfc, sparse = FALSE)[, 1] What about features that do not touch, but *almost touch* the selection object? These can be selected using `st_is_within_distance()`, which has an additional `dist` argument. It can be used to set how close target objects need to be before they are selected. -Note that although point 4 is one unit of distance from the nearest node of `a` (at point 2 in Figure \@ref(fig:relation-objects)), it is still selected when the distance is set to 0.9. -This is illustrated in the code chunk below, the second line of which converts the lengthy list output into a `logical` object: +Note that although point 4 is one unit of distance from the nearest node of `polygon_sfc` (at point 2 in Figure \@ref(fig:relation-objects)), it is still selected when the distance is set to 0.9. +This is illustrated in the code chunk below, the second line of which uses the function `lengths()` to convert the lengthy list output into a `logical` object: ```{r 04-spatial-operations-14} -sel = st_is_within_distance(point_sf, polygon_sfc, dist = 0.9) # can only return a sparse matrix +# can only return a sparse matrix +sel = st_is_within_distance(point_sf, polygon_sfc, dist = 0.1) lengths(sel) > 0 ``` From 876e4ac70a793f77ebe9c6894743899440bd0e0f Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 6 Nov 2021 23:03:25 +0000 Subject: [PATCH 042/104] Split-out pkg loading and reading for #669 --- 04-spatial-operations.Rmd | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index b28a364b8..509aad0a4 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -9,10 +9,16 @@ library(sf) library(terra) library(dplyr) library(spData) +``` + +- You also need to read in a couple of datasets as follows for Section \@ref(spatial-ras) + +```{r 04-spatial-operations-1-1} elev = rast(system.file("raster/elev.tif", package = "spData")) grain = rast(system.file("raster/grain.tif", package = "spData")) ``` + ## Introduction Spatial operations are a vital part of geocomputation\index{geocomputation}. From 490e37a0b62819aaa8e26c3c1bc52cfa5361739b Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 8 Nov 2021 11:55:44 +0000 Subject: [PATCH 043/104] Refactor index.Rmd text about bookdown, update links --- index.Rmd | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/index.Rmd b/index.Rmd index 1b75df6a4..47dcffba2 100644 --- a/index.Rmd +++ b/index.Rmd @@ -34,12 +34,11 @@ This is the online home of *Geocomputation with R*, a book on geographic data an The geocompr book cover **Note**: The first edition of the book has been published by CRC Press in the [R Series](https://www.routledge.com/Chapman--HallCRC-The-R-Series/book-series/CRCTHERSER). -You can buy the book from [CRC Press](https://www.routledge.com/9781138304512), or [Amazon](https://www.amazon.com/Geocomputation-R-Robin-Lovelace-dp-0367670577/dp/0367670577/), and see the archived first edition on the open source book platform [bookdown.org](https://bookdown.org/robinlovelace/geocompr/). +You can buy the book from [CRC Press](https://www.routledge.com/9781138304512), or [Amazon](https://www.amazon.com/Geocomputation-R-Robin-Lovelace-dp-0367670577/dp/0367670577/), and see the archived **First Edition** hosted on [bookdown.org](https://bookdown.org/robinlovelace/geocompr/). -Inspired by [**bookdown**](https://github.com/rstudio/bookdown) and the Free and Open Source Software for Geospatial ([FOSS4G](https://foss4g.org/)) movement, this book is open source. -This ensures its contents are reproducible and publicly accessible for people worldwide. - -The online version of the book is hosted at [geocompr.robinlovelace.net](https://geocompr.robinlovelace.net) and kept up-to-date by [GitHub Actions](https://github.com/Robinlovelace/geocompr/actions), which provides information on its 'build status' as follows: +Inspired by [**bookdown**](https://bookdown.org/) and the Free and Open Source Software for Geospatial ([FOSS4G](https://foss4g.org/)) movement, the code and prose underlying this book is open source, ensuring it's reproducible, accessible and modifiable (e.g. in case you find an inevitable typo) for the benefit of people worldwide. +The online version of the book is hosted at [geocompr.robinlovelace.net](https://geocompr.robinlovelace.net) and kept up-to-date by [GitHub Actions](https://github.com/Robinlovelace/geocompr/actions). +Its current 'build status' as follows: [![Actions](https://github.com/Robinlovelace/geocompr/workflows/Render/badge.svg)](https://github.com/Robinlovelace/geocompr/actions) ``` From 5e7611db3c5c005999f2ffa0f8cdb177ffa3338b Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 11 Nov 2021 18:04:09 +0000 Subject: [PATCH 044/104] Define spatial ops upfront --- 04-spatial-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 509aad0a4..739a23c49 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -21,7 +21,7 @@ grain = rast(system.file("raster/grain.tif", package = "spData")) ## Introduction -Spatial operations are a vital part of geocomputation\index{geocomputation}. +Spatial operations, including spatial joins between vector datasets and local and focal operations on raster datasets, are a vital part of geocomputation\index{geocomputation}. This chapter shows how spatial objects can be modified in a multitude of ways based on their location and shape. The content builds on the previous chapter because many spatial operations have a non-spatial (attribute) equivalent. This is especially true for *vector* operations: Section \@ref(vector-attribute-manipulation) on vector attribute manipulation provides the basis for understanding its spatial counterpart, namely spatial subsetting (covered in Section \@ref(spatial-subsetting)). From 70140ccb9199bebc6f7e55a720147b6939a39ce3 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 11 Nov 2021 18:41:39 +0000 Subject: [PATCH 045/104] Update prose in Section 4.1 --- 04-spatial-operations.Rmd | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 739a23c49..40ef57f1b 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -23,27 +23,25 @@ grain = rast(system.file("raster/grain.tif", package = "spData")) Spatial operations, including spatial joins between vector datasets and local and focal operations on raster datasets, are a vital part of geocomputation\index{geocomputation}. This chapter shows how spatial objects can be modified in a multitude of ways based on their location and shape. -The content builds on the previous chapter because many spatial operations have a non-spatial (attribute) equivalent. +Many spatial operations have a non-spatial (attribute) equivalent, so concepts such as subsetting and joining datasets demonstrated in the previous chapter are applicable here. This is especially true for *vector* operations: Section \@ref(vector-attribute-manipulation) on vector attribute manipulation provides the basis for understanding its spatial counterpart, namely spatial subsetting (covered in Section \@ref(spatial-subsetting)). Spatial joining (Section \@ref(spatial-joining)) and aggregation (Section \@ref(spatial-aggr)) also have non-spatial counterparts, covered in the previous chapter. Spatial operations differ from non-spatial operations in some ways, however. To illustrate the point, imagine you are researching road safety. Spatial joins can be used to find road speed limits related with administrative zones, even when no zone ID is provided. -But this raises the question: should the road completely fall inside a zone for its values to be joined? -Or is simply crossing or being within a certain distance sufficient? -When posing such questions, it becomes apparent that spatial operations differ substantially from attribute operations on data frames: +But this raises the question: should the road completely fall inside a zone for its values to be joined, or is simply crossing or being within a certain distance sufficient? +When posing such questions, it becomes apparent that spatial operations differ substantially from attribute operations on data frames and that spatial operations are often more complex: the *type* of spatial relationship between objects must be considered. -These are covered in Section \@ref(topological-relations), on topological relations. - +Various types of spatial operations are covered in Section \@ref(topological-relations), on topological relations between vector features. \index{spatial operations} Another unique aspect of spatial objects is distance. -All spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship. These are covered in Section \@ref(distance-relations). +All spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship, as described in the context of vector data in Section \@ref(distance-relations). -Spatial operations also apply to raster objects. -Spatial subsetting of raster objects is covered in Section \@ref(spatial-raster-subsetting); merging several raster 'tiles' into a single object is covered in Section \@ref(merging-rasters). -For many applications, the most important spatial operation on raster objects is *map algebra*, as we will see in Sections \@ref(map-algebra) to \@ref(global-operations-and-distances). -Map algebra is also the prerequisite for distance calculations on rasters, a technique which is covered in Section \@ref(global-operations-and-distances). +Spatial operations on raster objects include subsetting --- covered in Section \@ref(spatial-raster-subsetting) --- and merging several raster 'tiles' into a single object, as demonstrated in Section \@ref(merging-rasters). +*Map algebra* cover a range of operations that modify raster cell values, with or without reference to surrounding cell values, is is vital for many applications. +The concept of map algebra is introduced in Section \@ref(map-algebra); local, focal and zonal map algebra operations are covered in sections \@ref(local-operations), \@ref(focal-operations), and \@ref(zonal-operations), respectively. Global map algebra operations, which generate summary statistics representing an entire raster dataset, and distance calculations on rasters, are discussed in Section \@ref(global-operations-and-distances). +In the final section before the exercises (\@ref(merging-rasets)) the process of merging two raster datasets is discussed and demonstrated with reference to a reproducible example. ```{block2 04-spatial-operations-2, type='rmdnote'} It is important to note that spatial operations that use two spatial objects rely on both objects having the same coordinate reference system, a topic that was introduced in Section \@ref(crs-intro) and which will be covered in more depth in Chapter \@ref(reproj-geo-data). @@ -607,7 +605,7 @@ The next subsection explores these and related operations in more detail. ### Map algebra \index{map algebra} -Map algebra makes raster processing really fast. +Map algebra operations tend to be fast. This is because raster datasets only implicitly store coordinates. To derive the coordinate of a specific cell, we have to calculate it using its matrix position and the raster resolution and origin. For the processing, however, the geographic position of a cell is barely relevant as long as we make sure that the cell position is still the same after the processing. From 8f8dfe0736d948f9652ab36f8980b57755969948 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 08:21:02 +0000 Subject: [PATCH 046/104] Define map algebra --- 04-spatial-operations.Rmd | 3 +++ 1 file changed, 3 insertions(+) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 40ef57f1b..a5d41ae27 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -605,6 +605,9 @@ The next subsection explores these and related operations in more detail. ### Map algebra \index{map algebra} +The term 'map algebra' was coined in the late 1970s to describe a "set of conventions, capabilities, and techniques" for the analysis of geographic data [@tomlin_1994_map]. + + Map algebra operations tend to be fast. This is because raster datasets only implicitly store coordinates. To derive the coordinate of a specific cell, we have to calculate it using its matrix position and the raster resolution and origin. From 13272f5ab212c1dd5716efccb9541bc9e7652012 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 08:59:03 +0000 Subject: [PATCH 047/104] Update map algebra section --- 04-spatial-operations.Rmd | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index a5d41ae27..ab433c076 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -39,7 +39,7 @@ Another unique aspect of spatial objects is distance. All spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship, as described in the context of vector data in Section \@ref(distance-relations). Spatial operations on raster objects include subsetting --- covered in Section \@ref(spatial-raster-subsetting) --- and merging several raster 'tiles' into a single object, as demonstrated in Section \@ref(merging-rasters). -*Map algebra* cover a range of operations that modify raster cell values, with or without reference to surrounding cell values, is is vital for many applications. +*Map algebra* covers a range of operations that modify raster cell values, with or without reference to surrounding cell values, is is vital for many applications. The concept of map algebra is introduced in Section \@ref(map-algebra); local, focal and zonal map algebra operations are covered in sections \@ref(local-operations), \@ref(focal-operations), and \@ref(zonal-operations), respectively. Global map algebra operations, which generate summary statistics representing an entire raster dataset, and distance calculations on rasters, are discussed in Section \@ref(global-operations-and-distances). In the final section before the exercises (\@ref(merging-rasets)) the process of merging two raster datasets is discussed and demonstrated with reference to a reproducible example. @@ -605,20 +605,21 @@ The next subsection explores these and related operations in more detail. ### Map algebra \index{map algebra} -The term 'map algebra' was coined in the late 1970s to describe a "set of conventions, capabilities, and techniques" for the analysis of geographic data [@tomlin_1994_map]. +The term 'map algebra' was coined in the late 1970s to describe a "set of conventions, capabilities, and techniques" for the analysis of geographic raster *and* (although less prominently) vector data [@tomlin_1994_map]. +Although the concept never became widely adopted, the term usefully encapsulates and helps classify the range operations that can be undertaken on raster datasets. +In this context we define map algebra more narrowly, as operations that modify or summarise raster cell values, with reference to surrounding cells, zones, or statistical functions that apply to every cell. - -Map algebra operations tend to be fast. -This is because raster datasets only implicitly store coordinates. -To derive the coordinate of a specific cell, we have to calculate it using its matrix position and the raster resolution and origin. +Map algebra operations tend to be fast, because raster datasets only implicitly store coordinates (hence the oversimplifying phrase "raster is faster but vector is corrector"). +The location of cells in raster datasets can be calculated it using its matrix position and the resolution and origin of the dataset (stored in the header). For the processing, however, the geographic position of a cell is barely relevant as long as we make sure that the cell position is still the same after the processing. Additionally, if two or more raster datasets share the same extent, projection and resolution, one could treat them as matrices for the processing. -This is exactly what map algebra is doing in R. -First, the **terra** package checks the headers of the rasters on which to perform any algebraic operation, and only if they are correspondent to each other, the processing goes on. -And secondly, map algebra retains the so-called one-to-one locational correspondence. -This is where it substantially differs from matrix algebra which changes positions when for example multiplying or dividing matrices. -Map algebra (or cartographic modeling) divides raster operations into four subclasses [@tomlin_geographic_1990], with each working on one or several grids simultaneously: +This is the way that map algebra works with the **terra** package. +First, the headers of the raster datasets are queried and (in cases where map algebra operations work on more than one dataset) checked to ensure the datasets are compatible. +Second, map algebra retains the so-called one-to-one locational correspondence, meaning that cells cannot move. +This differs from matrix algebra, in which values change position, for example when multiplying or dividing matrices. + +Map algebra (or cartographic modeling with raster data) divides raster operations into four subclasses [@tomlin_geographic_1990], with each working on one or several grids simultaneously: 1. *Local* or per-cell operations 2. *Focal* or neighborhood operations. From a4f3a566a7e4bbed97c7e46a1e009ae5b6585c15 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 20:28:19 +0000 Subject: [PATCH 048/104] Shorten and simplify description of spatial vs attribute operations --- 04-spatial-operations.Rmd | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index ab433c076..af0c2fca4 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -27,14 +27,10 @@ Many spatial operations have a non-spatial (attribute) equivalent, so concepts s This is especially true for *vector* operations: Section \@ref(vector-attribute-manipulation) on vector attribute manipulation provides the basis for understanding its spatial counterpart, namely spatial subsetting (covered in Section \@ref(spatial-subsetting)). Spatial joining (Section \@ref(spatial-joining)) and aggregation (Section \@ref(spatial-aggr)) also have non-spatial counterparts, covered in the previous chapter. -Spatial operations differ from non-spatial operations in some ways, however. -To illustrate the point, imagine you are researching road safety. -Spatial joins can be used to find road speed limits related with administrative zones, even when no zone ID is provided. -But this raises the question: should the road completely fall inside a zone for its values to be joined, or is simply crossing or being within a certain distance sufficient? -When posing such questions, it becomes apparent that spatial operations differ substantially from attribute operations on data frames and that spatial operations are often more complex: -the *type* of spatial relationship between objects must be considered. -Various types of spatial operations are covered in Section \@ref(topological-relations), on topological relations between vector features. -\index{spatial operations} +Spatial operations differ from non-spatial operations in in a number of ways, however: +Spatial joins, for example, can be done in a number of ways --- including matching entities that intersect with or are within a certain distance of the target dataset --- while the attribution joins discussed in Section \@ref(vector-attribute-joining) in the previous chapter can only be done in one way (except when using fuzzy joins, as described in the documentation of the [**fuzzyjoin**](https://cran.r-project.org/package=fuzzyjoin) package). +The *type* of spatial relationship between objects must be considered when undertaking spatial operations, as described in Section \@ref(topological-relations), on topological relations between vector features. +\index{spatial operations}. Another unique aspect of spatial objects is distance. All spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship, as described in the context of vector data in Section \@ref(distance-relations). From 0061610e64e2f5f37c409dfdb88f9dc808c3f271 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 21:03:28 +0000 Subject: [PATCH 049/104] Minor tweak to c4 --- 04-spatial-operations.Rmd | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index af0c2fca4..8dced00db 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -31,8 +31,7 @@ Spatial operations differ from non-spatial operations in in a number of ways, ho Spatial joins, for example, can be done in a number of ways --- including matching entities that intersect with or are within a certain distance of the target dataset --- while the attribution joins discussed in Section \@ref(vector-attribute-joining) in the previous chapter can only be done in one way (except when using fuzzy joins, as described in the documentation of the [**fuzzyjoin**](https://cran.r-project.org/package=fuzzyjoin) package). The *type* of spatial relationship between objects must be considered when undertaking spatial operations, as described in Section \@ref(topological-relations), on topological relations between vector features. \index{spatial operations}. -Another unique aspect of spatial objects is distance. -All spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship, as described in the context of vector data in Section \@ref(distance-relations). +Another unique aspect of spatial objects is distance: all spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship, as described in the context of vector data in Section \@ref(distance-relations). Spatial operations on raster objects include subsetting --- covered in Section \@ref(spatial-raster-subsetting) --- and merging several raster 'tiles' into a single object, as demonstrated in Section \@ref(merging-rasters). *Map algebra* covers a range of operations that modify raster cell values, with or without reference to surrounding cell values, is is vital for many applications. From 40d7f491a72c8561700cc5c2eed5dcc2d58f2ee1 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 21:05:09 +0000 Subject: [PATCH 050/104] It uses the terra not raster package --- 04-spatial-operations.Rmd | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 8dced00db..fe956b937 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -44,7 +44,8 @@ It is important to note that spatial operations that use two spatial objects rel ## Spatial operations on vector data {#spatial-vec} -This section provides an overview of spatial operations on vector geographic data represented as simple features in the **sf** package before Section \@ref(spatial-ras), which presents spatial methods using the **raster** package. +This section provides an overview of spatial operations on vector geographic data represented as simple features in the **sf** package. +Section \@ref(spatial-ras) presents spatial operations on raster datasets using functions from the **terra** package. ### Spatial subsetting From cebae2c1be44dcab174e79dc24a2a272fc7d446f Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 21:22:26 +0000 Subject: [PATCH 051/104] Introduce terra pkg in context of raster in c2 --- 02-spatial-data.Rmd | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/02-spatial-data.Rmd b/02-spatial-data.Rmd index 7e8c08615..51160cda4 100644 --- a/02-spatial-data.Rmd +++ b/02-spatial-data.Rmd @@ -723,11 +723,23 @@ rast_table = tibble::tribble( ### An introduction to terra -The **terra** package supports raster objects in R. -It provides an extensive set of functions to create, read, export, manipulate and process raster datasets. +The **terra** package supports raster objects in R. +Like its predecessor **raster** (created by the same developer, Robert Hijmans), it provides an extensive set of functions to create, read, export, manipulate and process raster datasets. +**terra**'s functionality is largely the same as the more mature **raster** package, but there are some differences: **terra** functions are usually more computationally efficient than **raster** equivalents. + +On the other hand, the **raster** class system is popular and used by many other packages; as with **sf** and **sp**, the good news is you can seamlessly translate between the two types of object to ensure backwards compatibility with older scripts and packages, for example, with the functions `as.raster(my_rast)` . + +```{r, echo=FALSE, eval=FALSE} +# test raster/terra conversions +my_raster = as.raster(my_rast) +library(raster) +my_rast2 = rast(my_raster) # how to convert back to terra? fails for me (RL 2021-11) +identical(my_rast, my_rast2) +``` + Aside from general raster data manipulation, **terra** provides many low-level functions that can form the basis to develop more advanced raster functionality. \index{terra (package)|see {terra}} -**terra** also lets you work on large raster datasets that are too large to fit into the main memory. +**terra** also lets you work on large raster datasets that are too large to fit into the main memory. In this case, **terra** provides the possibility to divide the raster into smaller chunks, and processes these iteratively instead of loading the whole raster file into RAM. For the illustration of **terra** concepts, we will use datasets from the **spDataLarge**. @@ -738,6 +750,7 @@ First, let's create a `SpatRaster` object named `my_rast`: ```{r 02-spatial-data-37, message=FALSE} raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge") my_rast = rast(raster_filepath) +class(my_rast) ``` Typing the name of the raster into the console, will print out the raster header (dimensions, resolution, extent, CRS) and some additional information (class, data source, summary of the raster values): From 26bc96155b7dc1c95b9a6dce427ca5d080730cda Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 12 Nov 2021 21:25:48 +0000 Subject: [PATCH 052/104] Cross reference c1 --- 02-spatial-data.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/02-spatial-data.Rmd b/02-spatial-data.Rmd index 51160cda4..258ac9614 100644 --- a/02-spatial-data.Rmd +++ b/02-spatial-data.Rmd @@ -727,7 +727,7 @@ The **terra** package supports raster objects in R. Like its predecessor **raster** (created by the same developer, Robert Hijmans), it provides an extensive set of functions to create, read, export, manipulate and process raster datasets. **terra**'s functionality is largely the same as the more mature **raster** package, but there are some differences: **terra** functions are usually more computationally efficient than **raster** equivalents. -On the other hand, the **raster** class system is popular and used by many other packages; as with **sf** and **sp**, the good news is you can seamlessly translate between the two types of object to ensure backwards compatibility with older scripts and packages, for example, with the functions `as.raster(my_rast)` . +On the other hand, the **raster** class system is popular and used by many other packages; as with **sf** and **sp**, the good news is you can seamlessly translate between the two types of object to ensure backwards compatibility with older scripts and packages, for example, with the functions `as.raster(my_rast)` (see the previous chapter for more on the evolution of R packages for working with geographic data). ```{r, echo=FALSE, eval=FALSE} # test raster/terra conversions From 627214c1fc610b33aaaa4942129d94827285a7a3 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 00:48:59 +0000 Subject: [PATCH 053/104] Reword 4.2.1, make it clearer (hopefully) --- 04-spatial-operations.Rmd | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index fe956b937..acc5821e3 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -45,21 +45,22 @@ It is important to note that spatial operations that use two spatial objects rel ## Spatial operations on vector data {#spatial-vec} This section provides an overview of spatial operations on vector geographic data represented as simple features in the **sf** package. -Section \@ref(spatial-ras) presents spatial operations on raster datasets using functions from the **terra** package. +Section \@ref(spatial-ras) presents spatial operations on raster datasets using classes and functions from the **terra** package. ### Spatial subsetting -Spatial subsetting is the process of selecting features of a spatial object based on whether or not they in some way *relate* in space to another object. -It is analogous to *attribute subsetting* (covered in Section \@ref(vector-attribute-subsetting)) and can be done with the base R square bracket (`[`) operator or with the `filter()` function from the **tidyverse**\index{tidyverse (package)}. +Spatial subsetting is the process of taking a spatial object and returning a new object containing only features that *relate* in space to another object. +Analogous to *attribute subsetting* (covered in Section \@ref(vector-attribute-subsetting)), subsets of `sf` data frames can be created with square bracket (`[`) operator using the syntax `x[y, , op = st_intersects]`, where `x` is an `sf` object from which a subset of rows will be returned, `y` is the 'subsetting object' and `, op = st_intersects` is an optional argument that specifies the topological relation (also known as the binary predicate) used to do the subsetting. +The default topological relation used when an `op` argument is not provided is `st_intersects()`: the command `x[y, ]` is identical to `x[y, , op = st_intersects]` shown above but not `x[y, , op = st_disjoint]` (the meaning of these and other topological relations is described in the next section). +The `filter()` function from the **tidyverse**\index{tidyverse (package)} can also be used but this approach more verbose, as we will see in the examples below. \index{vector!subsetting} \index{spatial!subsetting} -An example of spatial subsetting is provided by the `nz` and `nz_height` datasets in **spData**. -These contain projected data on the 16 main regions and 101 highest points in New Zealand, respectively (Figure \@ref(fig:nz-subset)). -The following code chunk first creates an object representing Canterbury, then uses spatial subsetting to return all high points in the region: +To demonstrate spatial subsetting, we will use the `nz` and `nz_height` datasets in the **spData** package, which contain geographic data on the 16 main regions and 101 highest points in New Zealand, respectively (Figure \@ref(fig:nz-subset)), in a projected coordinate system. +The following code chunk creates an object representing Canterbury, then uses spatial subsetting to return all high points in the region: ```{r 04-spatial-operations-3} canterbury = nz %>% filter(Name == "Canterbury") @@ -80,8 +81,9 @@ p_hpnz2 = tm_shape(nz) + tm_polygons(col = "white") + tmap_arrange(p_hpnz1, p_hpnz2, ncol = 2) ``` -Like attribute subsetting `x[y, ]` subsets features of a *target* `x` using the contents of a *source* object `y`. -Instead of `y` being of class `logical` or `integer` --- a vector of `TRUE` and `FALSE` values or whole numbers --- for spatial subsetting it is another spatial (`sf`) object. +Like attribute subsetting, the command `x[y, ]` subsets features of a *target* `x` using the contents of a *source* object `y`. +Instead of `y` being a vector of class `logical` or `integer`, however, for spatial subsetting both `x` and `y` must be geographic objects. +Specifically, objects used for spatial subsetting in this way must have the class `sf` or `sfc`. Various *topological relations* can be used for spatial subsetting. These determine the type of spatial relationship that features in the target object must have with the subsetting object to be selected, including *touches*, *crosses* or *within* (see Section \@ref(topological-relations)). From d00f11638697cf6157d6e23b2cca5eb6081a499e Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 01:14:47 +0000 Subject: [PATCH 054/104] Complete review of Section 4.2.1 --- 04-spatial-operations.Rmd | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index acc5821e3..552df3e73 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -81,15 +81,14 @@ p_hpnz2 = tm_shape(nz) + tm_polygons(col = "white") + tmap_arrange(p_hpnz1, p_hpnz2, ncol = 2) ``` -Like attribute subsetting, the command `x[y, ]` subsets features of a *target* `x` using the contents of a *source* object `y`. +Like attribute subsetting, the command `x[y, ]` (equivalent to `nz_height[canterbury, ]`) subsets features of a *target* `x` using the contents of a *source* object `y`. Instead of `y` being a vector of class `logical` or `integer`, however, for spatial subsetting both `x` and `y` must be geographic objects. -Specifically, objects used for spatial subsetting in this way must have the class `sf` or `sfc`. +Specifically, objects used for spatial subsetting in this way must have the class `sf` or `sfc`: both `nz` and `nz_height` are geographic vector data frames and have the class `sf`, and the result of the operation returns another `sf` object representing the features in the target `nz_height` object that intersect with (in this case high points that are located within) the `canterbury` region. -Various *topological relations* can be used for spatial subsetting. -These determine the type of spatial relationship that features in the target object must have with the subsetting object to be selected, including *touches*, *crosses* or *within* (see Section \@ref(topological-relations)). -*Intersects* is the default spatial subsetting operator, a default that returns `TRUE` for many types of spatial relations, including *touches*, *crosses* and *is within*. -These alternative spatial operators can be specified with the `op =` argument, a third argument that can be passed to the `[` operator for `sf` objects. -This is demonstrated in the following command which returns the opposite of `st_intersects()`, points that do not intersect with Canterbury (see Section \@ref(topological-relations)): +Various *topological relations* can be used for spatial subsetting which determine the type of spatial relationship that features in the target object must have with the subsetting object to be selected. +These include *touches*, *crosses* or *within*, as we will see shortly in Section \@ref(topological-relations). +The default setting `st_intersects` is a 'catch all' topological relation that will return features in the target that *touch*, *cross* or are *within* the source 'subsetting' object. +As indicated above, alternative spatial operators can be specified with the `op =` argument, as demonstrated in the following command which returns the opposite of `st_intersects()`, points that do not intersect with Canterbury (see Section \@ref(topological-relations)): ```{r 04-spatial-operations-4, eval=FALSE} nz_height[canterbury, , op = st_disjoint] @@ -101,16 +100,17 @@ One can use this to change the subsetting operation in many ways. `nz_height[canterbury, 2, op = st_disjoint]`, for example, returns the same rows but only includes the second attribute column (see `` sf:::`[.sf` `` and the `?sf` for details). ``` -For many applications, this is all you'll need to know about spatial subsetting for vector data. -In this case, you can safely skip to Section \@ref(topological-relations). - +For many applications, this is all you'll need to know about spatial subsetting for vector data: it just works. +If you are impatient to learn about more topological relations, beyond `st_intersects()` and `st_disjoint()`, skip to the next section (\@ref(topological-relations)). If you're interested in the details, including other ways of subsetting, read on. -Another way of doing spatial subsetting uses objects returned by *topological operators*. -This is demonstrated in the first command below: + +Another way of doing spatial subsetting uses objects returned by topological operators. +These objects can be useful in their own right, for example when exploring the graph network of relationships between contiguous regions, but they can also be used for subsetting, as demonstrated in the code chunk below: ```{r 04-spatial-operations-6} sel_sgbp = st_intersects(x = nz_height, y = canterbury) class(sel_sgbp) +sel_sgbp sel_logical = lengths(sel_sgbp) > 0 canterbury_height2 = nz_height[sel_logical, ] ``` @@ -125,7 +125,7 @@ Note: another way to return a logical output is by setting `sparse = FALSE` (mea Note: the solution involving `sgbp` objects is more generalisable though, as it works for many-to-many operations and has lower memory requirements. ``` -It should be noted that a logical can also be used with `filter()` as follows (`sparse = FALSE` is explained in Section \@ref(topological-relations)): +It should be noted that a logical can also be used with `filter()` as follows (`sparse = FALSE` is explained in Section \@ref(topological-relations)): ```{r 04-spatial-operations-7b} canterbury_height3 = nz_height %>% From f4974ec3b16ff9ca9fa46b1dda10320413939684 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 07:57:33 +0000 Subject: [PATCH 055/104] Improve segway to topological relations section --- 04-spatial-operations.Rmd | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 552df3e73..324e89cc1 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -133,7 +133,8 @@ canterbury_height3 = nz_height %>% ``` At this point, there are three versions of `canterbury_height`, one created with spatial subsetting directly and the other two via intermediary selection objects. -To explore these objects and spatial subsetting in more detail, see the supplementary vignettes on `subsetting` and [`tidyverse-pitfalls`](https://geocompr.github.io/geocompkg/articles/). +To explore spatial subsetting in more detail, see the supplementary vignettes on `subsetting` and [`tidyverse-pitfalls`](https://geocompr.github.io/geocompkg/articles/) on the [geocompkg website](https://geocompr.github.io/geocompkg/articles/). +The next section explores different types of spatial relation, also known as binary predicates, that can be used to identify whether or not two features are spatially related or not. ### Topological relations From 788081a16d3398d462fcca1017dced6dc87f8843 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 08:39:19 +0000 Subject: [PATCH 056/104] Update with reference for #677 --- 04-spatial-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 324e89cc1..774ec9f24 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -138,7 +138,7 @@ The next section explores different types of spatial relation, also known as bin ### Topological relations -Topological relations describe the spatial relationships between objects. +Topological relations describe the spatial relationships between objects [@egenhofer_1990_mathematical]. To understand them, it helps to have some simple test data to work with. Figure \@ref(fig:relation-objects) contains a polygon (`a`), a line (`l`) and some points (`p`), which are created in the code below. \index{topological relations} From dead3a1872660fede3991a65f24a83c543a46dee Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 13:34:57 +0000 Subject: [PATCH 057/104] Add new refs and update bib --- geocompr.bib | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/geocompr.bib b/geocompr.bib index c5c9b98b6..e7a9302d4 100644 --- a/geocompr.bib +++ b/geocompr.bib @@ -529,6 +529,14 @@ @article{eddelbuettel_extending_2018 keywords = {\#nosource} } +@inproceedings{egenhofer_mathematical_1990, + title = {A Mathematical Framework for the Definition of Topological Relations}, + booktitle = {Proc. the Fourth International Symposium on Spatial Data Handing}, + author = {Egenhofer, Max and Herring, John}, + date = {1990}, + pages = {803--813} +} + @article{essd-11-647-2019, title = {{{ICGEM}} – 15 Years of Successful Collection and Distribution of Global Gravitational Models, Associated Services, and Future Plans}, author = {Ince, E. S. and Barthelmes, F. and Reißland, S. and Elger, K. and Förste, C. and Flechtner, F. and Schuh, H.}, @@ -676,7 +684,7 @@ @article{hengl_random_2018 } @article{hesselbarth_opensource_2021, - title = {Open-Source Tools in r for Landscape Ecology}, + title = {Open-Source Tools in {{R}} for Landscape Ecology}, author = {Hesselbarth, Maximillian H.K. and Nowosad, Jakub and Signer, Johannes and Graham, Laura J.}, date = {2021-06}, volume = {6}, @@ -1434,7 +1442,8 @@ @article{pebesma_simple_2018 @book{pebesma_spatial_2022, title = {Spatial {{Data Science}} with Applications in {{R}}}, author = {Pebesma, Edzer and Bivand, Roger}, - date = {2022} + date = {2022}, + url = {https://r-spatial.org/book} } @book{perpinan_rastervis_2016, @@ -1700,6 +1709,19 @@ @book{tomlin_geographic_1990 keywords = {\#nosource,Cartography,Data processing,Geographic information systems} } +@article{tomlin_map_1994, + title = {Map Algebra: One Perspective}, + shorttitle = {Map Algebra}, + author = {Tomlin, C. Dana}, + date = {1994}, + journaltitle = {Landscape and Urban Planning}, + volume = {30}, + number = {1-2}, + pages = {3--12}, + publisher = {{Elsevier}}, + doi = {10/dm2qm2} +} + @book{usgs_geological_2016, title = {U.{{S}}. {{Geological Survey}} ({{USGS}}) {{Earth Resources Observation}} and {{Science}} ({{EROS}}) {{Center}}}, author = {{USGS}}, From bc0288bcae65c5154784394f609f7218c2e581d2 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 14:37:36 +0000 Subject: [PATCH 058/104] Add new references on algebraic topology --- geocompr.bib | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/geocompr.bib b/geocompr.bib index e7a9302d4..dd9162031 100644 --- a/geocompr.bib +++ b/geocompr.bib @@ -472,6 +472,20 @@ @article{coppock_history_1991 keywords = {\#nosource,⛔ No DOI found,History of GIS} } +@book{dieck_algebraic_2008, + title = {Algebraic Topology}, + author = {tom Dieck, Tammo}, + date = {2008}, + series = {{{EMS}} Textbooks in Mathematics}, + publisher = {{European Mathematical Society}}, + location = {{Zürich}}, + url = {https://www.maths.ed.ac.uk/~v1ranick/papers/diecktop.pdf}, + isbn = {978-3-03719-048-7}, + pagetotal = {567}, + keywords = {Algebraic topology,Homology theory,Homotopy theory}, + annotation = {OCLC: ocn261176011} +} + @book{diggle_modelbased_2007, title = {Model-Based Geostatistics}, author = {Diggle, Peter and Ribeiro, Paulo Justiniano}, @@ -1594,6 +1608,17 @@ @book{sherman_desktop_2008 keywords = {\#nosource} } +@book{spanier_algebraic_1995, + title = {Algebraic Topology}, + author = {Spanier, Edwin Henry}, + date = {1995}, + edition = {1. corr. Springer ed}, + publisher = {{Springer}}, + location = {{New York Berlin Barcelona Budapest}}, + isbn = {978-0-387-94426-5 978-3-540-90646-9 978-0-387-90646-1}, + pagetotal = {528} +} + @book{talbert_ancient_2014, title = {Ancient {{Perspectives}}: Maps and {{Their Place}} in {{Mesopotamia}}, {{Egypt}}, {{Greece}}, and {{Rome}}}, shorttitle = {Ancient {{Perspectives}}}, From 7546fdceb21be49450e017a3005e44c6da149efa Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 14:40:02 +0000 Subject: [PATCH 059/104] Update text with refs for #677 --- 04-spatial-operations.Rmd | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 774ec9f24..0f0a2f311 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -138,7 +138,14 @@ The next section explores different types of spatial relation, also known as bin ### Topological relations -Topological relations describe the spatial relationships between objects [@egenhofer_1990_mathematical]. +Topological relations describe the spatial relationships between objects. +Specifically, "binary topological relationships" are logical statements (in that the answer can only be `TRUE` or `FALSE`) about the way that discrete regions of space defined by ordered sets of points (typically forming points, lines and polygons) relate to each other in two or more dimensions [@egenhofer_1990_mathematical]. +That may sound rather abstract and, indeed, the definition and classification of topological relations is based on earlier mathematical foundations first published as a book in 1966 [@spanier_algebraic_1995].^[ +See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. +] +However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for different types of topological relationships, as illustrated in Figure \@ref(fig:relation-objects). +In `sf`, functions testing for different types of topological relations are as 'binary predicates', as described in the vignette *Manipulating Simple Feature Geometries*, which can be viewed with the command [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), and in the help page [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) [@pebesma_simple_2018]. + To understand them, it helps to have some simple test data to work with. Figure \@ref(fig:relation-objects) contains a polygon (`a`), a line (`l`) and some points (`p`), which are created in the code below. \index{topological relations} From 0548e7eb7509ee08ae5a96fc67b11081ad1f5d6f Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 21:02:33 +0000 Subject: [PATCH 060/104] Test ggplot2 code for #677 --- 04-spatial-operations.Rmd | 49 ++++++++++++++++++++++++++++++++++----- 1 file changed, 43 insertions(+), 6 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 0f0a2f311..9ba950a24 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -146,16 +146,53 @@ See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for different types of topological relationships, as illustrated in Figure \@ref(fig:relation-objects). In `sf`, functions testing for different types of topological relations are as 'binary predicates', as described in the vignette *Manipulating Simple Feature Geometries*, which can be viewed with the command [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), and in the help page [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) [@pebesma_simple_2018]. -To understand them, it helps to have some simple test data to work with. -Figure \@ref(fig:relation-objects) contains a polygon (`a`), a line (`l`) and some points (`p`), which are created in the code below. +Figure \@ref(fig:relation-objects) + \index{topological relations} -```{r 04-spatial-operations-8} +```{r 04-spatial-operations-8, echo=FALSE} +library(ggplot2) # create a polygon -a_poly = st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1)))) -a = st_sfc(a_poly) +p1 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) +p2 = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) +p3 = st_sfc(st_polygon(list(rbind(c(0.7, 0.3), c(0.7, 0.95), c(0.9, 0.95), c(0.7, 0.3))))) +p4 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(0.9, 0.5), c(0.6, 0.1))))) +p5 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) + +# convert to sf objects +psf1 = st_sf(data.frame(Object = "Polygon1"), geometry = p1) +psf2 = st_sf(data.frame(Object = "Polygon2"), geometry = p2) +psf3 = st_sf(data.frame(Object = "Polygon3"), geometry = p3) +psf4 = st_sf(data.frame(Object = "Polygon4"), geometry = p4) +psf5 = st_sf(data.frame(Object = "Polygon5"), geometry = p5) + +# plot them +ps1 = rbind(psf1, psf2) +ps2 = rbind(psf1, psf3) +# sf::st_intersects(psf1, psf2, sparse = FALSE) +# sf::st_disjoint(psf1, psf2, sparse = FALSE) +# sf::st_overlaps(psf1, psf2, sparse = FALSE) +# sf::st_relate(psf1, psf2, sparse = FALSE) +# sf::st_relate(psf1, psf3) +# tm_shape(ps1) + tm_polygons(col = "Object") + tm_layout(legend.position = c("left", "top")) + +theme_set(new = theme_void()) +g1 = ggplot(ps1) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) +# g1 + annotate("text", x = 0.3, y = 0.9, label = "st_intersects(Polygon1, Polygon2)") +g1 + annotate("text", x = 0.1, y = 0.95, label = "intersects TRUE\ndisjoint FALSE\ntouches TRUE\n", hjust = "left", vjust = "top") +# Try annotating only which type of relations apply +# g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches", hjust = "left", vjust = "top") +g1an = g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches\nDE-9IM: FF2F11212", hjust = "left", vjust = "top") + +g2 = ggplot(ps2) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) +g2an = g2 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects,\ntouches, overlaps\nDE-9IM: 212101212", hjust = "left", vjust = "top") + +library(patchwork) +g1an + g2an + + # create a line -l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 2)) +l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 3)) l = st_sfc(l_line) # create points p_matrix = matrix(c(0.5, 1, -1, 0, 0, 1, 0.5, 1), ncol = 2) From 2c72f4ac4b3d7ab9a3a2663c2289352ade488ed3 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sat, 13 Nov 2021 22:15:34 +0000 Subject: [PATCH 061/104] Add placeholder section on DE-9IM #677 --- 04-spatial-operations.Rmd | 18 ++++++++++++++++-- geocompr.bib | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 48 insertions(+), 2 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 9ba950a24..c1ba199be 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -139,11 +139,13 @@ The next section explores different types of spatial relation, also known as bin ### Topological relations Topological relations describe the spatial relationships between objects. -Specifically, "binary topological relationships" are logical statements (in that the answer can only be `TRUE` or `FALSE`) about the way that discrete regions of space defined by ordered sets of points (typically forming points, lines and polygons) relate to each other in two or more dimensions [@egenhofer_1990_mathematical]. +"Binary topological relationships", to give them their full name, are logical statements (in that the answer can only be `TRUE` or `FALSE`) about the spatial relationships between two objects defined by ordered sets of points (typically forming points, lines and polygons) in two or more dimensions [@egenhofer_1990_mathematical]. That may sound rather abstract and, indeed, the definition and classification of topological relations is based on earlier mathematical foundations first published as a book in 1966 [@spanier_algebraic_1995].^[ See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. ] -However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for different types of topological relationships, as illustrated in Figure \@ref(fig:relation-objects). +However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships defined in English (and other human) languages as illustrated in Figure \@ref(fig:relations). + + In `sf`, functions testing for different types of topological relations are as 'binary predicates', as described in the vignette *Manipulating Simple Feature Geometries*, which can be viewed with the command [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), and in the help page [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) [@pebesma_simple_2018]. Figure \@ref(fig:relation-objects) @@ -310,6 +312,18 @@ plot(l, add = TRUE) plot(p, add = TRUE) ``` +### DE-9IM strings + +Underlying the binary predicates demonstrated in the previous section is the Dimensionally Extended 9-Intersection Model (DE-9IM). +The model was originally labelled "DE + 9IM" by its inventors, referring to the "dimension of the intersections of' boundaries, interiors, and exteriors of two features" [@clementini_1995_comparison] + +The way the model works can be demonstrated with reference to two intersecting polygons, as illustrated in Figure \@ref(fig:de-9im). + +```{r} + +``` + + ### Spatial joining Joining two non-spatial datasets relies on a shared 'key' variable, as described in Section \@ref(vector-attribute-joining). diff --git a/geocompr.bib b/geocompr.bib index dd9162031..1f6a7e1d9 100644 --- a/geocompr.bib +++ b/geocompr.bib @@ -429,6 +429,22 @@ @incollection{cheshire_spatial_2015 keywords = {\#nosource} } +@article{clementini_comparison_1995a, + title = {A Comparison of Methods for Representing Topological Relationships}, + author = {Clementini, Eliseo and Di Felice, Paolino}, + date = {1995-05-01}, + journaltitle = {Information Sciences - Applications}, + shortjournal = {Information Sciences - Applications}, + volume = {3}, + number = {3}, + pages = {149--178}, + issn = {1069-0115}, + doi = {10.1016/1069-0115(94)00033-X}, + url = {https://www.sciencedirect.com/science/article/pii/106901159400033X}, + urldate = {2021-11-13}, + langid = {english} +} + @article{conrad_system_2015, title = {System for {{Automated Geoscientific Analyses}} ({{SAGA}}) v. 2.1.4}, author = {Conrad, O. and Bechtel, B. and Bock, M. and Dietrich, H. and Fischer, E. and Gerlitz, L. and Wehberg, J. and Wichmann, V. and Böhner, J.}, @@ -1600,6 +1616,22 @@ @article{schratz_performance_nodate keywords = {\#nosource,⛔ No DOI found,Computer Science - Machine Learning,Statistics - Machine Learning,Statistics - Methodology} } +@article{shen_classification_2018, + title = {Classification of Topological Relations between Spatial Objects in Two-Dimensional Space within the Dimensionally Extended 9-Intersection Model}, + author = {Shen, Jingwei and Chen, Min and Liu, Xintao}, + date = {2018}, + journaltitle = {Transactions in GIS}, + volume = {22}, + number = {2}, + pages = {514--541}, + issn = {1467-9671}, + doi = {10.1111/tgis.12328}, + url = {https://onlinelibrary.wiley.com/doi/abs/10.1111/tgis.12328}, + urldate = {2021-11-13}, + langid = {english}, + annotation = {\_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/tgis.12328} +} + @book{sherman_desktop_2008, title = {Desktop {{GIS}}: Mapping the {{Planet}} with {{Open Source Tools}}}, author = {Sherman, Gary}, From 209556413cf2d0047c7bf15b926565785cb9a780 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sun, 14 Nov 2021 11:09:12 +0000 Subject: [PATCH 062/104] Tweak comment and replace code with link to terra issue in c2 --- 02-spatial-data.Rmd | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/02-spatial-data.Rmd b/02-spatial-data.Rmd index 258ac9614..e3fd940f5 100644 --- a/02-spatial-data.Rmd +++ b/02-spatial-data.Rmd @@ -730,14 +730,11 @@ Like its predecessor **raster** (created by the same developer, Robert Hijmans), On the other hand, the **raster** class system is popular and used by many other packages; as with **sf** and **sp**, the good news is you can seamlessly translate between the two types of object to ensure backwards compatibility with older scripts and packages, for example, with the functions `as.raster(my_rast)` (see the previous chapter for more on the evolution of R packages for working with geographic data). ```{r, echo=FALSE, eval=FALSE} -# test raster/terra conversions -my_raster = as.raster(my_rast) -library(raster) -my_rast2 = rast(my_raster) # how to convert back to terra? fails for me (RL 2021-11) -identical(my_rast, my_rast2) +# # test raster/terra conversions +# See https://github.com/rspatial/terra/issues/399 ``` -Aside from general raster data manipulation, **terra** provides many low-level functions that can form the basis to develop more advanced raster functionality. +In addition to functions for raster data manipulation, **terra** provides many low-level functions that can form a foundation for developing new tools for working with raster datasets. \index{terra (package)|see {terra}} **terra** also lets you work on large raster datasets that are too large to fit into the main memory. In this case, **terra** provides the possibility to divide the raster into smaller chunks, and processes these iteratively instead of loading the whole raster file into RAM. From 886f98ae2f6aaec79d4147a7e127af132a4aacd6 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Sun, 14 Nov 2021 11:09:36 +0000 Subject: [PATCH 063/104] Test code for #677 --- 04-spatial-operations.Rmd | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index c1ba199be..aaa7369db 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -320,7 +320,35 @@ The model was originally labelled "DE + 9IM" by its inventors, referring to the The way the model works can be demonstrated with reference to two intersecting polygons, as illustrated in Figure \@ref(fig:de-9im). ```{r} +b = st_sfc(st_point(c(0, 1)), st_point(c(1, 1))) # create 2 points +b = st_buffer(b, dist = 1) # convert points to circles +bsf = sf::st_sf(data.frame(Object = c("a", "b")), geometry = b) +b9 = replicate(bsf, n = 9, simplify = FALSE) +b9sf = do.call(rbind, b9) +domains = c("Interior", "Boundary", "Exterior") +b9sf$domain_a = rep(rep(domains, 3), each = 2) +b9sf$domain_b = rep(rep(domains, each = 3), each = 2) +library(ggplot2) +ggplot(b9sf) + + geom_sf() + + facet_grid(domain_a ~ domain_b) + +plot(b9sf) +tmap_arrange( + tm_shape(b) + tm_polygons(alpha = 0.5) + tm_layout(title = "Interior-Interior"), + tm_shape(b) + tm_polygons(alpha = 0.5) + tm_layout(title = "Interior-Boundary"), + tm_shape(b) + tm_polygons(alpha = 0.5), + tm_shape(b) + tm_polygons(alpha = 0.5), + tm_shape(b) + tm_polygons(alpha = 0.5), + tm_shape(b) + tm_polygons(alpha = 0.5), + tm_shape(b) + tm_polygons(alpha = 0.5), + tm_shape(b) + tm_polygons(alpha = 0.5), + tm_shape(b) + tm_polygons(alpha = 0.5), + nrow = 3 +) +plot(b) +text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y")) # add text ``` From ee22ec24538acfc933402422552d456f45d56af7 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 20:02:35 +0000 Subject: [PATCH 064/104] Add start fun for #677 --- code/de_9im.R | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 code/de_9im.R diff --git a/code/de_9im.R b/code/de_9im.R new file mode 100644 index 000000000..7881dd6f3 --- /dev/null +++ b/code/de_9im.R @@ -0,0 +1,18 @@ +#' This function visualises sf objects and returns info on the +#' types of spatial relationship there is between them +#' +#' Context: [robinlovelace/geocompr#677](https://github.com/Robinlovelace/geocompr/issues/677) +#' +#' @examples +#' library(sf) +#' x = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) +#' y = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) +#' de_9im +de_9im = function(x, y, object_names = c("x", "y")) { + if(is(x, "sfc") && is(y, "sfc")) { + x = st_sf(data.frame(Object = object_names[1]), geometry = x) + y = st_sf(data.frame(Object = object_names[2]), geometry = y) + } + browser() + ps1 = rbind(psf1, psf2) +} \ No newline at end of file From 3f3206ec6eb25909949d17431c43773528c36184 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 20:16:45 +0000 Subject: [PATCH 065/104] MVP for de_9im fun --- code/de_9im.R | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/code/de_9im.R b/code/de_9im.R index 7881dd6f3..a1b491e6c 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -7,12 +7,17 @@ #' library(sf) #' x = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) #' y = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) -#' de_9im -de_9im = function(x, y, object_names = c("x", "y")) { +#' de_9im(x, y) +de_9im = function(x, y, object_names = c("x", "y"), funs = list("st_intersects", "st_disjoint"), sparse = FALSE) { + requireNamespace("sf", quietly = TRUE) if(is(x, "sfc") && is(y, "sfc")) { x = st_sf(data.frame(Object = object_names[1]), geometry = x) y = st_sf(data.frame(Object = object_names[2]), geometry = y) } - browser() - ps1 = rbind(psf1, psf2) -} \ No newline at end of file + xy = rbind(x, y) + funs_matched = lapply(funs, match.fun) + res = lapply(seq(length(funs)), function(i) { + funs_matched[[i]](x, y, sparse = sparse) + }) + unlist(res) +} From 596aa447465afee5ac97e3cea643258f43672ca0 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 20:18:59 +0000 Subject: [PATCH 066/104] Move test code to function script, fix build --- 04-spatial-operations.Rmd | 33 --------------------------------- code/de_9im.R | 12 ++++++++++++ 2 files changed, 12 insertions(+), 33 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index aaa7369db..2a49736cf 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -153,7 +153,6 @@ Figure \@ref(fig:relation-objects) \index{topological relations} ```{r 04-spatial-operations-8, echo=FALSE} -library(ggplot2) # create a polygon p1 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) p2 = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) @@ -161,38 +160,6 @@ p3 = st_sfc(st_polygon(list(rbind(c(0.7, 0.3), c(0.7, 0.95), c(0.9, 0.95), c(0.7 p4 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(0.9, 0.5), c(0.6, 0.1))))) p5 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) -# convert to sf objects -psf1 = st_sf(data.frame(Object = "Polygon1"), geometry = p1) -psf2 = st_sf(data.frame(Object = "Polygon2"), geometry = p2) -psf3 = st_sf(data.frame(Object = "Polygon3"), geometry = p3) -psf4 = st_sf(data.frame(Object = "Polygon4"), geometry = p4) -psf5 = st_sf(data.frame(Object = "Polygon5"), geometry = p5) - -# plot them -ps1 = rbind(psf1, psf2) -ps2 = rbind(psf1, psf3) -# sf::st_intersects(psf1, psf2, sparse = FALSE) -# sf::st_disjoint(psf1, psf2, sparse = FALSE) -# sf::st_overlaps(psf1, psf2, sparse = FALSE) -# sf::st_relate(psf1, psf2, sparse = FALSE) -# sf::st_relate(psf1, psf3) -# tm_shape(ps1) + tm_polygons(col = "Object") + tm_layout(legend.position = c("left", "top")) - -theme_set(new = theme_void()) -g1 = ggplot(ps1) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) -# g1 + annotate("text", x = 0.3, y = 0.9, label = "st_intersects(Polygon1, Polygon2)") -g1 + annotate("text", x = 0.1, y = 0.95, label = "intersects TRUE\ndisjoint FALSE\ntouches TRUE\n", hjust = "left", vjust = "top") -# Try annotating only which type of relations apply -# g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches", hjust = "left", vjust = "top") -g1an = g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches\nDE-9IM: FF2F11212", hjust = "left", vjust = "top") - -g2 = ggplot(ps2) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) -g2an = g2 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects,\ntouches, overlaps\nDE-9IM: 212101212", hjust = "left", vjust = "top") - -library(patchwork) -g1an + g2an - - # create a line l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 3)) l = st_sfc(l_line) diff --git a/code/de_9im.R b/code/de_9im.R index a1b491e6c..9c15ae0bb 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -21,3 +21,15 @@ de_9im = function(x, y, object_names = c("x", "y"), funs = list("st_intersects", }) unlist(res) } + +# # Test code to functionalize: +# theme_set(new = theme_void()) +# g1 = ggplot(ps1) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) +# # g1 + annotate("text", x = 0.3, y = 0.9, label = "st_intersects(Polygon1, Polygon2)") +# g1 + annotate("text", x = 0.1, y = 0.95, label = "intersects TRUE\ndisjoint FALSE\ntouches TRUE\n", hjust = "left", vjust = "top") +# # Try annotating only which type of relations apply +# # g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches", hjust = "left", vjust = "top") +# g1an = g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches\nDE-9IM: FF2F11212", hjust = "left", vjust = "top") +# +# g2 = ggplot(ps2) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) +# g2an = g2 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects,\ntouches, overlaps\nDE-9IM: 212101212", hjust = "left", vjust = "top") From 8f1392166c460f144f834857e3df18f720f9840e Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 20:28:13 +0000 Subject: [PATCH 067/104] Working de_9im function --- code/de_9im.R | 28 +++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/code/de_9im.R b/code/de_9im.R index 9c15ae0bb..d05c425ee 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -1,16 +1,22 @@ -#' This function visualises sf objects and returns info on the +#' This function visualises sf objects and returns info on the #' types of spatial relationship there is between them -#' +#' #' Context: [robinlovelace/geocompr#677](https://github.com/Robinlovelace/geocompr/issues/677) -#' -#' @examples +#' +#' @examples #' library(sf) #' x = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) #' y = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) #' de_9im(x, y) -de_9im = function(x, y, object_names = c("x", "y"), funs = list("st_intersects", "st_disjoint"), sparse = FALSE) { +de_9im = function(x, + y, + object_names = c("x", "y"), + funs = list("st_intersects", "st_disjoint"), + sparse = FALSE, + output = "character" + ) { requireNamespace("sf", quietly = TRUE) - if(is(x, "sfc") && is(y, "sfc")) { + if (is(x, "sfc") && is(y, "sfc")) { x = st_sf(data.frame(Object = object_names[1]), geometry = x) y = st_sf(data.frame(Object = object_names[2]), geometry = y) } @@ -19,17 +25,21 @@ de_9im = function(x, y, object_names = c("x", "y"), funs = list("st_intersects", res = lapply(seq(length(funs)), function(i) { funs_matched[[i]](x, y, sparse = sparse) }) - unlist(res) + res = unlist(res) + if(output == "character") { + res = unlist(funs)[res] + } + res } # # Test code to functionalize: # theme_set(new = theme_void()) # g1 = ggplot(ps1) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) # # g1 + annotate("text", x = 0.3, y = 0.9, label = "st_intersects(Polygon1, Polygon2)") -# g1 + annotate("text", x = 0.1, y = 0.95, label = "intersects TRUE\ndisjoint FALSE\ntouches TRUE\n", hjust = "left", vjust = "top") +# g1 + annotate("text", x = 0.1, y = 0.95, label = "intersects TRUE\ndisjoint FALSE\ntouches TRUE\n", hjust = "left", vjust = "top") # # Try annotating only which type of relations apply # # g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches", hjust = "left", vjust = "top") # g1an = g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches\nDE-9IM: FF2F11212", hjust = "left", vjust = "top") -# +# # g2 = ggplot(ps2) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) # g2an = g2 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects,\ntouches, overlaps\nDE-9IM: 212101212", hjust = "left", vjust = "top") From 4569a2aa0e3c900ef4280f03098e6cc9973a625c Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 21:20:13 +0000 Subject: [PATCH 068/104] Update code, fix build mkII --- 04-spatial-operations.Rmd | 13 ++++++++++++- code/de_9im.R | 16 +++++++++++++++- 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 2a49736cf..6e8b43ee9 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -153,12 +153,16 @@ Figure \@ref(fig:relation-objects) \index{topological relations} ```{r 04-spatial-operations-8, echo=FALSE} -# create a polygon +source("code/de_9im.R") p1 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) p2 = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) p3 = st_sfc(st_polygon(list(rbind(c(0.7, 0.3), c(0.7, 0.95), c(0.9, 0.95), c(0.7, 0.3))))) p4 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(0.9, 0.5), c(0.6, 0.1))))) p5 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) +de_9im(p1, p2) +de_9im(p1, p3) +de_9im(p1, p4) +de_9im(p1, p5) # create a line l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 3)) @@ -169,6 +173,13 @@ p_multi = st_multipoint(x = p_matrix) p = st_cast(st_sfc(p_multi), "POINT") ``` +```{r} +a_poly = st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1)))) +a = st_sfc(a_poly) +l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 2)) +``` + + ```{r relation-objects, echo=FALSE, fig.cap="Points (p 1 to 4), line and polygon objects arranged to illustrate topological relations.", fig.asp=1, out.width="50%", fig.scap="Demonstration of topological relations."} par(pty = "s") plot(a, border = "red", col = "gray", axes = TRUE) diff --git a/code/de_9im.R b/code/de_9im.R index d05c425ee..f7760a39c 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -11,7 +11,21 @@ de_9im = function(x, y, object_names = c("x", "y"), - funs = list("st_intersects", "st_disjoint"), + funs = list( + "st_intersects", + "st_disjoint", + "st_touches", + "st_crosses", + "st_within", + "st_contains", + "st_contains_properly", + "st_overlaps", + "st_equals", + "st_covers", + "st_covered_by" + # , + # "st_equals_exact" # requuires par argument + ), sparse = FALSE, output = "character" ) { From 2fe134547a6e2a83ba80588b1308c97373d1aad6 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 21:51:52 +0000 Subject: [PATCH 069/104] Progress on #677 --- 04-spatial-operations.Rmd | 1 + code/de_9im.R | 25 +++++++++++++++++++++---- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 6e8b43ee9..7b79cc4dd 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -162,6 +162,7 @@ p5 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) de_9im(p1, p2) de_9im(p1, p3) de_9im(p1, p4) +de_9im(p4, p1) de_9im(p1, p5) # create a line diff --git a/code/de_9im.R b/code/de_9im.R index f7760a39c..c1c9a9d92 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -8,9 +8,14 @@ #' x = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) #' y = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) #' de_9im(x, y) +#' p3 = st_sfc(st_polygon(list(rbind(c(0.7, 0.3), c(0.7, 0.95), c(0.9, 0.95), c(0.7, 0.3))))) +#' p4 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(0.9, 0.5), c(0.6, 0.1))))) +#' p5 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) +#' de_9im(x, p3) de_9im = function(x, y, object_names = c("x", "y"), + plot = TRUE, funs = list( "st_intersects", "st_disjoint", @@ -27,7 +32,8 @@ de_9im = function(x, # "st_equals_exact" # requuires par argument ), sparse = FALSE, - output = "character" + output = "character", + collapse = "\n" ) { requireNamespace("sf", quietly = TRUE) if (is(x, "sfc") && is(y, "sfc")) { @@ -43,9 +49,22 @@ de_9im = function(x, if(output == "character") { res = unlist(funs)[res] } + if(plot) { + res_text = paste(res, collapse = collapse) + message("Object x has the following spatial relations to y: ", res_text) + res = de_9im_plot(xy, label = res_text) + } res } +de_9im_plot = function(xy, label = "test", alpha = 0.5, show.legend = FALSE, x = 0.1, y = 0.95) { + require("ggplot2", quietly = TRUE) + # browser() + ggplot(xy) + geom_sf(aes(fill = Object), alpha = alpha, show.legend = show.legend) + + annotate("text", x = 0.1, y = 0.95, label = label, hjust = "left", vjust = "top") + + ggplot2::theme_void() +} + # # Test code to functionalize: # theme_set(new = theme_void()) # g1 = ggplot(ps1) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) @@ -53,7 +72,5 @@ de_9im = function(x, # g1 + annotate("text", x = 0.1, y = 0.95, label = "intersects TRUE\ndisjoint FALSE\ntouches TRUE\n", hjust = "left", vjust = "top") # # Try annotating only which type of relations apply # # g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches", hjust = "left", vjust = "top") -# g1an = g1 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects, touches\nDE-9IM: FF2F11212", hjust = "left", vjust = "top") +# g1an = g1 + # -# g2 = ggplot(ps2) + geom_sf(aes(fill = Object), alpha = 0.5, show.legend = FALSE) -# g2an = g2 + annotate("text", x = 0.1, y = 0.95, label = "Relations: intersects,\ntouches, overlaps\nDE-9IM: 212101212", hjust = "left", vjust = "top") From 29e45872a4c210c517c020a0dc8bbf0256204868 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 21:59:00 +0000 Subject: [PATCH 070/104] Require sf --- code/de_9im.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/code/de_9im.R b/code/de_9im.R index c1c9a9d92..e3b553ce5 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -35,7 +35,7 @@ de_9im = function(x, output = "character", collapse = "\n" ) { - requireNamespace("sf", quietly = TRUE) + require("sf", quietly = TRUE) if (is(x, "sfc") && is(y, "sfc")) { x = st_sf(data.frame(Object = object_names[1]), geometry = x) y = st_sf(data.frame(Object = object_names[2]), geometry = y) From 20c292d27d558e34df4e022d98f2b3eceb744c92 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 22:18:24 +0000 Subject: [PATCH 071/104] New fig 4.2 for #677 --- 04-spatial-operations.Rmd | 12 +++++++++--- code/de_9im.R | 5 +++-- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 7b79cc4dd..0fcd56e58 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -152,19 +152,25 @@ Figure \@ref(fig:relation-objects) \index{topological relations} -```{r 04-spatial-operations-8, echo=FALSE} -source("code/de_9im.R") +```{r 04-spatial-operations-8, echo=FALSE, fig.cap="Topological relations between vector geometries", fig.show='hold', out.width="33%", message=FALSE, fig.width=3} +source("https://github.com/Robinlovelace/geocompr/raw/c4-v2-updates-rl/code/de_9im.R") +library(sf) p1 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) p2 = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) p3 = st_sfc(st_polygon(list(rbind(c(0.7, 0.3), c(0.7, 0.95), c(0.9, 0.95), c(0.7, 0.3))))) p4 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(0.9, 0.5), c(0.6, 0.1))))) -p5 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) +p5 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(1, 0.5), c(0.6, 0.1))))) +p6 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) de_9im(p1, p2) de_9im(p1, p3) de_9im(p1, p4) de_9im(p4, p1) de_9im(p1, p5) +de_9im(p1, p6) +``` + +```{r} # create a line l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 3)) l = st_sfc(l_line) diff --git a/code/de_9im.R b/code/de_9im.R index e3b553ce5..f9a61f1b7 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -27,7 +27,8 @@ de_9im = function(x, "st_overlaps", "st_equals", "st_covers", - "st_covered_by" + "st_covered_by", + ... # , # "st_equals_exact" # requuires par argument ), @@ -35,7 +36,7 @@ de_9im = function(x, output = "character", collapse = "\n" ) { - require("sf", quietly = TRUE) + require("sf") if (is(x, "sfc") && is(y, "sfc")) { x = st_sf(data.frame(Object = object_names[1]), geometry = x) y = st_sf(data.frame(Object = object_names[2]), geometry = y) From 90f5f2b2ad7528d15c0160c54f1f201842806df8 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 23:04:52 +0000 Subject: [PATCH 072/104] Fix polygons, improve vis for #677 --- 04-spatial-operations.Rmd | 21 ++++++++++++--------- code/de_9im.R | 15 ++++++++++----- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 0fcd56e58..fbfc8068f 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -152,19 +152,22 @@ Figure \@ref(fig:relation-objects) \index{topological relations} -```{r 04-spatial-operations-8, echo=FALSE, fig.cap="Topological relations between vector geometries", fig.show='hold', out.width="33%", message=FALSE, fig.width=3} -source("https://github.com/Robinlovelace/geocompr/raw/c4-v2-updates-rl/code/de_9im.R") +```{r 04-spatial-operations-8, echo=FALSE, fig.cap="Topological relations between vector geometries", fig.show='hold', out.width="33%", message=FALSE, fig.width=3, fig.height=3} +# source("https://github.com/Robinlovelace/geocompr/raw/c4-v2-updates-rl/code/de_9im.R") +source("code/de_9im.R") library(sf) -p1 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 0))))) -p2 = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(0, 0))))) -p3 = st_sfc(st_polygon(list(rbind(c(0.7, 0.3), c(0.7, 0.95), c(0.9, 0.95), c(0.7, 0.3))))) -p4 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(0.9, 0.5), c(0.6, 0.1))))) -p5 = st_sfc(st_polygon(list(rbind(c(0.6, 0.1), c(0.7, 0.5), c(1, 0.5), c(0.6, 0.1))))) -p6 = st_sfc(st_polygon(list(rbind(c(0, 0.2), c(0, 1), c(0.9, 1), c(0, 0.2))))) +p1 = st_sfc(st_polygon(list(rbind(c(0, 0), c(0, 1), c(1, 1), c(1, 0.5), c(0, 0))))) +p2 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 0.5), c(0, 0))))) +p3 = st_sfc(st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 0.7), c(0, 0))))) +p4 = st_sfc(st_polygon(list(rbind(c(0.7, 0.8), c(0.7, 0.5), c(0.9, 0.5), c(0.7, 0.8))))) +p5 = st_sfc(st_polygon(list(rbind(c(0.6, 0.7), c(0.7, 0.5), c(1, 0.5), c(0.6, 0.7))))) +p6 = st_sfc(st_polygon(list(rbind(c(0.1, 0), c(1, 0), c(1, 0.3), c(0.1, 0))))) +p7 = st_sfc(st_polygon(list(rbind(c(0.05, 0.4), c(0.05, 0.97), c(0.6, 0.97), c(0.5, 0.4), c(0.05, 0.4))))) + de_9im(p1, p2) de_9im(p1, p3) de_9im(p1, p4) -de_9im(p4, p1) +de_9im(p7, p1) de_9im(p1, p5) de_9im(p1, p6) ``` diff --git a/code/de_9im.R b/code/de_9im.R index f9a61f1b7..d27f52aad 100644 --- a/code/de_9im.R +++ b/code/de_9im.R @@ -27,14 +27,14 @@ de_9im = function(x, "st_overlaps", "st_equals", "st_covers", - "st_covered_by", - ... + "st_covered_by" # , # "st_equals_exact" # requuires par argument ), + include_relate = TRUE, sparse = FALSE, output = "character", - collapse = "\n" + collapse = " ✓\n" ) { require("sf") if (is(x, "sfc") && is(y, "sfc")) { @@ -50,6 +50,11 @@ de_9im = function(x, if(output == "character") { res = unlist(funs)[res] } + if(include_relate) { + relation = sf::st_relate(x, y) + relate_text = paste0(" \nDE-9IM string: \n", relation) + res = c(res, relate_text) + } if(plot) { res_text = paste(res, collapse = collapse) message("Object x has the following spatial relations to y: ", res_text) @@ -58,12 +63,12 @@ de_9im = function(x, res } -de_9im_plot = function(xy, label = "test", alpha = 0.5, show.legend = FALSE, x = 0.1, y = 0.95) { +de_9im_plot = function(xy, label = "test", alpha = 0.5, show.legend = FALSE, x = 0.1, y = 0.95, theme = ggplot2::theme_void()) { require("ggplot2", quietly = TRUE) # browser() ggplot(xy) + geom_sf(aes(fill = Object), alpha = alpha, show.legend = show.legend) + annotate("text", x = 0.1, y = 0.95, label = label, hjust = "left", vjust = "top") + - ggplot2::theme_void() + theme } # # Test code to functionalize: From 60481d723a04d9fb0c3c5acfd6c72553b0325ef8 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 15 Nov 2021 23:58:38 +0000 Subject: [PATCH 073/104] Finish updating s. 4.2.2, close #677 --- 04-spatial-operations.Rmd | 96 +++++++++++++++++++++++---------------- 1 file changed, 57 insertions(+), 39 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index fbfc8068f..09761b643 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -143,16 +143,11 @@ Topological relations describe the spatial relationships between objects. That may sound rather abstract and, indeed, the definition and classification of topological relations is based on earlier mathematical foundations first published as a book in 1966 [@spanier_algebraic_1995].^[ See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. ] -However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships defined in English (and other human) languages as illustrated in Figure \@ref(fig:relations). - - -In `sf`, functions testing for different types of topological relations are as 'binary predicates', as described in the vignette *Manipulating Simple Feature Geometries*, which can be viewed with the command [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), and in the help page [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) [@pebesma_simple_2018]. - -Figure \@ref(fig:relation-objects) - +However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships defined in English (and other human languages) as illustrated in Figure \@ref(fig:relations), which shows a variety of geometry pairs and their associated relations. +The third and fourth pairs in Figure \@ref(fig:relations) (from left to right and then down) demonstrate that, for some relations, order is important: while the relations *equals*, *intersects*, *crosses*, *touches* and *overlaps* are symmetrical, meaning that if `function(x, y)` is true, `function(y, x)` will also by true, relations in which the order of the geomtries are important such as *contains* and *within* are not. \index{topological relations} -```{r 04-spatial-operations-8, echo=FALSE, fig.cap="Topological relations between vector geometries", fig.show='hold', out.width="33%", message=FALSE, fig.width=3, fig.height=3} +```{r relations, echo=FALSE, fig.cap="Topological relations between vector geometries", fig.show='hold', out.width="33%", message=FALSE, fig.width=3, fig.height=3} # source("https://github.com/Robinlovelace/geocompr/raw/c4-v2-updates-rl/code/de_9im.R") source("code/de_9im.R") library(sf) @@ -170,56 +165,79 @@ de_9im(p1, p4) de_9im(p7, p1) de_9im(p1, p5) de_9im(p1, p6) +# todo: add 3 more with line/point relations? ``` +In `sf`, functions testing for different types of topological relations are called 'binary predicates', as described in the vignette *Manipulating Simple Feature Geometries*, which can be viewed with the command [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), and in the help page [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) [@pebesma_simple_2018]. +To see how topological relations work in practice, let's create a simple reproducible example, building on the relations illustrated in Figure \@ref(fig:relations). -```{r} -# create a line -l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 3)) -l = st_sfc(l_line) -# create points -p_matrix = matrix(c(0.5, 1, -1, 0, 0, 1, 0.5, 1), ncol = 2) -p_multi = st_multipoint(x = p_matrix) -p = st_cast(st_sfc(p_multi), "POINT") +```{r, echo=FALSE} +# alternative way of getting reader to data, less clear... +# polygon_matrix = matrix(c( +# 0, 0.0, +# 0, 1.0, +# 1, 1.0, +# 1, 0.5, +# 0, 0.0 +# ), +# ncol = 2, +# byrow = TRUE +# ) ``` + ```{r} -a_poly = st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1)))) -a = st_sfc(a_poly) -l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 2)) +polygon_df = read.csv(text = "x, y +0, 0.0 +0, 1.0 +1, 1.0 +1, 0.5 +0, 0.0") +polygon_matrix = as.matrix(polygon_df) +polygon = st_polygon(list(polygon_matrix)) +polygon_sfc = st_sfc(polygon) +line_df = read.csv(text = "x,y +0.1, 0 +1, 0.1") +line = st_linestring(x = as.matrix(line_df)) +line_sfc = st_sfc(line) +# create points +point_df = read.csv(text = "x,y +0.1,0 +0.7,0.2 +0.4,0.8") +point_sf = st_as_sf(point_df, coords = c("x", "y")) ``` - ```{r relation-objects, echo=FALSE, fig.cap="Points (p 1 to 4), line and polygon objects arranged to illustrate topological relations.", fig.asp=1, out.width="50%", fig.scap="Demonstration of topological relations."} par(pty = "s") -plot(a, border = "red", col = "gray", axes = TRUE) -plot(l, add = TRUE) -plot(p, add = TRUE, lab = 1:4) -text(p_matrix[, 1] + 0.04, p_matrix[, 2] - 0.06, 1:4, cex = 1.3) +plot(polygon_sfc, border = "red", col = "gray", axes = TRUE) +plot(line_sfc, add = TRUE) +plot(point_sf, add = TRUE, lab = 1:4) +text(point_df[, 1] + 0.02, point_df[, 2] + 0.03, 1:3, cex = 1.3) ``` -A simple query is: which of the points in `p` intersect in some way with polygon `a`? +A simple query is: which of the points in `point_sf` intersect in some way with polygon `polygon_sfc`? The question can be answered by inspection (points 1 and 2 are over or touch the triangle). It can also be answered by using a *spatial predicate* such as *do the objects intersect*? This is implemented in **sf** as follows: ```{r 04-spatial-operations-9, eval=FALSE} -st_intersects(p, a) +st_intersects(point_sf, polygon_sfc) #> Sparse geometry binary ..., where the predicate was `intersects' -#> 1: 1 -#> 2: 1 -#> 3: (empty) -#> 4: (empty) +#> 1: (empty) +#> 2: (empty) +#> 3: 1 ``` The contents of the result should be as you expected: -the function returns a positive (`1`) result for the first two points, and a negative result (represented by an empty vector) for the last two. +the function returns a positive (`1`) for the third point, and a negative result (represented by an empty vector) for the first two which are outside the polygon's border. What may be unexpected is that the result comes in the form of a list of vectors. This *sparse matrix* output only registers a relation if one exists, reducing the memory requirements of topological operations on multi-feature objects. As we saw in the previous section, a *dense matrix* consisting of `TRUE` or `FALSE` values for each combination of features can also be returned when `sparse = FALSE`: ```{r 04-spatial-operations-10} -st_intersects(p, a, sparse = FALSE) +st_intersects(point_sf, polygon_sfc, sparse = FALSE) ``` The output is a matrix in which each row represents a feature in the target object and each column represents a feature in the selecting object. @@ -231,21 +249,21 @@ Note that `st_intersects()` returns `TRUE` for the second feature in the object The opposite of `st_intersects()` is `st_disjoint()`, which returns only objects that do not spatially relate in any way to the selecting object (note `[, 1]` converts the result into a vector): ```{r 04-spatial-operations-11} -st_disjoint(p, a, sparse = FALSE)[, 1] +st_disjoint(point_sf, polygon_sfc, sparse = FALSE)[, 1] ``` `st_within()` returns `TRUE` only for objects that are completely within the selecting object. This applies only to the first object, which is inside the triangular polygon, as illustrated below: ```{r 04-spatial-operations-12} -st_within(p, a, sparse = FALSE)[, 1] +st_within(point_sf, polygon_sfc, sparse = FALSE)[, 1] ``` Note that although the first point is *within* the triangle, it does not *touch* any part of its border. For this reason `st_touches()` only returns `TRUE` for the second point: ```{r 04-spatial-operations-13} -st_touches(p, a, sparse = FALSE)[, 1] +st_touches(point_sf, polygon_sfc, sparse = FALSE)[, 1] ``` What about features that do not touch, but *almost touch* the selection object? @@ -255,7 +273,7 @@ Note that although point 4 is one unit of distance from the nearest node of `a` This is illustrated in the code chunk below, the second line of which converts the lengthy list output into a `logical` object: ```{r 04-spatial-operations-14} -sel = st_is_within_distance(p, a, dist = 0.9) # can only return a sparse matrix +sel = st_is_within_distance(point_sf, polygon_sfc, dist = 0.9) # can only return a sparse matrix lengths(sel) > 0 ``` @@ -269,9 +287,9 @@ You can learn more at https://www.r-spatial.org/r/2017/06/22/spatial-index.html. ```{r 04-spatial-operations-16, eval=FALSE, echo=FALSE} # other tests -st_overlaps(p, a, sparse = FALSE) -st_covers(p, a, sparse = FALSE) -st_covered_by(p, a, sparse = FALSE) +st_overlaps(point_sf, polygon_sfc, sparse = FALSE) +st_covers(point_sf, polygon_sfc, sparse = FALSE) +st_covered_by(point_sf, polygon_sfc, sparse = FALSE) ``` ```{r 04-spatial-operations-17, eval=FALSE, echo=FALSE} From 2d21ff82618dd0fb08f3ce19c3d56aa2f7fafba2 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 00:39:26 +0000 Subject: [PATCH 074/104] Add chunk to get around failing terra code --- 04-spatial-operations.Rmd | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 09761b643..26e68c01a 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -611,6 +611,11 @@ For the reader's convenience, these datasets can be also found in the **spData** ### Spatial subsetting {#spatial-raster-subsetting} +```{r} +knitr::opts_chunk$set(eval = FALSE) # terra code failing locally +``` + + The previous chapter (Section \@ref(manipulating-raster-objects)) demonstrated how to retrieve values associated with specific cell IDs or row and column combinations. Raster objects can also be extracted by location (coordinates) and other spatial objects. To use coordinates for subsetting, one can 'translate' the coordinates into a cell ID with the **terra** function `cellFromXY()`. From d5cd6c9741e1cbda7eaec8c4ed40e29da7465181 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 00:49:02 +0000 Subject: [PATCH 075/104] Fix references in c4 --- 04-spatial-operations.Rmd | 15 +++++---------- geocompr.bib | 39 ++++++++++++++++----------------------- 2 files changed, 21 insertions(+), 33 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 26e68c01a..3aeab48fe 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -139,7 +139,7 @@ The next section explores different types of spatial relation, also known as bin ### Topological relations Topological relations describe the spatial relationships between objects. -"Binary topological relationships", to give them their full name, are logical statements (in that the answer can only be `TRUE` or `FALSE`) about the spatial relationships between two objects defined by ordered sets of points (typically forming points, lines and polygons) in two or more dimensions [@egenhofer_1990_mathematical]. +"Binary topological relationships", to give them their full name, are logical statements (in that the answer can only be `TRUE` or `FALSE`) about the spatial relationships between two objects defined by ordered sets of points (typically forming points, lines and polygons) in two or more dimensions [@egenhofer_mathematical_1990]. That may sound rather abstract and, indeed, the definition and classification of topological relations is based on earlier mathematical foundations first published as a book in 1966 [@spanier_algebraic_1995].^[ See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. ] @@ -321,7 +321,7 @@ plot(p, add = TRUE) ### DE-9IM strings Underlying the binary predicates demonstrated in the previous section is the Dimensionally Extended 9-Intersection Model (DE-9IM). -The model was originally labelled "DE + 9IM" by its inventors, referring to the "dimension of the intersections of' boundaries, interiors, and exteriors of two features" [@clementini_1995_comparison] +The model was originally labelled "DE + 9IM" by its inventors, referring to the "dimension of the intersections of' boundaries, interiors, and exteriors of two features" [@clementini_comparison_1995] The way the model works can be demonstrated with reference to two intersecting polygons, as illustrated in Figure \@ref(fig:de-9im). @@ -611,11 +611,6 @@ For the reader's convenience, these datasets can be also found in the **spData** ### Spatial subsetting {#spatial-raster-subsetting} -```{r} -knitr::opts_chunk$set(eval = FALSE) # terra code failing locally -``` - - The previous chapter (Section \@ref(manipulating-raster-objects)) demonstrated how to retrieve values associated with specific cell IDs or row and column combinations. Raster objects can also be extracted by location (coordinates) and other spatial objects. To use coordinates for subsetting, one can 'translate' the coordinates into a cell ID with the **terra** function `cellFromXY()`. @@ -636,7 +631,7 @@ terra::extract(elev, matrix(c(0.1, 0.1), ncol = 2)) Raster objects can also be subset with another raster object, as demonstrated in the code chunk below: -```{r 04-spatial-operations-35} +```{r 04-spatial-operations-35, eval=FALSE} clip = rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45, resolution = 0.3, vals = rep(1, 9)) elev[clip] @@ -701,7 +696,7 @@ The next subsection explores these and related operations in more detail. ### Map algebra \index{map algebra} -The term 'map algebra' was coined in the late 1970s to describe a "set of conventions, capabilities, and techniques" for the analysis of geographic raster *and* (although less prominently) vector data [@tomlin_1994_map]. +The term 'map algebra' was coined in the late 1970s to describe a "set of conventions, capabilities, and techniques" for the analysis of geographic raster *and* (although less prominently) vector data [@tomlin_map_1994]. Although the concept never became widely adopted, the term usefully encapsulates and helps classify the range operations that can be undertaken on raster datasets. In this context we define map algebra more narrowly, as operations that modify or summarise raster cell values, with reference to surrounding cells, zones, or statistical functions that apply to every cell. @@ -872,7 +867,7 @@ This is in contrast to focal operations which return a raster object (see previo For example, to find the mean elevation for each grain size class (Figure \@ref(fig:cont-raster)), we use the `zonal()` function. -```{r 04-spatial-operations-43} +```{r 04-spatial-operations-43, eval=FALSE} z = zonal(elev, grain, fun = "mean") z ``` diff --git a/geocompr.bib b/geocompr.bib index 1f6a7e1d9..be2214d87 100644 --- a/geocompr.bib +++ b/geocompr.bib @@ -156,29 +156,18 @@ @inproceedings{bivand_open_2000 keywords = {\#nosource,⛔ No DOI found} } -@article{bivand_progress_2020a, +@article{bivand_progress_2020, + ids = {bivand_progress_2020a,bivand_progress_2021}, title = {Progress in the {{R}} Ecosystem for Representing and Handling Spatial Data}, author = {Bivand, Roger S.}, - date = {2020-10}, - journaltitle = {Journal of Geographical Systems}, - publisher = {{Springer Science and Business Media LLC}}, - doi = {10/ghnwg3}, - url = {https://doi.org/10.1007/s10109-020-00336-0} -} - -@article{bivand_progress_2021, - title = {Progress in the {{R}} Ecosystem for Representing and Handling Spatial Data}, - author = {Bivand, Roger S.}, - date = {2021-10-01}, + date = {2020-10-16}, journaltitle = {Journal of Geographical Systems}, shortjournal = {J Geogr Syst}, - volume = {23}, - number = {4}, - pages = {515--546}, + publisher = {{Springer Science and Business Media LLC}}, issn = {1435-5949}, doi = {10/ghnwg3}, url = {https://doi.org/10.1007/s10109-020-00336-0}, - urldate = {2021-11-04}, + urldate = {2020-11-08}, langid = {english} } @@ -429,7 +418,7 @@ @incollection{cheshire_spatial_2015 keywords = {\#nosource} } -@article{clementini_comparison_1995a, +@article{clementini_comparison_1995, title = {A Comparison of Methods for Representing Topological Relationships}, author = {Clementini, Eliseo and Di Felice, Paolino}, date = {1995-05-01}, @@ -1049,6 +1038,7 @@ @article{loecher_rgooglemaps_2015 } @article{loidl_spatial_2016, + ids = {loidl_spatial_2016a}, title = {Spatial Patterns and Temporal Dynamics of Urban Bicycle Crashes—{{A}} Case Study from {{Salzburg}} ({{Austria}})}, author = {Loidl, Martin and Traun, Christoph and Wallentin, Gudrun}, date = {2016-04-01}, @@ -1106,16 +1096,19 @@ @book{lovelace_geocomputation_2019 isbn = {1-138-30451-4} } -@article{lovelace_open_2021a, +@article{lovelace_open_2021, + ids = {lovelace_open_2021a}, title = {Open Source Tools for Geographic Analysis in Transport Planning}, author = {Lovelace, Robin}, - date = {2021-01}, - volume = {23}, - number = {4}, - pages = {547--578}, + date = {2021-01-16}, + journaltitle = {Journal of Geographical Systems}, + shortjournal = {J Geogr Syst}, publisher = {{Springer Science and Business Media LLC}}, + issn = {1435-5949}, doi = {10/ghtnrp}, - url = {https://doi.org/10.1007/s10109-020-00342-2} + url = {https://doi.org/10.1007/s10109-020-00342-2}, + urldate = {2021-01-17}, + langid = {english} } @article{lovelace_propensity_2017, From d305f94bf9469b178549c4fda883a597cd52c2c1 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 00:51:57 +0000 Subject: [PATCH 076/104] Comment out experimental new de-9im section for now... --- 04-spatial-operations.Rmd | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 3aeab48fe..0d11bcc03 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -318,14 +318,15 @@ plot(l, add = TRUE) plot(p, add = TRUE) ``` -### DE-9IM strings + + -Underlying the binary predicates demonstrated in the previous section is the Dimensionally Extended 9-Intersection Model (DE-9IM). -The model was originally labelled "DE + 9IM" by its inventors, referring to the "dimension of the intersections of' boundaries, interiors, and exteriors of two features" [@clementini_comparison_1995] + + -The way the model works can be demonstrated with reference to two intersecting polygons, as illustrated in Figure \@ref(fig:de-9im). + -```{r} +```{r, echo=FALSE, eval=FALSE} b = st_sfc(st_point(c(0, 1)), st_point(c(1, 1))) # create 2 points b = st_buffer(b, dist = 1) # convert points to circles bsf = sf::st_sf(data.frame(Object = c("a", "b")), geometry = b) @@ -357,7 +358,6 @@ plot(b) text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y")) # add text ``` - ### Spatial joining Joining two non-spatial datasets relies on a shared 'key' variable, as described in Section \@ref(vector-attribute-joining). From 89503062471bc0a656009189f2e1baa015e9d619 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 00:55:03 +0000 Subject: [PATCH 077/104] Split paragraph in s. 4.2.2 in 2 --- 04-spatial-operations.Rmd | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 0d11bcc03..031278611 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -143,7 +143,8 @@ Topological relations describe the spatial relationships between objects. That may sound rather abstract and, indeed, the definition and classification of topological relations is based on earlier mathematical foundations first published as a book in 1966 [@spanier_algebraic_1995].^[ See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. ] -However, topological relations can be clearly understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships defined in English (and other human languages) as illustrated in Figure \@ref(fig:relations), which shows a variety of geometry pairs and their associated relations. + +Despite their mathematical origins, topological relations can be understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships defined in English (and other human languages) as illustrated in Figure \@ref(fig:relations), which shows a variety of geometry pairs and their associated relations. The third and fourth pairs in Figure \@ref(fig:relations) (from left to right and then down) demonstrate that, for some relations, order is important: while the relations *equals*, *intersects*, *crosses*, *touches* and *overlaps* are symmetrical, meaning that if `function(x, y)` is true, `function(y, x)` will also by true, relations in which the order of the geomtries are important such as *contains* and *within* are not. \index{topological relations} From 6c58bd9a02ef58067e0c99341db5061bf9e75f92 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 00:59:57 +0000 Subject: [PATCH 078/104] Update text in 4.2.2, improve figure 4.2 caption --- 04-spatial-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 031278611..d949c5523 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -148,7 +148,7 @@ Despite their mathematical origins, topological relations can be understood intu The third and fourth pairs in Figure \@ref(fig:relations) (from left to right and then down) demonstrate that, for some relations, order is important: while the relations *equals*, *intersects*, *crosses*, *touches* and *overlaps* are symmetrical, meaning that if `function(x, y)` is true, `function(y, x)` will also by true, relations in which the order of the geomtries are important such as *contains* and *within* are not. \index{topological relations} -```{r relations, echo=FALSE, fig.cap="Topological relations between vector geometries", fig.show='hold', out.width="33%", message=FALSE, fig.width=3, fig.height=3} +```{r relations, echo=FALSE, fig.cap="Topological relations between vector geometries, inspired by Figures 1 and 2 in Egenhofer and Herring (1990). The relations for which the function(x, y) is true are printed for each geometry pair, with x represented in pink and y represented in blue.", fig.show='hold', out.width="33%", message=FALSE, fig.width=3, fig.height=3} # source("https://github.com/Robinlovelace/geocompr/raw/c4-v2-updates-rl/code/de_9im.R") source("code/de_9im.R") library(sf) From 262c1589ccb75eeb082cb7581af5be4f8ce6f939 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 01:08:53 +0000 Subject: [PATCH 079/104] Break-up geometry creation code chunk --- 04-spatial-operations.Rmd | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index d949c5523..a4dfce107 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -170,7 +170,8 @@ de_9im(p1, p6) ``` In `sf`, functions testing for different types of topological relations are called 'binary predicates', as described in the vignette *Manipulating Simple Feature Geometries*, which can be viewed with the command [`vignette("sf3")`](https://r-spatial.github.io/sf/articles/sf3.html), and in the help page [`?geos_binary_pred`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) [@pebesma_simple_2018]. -To see how topological relations work in practice, let's create a simple reproducible example, building on the relations illustrated in Figure \@ref(fig:relations). +To see how topological relations work in practice, let's create a simple reproducible example, building on the relations illustrated in Figure \@ref(fig:relations) and consolidating knowledge of how vector geometries are represented from a previous chapter (Section \@ref(geometry)). +Note that to create tabular data representing coordinates (x and y) of the polygon vertices, we use the base R function `read.csv()` and then convert the result into a matrix, a `POLYGON` and finally an `sfc` object: ```{r, echo=FALSE} # alternative way of getting reader to data, less clear... @@ -197,6 +198,12 @@ polygon_df = read.csv(text = "x, y polygon_matrix = as.matrix(polygon_df) polygon = st_polygon(list(polygon_matrix)) polygon_sfc = st_sfc(polygon) +``` + +We will create additional geometries to demonstrate spatial relations with the following commands. +Note the use of the function `st_as_sf()` and the argument `coords` to efficiently convert from a data frame containing columns representing coordinates to an `sf` object containing points: + +```{r} line_df = read.csv(text = "x,y 0.1, 0 1, 0.1") From 4fe30d5762ba3da542764cd6f74c0697f99b54e3 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Tue, 16 Nov 2021 01:20:44 +0000 Subject: [PATCH 080/104] Finish review of Section 4.2.2 --- 04-spatial-operations.Rmd | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index a4dfce107..b1f16d1c6 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -200,7 +200,7 @@ polygon = st_polygon(list(polygon_matrix)) polygon_sfc = st_sfc(polygon) ``` -We will create additional geometries to demonstrate spatial relations with the following commands. +We will create additional geometries to demonstrate spatial relations with the following commands which, when plotted on top of the polygon created above, relate in space to one another, as shown in Figure \@ref(fig:relation-objects). Note the use of the function `st_as_sf()` and the argument `coords` to efficiently convert from a data frame containing columns representing coordinates to an `sf` object containing points: ```{r} @@ -217,7 +217,7 @@ point_df = read.csv(text = "x,y point_sf = st_as_sf(point_df, coords = c("x", "y")) ``` -```{r relation-objects, echo=FALSE, fig.cap="Points (p 1 to 4), line and polygon objects arranged to illustrate topological relations.", fig.asp=1, out.width="50%", fig.scap="Demonstration of topological relations."} +```{r relation-objects, echo=FALSE, fig.cap="Points (`point_df` 1 to 3), line and polygon objects arranged to illustrate topological relations.", fig.asp=1, out.width="50%", fig.scap="Demonstration of topological relations."} par(pty = "s") plot(polygon_sfc, border = "red", col = "gray", axes = TRUE) plot(line_sfc, add = TRUE) @@ -249,11 +249,11 @@ st_intersects(point_sf, polygon_sfc, sparse = FALSE) ``` The output is a matrix in which each row represents a feature in the target object and each column represents a feature in the selecting object. -In this case, only the first two features in `p` intersect with `a` and there is only one feature in `a` so the result has only one column. +In this case, only the third feature in `point_sf` intersects with `polygon_sfc`. +There is only one feature in `polygon_sfc` so the result has only one column. The result can be used for subsetting as we saw in Section \@ref(spatial-subsetting). -Note that `st_intersects()` returns `TRUE` for the second feature in the object `p` even though it just touches the polygon `a`: *intersects* is a 'catch-all' topological operation which identifies many types of spatial relation. - +Note: `st_intersects()` returns `TRUE` even in cases where the features just touch: *intersects* is a 'catch-all' topological operation which identifies many types of spatial relation, as illustrated in Figure \@ref(fig:relate). The opposite of `st_intersects()` is `st_disjoint()`, which returns only objects that do not spatially relate in any way to the selecting object (note `[, 1]` converts the result into a vector): ```{r 04-spatial-operations-11} @@ -261,14 +261,14 @@ st_disjoint(point_sf, polygon_sfc, sparse = FALSE)[, 1] ``` `st_within()` returns `TRUE` only for objects that are completely within the selecting object. -This applies only to the first object, which is inside the triangular polygon, as illustrated below: +This also only applies to the third point, which is inside the polygon, as illustrated below: ```{r 04-spatial-operations-12} st_within(point_sf, polygon_sfc, sparse = FALSE)[, 1] ``` Note that although the first point is *within* the triangle, it does not *touch* any part of its border. -For this reason `st_touches()` only returns `TRUE` for the second point: +For this reason `st_touches()` only returns `FALSE` for all points: ```{r 04-spatial-operations-13} st_touches(point_sf, polygon_sfc, sparse = FALSE)[, 1] @@ -277,11 +277,12 @@ st_touches(point_sf, polygon_sfc, sparse = FALSE)[, 1] What about features that do not touch, but *almost touch* the selection object? These can be selected using `st_is_within_distance()`, which has an additional `dist` argument. It can be used to set how close target objects need to be before they are selected. -Note that although point 4 is one unit of distance from the nearest node of `a` (at point 2 in Figure \@ref(fig:relation-objects)), it is still selected when the distance is set to 0.9. -This is illustrated in the code chunk below, the second line of which converts the lengthy list output into a `logical` object: +Note that although point 4 is one unit of distance from the nearest node of `polygon_sfc` (at point 2 in Figure \@ref(fig:relation-objects)), it is still selected when the distance is set to 0.9. +This is illustrated in the code chunk below, the second line of which uses the function `lengths()` to convert the lengthy list output into a `logical` object: ```{r 04-spatial-operations-14} -sel = st_is_within_distance(point_sf, polygon_sfc, dist = 0.9) # can only return a sparse matrix +# can only return a sparse matrix +sel = st_is_within_distance(point_sf, polygon_sfc, dist = 0.1) lengths(sel) > 0 ``` From a93c16babefe1b294ad0de2b77893b3b355a80f0 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 3 Dec 2021 04:49:10 +0000 Subject: [PATCH 081/104] Update bib --- geocompr.bib | 68 +++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 65 insertions(+), 3 deletions(-) diff --git a/geocompr.bib b/geocompr.bib index be2214d87..c333f14fd 100644 --- a/geocompr.bib +++ b/geocompr.bib @@ -428,7 +428,7 @@ @article{clementini_comparison_1995 number = {3}, pages = {149--178}, issn = {1069-0115}, - doi = {10.1016/1069-0115(94)00033-X}, + doi = {10/ddtnhx}, url = {https://www.sciencedirect.com/science/article/pii/106901159400033X}, urldate = {2021-11-13}, langid = {english} @@ -531,6 +531,14 @@ @article{douglas_algorithms_1973 keywords = {\#nosource} } +@manual{dunnington_ggspatial_2021, + type = {manual}, + title = {Ggspatial: Spatial Data Framework for Ggplot2}, + author = {Dunnington, Dewey}, + date = {2021}, + url = {https://CRAN.R-project.org/package=ggspatial} +} + @article{eddelbuettel_extending_2018, title = {Extending {{R}} with {{C}}++: A {{Brief Introduction}} to {{Rcpp}}}, shorttitle = {Extending {{R}} with {{C}}++}, @@ -553,7 +561,8 @@ @inproceedings{egenhofer_mathematical_1990 booktitle = {Proc. the Fourth International Symposium on Spatial Data Handing}, author = {Egenhofer, Max and Herring, John}, date = {1990}, - pages = {803--813} + pages = {803--813}, + keywords = {⛔ No DOI found} } @article{essd-11-647-2019, @@ -615,6 +624,14 @@ @book{gillespie_efficient_2016 keywords = {\#nosource} } +@manual{giraud_mapsf_2021, + type = {manual}, + title = {Mapsf: Thematic Cartography}, + author = {Giraud, Timothée}, + date = {2021}, + url = {https://CRAN.R-project.org/package=mapsf} +} + @article{goetz_evaluating_2015, title = {Evaluating Machine Learning and Statistical Prediction Techniques for Landslide Susceptibility Modeling}, author = {Goetz, J.N. and Brenning, A. and Petschko, H. and Leopold, P.}, @@ -740,6 +757,14 @@ @book{hijmans_geosphere_2016 keywords = {\#nosource} } +@manual{hijmans_terra_2021, + type = {manual}, + title = {Terra: Spatial Data Analysis}, + author = {Hijmans, Robert J.}, + date = {2021}, + url = {https://rspatial.org/terra/} +} + @book{hollander_transport_2016, title = {Transport {{Modelling}} for a {{Complete Beginner}}}, author = {Hollander, Yaron}, @@ -1226,6 +1251,14 @@ @article{moreno-monroy_public_2017 keywords = {\#nosource,Accessibility,Inequality,Latin America,Public transport,Schools} } +@manual{morganwall_rayshader_2021, + type = {manual}, + title = {Rayshader: Create Maps and Visualize Data in {{2D}} and {{3D}}}, + author = {Morgan-Wall, Tyler}, + date = {2021}, + url = {https://CRAN.R-project.org/package=rayshader} +} + @article{muenchow_geomorphic_2012, title = {Geomorphic Process Rates of Landslides along a Humidity Gradient in the Tropical {{Andes}}}, author = {Muenchow, Jannes and Brenning, Alexander and Richter, Michael}, @@ -1469,6 +1502,22 @@ @book{pebesma_spatial_2022 url = {https://r-spatial.org/book} } +@manual{pebesma_stars_2021, + type = {manual}, + title = {Stars: Spatiotemporal Arrays, Raster and Vector Data Cubes}, + author = {Pebesma, Edzer}, + date = {2021}, + url = {https://CRAN.R-project.org/package=stars} +} + +@manual{pedersen_gganimate_2020, + type = {manual}, + title = {Gganimate: A Grammar of Animated Graphics}, + author = {Pedersen, Thomas Lin and Robinson, David}, + date = {2020}, + url = {https://CRAN.R-project.org/package=gganimate} +} + @book{perpinan_rastervis_2016, title = {{{rasterVis}}}, author = {Perpiñán, Oscar and Hijmans, Robert}, @@ -1553,6 +1602,19 @@ @book{rodrigue_geography_2013 pagetotal = {432} } +@article{Roussel2020, + title = {{{lidR}}: An r Package for Analysis of Airborne Laser Scanning ({{ALS}}) Data}, + author = {Roussel, Jean-Romain and Auty, David and Coops, Nicholas C. and Tompalski, Piotr and Goodbody, Tristan R.H. and Meador, Andrew Sánchez and Bourdon, Jean-François and de Boissieu, Florian and Achim, Alexis}, + options = {useprefix=true}, + date = {2020-12}, + journaltitle = {Remote Sensing of Environment}, + volume = {251}, + pages = {112061}, + publisher = {{Elsevier BV}}, + doi = {10/ghddxb}, + url = {https://doi.org/10.1016/j.rse.2020.112061} +} + @inproceedings{rowlingson_rasp:_2003, title = {Rasp: A {{Package}} for {{Spatial Statistics}}}, booktitle = {Proceedings of the 3rd {{International Workshop}} on {{Distributed Statistical Computing}}}, @@ -1618,7 +1680,7 @@ @article{shen_classification_2018 number = {2}, pages = {514--541}, issn = {1467-9671}, - doi = {10.1111/tgis.12328}, + doi = {10/gnhcx9}, url = {https://onlinelibrary.wiley.com/doi/abs/10.1111/tgis.12328}, urldate = {2021-11-13}, langid = {english}, From 47601955d6aeb69bccea3691a45753ae7e70396c Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Fri, 3 Dec 2021 04:57:52 +0000 Subject: [PATCH 082/104] Tweak prose, see https://github.com/Robinlovelace/geocompr/pull/686/commits/df336172a30245fff8c03aa83cc96eaca79feaa2 --- 07-read-write-plot.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/07-read-write-plot.Rmd b/07-read-write-plot.Rmd index e06370b2b..f7ac63966 100644 --- a/07-read-write-plot.Rmd +++ b/07-read-write-plot.Rmd @@ -47,7 +47,7 @@ Various 'geoportals' (web services providing geospatial datasets such as [Data.g \index{geoportals} Some global geoportals overcome this issue. The [GEOSS portal](http://www.geoportal.org/) and the [Copernicus Open Access Hub](https://scihub.copernicus.eu/), for example, contain many raster datasets with global coverage. -A wealth of vector datasets can be accessed from the National Aeronautics and Space Administration agency's (NASA) [SEDAC](http://sedac.ciesin.columbia.edu/) portal and the European Union's [INSPIRE geoportal](http://inspire-geoportal.ec.europa.eu/), with global and regional coverage. +A wealth of vector datasets can be accessed from the [SEDAC](http://sedac.ciesin.columbia.edu/) portal run by the National Aeronautics and Space Administration (NASA) and the European Union's [INSPIRE geoportal](http://inspire-geoportal.ec.europa.eu/), with global and regional coverage. Most geoportals provide a graphical interface allowing datasets to be queried based on characteristics such as spatial and temporal extent, the United States Geological Services' [EarthExplorer](https://earthexplorer.usgs.gov/) being a prime example. *Exploring* datasets interactively on a browser is an effective way of understanding available layers. From db7daee1974756d73c08fe66dcc4de869758fe35 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 20:19:17 +0000 Subject: [PATCH 083/104] Refactor code + prose in 4.2.3 --- 04-spatial-operations.Rmd | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 50bc73816..5395acc11 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -370,31 +370,30 @@ text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y")) # add text ### Spatial joining Joining two non-spatial datasets relies on a shared 'key' variable, as described in Section \@ref(vector-attribute-joining). -Spatial data joining applies the same concept, but instead relies on shared areas of geographic space (it is also known as spatial overlay). -As with attribute data, joining adds a new column to the target object (the argument `x` in joining functions), from a source object (`y`). +Spatial data joining applies the same concept, but instead relies on spatial relations, described in the previous section. +As with attribute data, joining adds new columns to the target object (the argument `x` in joining functions), from a source object (`y`). \index{join!spatial} \index{spatial!join} -The process can be illustrated by an example. -Imagine you have ten points randomly distributed across the Earth's surface. -Of the points that are on land, which countries are they in? -Random points to demonstrate spatial joining are created as follows: +The process is illustrated by the following example: imagine you have ten points randomly distributed across the Earth's surface and you ask, for the points that are on land, which countries are they in? +Implementing this idea in a [reproducible example](https://github.com/Robinlovelace/geocompr/blob/main/code/04-spatial-join.R) will build your geographic data handling skills and show how spatial joins work. +The starting point is to create points that are randomly scattered over the Earth's surface: ```{r 04-spatial-operations-19} set.seed(2018) # set seed for reproducibility -(bb_world = st_bbox(world)) # the world's bounds +(bb = st_bbox(world)) # the world's bounds random_df = tibble( - x = runif(n = 10, min = bb_world[1], max = bb_world[3]), - y = runif(n = 10, min = bb_world[2], max = bb_world[4]) + x = runif(n = 10, min = bb[1], max = bb[3]), + y = runif(n = 10, min = bb[2], max = bb[4]) ) random_points = random_df %>% st_as_sf(coords = c("x", "y")) %>% # set coordinates st_set_crs(4326) # set geographic CRS ``` -The scenario is illustrated in Figure \@ref(fig:spatial-join). -The `random_points` object (top left) has no attribute data, while the `world` (top right) does. -The spatial join operation is done by `st_join()`, which adds the `name_long` variable to the points, resulting in `random_joined` which is illustrated in Figure \@ref(fig:spatial-join) (bottom left --- see [`04-spatial-join.R`](https://github.com/Robinlovelace/geocompr/blob/main/code/04-spatial-join.R)). +The scenario illustrated in Figure \@ref(fig:spatial-join) shows that the `random_points` object (top left) lacks attribute data, while the `world` (top right) has attributes, including country names shown for a sample of countries in the legend. +Spatial joins are implemented with `st_join()`, as illustrated in the code chunk below. +The output is the `random_joined` object which is illustrated in Figure \@ref(fig:spatial-join) (bottom left). Before creating the joined dataset, we use spatial subsetting to create `world_random`, which contains only countries that contain random points, to verify the number of country names returned in the joined dataset should be four (see the top right panel of Figure \@ref(fig:spatial-join)). ```{r 04-spatial-operations-20, message=FALSE} From d8aad30e9ba48dcb1866824008f7b49bb1007d8f Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 20:27:30 +0000 Subject: [PATCH 084/104] Complete review of 4.2.3 --- 04-spatial-operations.Rmd | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 5395acc11..5d2be4780 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -407,13 +407,10 @@ source("https://github.com/Robinlovelace/geocompr/raw/main/code/04-spatial-join. tmap_arrange(jm1, jm2, jm3, jm4, nrow = 2, ncol = 2) ``` -By default, `st_join()` performs a left join (see Section \@ref(vector-attribute-joining)), but it can also do inner joins by setting the argument `left = FALSE`. -Like spatial subsetting, the default topological operator used by `st_join()` is `st_intersects()`. -This can be changed with the `join` argument (see `?st_join` for details). -In the example above, we have added features of a polygon layer to a point layer. -In other cases, we might want to join point attributes to a polygon layer. -There might be occasions where more than one point falls inside one polygon. -In such a case `st_join()` duplicates the polygon feature: it creates a new row for each match. +By default, `st_join()` performs a left join, meaning that the result is an object containing all rows from `x` including rows with no match in `y` (see Section \@ref(vector-attribute-joining)), but it can also do inner joins by setting the argument `left = FALSE`. +Like spatial subsetting, the default topological operator used by `st_join()` is `st_intersects()`, which can be changed by setting the `join` argument (see `?st_join` for details). +The example above demonstrates the addition of a column from a polygon layer to a point layer, but same approach works regardless of geometry types. +In such cases, for example when `x` contains polygons, each of which match multiple objects in `y`, spatial joins will result in duplicate features, creates a new row for each match in `y`. ### Non-overlapping joins From 471defa97bb91e09f74d81431434fb60c41d80e1 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 21:05:42 +0000 Subject: [PATCH 085/104] Add idea --- 04-spatial-operations.Rmd | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 5d2be4780..4763339a1 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -412,9 +412,11 @@ Like spatial subsetting, the default topological operator used by `st_join()` is The example above demonstrates the addition of a column from a polygon layer to a point layer, but same approach works regardless of geometry types. In such cases, for example when `x` contains polygons, each of which match multiple objects in `y`, spatial joins will result in duplicate features, creates a new row for each match in `y`. + + ### Non-overlapping joins -Sometimes two geographic datasets do not touch but still have a strong geographic relationship enabling joins. +Sometimes two geographic datasets do not touch but still have a strong geographic relationship. The datasets `cycle_hire` and `cycle_hire_osm`, already attached in the **spData** package, provide a good example. Plotting them shows that they are often closely related but they do not touch, as shown in Figure \@ref(fig:cycle-hire), a base version of which is created with the following code below: \index{join!non-overlapping} From 7ecb0de8650005c7d95d662c274e86fb19e66001 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 21:36:10 +0000 Subject: [PATCH 086/104] Complete review of Section 4.2.4 --- 04-spatial-operations.Rmd | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 4763339a1..a1f229c40 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -460,30 +460,43 @@ if (knitr::is_latex_output()){ Imagine that we need to join the `capacity` variable in `cycle_hire_osm` onto the official 'target' data contained in `cycle_hire`. This is when a non-overlapping join is needed. -The simplest method is to use the topological operator `st_is_within_distance()` shown in Section \@ref(topological-relations), using a threshold distance of 20 m. -Note that, before performing the relation, both objects are transformed into a projected CRS. -These projected objects are created below (note the affix `_P`, short for projected): +The simplest method is to use the topological operator `st_is_within_distance()`, as demonstrated below using a threshold distance of 20 m (note that this works with projected and unprojected data). ```{r 04-spatial-operations-24} +sel = st_is_within_distance(cycle_hire, cycle_hire_osm, dist = 20) +summary(lengths(sel) > 0) +``` + +```{r 04-spatial-operations-24-without-s2-test, eval=FALSE, echo=FALSE} +sf::sf_use_s2(FALSE) +sel = st_is_within_distance(cycle_hire, cycle_hire_osm, dist = 20) +summary(lengths(sel) > 0) +# still works: must be lwgeom or some other magic! +``` + + +```{r 04-spatial-operations-24-projected, eval=FALSE, echo=FALSE} +# This chunk contains the non-overlapping join on projected data, a step that is no longer needed: +# Note that, before performing the relation, both objects are transformed into a projected CRS. +# These projected objects are created below (note the affix `_P`, short for projected): cycle_hire_P = st_transform(cycle_hire, 27700) cycle_hire_osm_P = st_transform(cycle_hire_osm, 27700) sel = st_is_within_distance(cycle_hire_P, cycle_hire_osm_P, dist = 20) summary(lengths(sel) > 0) ``` -This shows that there are `r sum(lengths(sel) > 0)` points in the target object `cycle_hire_P` within the threshold distance of `cycle_hire_osm_P`. -How to retrieve the *values* associated with the respective `cycle_hire_osm_P` points? +This shows that there are `r sum(lengths(sel) > 0)` points in the target object `cycle_hire` within the threshold distance of `cycle_hire_osm`. +How to retrieve the *values* associated with the respective `cycle_hire_osm` points? The solution is again with `st_join()`, but with an addition `dist` argument (set to 20 m below): ```{r 04-spatial-operations-25} -z = st_join(cycle_hire_P, cycle_hire_osm_P, - join = st_is_within_distance, dist = 20) +z = st_join(cycle_hire, cycle_hire_osm, st_is_within_distance, dist = 20) nrow(cycle_hire) nrow(z) ``` Note that the number of rows in the joined result is greater than the target. -This is because some cycle hire stations in `cycle_hire_P` have multiple matches in `cycle_hire_osm_P`. +This is because some cycle hire stations in `cycle_hire` have multiple matches in `cycle_hire_osm`. To aggregate the values for the overlapping points and return the mean, we can use the aggregation methods learned in Chapter \@ref(attr), resulting in an object with the same number of rows as the target: ```{r 04-spatial-operations-26} @@ -500,7 +513,7 @@ plot(cycle_hire_osm["capacity"]) plot(z["capacity"]) ``` -The result of this join has used a spatial operation to change the attribute data associated with simple features; the geometry associated with each feature has remained unchanged. +The result of this join has used a spatial operation to change the attribute data associated with simple features; the geometry associated with each feature has remained unchanged. ### Spatial data aggregation {#spatial-aggr} From bd857de83c8ab5f4aae096b01053f706a767179e Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 22:20:34 +0000 Subject: [PATCH 087/104] Review section on spatial aggregation --- 04-spatial-operations.Rmd | 37 +++++++++++++++++++------------------ 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index a1f229c40..5bacf1e42 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -517,42 +517,43 @@ The result of this join has used a spatial operation to change the attribute dat ### Spatial data aggregation {#spatial-aggr} -Like attribute data aggregation, covered in Section \@ref(vector-attribute-aggregation), spatial data aggregation can be a way of *condensing* data. -Aggregated data show some statistics\index{statistics} about a variable (typically average or total) in relation to some kind of *grouping variable*. -Section \@ref(vector-attribute-aggregation) demonstrated how `aggregate()` and `group_by() %>% summarize()` condense data based on attribute variables. -This section demonstrates how the same functions work using spatial grouping variables. +As with attribute data aggregation, spatial data aggregation *condenses* data: aggregated outputs have fewer rows than non-aggregated inputs. +Statistical *aggregating functions*, such as mean average or sum, summarise multiple values \index{statistics} of a variable, and return a single value per *grouping variable*. +Section \@ref(vector-attribute-aggregation) demonstrated how `aggregate()` and `group_by() %>% summarize()` condense data based on attribute variables, this section shows how the same functions work with spatial objects. \index{aggregation!spatial} -Returning to the example of New Zealand, imagine you want to find out the average height of high points in each region. -This is a good example of spatial aggregation: it is the geometry of the source (`y` or `nz` in this case) that defines how values in the target object (`x` or `nz_height`) are grouped. -This is illustrated using the base `aggregate()` function below: +Returning to the example of New Zealand, imagine you want to find out the average height of high points in each region: it is the geometry of the source (`y` or `nz` in this case) that defines how values in the target object (`x` or `nz_height`) are grouped. +This can be done in a single line of code with base R's `aggregate()` method: ```{r 04-spatial-operations-28} -nz_avheight = aggregate(x = nz_height, by = nz, FUN = mean) +nz_agg = aggregate(x = nz_height, by = nz, FUN = mean) ``` -The result of the previous command is an `sf` object with the same geometry as the (spatial) aggregating object (`nz`).^[ -This can be verified with `identical(st_geometry(nz), st_geometry(nz_avheight))`. -] -The result of the previous operation is illustrated in Figure \@ref(fig:spatial-aggregation). -The same result can also be generated using the 'tidy' functions `group_by()` and `summarize()` (used in combination with `st_join()`): +The result of the previous command is an `sf` object with the same geometry as the (spatial) aggregating object (`nz`), which you can verify with the command `identical(st_geometry(nz), st_geometry(nz_agg))`. +The result of the previous operation is illustrated in Figure \@ref(fig:spatial-aggregation), which shows the average value of features in `nz_height` within each of New Zealand's 16 regions. +The same result can also be generated by piping the output from `st_join()` into the 'tidy' functions `group_by()` and `summarize()` as follows: ```{r spatial-aggregation, echo=FALSE, fig.cap="Average height of the top 101 high points across the regions of New Zealand.", fig.asp=1, message=FALSE, out.width="50%"} library(tmap) -tm_shape(nz_avheight) + +tm_shape(nz_agg) + tm_fill("elevation", breaks = seq(27, 30, by = 0.5) * 1e2) + tm_borders() ``` ```{r 04-spatial-operations-29} -nz_avheight2 = nz %>% - st_join(nz_height) %>% +nz_agg2 = st_join(x = nz, y = nz_height) %>% group_by(Name) %>% summarize(elevation = mean(elevation, na.rm = TRUE)) ``` -The resulting `nz_avheight` objects have the same geometry as the aggregating object `nz` but with a new column representing the mean average height of points within each region of New Zealand (other summary functions such as `median()` and `sd()` can be used in place of `mean()`). -Note that regions containing no points have an associated `elevation` value of `NA`. +```{r test-tidy-spatial-join, eval=FALSE, echo=FALSE} +plot(nz_agg) +plot(nz_agg2) +# aggregate looses the name of aggregating objects +``` + +The resulting `nz_agg` objects have the same geometry as the aggregating object `nz` but with a new column summarising the values of `x` in each region using the function `mean()` (which cold, of course, be replaced by `median()`, `sd()` and other functions that return a single value). +Note: one difference between the `aggregate()` and `group_by() %>% summarize()` approaches is that the former results in `NA` values for unmatching region names, while the latter preserves region names and is more flexible in terms of aggregating functions and the column names of the results. For aggregating operations which also create new geometries, see Section \@ref(geometry-unions). Spatial congruence\index{spatial congruence} is an important concept related to spatial aggregation. From 59639bc66b10e30d8f06d41d36ae9614a12ed45e Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 22:21:05 +0000 Subject: [PATCH 088/104] Rename spatia-aggr section --- 04-spatial-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 5bacf1e42..c0f590467 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -515,7 +515,7 @@ plot(z["capacity"]) The result of this join has used a spatial operation to change the attribute data associated with simple features; the geometry associated with each feature has remained unchanged. -### Spatial data aggregation {#spatial-aggr} +### Spatial aggregation {#spatial-aggr} As with attribute data aggregation, spatial data aggregation *condenses* data: aggregated outputs have fewer rows than non-aggregated inputs. Statistical *aggregating functions*, such as mean average or sum, summarise multiple values \index{statistics} of a variable, and return a single value per *grouping variable*. From 1d28020d8c19fec6d5c1d414c48f6ca95f623ecc Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 22:25:20 +0000 Subject: [PATCH 089/104] Create new section on incongruent zones, they deserve one! --- 04-spatial-operations.Rmd | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index c0f590467..191e8ab89 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -556,14 +556,15 @@ The resulting `nz_agg` objects have the same geometry as the aggregating object Note: one difference between the `aggregate()` and `group_by() %>% summarize()` approaches is that the former results in `NA` values for unmatching region names, while the latter preserves region names and is more flexible in terms of aggregating functions and the column names of the results. For aggregating operations which also create new geometries, see Section \@ref(geometry-unions). +### Aggregating incongruent zones + Spatial congruence\index{spatial congruence} is an important concept related to spatial aggregation. An *aggregating object* (which we will refer to as `y`) is *congruent* with the target object (`x`) if the two objects have shared borders. Often this is the case for administrative boundary data, whereby larger units --- such as Middle Layer Super Output Areas ([MSOAs](https://www.ons.gov.uk/methodology/geography/ukgeographies/censusgeography)) in the UK or districts in many other European countries --- are composed of many smaller units. *Incongruent* aggregating objects, by contrast, do not share common borders with the target [@qiu_development_2012]. -This is problematic for spatial aggregation (and other spatial operations) illustrated in Figure \@ref(fig:areal-example). -Areal interpolation overcomes this issue by transferring values from one set of areal units to another. -Algorithms developed for this task include area weighted and 'pycnophylactic' areal interpolation methods [@tobler_smooth_1979]. +This is problematic for spatial aggregation (and other spatial operations) illustrated in Figure \@ref(fig:areal-example): aggregating the centroid of each sub-zone will not return accurate results. +Areal interpolation overcomes this issue by transferring values from one set of areal units to another, using a range of algorithms including simple area weighted approaches and more sophisticated approaches such as 'pycnophylactic' methods [@tobler_smooth_1979]. ```{r areal-example, echo=FALSE, fig.cap="Illustration of congruent (left) and incongruent (right) areal units with respect to larger aggregating zones (translucent blue borders).", fig.asp=0.2, fig.scap="Illustration of congruent and incongruent areal units."} source("https://github.com/Robinlovelace/geocompr/raw/main/code/04-areal-example.R", print.eval = TRUE) From 6e7f825d4c2fd239a87be3d1e6881eeb880dc7e1 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 22:28:27 +0000 Subject: [PATCH 090/104] Remove unclear explanation of st_interpolate_aw() --- 04-spatial-operations.Rmd | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 191e8ab89..272616f38 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -574,9 +574,7 @@ The **spData** package contains a dataset named `incongruent` (colored polygons Let us assume that the `value` column of `incongruent` refers to the total regional income in million Euros. How can we transfer the values of the underlying nine spatial polygons into the two polygons of `aggregating_zones`? -The simplest useful method for this is *area weighted* spatial interpolation. -In this case values from the `incongruent` object are allocated to the `aggregating_zones` in proportion to area; the larger the spatial intersection between input and output features, the larger the corresponding value. -For instance, if one intersection of `incongruent` and `aggregating_zones` is 1.5 km^2^ but the whole incongruent polygon in question has 2 km^2^ and a total income of 4 million Euros, then the target aggregating zone will obtain three quarters of the income, in this case 3 million Euros. +The simplest useful method for this is *area weighted* spatial interpolation, which transfers values from the `incongruent` object to a new column in `aggregating_zones` in proportion with the area of overlap: the larger the spatial intersection between input and output features, the larger the corresponding value. This is implemented in `st_interpolate_aw()`, as demonstrated in the code chunk below. ```{r 04-spatial-operations-30} From 1232ac2370f571711eff295aeba19110ed66d7cb Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 22:38:00 +0000 Subject: [PATCH 091/104] Complete review of new section on incongruent data --- 04-spatial-operations.Rmd | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 272616f38..a62114b2e 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -578,16 +578,14 @@ The simplest useful method for this is *area weighted* spatial interpolation, wh This is implemented in `st_interpolate_aw()`, as demonstrated in the code chunk below. ```{r 04-spatial-operations-30} -agg_aw = st_interpolate_aw(incongruent[, "value"], aggregating_zones, - extensive = TRUE) -# show the aggregated result +iv = incongruent["value"] # keep only the values to be transferred +agg_aw = st_interpolate_aw(iv, aggregating_zones, ext = TRUE) agg_aw$value ``` -In our case it is meaningful to sum up the values of the intersections falling into the aggregating zones since total income is a so-called spatially extensive variable. -This would be different for spatially intensive variables, which are independent of the spatial units used, such as income per head or [percentages](http://ibis.geog.ubc.ca/courses/geob370/notes/intensive_extensive.htm). -In this case it is more meaningful to apply an average function when doing the aggregation instead of a sum function. -To do so, one would only have to set the `extensive` parameter to `FALSE`. +In our case it is meaningful to sum up the values of the intersections falling into the aggregating zones since total income is a so-called spatially extensive variable (which increases with area), assuming income is evenly distributed across the smaller zones (hence the warning message above). +This would be different for spatially [intensive](http://ibis.geog.ubc.ca/courses/geob370/notes/intensive_extensive.htm) variables such as *average* income or percentages, which do not increase as the area increases. +`st_interpolate_aw()` works equally with spatially intensive variables: set the `extensive` parameter to `FALSE` and it will use an average rather than a sum function when doing the aggregation. ### Distance relations From ca834351456a05e16f5fa6af7a12d831b96b8348 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 22:39:19 +0000 Subject: [PATCH 092/104] Rename + label incongruent section --- 04-spatial-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index a62114b2e..15673a8f0 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -556,7 +556,7 @@ The resulting `nz_agg` objects have the same geometry as the aggregating object Note: one difference between the `aggregate()` and `group_by() %>% summarize()` approaches is that the former results in `NA` values for unmatching region names, while the latter preserves region names and is more flexible in terms of aggregating functions and the column names of the results. For aggregating operations which also create new geometries, see Section \@ref(geometry-unions). -### Aggregating incongruent zones +### Joining incongruent layers {#incongruent} Spatial congruence\index{spatial congruence} is an important concept related to spatial aggregation. An *aggregating object* (which we will refer to as `y`) is *congruent* with the target object (`x`) if the two objects have shared borders. From 8618b4ca36989b2e945dc39ac8ab21fad05ea854 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 23:13:58 +0000 Subject: [PATCH 093/104] Refactor prose in raster section --- 04-spatial-operations.Rmd | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 15673a8f0..f7ff81d00 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -657,21 +657,18 @@ elev[clip] # terra::extract(elev, ext(clip)) ``` -Basically, this amounts to retrieving the values of the first raster (here: `elev`) falling within the extent of a second raster (here: `clip`). +This amounts to retrieving the values of the first raster object (in this case `elev`) that fall within the extent of a second raster (here: `clip`), as illustrated in Figure \@ref(fig:raster-subset). ```{r raster-subset, echo = FALSE, fig.cap = "Original raster (left). Raster mask (middle). Output of masking a raster (right).", fig.scap="Subsetting raster values."} knitr::include_graphics("figures/04_raster_subset.png") ``` -So far, the subsetting returned the values of specific cells, however, when doing spatial subsetting, one often also expects a spatial object as an output. -To do this, we can use again the `[` when we additionally set the `drop` parameter to `FALSE`. -Let's illustrate this by retrieving the first two cells of `elev` as an individual raster object. -As mentioned in Section \@ref(manipulating-raster-objects), the `[` operator accepts various inputs to subset rasters and returns a raster object when `drop = FALSE`. -The code chunk below subsets the `elev` raster by cell ID and row-column index with identical results: the first two cells on the top row (only the first 2 lines of the output is shown): +The example above returned the values of specific cells, but in many cases spatial outputs from subsetting operations on raster datasets are needed. +This can be done using the `[` operator, with `drop = FALSE`, as outlined in Section \@ref(manipulating-raster-objects), which also shows how raster objects can be subsetted by various objects. +This is demonstrated in the code below, which returns the first two cells of `elev` as a raster object the first two cells on the top row (only the first 2 lines of the output is shown): ```{r 04-spatial-operations-36, eval=FALSE} elev[1:2, drop = FALSE] # spatial subsetting with cell IDs -elev[1, 1:2, drop = FALSE] # spatial subsetting by row,column indices #> class : SpatRaster #> dimensions : 1, 2, 1 (nrow, ncol, nlyr) #> ... @@ -876,20 +873,22 @@ Chapter \@ref(gis) shows how to access such GIS functionality from within R. \index{map algebra!zonal operations} Just like focal operations, *zonal* operations apply an aggregation function to multiple raster cells. -However, a second raster, usually a categorical raster, defines the *zonal filters* (or 'zones') in the case of zonal operations as opposed to a predefined neighborhood window in the case of focal operations (see previous Section). -Consequently, the raster cells defining the zonal filter do not necessarily have to be neighbors. -Our grain size raster is a good example (right panel of Figure \@ref(fig:cont-raster)) because the different grain sizes are spread in an irregular fashion throughout the raster. +However, a second raster, usually with categorical values, defines the *zonal filters* (or 'zones') in the case of zonal operations, as opposed to a predefined neighborhood window in the case of focal operation presented in the previous section. +Consequently, raster cells defining the zonal filter do not necessarily have to be neighbors. +Our grain size raster is a good example, as illustrated in the right panel of Figure \@ref(fig:cont-raster)): different grain sizes are spread irregularly throughout the raster. Finally, the result of a zonal operation is a summary table grouped by zone which is why this operation is also known as *zonal statistics* in the GIS world\index{GIS}. -This is in contrast to focal operations which return a raster object (see previous Section). +This is in contrast to focal operations which return a raster object. -For example, to find the mean elevation for each grain size class (Figure \@ref(fig:cont-raster)), we use the `zonal()` function. +The following code chunk uses the `zonal()` function to calculate the mean elevation associated with each grain size class, for example. +The output is shown in Figure \@ref(fig:cont-raster)). ```{r 04-spatial-operations-43, eval=FALSE} z = zonal(elev, grain, fun = "mean") z ``` -This returns the statistics\index{statistics} for each category, here the mean altitude for each grain size class.^[It is also possible to get a raster with calculated statistics for each zone by setting the `as.raster` argument to `TRUE`.] +This returns the statistics\index{statistics} for each category, here the mean altitude for each grain size class. +Note: it is also possible to get a raster with calculated statistics for each zone by setting the `as.raster` argument to `TRUE`. ### Global operations and distances From 2f256c16d1fc284c3d078274768850807e4fdb2a Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 23:16:41 +0000 Subject: [PATCH 094/104] Evaluate zonal fun --- 04-spatial-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index f7ff81d00..49ff8631c 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -882,7 +882,7 @@ This is in contrast to focal operations which return a raster object. The following code chunk uses the `zonal()` function to calculate the mean elevation associated with each grain size class, for example. The output is shown in Figure \@ref(fig:cont-raster)). -```{r 04-spatial-operations-43, eval=FALSE} +```{r 04-spatial-operations-43} z = zonal(elev, grain, fun = "mean") z ``` From 887d73d37349d065d4b6b0744b53b008cd8e3ac5 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 23:20:04 +0000 Subject: [PATCH 095/104] Global object rename w. fd: sd 'bb_world' 'bb' 01-introduction.Rmd 02-spatial-data.Rmd 03-attribute-operations.Rmd 04-spatial-operations.Rmd 05-geometry-operations.Rmd 06-reproj.Rmd 07-read-write-plot.Rmd 08-mapping.Rmd 09-gis.Rmd 10-algorithms.Rmd 11-spatial-cv.Rmd 12-transport.Rmd 13-location.Rmd 14-eco.Rmd 15-synthesis.Rmd DESCRIPTION README.Rmd README.md _01-ex.Rmd _02-ex.Rmd _03-ex.Rmd _04-ex.Rmd _05-ex.Rmd _06-ex.Rmd _07-ex.Rmd _404.Rmd apps/CycleHireApp/app.R apps/coffeeApp/app.R code/01-cranlogs.R code/02-contpop.R code/02-datum-fig.R code/02-raster-crs.R code/02-raster-intro-plot.R code/02-raster-intro-plot2.R code/02-vector-crs.R code/02-vectorplots.R code/03-cont-raster-plot.R code/04-areal-example.R code/04-focal-example.R code/04-local-operations.R code/04-ndvi.R code/04-raster-subset.R code/04-spatial-join.R code/05-bilinear.R code/05-contour-tmap.R code/05-extend-example.R code/05-pointextr.R code/05-raster-vectorization1.R code/05-raster-vectorization2.R code/05-us-regions.R code/05-vector-rasterization1.R code/05-vector-rasterization2.R code/05-venn-clip.R code/08-break-styles.R code/08-layout1.R code/08-layout2.R code/08-layout3.R code/08-map-pkgs.R code/08-tmpal.R code/08-tmshape.R code/08-tmstyles.R code/08-urban-animation.R code/08-usboundaries.R code/09-saga-wetness.R code/09-sliver.R code/09-tsp.R code/10-centroid-alg.R code/10-centroid-setup.R code/10-hello.R code/10-polycent.R code/11-partitioning.R code/11-spatial-cv-jm.R code/11-svm.R code/12-cycleways.R code/12-desire.R code/12-transport-data-gen.R code/12-zones.R code/13-location-jm.R code/14-eco.R code/14-rf.R code/add-impact.R code/before_script.R code/chapters/01-introduction.R code/chapters/02-spatial-data.R code/chapters/03-attribute-operations.R code/chapters/04-spatial-operations.R code/chapters/05-geometry-operations.R code/chapters/06-reproj.R code/chapters/07-read-write-plot.R code/chapters/08-mapping.R code/chapters/09-gis.R code/chapters/10-algorithms.R code/chapters/11-spatial-cv.R code/chapters/12-transport.R code/chapters/13-location.R code/chapters/14-eco.R code/chapters/15-synthesis.R code/chapters/README.R code/chapters/index.R code/chapters/references.R code/de_9im.R code/extra-pkgs.R code/frontcover.R code/generate-chapter-code.R code/hex_sticker.R code/list-contributors.R code/old-to-future-remove/04_spatial_incongruence.R code/old-to-future-remove/06_raster_reprojection_tests.R code/old-to-future-remove/08-uscolonize.R code/old-to-future-remove/10-centroid.R code/old-to-future-remove/10-earthquakes.R code/old-to-future-remove/12-code-extension.R code/old-to-future-remove/12-desire-front.R code/old-to-future-remove/globe.R code/old-to-future-remove/sfr-class-diagram-gen.R code/old-to-future-remove/spData.R code/sf-classes.R code/sf-revdep.R extdata/postgis_data.Rdata geocompr.Rproj index.Rmd references.Rmd --- code/04-spatial-join.R | 14 +++++++------- code/chapters/04-spatial-operations.R | 6 +++--- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/code/04-spatial-join.R b/code/04-spatial-join.R index de3ebb080..2c3857ad5 100644 --- a/code/04-spatial-join.R +++ b/code/04-spatial-join.R @@ -10,10 +10,10 @@ library(tmap) if(!exists("random_joined")) { set.seed(2018) - bb_world = st_bbox(world) + bb = st_bbox(world) random_df = tibble::tibble( - x = runif(n = 10, min = bb_world[1], max = bb_world[3]), - y = runif(n = 10, min = bb_world[2], max = bb_world[4]) + x = runif(n = 10, min = bb[1], max = bb[3]), + y = runif(n = 10, min = bb[2], max = bb[4]) ) random_points = st_as_sf(random_df, coords = c("x", "y")) %>% st_set_crs(4326) @@ -29,26 +29,26 @@ random_joined$name_long = as.character(random_joined$name_long) jm0 = tm_shape(world) + tm_borders(lwd = 0.2) + tm_format("World") jm1 = jm0 + - tm_shape(shp = random_points, bbox = bb_world) + + tm_shape(shp = random_points, bbox = bb) + tm_symbols(col = "black", shape = 4, border.lwd = 2) + tm_layout(scale = 1, legend.bg.color = "white", legend.bg.alpha = 0.3, legend.position = c("right", "bottom")) jm2 = jm0 + - tm_shape(world_random, bbox = bb_world) + + tm_shape(world_random, bbox = bb) + tm_fill(col = "name_long", palette = "Dark2") + tm_layout(legend.show = FALSE) #tm_borders(col = "name_long", lwd = 4) + # issue with tmap: no variable border # tm_layout(scale = 1, legend.bg.color = "white", legend.bg.alpha = 0.3, legend.position = c("right", "bottom")) jm3 = jm0 + - tm_shape(shp = random_joined, bbox = bb_world) + + tm_shape(shp = random_joined, bbox = bb) + tm_symbols(col = "name_long", shape = 4, border.lwd = 2, palette = "Dark2") + tm_layout(legend.show = FALSE) # + # tm_layout(scale = 1, legend.bg.color = "white", legend.bg.alpha = 0.3, legend.position = c("right", "bottom")) jm4 = jm0 + - tm_shape(shp = random_joined, bbox = bb_world) + + tm_shape(shp = random_joined, bbox = bb) + tm_symbols(col = "name_long", shape = 4, border.lwd = 2, palette = "Dark2") + tm_layout(legend.only = TRUE) diff --git a/code/chapters/04-spatial-operations.R b/code/chapters/04-spatial-operations.R index 7b515255f..d9c2373ac 100644 --- a/code/chapters/04-spatial-operations.R +++ b/code/chapters/04-spatial-operations.R @@ -117,10 +117,10 @@ lengths(sel) > 0 ## ----04-spatial-operations-19-------------------------------------------- set.seed(2018) # set seed for reproducibility -(bb_world = st_bbox(world)) # the world's bounds +(bb = st_bbox(world)) # the world's bounds random_df = tibble( - x = runif(n = 10, min = bb_world[1], max = bb_world[3]), - y = runif(n = 10, min = bb_world[2], max = bb_world[4]) + x = runif(n = 10, min = bb[1], max = bb[3]), + y = runif(n = 10, min = bb[2], max = bb[4]) ) random_points = random_df %>% st_as_sf(coords = c("x", "y")) %>% # set coordinates From 8ca4363b2197d321e60e18352927a725458463a9 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Thu, 9 Dec 2021 23:26:30 +0000 Subject: [PATCH 096/104] Source local script for CI --- 04-spatial-operations.Rmd | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 49ff8631c..f8128db03 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -403,7 +403,8 @@ random_joined = st_join(random_points, world["name_long"]) ``` ```{r spatial-join, echo=FALSE, fig.cap="Illustration of a spatial join. A new attribute variable is added to random points (top left) from source world object (top right) resulting in the data represented in the final panel.", fig.asp=0.5, warning=FALSE, message=FALSE, out.width="100%", fig.scap="Illustration of a spatial join."} -source("https://github.com/Robinlovelace/geocompr/raw/main/code/04-spatial-join.R") +# source("https://github.com/Robinlovelace/geocompr/raw/main/code/04-spatial-join.R") +source("code/04-spatial-join.R") tmap_arrange(jm1, jm2, jm3, jm4, nrow = 2, ncol = 2) ``` From 92fc311c480ed6689852c588a513036b408a01e9 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 13 Dec 2021 15:36:17 +0000 Subject: [PATCH 097/104] Demonstrate st_filter, close #688 --- 04-spatial-operations.Rmd | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index f8128db03..48f40726c 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -49,9 +49,6 @@ Section \@ref(spatial-ras) presents spatial operations on raster datasets using ### Spatial subsetting - - - Spatial subsetting is the process of taking a spatial object and returning a new object containing only features that *relate* in space to another object. Analogous to *attribute subsetting* (covered in Section \@ref(vector-attribute-subsetting)), subsets of `sf` data frames can be created with square bracket (`[`) operator using the syntax `x[y, , op = st_intersects]`, where `x` is an `sf` object from which a subset of rows will be returned, `y` is the 'subsetting object' and `, op = st_intersects` is an optional argument that specifies the topological relation (also known as the binary predicate) used to do the subsetting. The default topological relation used when an `op` argument is not provided is `st_intersects()`: the command `x[y, ]` is identical to `x[y, , op = st_intersects]` shown above but not `x[y, , op = st_disjoint]` (the meaning of these and other topological relations is described in the next section). @@ -115,7 +112,7 @@ sel_logical = lengths(sel_sgbp) > 0 canterbury_height2 = nz_height[sel_logical, ] ``` -In the above code chunk, an object of class `sgbp` (a sparse geometry binary predicate, a list of length `x` in the spatial operation) is created and then converted into a logical vector `sel_logical` (containing only `TRUE` and `FALSE` values). +The above code chunk creates an object of class `sgbp` (a sparse geometry binary predicate, a list of length `x` in the spatial operation) and then converts it into a logical vector `sel_logical` (containing only `TRUE` and `FALSE` values, something that can also be used by **dplyr**'s filter function). \index{binary predicate|seealso {topological relations}} The function `lengths()` identifies which features in `nz_height` intersect with *any* objects in `y`. In this case 1 is the greatest possible value but for more complex operations one could use the method to subset only features that intersect with, for example, 2 or more features from the source object. @@ -125,14 +122,28 @@ Note: another way to return a logical output is by setting `sparse = FALSE` (mea Note: the solution involving `sgbp` objects is more generalisable though, as it works for many-to-many operations and has lower memory requirements. ``` -It should be noted that a logical can also be used with `filter()` as follows (`sparse = FALSE` is explained in Section \@ref(topological-relations)): +Note: since late 2019, the same result can be achieved with the function `st_filter()`: -```{r 04-spatial-operations-7b} +```{r} canterbury_height3 = nz_height %>% + st_filter(y = canterbury, .predicate = st_intersects) +``` + + +```{r 04-spatial-operations-7b-old, eval=FALSE, echo=FALSE} +# Additional tests of subsetting +canterbury_height4 = nz_height %>% filter(st_intersects(x = ., y = canterbury, sparse = FALSE)) +canterbury_height5 = nz_height %>% + filter(sel_logical) +identical(canterbury_height3, canterbury_height4) +identical(canterbury_height3, canterbury_height5) +identical(canterbury_height2, canterbury_height4) +identical(canterbury_height, canterbury_height4) +waldo::compare(canterbury_height2, canterbury_height4) ``` -At this point, there are three versions of `canterbury_height`, one created with spatial subsetting directly and the other two via intermediary selection objects. +At this point, there are three identical (in all but row names) versions of `canterbury_height`, one created using the `[` operator, one created via an intermediary selection object, and another using **sf**'s convenience function `st_filter()`. To explore spatial subsetting in more detail, see the supplementary vignettes on `subsetting` and [`tidyverse-pitfalls`](https://geocompr.github.io/geocompkg/articles/) on the [geocompkg website](https://geocompr.github.io/geocompkg/articles/). The next section explores different types of spatial relation, also known as binary predicates, that can be used to identify whether or not two features are spatially related or not. From a50f4795609b1997dc4e0d8f20e1ae05d0bd083b Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 13 Dec 2021 15:38:16 +0000 Subject: [PATCH 098/104] Improve description of st_filter --- 04-spatial-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 48f40726c..871a45272 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -122,7 +122,7 @@ Note: another way to return a logical output is by setting `sparse = FALSE` (mea Note: the solution involving `sgbp` objects is more generalisable though, as it works for many-to-many operations and has lower memory requirements. ``` -Note: since late 2019, the same result can be achieved with the function `st_filter()`: +The same result can be achieved with the **sf** function `st_filter()` which was [created](https://github.com/r-spatial/sf/issues/1148) to increase compatibility between `sf` objects and **dplyr** data manipulation code: ```{r} canterbury_height3 = nz_height %>% From dc2725e7f9fc88959c92069503a6ea7e5f5bab38 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 13 Dec 2021 15:40:12 +0000 Subject: [PATCH 099/104] Comment out link to tidyverse-pitfalls vignette --- 04-spatial-operations.Rmd | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 871a45272..fde8902f4 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -144,7 +144,8 @@ waldo::compare(canterbury_height2, canterbury_height4) ``` At this point, there are three identical (in all but row names) versions of `canterbury_height`, one created using the `[` operator, one created via an intermediary selection object, and another using **sf**'s convenience function `st_filter()`. -To explore spatial subsetting in more detail, see the supplementary vignettes on `subsetting` and [`tidyverse-pitfalls`](https://geocompr.github.io/geocompkg/articles/) on the [geocompkg website](https://geocompr.github.io/geocompkg/articles/). + + The next section explores different types of spatial relation, also known as binary predicates, that can be used to identify whether or not two features are spatially related or not. ### Topological relations From ce9cb74c418ce662f88fb25b3267597cbd9cf773 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 13 Dec 2021 15:41:39 +0000 Subject: [PATCH 100/104] Shorten sentence --- 04-spatial-operations.Rmd | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index fde8902f4..3174c51cc 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -156,7 +156,8 @@ That may sound rather abstract and, indeed, the definition and classification of See @dieck_algebraic_2008 for an updated textbook teaching algebraic topology. ] -Despite their mathematical origins, topological relations can be understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships defined in English (and other human languages) as illustrated in Figure \@ref(fig:relations), which shows a variety of geometry pairs and their associated relations. +Despite their mathematical origins, topological relations can be understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships. +Figure \@ref(fig:relations) shows a variety of geometry pairs and their associated relations. The third and fourth pairs in Figure \@ref(fig:relations) (from left to right and then down) demonstrate that, for some relations, order is important: while the relations *equals*, *intersects*, *crosses*, *touches* and *overlaps* are symmetrical, meaning that if `function(x, y)` is true, `function(y, x)` will also by true, relations in which the order of the geomtries are important such as *contains* and *within* are not. \index{topological relations} From 5ff0674ab14eff8317c86ab974b1fbca3fbcee57 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 13 Dec 2021 15:42:35 +0000 Subject: [PATCH 101/104] Use full epsg code --- 04-spatial-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 3174c51cc..0161637bc 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -401,7 +401,7 @@ random_df = tibble( ) random_points = random_df %>% st_as_sf(coords = c("x", "y")) %>% # set coordinates - st_set_crs(4326) # set geographic CRS + st_set_crs("EPSG:4326") # set geographic CRS ``` The scenario illustrated in Figure \@ref(fig:spatial-join) shows that the `random_points` object (top left) lacks attribute data, while the `world` (top right) has attributes, including country names shown for a sample of countries in the legend. From bf7572c9a5fc9b905df16fd961f1d711a3ef3817 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 13 Dec 2021 15:43:21 +0000 Subject: [PATCH 102/104] Drop sentence on history of map algebra (!) --- 04-spatial-operations.Rmd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 0161637bc..7fe1399d1 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -726,8 +726,8 @@ The next subsection explores these and related operations in more detail. \index{map algebra} The term 'map algebra' was coined in the late 1970s to describe a "set of conventions, capabilities, and techniques" for the analysis of geographic raster *and* (although less prominently) vector data [@tomlin_map_1994]. -Although the concept never became widely adopted, the term usefully encapsulates and helps classify the range operations that can be undertaken on raster datasets. -In this context we define map algebra more narrowly, as operations that modify or summarise raster cell values, with reference to surrounding cells, zones, or statistical functions that apply to every cell. + +In this context, we define map algebra more narrowly, as operations that modify or summarise raster cell values, with reference to surrounding cells, zones, or statistical functions that apply to every cell. Map algebra operations tend to be fast, because raster datasets only implicitly store coordinates (hence the oversimplifying phrase "raster is faster but vector is corrector"). The location of cells in raster datasets can be calculated it using its matrix position and the resolution and origin of the dataset (stored in the header). From 9ff09a16f63c0dff4e0333fbcdfa3edc3b999681 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 13 Dec 2021 15:45:32 +0000 Subject: [PATCH 103/104] Remove double is --- 04-spatial-operations.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/04-spatial-operations.Rmd b/04-spatial-operations.Rmd index 7fe1399d1..65c5a1a2e 100644 --- a/04-spatial-operations.Rmd +++ b/04-spatial-operations.Rmd @@ -34,7 +34,7 @@ The *type* of spatial relationship between objects must be considered when under Another unique aspect of spatial objects is distance: all spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship, as described in the context of vector data in Section \@ref(distance-relations). Spatial operations on raster objects include subsetting --- covered in Section \@ref(spatial-raster-subsetting) --- and merging several raster 'tiles' into a single object, as demonstrated in Section \@ref(merging-rasters). -*Map algebra* covers a range of operations that modify raster cell values, with or without reference to surrounding cell values, is is vital for many applications. +*Map algebra* covers a range of operations that modify raster cell values, with or without reference to surrounding cell values, is vital for many applications. The concept of map algebra is introduced in Section \@ref(map-algebra); local, focal and zonal map algebra operations are covered in sections \@ref(local-operations), \@ref(focal-operations), and \@ref(zonal-operations), respectively. Global map algebra operations, which generate summary statistics representing an entire raster dataset, and distance calculations on rasters, are discussed in Section \@ref(global-operations-and-distances). In the final section before the exercises (\@ref(merging-rasets)) the process of merging two raster datasets is discussed and demonstrated with reference to a reproducible example. From 70104db039602d39daffd56946f31aa4f16e3683 Mon Sep 17 00:00:00 2001 From: Robin Lovelace Date: Mon, 13 Dec 2021 15:47:04 +0000 Subject: [PATCH 104/104] Replace as.raster() with raster() etc --- 02-spatial-data.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/02-spatial-data.Rmd b/02-spatial-data.Rmd index a9f66583d..a917e489a 100644 --- a/02-spatial-data.Rmd +++ b/02-spatial-data.Rmd @@ -721,7 +721,7 @@ The **terra** package supports raster objects in R. Like its predecessor **raster** (created by the same developer, Robert Hijmans), it provides an extensive set of functions to create, read, export, manipulate and process raster datasets. **terra**'s functionality is largely the same as the more mature **raster** package, but there are some differences: **terra** functions are usually more computationally efficient than **raster** equivalents. -On the other hand, the **raster** class system is popular and used by many other packages; as with **sf** and **sp**, the good news is you can seamlessly translate between the two types of object to ensure backwards compatibility with older scripts and packages, for example, with the functions `as.raster(my_rast)` (see the previous chapter for more on the evolution of R packages for working with geographic data). +On the other hand, the **raster** class system is popular and used by many other packages; as with **sf** and **sp**, the good news is you can seamlessly translate between the two types of object to ensure backwards compatibility with older scripts and packages, for example, with the functions `raster()`, `stack()`, or `brick()` (see the previous chapter for more on the evolution of R packages for working with geographic data). ```{r, echo=FALSE, eval=FALSE} # # test raster/terra conversions