--- 
title: 'Introduction to open source GIS '
author: "Duncan Golicher"
date: "`r Sys.Date()`"
bibliography:
- book.bib
- packages.bib
- refs2.bib
biblio-style: apalike
description: 'Introduction to open source GIS:'
documentclass: book
link-citations: yes
site: bookdown::bookdown_site
always_allow_html: yes
---

# Introduction to the course

Placeholder


## Concepts: What is a GIS?
### Desktop GIS software
### Spatial data bases
### Analysing Geographical data using scripts
### Some key GIS concepts
#### Coordinate reference systems and projections.
#### Raster and vector data
#### Map making and cartography
### Exercise

<!--chapter:end:index.Rmd-->


# Getting started: Projects and basemaps

Placeholder


## Concepts
## QGIS: The GUI interface
### Starting a QGIS project
### Adding base maps
#### Exercise
## Getting started in RStudio
### Using RStudio
### Base maps in R {#basemaps}
#### Exercise
### Screenshots of a new user starting their first web map in R

<!--chapter:end:011_starting_a_project.Rmd-->


# Quick tour of giscourse R functions

Placeholder


## Introduction
## Quick map
## Adding a WMS
## WMS with query
## Deforestation mapping
## Capturing a polygon or set of points
## More examples from the package

<!--chapter:end:011a_giscourse-quick-start.Rmd-->


# Spatial data bases {#postgis}

Placeholder


## Concepts
##### WMS visualisation of all priority habitats and SSSIs in South of England
### Data base work flow
### The sf package in R
## QGIS:Connecting to postgis
### Providing username and password
### Connecting and adding layers
### Inspecting the tables
### Structure of the tables
## RStudio:Connecting to the data base {#conr}
### Using the sf package
### Loading only geometries
### Disconnect from the data base

<!--chapter:end:012_connecting_to_the_server.Rmd-->


# Coordinate reference systems {#crs}

Placeholder


## Concepts
## QGIS: Projection
### Checking the CRS of each layer 
### Reprojecting on the fly
### Reprojecting the layer itself
## R: Coordinate reference systems
### Finding the CRS
### Global projections using proj4
### Checking the coordinate reference system used in postgis
### Transforming
#### Transforming in postgis
### On the fly reprojection with mapview

<!--chapter:end:013_projections.Rmd-->


# Attribute queries {#attributes}

Placeholder


## Concepts
### Selecting subsets of data
### Working with single large tables in PostGIS
## QGIS
### Querying on attributes of loaded layers
## R: Querying attributes
### Listing the tables in R
### Example: Selecting the United Kingdom from a table of countries
### Another example
### More specific queries
### Querying numerical values
### Boolean operators
### Filtering directly in R
### Forming tables of results in the data base
### Easy testing app
### Conclusion

<!--chapter:end:020_attribute_sql.Rmd-->


# Spatial queries {#spatialqueries}

Placeholder


## Concepts
## QGIS
### Quick and easy spatial subsetting
### Using SQL
#### Breaking down the query {#sqlexp}
## Spatial queries in R
### Sharing the results of polygon capture in QGIS
### Joining using sprintf.
### Arne example
### Forming new geometries {#geoprocintro}
### Forming a new schema and adding tables
### Using OSM vector data
### Quick and easy spatial queries using Geocoding

<!--chapter:end:030_spatial_queries.Rmd-->


# Making maps; Digitising and cartography 

Placeholder


## QGIS
### Digitising in QGIS
### Cartography in QGIS
## R: Map making
### Web maps in R
### Saving web maps
### Base maps
### Geocoding {#geocode}
## Drawing and saving your own features in R
### Example for points
## Printable maps in R
### Using qtm
### Animated charts in R
### Animated bubble charts
#### Adonis blue animated map
### Animated barcharts
### Work in progress!

<!--chapter:end:031_making_maps.Rmd-->

# Geoprocessing {#geoprocessing1}

## Concepts

Geoprocessing is the term given to any operation that alters the geometry of spatial data. An example of geoprocessing has already been shown, in section \@ref(geoprocintro). There are many books and web bages that show details of the common vector geoprocesing operations, often in the form of venn diagrams. These won't be repeated here.  A good, comprehensive description of geoprocessing with R code to run it is provided by Lovelace's book.

https://geocompr.robinlovelace.net/geometric-operations.html#geo-vec

The wide range of possible geoprocessing operations might give the impression that geoprocessing provides lots of different answers in search of a problem. Geoprocessing is used very frequently when doing GIS. Sometimes it is a principal analytical method but it is also often needed in order to sort out rather mundane issues such as minor topological errors and generally put the geometry into shape ready for further analysis.


```{block2, type ="rmdnote",echo= TRUE}
Geoprocessing can be applied whenever a change to the geometry of the layers you are interested in can be achieved by applying a rule or a set of rules. 

The alternative to geoprocessing is to edit the layers manually. This may be necessary if you can't find a rule that does the job that you want. However it is usually much more time consuming.
```

The geoprocessing toolbox is the general problem solving workshop of vector GIS. It is quite difficult to teach the use of geoprocessing. Every problem that arises tends to be different in nature. Geoprocessing tools usually have to be used in a series of steps. As you learn what geoprocessing does it becomes increasingly easy to implement these steps and choose a method that will fix a problem you have.
Using R scripts or postgis makes the steps quick and reproducible. Postgis can be used directly for server side geoprocessing on very large data sets.

```{block2, type ="rmdwarn",echo= TRUE}
**Think carefully when running geoprocessing operations. It is important to be aware of the connection between the attribute table in the original layer(s) and the new geometry.** As geoprocessing works on the geometry the match between geometry and attributes often breaks down. In particular, areas and lengths usually have to be recalculated. 
```

## QGIS

Although geoprocessing is best achieved using code, the QGIS GUI produces visual output quickly and easily. It is very useful when learning geoprocessing techniques as the results can be seen quickly.


```{block2, type ="rmdnote",echo= TRUE}
Geoprocessing in QGIS became much easier when the program introduced temporary scratch pad layers. This is a great feature of QGIS. A series of operations can be built up using intermediate temporary layers without the need to save each one. This avoids filling a project folder with unecssary layers. However always remember that the layers are temporary and will be lost if not saved.
```



The geoprocessing tools are found as the first element below the vector heading on the top menu of QGIS (they can also be found in the toolbox if that is open).

If you are new to the concept it is well worth opening the menu and browsing the tools. Each one has a clear textual description of what it does along side. There is no point in repeating this text here.

Let's put geoprocessing to use.

**Action**

1. Load the Arne priority habitat layer from the Arne schema produced by the code in section \@ref(geoprocintro) from the postgis connection.
2. Right click on the layer to open the properties menu
3. Go to symbology and choose categorized (not shown in the figures below)
2. Shade the layer using the main habitat category.

```{r geoproc1, fig.cap="Shading the Arne priority habitats layer. Choose random colours based on the category. If the layer is going to be used for cartography the choice of colours would be important, but it is not so critical when you are just working on an analysis.",echo=FALSE}

if(knitr::is_html_output()) {
 knitr::include_graphics("gifs/geoproc_1.gif") 
}
```


```{r geoproc2, fig.cap="Add some basemaps and check the alignment of the priority habitats with features on the ground.",echo=FALSE}

if(knitr::is_html_output()) {
 knitr::include_graphics("gifs/geoproc_2.gif") 
}
```


```{r geoproc3, fig.cap="Zooming in to compare the priority habitat delimitation.",echo=FALSE}

if(knitr::is_html_output()) {
 knitr::include_graphics("gifs/geoproc_3.gif") 
}
```


**Action**

1. Inspect the layer with relation to some base maps. Zoom in and out. Notice that the priority habitats are divided into sites. For many uses this is **NOT** what we want. 

**A problem to solve**

For many purposes it would be preferable to join the elements together to form a continuous layer with no breaks. This is an issue that geoprocessing can solve. The idea is to **dissolve** the geometry in order to simplify it.

1. Click on the geoprocessing tools and find the **dissolve** option
2. Dissolve the geometry using the main habitat field as the input.
3. Click run
4. Right click the parent layer and find copy style. Copy the style to the clipboard then paste it into the new map style.
3. Look at the results

These steps are all shown in figure \@ref(fig:geoproc4)




```{block2, type ="rmdwarn",echo= TRUE}
Temporary scratch layers are used for each step of the analysis which are lost when the program is closed. This is very much a **feature**, not a bug. It is far simpler to discard all intermediate steps in an analysis than to fill a folder with saved layers.  
However always remember to save the final layer!

```



```{r geoproc4, fig.cap="Dissolving the priority habitats layer on the main habitat attribute. Go to the vevtor menu, choose dissolve then make sure that the main habitat attribute is selected. Click run. The figure also shows the style being copied to the new layer.",echo=FALSE}

if(knitr::is_html_output()) {
 knitr::include_graphics("gifs/geoproc_4.gif") 
}
```




```{block2, type ="rmdnote",echo= TRUE}
Notice that the style from the parent layer is being copied to each of the new layers. You can do this by right  clicking the layer, copying the style and then choosing paste style. This is a convenient way to make sure that the colour schemes of each of the layers match up so you can check the results..

```

In figure \@ref(fig:geoproc4) you can see the results of the dissolve operation. Most of the dividing lines between the habitats have been dissolved. The divisions based on paths are still present. However there are still a few odd slivers and unwanted lines.

```{r geoproc5, fig.cap="Pasting the style into the dissolved layer and looking at the results by switching between the original and the dissolved layer",echo=FALSE}

if(knitr::is_html_output()) {
 knitr::include_graphics("gifs/geoproc_5.gif") 
}
```

Although dissolving the layers has produced more or less the result we expected, there are still a few unwanted lines.  A trick to clean up these is to run a very narrow buffer. This is shown in figure \@ref(fig:geoproc6) and figure  \@ref(fig:geoproc7)

```{r geoproc6, fig.cap="Using a very small buffer to clean up a dissolved layer: This trick is often useful",echo=FALSE}

if(knitr::is_html_output()) {
 knitr::include_graphics("gifs/geoproc_6.gif") 
}
```


```{r geoproc7, fig.cap=".",echo=FALSE}

if(knitr::is_html_output()) {
 knitr::include_graphics("gifs/geoproc_7.gif") 
}
```

### Buffering and differencing

A common operation is to find all the features beyond a certain distance from others. For example all the habitat at Arne that is further than 50 meters from paths. This is a simple buffer operation followed by differencing.

**Action**

1. Load the paths layer from the Arne schema.
2. Run a 50m buffer around it. See figure  \@ref(fig:buffer1)
3. Take the difference between the habitat layer and the buffered paths See figure  \@ref(fig:buffer2)
4. Inspect the results figure \@ref(fig:buffer3)

```{block2, type ="rmdnote",echo= TRUE}
Notice that when running a buffer you often choose to dissolve the results. This is can be a good idea as it produces a simpler layer. In this case it doesn't alter the final result that is provided by the differencing operation. You may want to experiment with this.

```



```{r buffer1, fig.cap="Forming a 50m buffer around the paths",echo=FALSE}

if(knitr::is_html_output()) {
 knitr::include_graphics("gifs/buffer_1.gif") 
}
```



```{r buffer2, fig.cap="Taking the difference between the paths buffer and the ",echo=FALSE}

if(knitr::is_html_output()) {
 knitr::include_graphics("gifs/buffer_2.gif") 
}
```



```{r buffer3, fig.cap="Looking at the results of the buffering and differencing",echo=FALSE}

if(knitr::is_html_output()) {
 knitr::include_graphics("gifs/buffer_3.gif") 
}
```


## R: Geoprocessing using postgis and R

```{r}
knitr::opts_chunk$set(echo = TRUE,eval=TRUE)
```


```{r,message=FALSE,warning=FALSE}

library(giscourse)
conn <- connect()
```

Once you have learned enough about geoprocessing operations to know how they can be used, actually implementing geoprocessing can be much quicker using postgis than QGIS. 

The data structure and the functions in the sf (simple features) package in R match those in postgis in many respects.

```{r}
query1<-"select mhabitat,st_buffer(st_union(geom),0.1) geom 
from arne.phabitat group by mhabitat "
habitat<-st_read(conn,query=query1)

```

The burst option in mapviews can be useful when looking at a layer with a smallish number of categories, as this allows each one of the categories to be made visible or invisible separately. 

```{r}
mapview(habitat,zcol="mhabitat",burst=TRUE, legend=FALSE,zoom=FALSE)
```

The paths buffer query can also be run in postgis.

```{r}
library(tmap)
query2<-"select st_buffer(st_union(geom),50) geom 
from arne.paths"
paths<-st_read(conn,query=query2)
qtm(paths)

```

### Joining subqueries using sprintf

The sprintf trick is particularly useful when building a chain of geoprocessing steps into a single query.


```{r}

query3<-sprintf("select mhabitat,st_difference(h.geom,p.geom) 
    geom from (%s) h, (%s) p",query1,query2)


```

This trick saves having to type out all the text of each subquery again, which could easily result in a syntax error. Notice the way sprintf uses %s to represent strings that are replaced by the text in query1 and query2. This is the resulting query. It looks quite complex when all the text is placed together, but effectively this conducts all the operations that were demonstrated in QGIS in one shot in postgis.  

**`r query3`**

One way to learn how to construct queries like this is to experiment with small layers such as those in the Arne schema in the db manager of QGIS.


```{r,message=FALSE}

dif<-st_read(conn,query=query3)
qtm(dif,"mhabitat") +tm_style("classic")
```



```{block2, type ="rmdwarn",echo= TRUE}
Do not attempt to run geoprocessing queries or any other complex operations in postgis on large data tables such as the original priority habitat layers until you are more experienced. A key advantage of postgis is that it is possible to run analysis on very large data sets. In theory all the priority habitat patches further than fifty meters from all the paths in Dorset could be extracted from the data with one query. In practice this query could take hours to run through, so techniques involving prototyping are needed for big data operations. See section x on big data operations in postgis (which I have not written yet!)
```


## Calculating areas after geoprocessing

To summarise the series of operations.

1. The priority habitat layer was dissolved to unite all the areas with the same main habitat classification.
2. A buffer of 50m was produced around the paths in the study area.
3. This buffer was subtracted from the unified priority habitat map.

So the difference layer holds all the habitat patches that are further than 50m from a path. We can now easily find out the total area of these classes of habitat both before and after the differencing. 

```{r,message=FALSE,warning=FALSE}
library(dplyr)
habitat$total_area_ha<-round(st_area(habitat)/10000,1)
dif$core_area_ha<-round(st_area(dif)/10000,1)

habitat %>% st_set_geometry(NULL) -> d1
dif %>% st_set_geometry(NULL) -> d2

left_join(d1,d2, by = "mhabitat") %>% DT::datatable()
```

### Multi polygons to single polygons

The table pf results had one row for each habitat type. This was because the SQL geoprocessing query had grouped the results by main habitat class leading to a multi geometry column. If we wanted to know the area of each polygon we need to cast the multipolygon to single polygons. 

```{r, message=FALSE,warning=FALSE}
dif %>% st_cast("POLYGON") -> dif2
dif2$area_ha<-round(st_area(dif2)/10000,1)

mapview(dif2, zcol = "mhabitat", burst=TRUE,legend=FALSE)

```


### Reproducing the same analysis for a different site

The big advantage of using R lies in the reproducibility of all the analyses. Although it may have taken some time to develop an analysis for Arne it can now easily be used in a different setting. 

What if we want to move the focus to Hengistbury head? This is a local nature reserve with a polygon delimiting the boundary in the local_reserves table. So the query below which just drops in the new area will do the trick to extract the priority habitat.

```{r,message =FALSE}

hh_query<-"select * from local_reserves where lnr_name = 'Hengistbury Head'"

st_read(conn,query=hh_query) ->hh
qtm(hh)
```


```{r}
hh_habitat<-sprintf("select main_habit mhabitat, 
s1habtype,s2habtype,s3habtype habitat,
st_intersection(h.geom,s.geom) geom from 
(%s) s,
ph_v2_1 h 
where st_intersects(h.geom,s.geom)",hh_query)

habitats<-st_read(conn, query=hh_habitat)
qtm(habitats, "mhabitat") +tm_style("classic")
```



The query below extracts the paths.

```{r}
hh_paths<-sprintf("select st_intersection(way,geom) geom from  (%s) s,
dorset_line l where st_intersects(way,geom)", hh_query)

paths<-st_read(conn, query=hh_paths)
qtm(paths)
```


Now to reproduce the buffering analysis for Hengistbury head we just need to drop in  queries that returned results for Hengistbury head instead of those which pulled the layers from the Arne schema. The sprintf trick is invaluable again here, as typing out the full query would be very tedious and error prone. As it is, all the code has been tested for Arne and all that is needed is to place the Hengistbury select queries within brackets.


```{r,message=FALSE}
query1<-sprintf("select mhabitat,st_buffer(st_union(geom),0.1) geom 
from (%s) s group by mhabitat ",hh_habitat)

query2<-sprintf("select st_buffer(st_union(geom),50) geom 
from (%s) s",hh_paths)

query3<-sprintf("select mhabitat,st_difference(h.geom,p.geom) 
    geom from (%s) h, (%s) p",query1,query2)

dif<-st_read(conn,query=query3)

qtm(dif,"mhabitat") +tm_style("classic")
```

### Exercises

We can see that there is much less undisturbed habitat of any type at Hengistbury head. This could be quantified and experiments run with different buffer widths.




```{r,message=FALSE,results='hide',warning=FALSE,echo=FALSE}
if(exists("conn"))dbDisconnect(conn)
rm(list=ls())
lapply(paste('package:',names(sessionInfo()$otherPkgs),sep=""),detach,character.only=TRUE,unload=TRUE)
```


<!--chapter:end:040_geoprocessing.Rmd-->


# Terrain analysis {#terrain}

Placeholder


## Concepts
## R
### Required packages
### Summaries of raster layers
### Visualising the layers
### Deriving new terrain layers 
### Saving the results and passing to QGIS through postgis
### Using SAGA commands from R
### Adding raster layers to postgis
### Getting a DTM for any place in the world
### Another volcano using the same code.
### Load the results to postgis for work in QGIS
## Using QGIS
### Using the DSM for three dimensional visualisation

<!--chapter:end:060_terrain_analysis.Rmd-->


# Combining raster and vector data {#extracting}

Placeholder


## Loading data
### Loading the Arne raster data
#### Vector data
### Producing a canopy height model
### Simulating quadrats
#### Extracting values
## Finding the trees
#### Forming a distance matrix
### Exercises

<!--chapter:end:061_terrain_analysis.Rmd-->


# Map algebra

Placeholder


## Concepts 
## Map algebra in R
### Getting Worldclim layers {#worldclim}
#### The world altitude layer
#### SRTM data for the UK
#### Cropping a raster layer in R {#crop}
### Raster statistics
### Multiple rasters held in a stack (or brick)
### Calculating the mean monthly temperature
## Exercise
## Extra bonus: Making climate diagrams for anywhere in the World

<!--chapter:end:070_map_algebra.Rmd-->


# Species distribution analysis {#sdm}

Placeholder


## Concepts
## Using R and the dismo package
### Getting records from gbif
#### Simplifying the data {#gbif}
#### Writing the records to postgis
### Looking at the records 
### Aggregating counts to a grid
### Visualising the gridded data
### Testing the association between species occurrences
#### Association with geology
### Modelling species distribution using maxent {#maxent}
#### Step 1. Get some occurrence data.
#### Step 2. Define the area of interest
#### Step 3. Choose the predictors

<!--chapter:end:080_species_distribution_modelling.Rmd-->


# Large raster layers {#geoserver}

Placeholder


## Concepts
### Visualising a large raster layer
### Using a geoserver to cut out site level data
## Connecting to the geoserver using QGIS
### WMS connection
#### Loading WMS layers
#### Setting up a WCS connection
#### Loading WCS layers {#wcsload}
#### Saving WCS layers
#### Cleaning up
## R: Processing large raster layers
#### Loading a single reserve polygon from wdpa {#wdpa}
### Terrain analysis
### Adding WMS tiles from the BU geoserver {#buwms}
### Accessing the raw data
### Building a VRT to reference all the data files
### Referencing the data as raster layers in R
### Cropping to the reserve
### Reprojecting raster layers
### Calculating the area lost
### Exercises

<!--chapter:end:090_using_geoserver.Rmd-->


# Using GRASS algorithms 

Placeholder


## Concepts
### Why use GRASS today?
## Using GRASS from QGIS
## Using GRASS from R
### Network analysis in GRASS
### Getting OSM data using the osmdata library
### Linking to GRASS
#### Cleaning the topology
####  Adding some points for the TSP problem
#### Finding the solution
### Watersheds
### Load a digital evevation model

<!--chapter:end:100_using_GRASS.Rmd-->


# Analysing road networks with pgrouting

Placeholder


## Introduction
## Downloading the data 
## Importing 
## Finding the nodes
## Round trip via Ferndown
## Driving distances 

<!--chapter:end:101_using_pg_routing.Rmd-->


# Appendix: Writing SQL: Some tips and tricks {#tips}

Placeholder


## General advice
## Find table names in R {#findtable}
## Find field names in R {#findfields}
## Count the number of records in a table {#countrecs}
## Find all possible values of a categorical field {#findvals}
## Set limits to queries when testing them {#safelimits}
## Add a unique identifier {#uniqueids}
## DbSendQuery, DbGetQuery and st_read {#dbqueries}
### DbSendQuery {#dbsendquery}
### st_read {#stread}
## dbGetQuery {#dbgetquery}
## Checking coordinate reference system 
## Using sprintf {#sprintf}
### Saving typing field names {#quickfields}
### Using sprintf to join sub queries {#subqueries}
## Using st_simplify {#stsimplify}

<!--chapter:end:190_SQL_tips.Rmd-->


# Appendix: Setting up the postgis database {#setupdb}

Placeholder


## Creating a data base
## Connecting to the server
## Loading the data files
## Loading Open Street Map data
## Loading directly
## National grid graticules
## Lidar data {#lidardata}
## Moving tables from other data bases into the gis_course data base
## Ordenance survey road links and waterways
### Save and restore 
### Restore

<!--chapter:end:200_set_up_Postgis.Rmd-->

# Mapdeck



```{r,eval=FALSE}
key<-"pk.eyJ1IjoiZGdvbGljaGVyIiwiYSI6ImNqdThpbHFobDF5dDQ0NHRhNmlyeGIxOGgifQ.Intjo4_tKHPfznud_gBEmw"

library(mapdeck)
set_token(key)
crash_data = read.csv("https://git.io/geocompr-mapdeck")
crash_data = na.omit(crash_data)
ms = mapdeck_style("light")
mapdeck(style = ms, pitch = 45, location = c(0, 52), zoom = 4) %>%
add_grid(data = crash_data, lat = "lat", lon = "lng", cell_size = 1000,
         elevation_scale = 50, layer_id = "grid_layer",
         colour_range = viridisLite::plasma(6))



```


<!--chapter:end:mapdeck.Rmd-->




<!--chapter:end:plr.Rmd-->

