07-gis-r-vectors.Rmd

---
title: "GIS R vectors"
output:
  html_notebook: default
  html_document: default
---

```{r setup, echo=FALSE, cache=FALSE}
knitr::read_chunk('./_solutions/_solutions_exercises.R')
```

# Setup

Make sure your have `rgdal`, `rgeos`, `raster`, `ggmap`, `leaflet` and `sp` installed:

```{r install, eval=FALSE, include=FALSE}
install.packages(c("rgdal", "rgeos", "raster", "ggmap", "leaflet", "sp"))
```

```{r load, message=FALSE, warning=FALSE}
library('rgdal')
library('sp')
library('raster')
library('rgeos')
library('leaflet')
library('ggmap')
```

# Spatial (vector) data objects: `sp`

Package `sp` is the central package supporting spatial data analysis in R. `sp` defines a set of classes to represent spatial data. A class defines a particular data type. The `data.frame` is an example of a class. Any particular data.frame you create is an **object** (instance) of that class.

In fact, the `sp` package does not provide many functions to modify or analyse spatial data; but the classes it defines are used in >100 other R packages that provide specific functionalities. The classes `sp` defines are mostly starting with `Spatial*`. For vector data, the basic types are the `SpatialPoints`, `SpatialLines`, and `SpatialPolygons`. These classes only represent **geometries** (i.e. the spatial information). To link the spatial information to a data-table, classes are available with these names plus `DataFrame`, for example, `SpatialPolygonsDataFrame` and `SpatialPointsDataFrame`.

## Geometries

Consider the following vectors with coordinate information:
```{r intro_data}
longitude <- c(4.4543, 5.02789, 4.94202, 4.49238, 4.49054, 4.54044, 5.95192, 
               4.49496, 4.4958, 5.68327, 4.49054, 4.49054, 4.4958, 4.48938)
latitude <- c(51.33847, 51.24824, 51.24325, 51.26218, 51.25701, 51.27518, 
              51.30803, 51.25803, 51.26106, 51.04185, 51.25701, 51.25701, 
              51.26106, 51.25955)
lonlat <- cbind(longitude, latitude)
```

Besides the numbers itself, the variable `lonlat` is not aware of any spatial context (it is just a `matrix` object). 

```{r intro_spatialload}
pts <- SpatialPoints(lonlat)
class(pts)
```

We converted the vector to an object of `SpatialPoints`, which contains, beside the matrix data, additional information:

```{r intro_show}
showDefault(pts)
```

So we see that the object has the coordinates we supplied, but also a `bbox`. This is a ‘bounding box’, or the ‘spatial extent’ that was computed from the coordinates. There is also a `proj4string`. This stores the coordinate reference system (CRS). We did not provide the crs so it is unknown (NA). That is not good, so let’s recreate the object, and now provide a crs. We'll come back to the projection definition later:

```{r intro_setcrs}
crs_wgs84 <- CRS("+init=epsg:4326")
pts <- SpatialPoints(lonlat, proj4string = crs_wgs84)
showDefault(pts)
```

Hence, the projection information is provided and our geometry is aware of the spatial context. 

**Remark**: Assigning a CRS is like labeling something. You need to provide the label that corresponds to the item. Not to what you would like it to be. The latter requires a transformation.

The definition of the geometry as an object of the `SpatialPoints` class, provides a number of spatial functionalities that are not directly available to *default* matrices. As said in the introduction, the `sp` class mainly focuses on the classes (objects) and not on the functionalities. However, an essential functionality provided by `sp` itself (apart from `over` and `aggregate`), is the **transformation of the coordinate reference system (CRS)**:

```{r intro_transform}
crs_lambert <- CRS("+init=epsg:31370") 
pts_lambert <- spTransform(pts, crs_lambert)
showDefault(pts_lambert)
```

Other (geometric) operations are provided by (a huge amount of) additional R Packages. The `rgeos` package is an important package to check for these functionalities. As an example: adding a buffer to our `SpatialPoints` object:
```{r intro_buffer}
# gbuffer from the rgeos package (see later)
pts_lambert_buffer <- gBuffer(pts_lambert, byid = TRUE, width = 2000) 
```

Hence, the result is again a `Spatial*` object, i.e. a `SpatialPolygons` object. A plot of the `pts_lambert` and the `pts_lambert_buffered` together looks as follows:

```{r intro_plotstatic}
plot(pts_lambert_buffer, border = 'blue', 
     col = 'yellow', lwd = 3, axes = TRUE)
plot(pts_lambert, add = TRUE, col = 'red', las = 1)
```

For quick data exploration of such a small data set, an interactive visualization provides better insight. The [leaflet](https://rstudio.github.io/leaflet/) package is well-documented and easy to use, certainly when the `dplyr` piping is familiar to use:

```{r intro_plotinteract}
wgs_84 <- CRS("+init=epsg:4326")
leaflet() %>%  
    addTiles() %>%  # provide a default openstreetmap background layer to the image
    addCircles(data = pts)  %>% # Add the points to the map
    addPolygons(data = spTransform(pts_lambert_buffer, wgs_84)) # Add the buffer polygons to the map
```


<br><div class="alert alert-success">
<b>EXERCISE</b>

Check the [documentation of the leaflet Package](https://rstudio.github.io/leaflet/) and adapt the following items to the leaflet map:<br>
(1) The coordinates should be replaced by an icon marker instead of just a circle<br>
(2) The buffers should not be filled and have a red color.
</div>

The output should look like:
```{r exercise1, echo=FALSE, eval=TRUE}
```

## Spatial***DataFrame

The geometry classes provided by `sp` contain the spatial information, but do not add data to these points (or lines/polygons). The data attribute table is typically provided as a `data.frame`. For this reason, `sp` also provides classes that link the geometries with DataFrames: `SpatialPointsDataFrame`, `SpatialLinesDataFrame` and `SpatialPolygonsDataFrame`. 

Let's assume that for each of the coordinates defined above, we also have information about an `identifier` and a `scientificName`:

```{r sp_create_data}
identifiers <- c('INBO:NBN:BFN00179000029DL', 'INBO:NBN:BFN001790000A584', 
                 'INBO:NBN:BFN001790000AACI', 'INBO:NBN:BFN0017900009DOC', 
                 'INBO:NBN:BFN001790000A4J0', 'INBO:NBN:BFN0017900009DO8', 
                 'INBO:NBN:BFN001790000AACF', 'INBO:NBN:BFN001790000A4JO', 
                 'INBO:NBN:BFN001790000A4K8', 'INBO:NBN:BFN001790000A51G', 
                 'INBO:NBN:BFN001790000A4J3', 'INBO:NBN:BFN001790000A4J5', 
                 'INBO:NBN:BFN001790000A4K6', 'INBO:NBN:BFN001790000A4L5')
scientific_names <- c('Lagarosiphon major', 'Lagarosiphon major', 
                      'Muntiacus reevesii', 'Muntiacus reevesii', 
                      'Muntiacus reevesii', 'Muntiacus reevesii', 
                      'Muntiacus reevesii', 'Muntiacus reevesii', 
                      'Muntiacus reevesii', 'Muntiacus reevesii', 
                      'Muntiacus reevesii', 'Muntiacus reevesii', 
                      'Muntiacus reevesii', 'Muntiacus reevesii')
recorded <- as.data.frame(cbind(identifiers, scientific_names))
head(recorded)
```

Hence, the `SpatialPointsDataFrame` combines the spatial information (`pts_lambert`) with the recorded data (`recorded`):

```{r sp_create_df}
pts_recorded <- SpatialPointsDataFrame(pts_lambert, data = recorded)
pts_recorded
```

For the moment, the take home message is that `Spatial***DataFrame` bring together the attribute data (`data.frame`), the coordinates (`geometry`) and the coordinate reference system (`crs`): 

```{r sp_data_attrib}
head(pts_recorded@data)
```

```{r sp_coord_attrib}
head(pts_recorded@coords)
```

```{r sp_proj_attrib}
pts_recorded@proj4string
```

## Projection

In the previous section, we already used the `CRS` class definition to specify the coordinate reference system. As many other systems, R uses the [PROJ.4](http://proj4.org/) notation to define the CRS. The number of parameters depends on the projection. However, the easiest way is (mostly) using the [`EPSG` code](http://spatialreference.org/), for example:

```{r proj_epsg, eval=FALSE, include=FALSE}
CRS("+init=epsg:4326")
CRS("+init=epsg:31370")
```

But other options are available as well:

```{r proj_utm}
CRS("+proj=utm +zone=32")
```

In general, the definition of the CRS (input argument) should ba a string setup up by one or more `+<arg>=<value>` combinations, each of them separated by a blank value: `+<arg1>=<value1> +<arg2>=<value2>`. The CRS-creation returns an object (with a variable name you can assign):
```{r proj_objects}
wgs_84 <- CRS("+init=epsg:4326")
wgs_lambert <- CRS("+init=epsg:31370")
wgs_84
wgs_lambert
```

which can be used in other functions/methods or used to **support the transformation of the coordinate reference system** of *spatial aware objects*. Hence, transformation won't work on a default vector or matrix object:

```{r proj_transf_fail, eval=FALSE, include=FALSE}
sp::spTransform(c(98710.32800000161, 162573.7030000016), wgs_84)
```

But after transformation of the single coordinate to a `SpatialPoints` object, the `spTransform` function can be used:
```{r proj_transf_point}
single_coordinate <- SpatialPoints(matrix(c(98710.32800000161, 
                                            162573.7030000016), 
                                          ncol = 2), 
                                   proj4string = wgs_lambert)
                                  
spTransform(single_coordinate, wgs_84)
```

The application is similar to other spatial objects, such as a `SpatialPointsDataFrame`:
```{r proj_transf_df}
pts_recorded_wgs84 <- spTransform(pts_recorded, wgs_84)
pts_recorded_wgs84
```

<div class="alert alert-success">
<b>EXERCISE</b>

As your regularly encounter point coordinates that need conversion to another projection system, you decide to make a functions for future use:<br>

For any `data.frame`(!) with coordinate information in a known coordinate reference system, update the columns with the x/y information to another coordinate system. The function returns the updated `data.frame`.

The function will probably look more or less like this:<br><br>

    reproject_points <- function(df, col_long, col_lat,
                                 project_input, project_output){
        # ...
        
        return(df)
</div>

```{r exercise2, include=FALSE}
```

Following test should be able to work with the `reproject_points` function:
```{r proj_ex_test, message=FALSE}
pts_ex2_test <- as.data.frame(pts_recorded)
ex2_test <- reproject_points(pts_ex2_test, "longitude", "latitude",
                             CRS("+init=epsg:31370"), 
                             CRS("+init=epsg:4326"))
head(ex2_test)
```

**Remark:** This function is part of an [`inbo-rutil` script,  ](https://github.com/inbo/inbo-rutils/blob/master/gis/guess_projection.R) to guess a given projection of a given data set, when no CRS information is provided. The `guess_projection` function creates a *leaflet* plot with the coordinates in a variety of predefined coordinate reference systems.

# Reading data: `rgdal`

Mostly, the starting point of a data analysis involving spatial data is a data set in any kind of GIS format. The functionality to read **vector** data is provide by the `readOGR` function or the `rgdal` package:

```{r read_shape}
deelbekkens <- readOGR("../data/deelbekkens/Deelbekken.shp")
```

```{r read_show_shape}
deelbekkens
```

```{r read_plot_object}
plot(deelbekkens)
```

Besides shapefiles, other vector data formats can be read as well. Certainly `geojson` is a increasingly used data format for feature (vector) data sets:

```{r read_geojson}
eu_10grid <- readOGR("../data/EUgrid10.geojson")
```
```{r read_show_geojson}
eu_10grid
```

Actually, the set of data formats you can interact with, *does not depend on R*, but is dependent on your GDAL installation. To check if the installation supports a specific **vector** data format, you can get an overview of them by the `ogrDrivers()` command:

```{r read_show_drivers}
head(ogrDrivers()) # Just the first 6 records are shown here
```

For data sets with multiple layers, the reading of a specific Layer is supported as well. For example, when the driver to read ESRI FileGeoDatabases is present, following commands will be useful as well:

```{r read_gdb_example, eval=FALSE}
# Information about the data (without actually reading the data itself)
ogrInfo(dsn = "name_filegeodatabase.gdb") 
# List the names of the layers in the dataset
ogrListLayers(dsn = "name_filegeodatabase.gdb")
# Reading a specific layer
readOGR(dsn = "name_filegeodatabase.gdb", layer = "layer_of_interest"
```

# Vector data manipulation

## DataFrame alike manipulation

Vector data represented as `Spatial***DataFrame` can be sub selected as a DataFrame. For example, selecting only the `deelbekkens` that are in the `Netebekken` (Boolean indexing):

```{r man_subset}
nete <- deelbekkens[deelbekkens$BEKNAAM == "Netebekken", ]
plot(nete)
```

The result of the selection is again a `SpatialPolygonsDataFrame`. When only interested in the data itself, the selection of the attribute data alone can be done by converting the data to a DataFrame. 

```{r man_subset_show}
head(as.data.frame(nete))
```

Hence, the geometry and projection data is no longer contained, only the data itself.

**Remark:** The `sp` data types do NOT support the `dplyr` package. If you want to use the `%>%` operator and the *verbs* such as `filter`, `mutate`,... as used in the `dplyr` package, you could check and install the [`spdplyr`](https://cran.r-project.org/web/packages/spdplyr/spdplyr.pdf) package. Another option is to start using the [`sf` package](https://cran.r-project.org/web/packages/sf/index.html), the newest package setup with the `tidyverse` environment in mind. [This introduction](http://strimas.com/r/tidy-sf/) and the [documentation of](https://github.com/edzer/sfr) of `sf` is good to start.

## `sp` overlay and aggregations

The `sp` package provides mostly the data types. However, the `spTransform` function is an essential function provided by the package itself. Two other functions of the `sp` package are worth mentioning:
*  `disaggregate`: groups of `Line` or `Polygon` (multi) are disaggregated to one `Line`/`Polygon` for each feature. When applied to the `deelbekkens` data set:

```{r sp_disaggreg}
length(deelbekkens)
deelbekkens_disaggr <- disaggregate(deelbekkens)
length(deelbekkens_disaggr)
```
The number of features increases, as each polygon in the data set is now defined as an individual feature of the data set.

* `over(x, y)` is used for spatial JOINs 

```{r over_join}
ids <- over(spTransform(eu_10grid, wgs_lambert), pts_recorded)
ids <- na.omit(ids)
ids
```
The result provides the indices of the `SpatialPolygonDataFrame` together with the `SpatialPointDataFrame` data entries.

**Remark** `Lines`, and `Polygon`-`Polygon` overlays require `rgeos`

<br><div class="alert alert-info">
Check the vignette (`vignette('over')` in console) to get more details about all the different spatial JOIN options.
</div>

## `rgeos` geometry operations

The `rgeos` package provide a number of geometric operations, such as `gArea()`, `gLength()`, `gDistance()`, `gBuffer()`,... (single geometry) or `gIntersection()`, `gUnion()`, `gContains`,... (combining geometries). 

<br><div class="alert alert-info">
As most commands of rgeos start with g***, use the power of the TAB-button to explore the available functions of the `rgeos` package: `rgeos::g` + TAB
</div>

These functions operate both on `Spatial*` objects (see above) as well as `Spatial***DataFrame` objects. Most of these functions speak for themselves (or read the manual). An important argument of these functions is `byid`. The `byid` argument determines if the function should be applied across sub geometries (`TRUE`) or the entire object (`FALSE`). Default is `FALSE`. For example, for the `area` calculation:

```{r man_area}
gArea(nete) # default is FALSE -> entire object
```
versus
```{r man_area_byid}
gArea(nete, byid = TRUE)
```

Comparison of geometries involves two geometries:

```{r man_contains}
vertical <- eu_10grid[grep('E400', eu_10grid$CellCode), ]
horizontal <- eu_10grid[grep('N311', eu_10grid$CellCode), ]

comparison <- gContains(vertical, horizontal, 
                        byid = TRUE, checkValidity = TRUE, 
                        returnDense = TRUE) # Test also with gIntersects
true_ids <- as.data.frame(which(comparison, arr.ind = TRUE)) # Find the TRUE elements
true_ids
```

```{r man_contains_plot}
# adjusting the colors transparency
blue <- adjustcolor("blue", alpha.f = 0.3) 
red <- adjustcolor("darkred", alpha.f = 0.3) 

# plot the cells selected by the geometry comparison:
plot(vertical)
plot(vertical[true_ids$col,], add = TRUE, col = blue)
plot(horizontal, add = TRUE)
plot(horizontal[true_ids$row,], add = TRUE, col = red)
```

<br><div class="alert alert-success">
<b>EXERCISE</b>

Make a `leaflet` plot (with open street map background) of the `nete` sub catchments, together with the `Centroid` for each of the sub geometries.

<b>Tip 1:</b> Looking for a command that matches the action `centroid` in R, is done by using `??centroid` in the Console.<br>
<b>Tip 2:</b> Remember that leaflet always expects WGS84 as CRS to make a plot(!)
</div>

The output should look like:
```{r exercise3, echo=FALSE, eval=TRUE}
```

# Mapping vector data

## `ggplot`

When working with DataFrame data, the application of `ggplot2` is typical. Hence, it would be convenient to use `ggplot2` (and the syntax) to plot `Spatial***DataFrame` as well. As `ggplot2` only works on `data.frame` objects, a translation need to be made between the `Spatial***DataFrame` and a default `data.frame`. 

For `SpatialPointsDataFrame` this is not really an issue, as the `longitude` and `latitude` columns provide the geometry information and can be used as such in `ggplot` by directly converting the data to a `data.fram`

```{r plot_point_df}
pts_recorded_df <- as.data.frame(pts_recorded)
pts_recorded_df
```

```{r plot_point_plot}
ggplot(pts_recorded_df, aes(x = longitude, y = latitude)) +
    geom_point() +
    geom_text(size = 4, aes(label = scientific_names), 
              check_overlap = TRUE, hjust = -0.1) +
    xlim(150000, 300000)
```

For `Lines` and `Polygons`, the translation requires an entire interpretation of the configuration towards a `data.frame`. This functionality is provided by `fortify`, which translates the **spatial** information:

```{r plot_df_geom, message=FALSE}
nete_df <- fortify(nete)
head(nete_df)
```

```{r plot_df_geom_show}
class(nete_df)
```

The individual columns of the resulting `data.frame` describe the following elements:

* `long` and `lat`: x and y coordinates
* `order`: sequence of the polygones or lines
* `hole`: identifies whether the feature the coordinates are part of define a hole
* `piece`: collects all the points that make a single feature (polygon/line), i.e. 1 attributes table line
• `group`: defines the subpolygons or sublines (in case of multipolygon/multiline features)
• `id`: row identifiers of the attribute table (data)

Still, the added value of using `ggplot2` is to use the attribute data to define colors,... As the attribute data table is available and we have an identifier `id`, we can JOIN the attribute data to the spatial information after retrieving the **attribute** data with the `as.data.frame` function:

```{r plot_df_attr}
nete_df <- merge(nete_df, as.data.frame(nete))
head(nete_df)
```

<br><div class="alert alert-info">
<b>REMEMBER</b>

In order to plot `SpatialLinesDataFrame` or `SpatialPolygonsDataFrame` data, the **spatial data** (with `fortify()`) need to be merged with the **attribute data** (with `as.data.frame`) to have a `ggplot2` compatible `data.frame`!
</div>

This seems to be a useful small tooling for a custom R function...

<br><div class="alert alert-success">
<b>EXERCISE</b>

Write a function called `prepare_spatial_for_ggplot` that takes a `SpatialLinesDataFrame` or `SpatialPolygonsDataFrame` as input and execute the necessary steps (described above: `fortify`/`as.data.frame`/`merge`) to enable the plotting with `ggplot2`. The output of the function is a `data.frame`.
</div>

```{r exercise4, include=FALSE}
```

Following test should be able to work with the `prepare_spatial_for_ggplot` function:
```{r plot_ex_test, message=FALSE}
ex4_test <- prepare_spatial_for_ggplot(nete)
head(ex4_test)
```



For plotting `Lines`, use the `ggplot` function `geom_path()` (and NOT `geom_line()`). For polygons, you can choose:
* `geom_path()`: only plot the borders
* `geom_polygon()`: usa a fill color for the polygons

```{r plot_polygon_plot}
ggplot(nete_df, aes(x = long, y = lat)) +
    geom_polygon(aes(group = group, fill = id)) +
    coord_fixed()
```

**Remark:** Adding `group = group` (grouping according to the by fortify derived groups) is required to plot the individual polygons separately(!)

## `ggmap`

The `ggmap` package is the spatial plotting extension to `ggplot2`. It provides the ability to integrate the graphs of `ggplot2` with spatial backgrounds coming from online resources such as Google maps, [stamen maps](http://maps.stamen.com/#terrain/12/37.7706/-122.3782) and OpenStreetMap. The most convenient way to setup these background maps is to use the `ggmap::get_***` functions. For example, using the `toner-background` of stamen:

```{r plot_ggmap_toner, message=FALSE}
extent <- c(left = 2.54, bottom = 49.46, right = 6.4, top = 51.51)
map <- get_stamenmap(extent, zoom = 8, maptype = "toner-background")
belgium <- ggmap(map)
belgium
```

Just as with the OpenStreetMap integration of the `leaflet` package, these spatial background coming from online resources use `WGS84` (decimal degrees). Make sure to convert the vector data to the WGS84 `CRS`:

```{r plot_ggmap_toner_nete, message=FALSE}
nete_wgs84 <- spTransform(nete, wgs_84)
nete_df_wgs84 <- fortify(nete_wgs84)
nete_df_wgs84 <- merge(nete_df_wgs84, as.data.frame(nete_wgs84))
    
belgium +
    geom_path(aes(x = long, y = lat, group = group), 
          data = nete_df_wgs84, colour = "red")
```

The `bbox` can be defined more specifically to improve the boundaries of the background image. Focusing on the Nete catchment:

```{r plot_ggmap_toner_zoom, message=FALSE}
extent <- c(left = 4.42, bottom = 50.9, right = 5.4, top = 51.4)
map <- get_stamenmap(bbox = extent, zoom = 10, maptype = "toner-hybrid")
ggmap(map) + 
    geom_path(aes(x = long, y = lat, group = group), 
              data = nete_df_wgs84, colour = "red")
```

**Remark:** each of these online providers have different `maptype` options. Check the help of these functions to see the options (e.g. `?get_stamenmap`).

The usage of Google maps enables the option to provide a *google maps search* as center of the image
```{r plot_ggmap_google, message=FALSE}
map <- get_googlemap(center = "Geel Belgium",  scale = 2, 
                     maptype = "terrain", zoom = 9)
ggmap(map) +
    geom_polygon(aes(x = long, y = lat, group = group, fill = id, alpha = 0.5), 
          data = nete_df_wgs84, colour = "white")
```

The [Github repository](https://github.com/dkahle/ggmap) and [the paper](https://journal.r-project.org/archive/2013-1/kahle-wickham.pdf) provide more introductionary material to get you started with `ggmap`. Another potentially useful plotting package is the [`tmap`](https://cran.r-project.org/web/packages/tmap/vignettes/tmap-nutshell.html) package.

<br><div class="alert alert-success">
<b>EXERCISE</b>

Pick **any spatial vector** data set (shapefile, geojson,...) currently on your computer:

(1) `read` the file into R as a Spatial***DataFrame
(2) make a subselection of the dataset using one of the data columns and a condition (e.g. `length  > 10000`)
(3) create a `ggplot` with the color of each spatial entity defined by a chosen data columns
</div>