Spatial Data Analysis in R

class: center, middle, inverse, title-slide

# Spatial Data Analysis in R
### Jakub Nowosad, Paweł Bogawski
### Workshop, 2019

---

class: inverse, left, bottom
# Ecosystem(s)

---
# Ecosystem(s)

<div class="figure" style="text-align: center">
<img src="figs/rb.png" alt="From Roger Bivand's slides 'A practical history of R'" width="80%" />
<p class="caption">From Roger Bivand's slides 'A practical history of R'</p>
</div>

---
# Ecosystem(s)

<div class="figure" style="text-align: center">
<img src="figs/rb2.png" alt="From Roger Bivand's slides 'A practical history of R-sig-geo'" width="80%" />
<p class="caption">From Roger Bivand's slides 'A practical history of R-sig-geo'</p>
</div>

---
# Ecosystem(s)

.pull-left[
R packages:

- **sf**, **raster** - spatial classes
- **dplyr**, **rmapshaper** - processing of atttribute tables/geometries 
- **rnaturalearth**, **osmdata**, **getlandsat** - spatial data download
- **rgrass7**, **RQGIS**, **RSAGA**, **link2GI** - connecting with GIS software
- **gstat**, **mlr**, **CAST** - spatial data modeling
- **rasterVis**, **tmap**, **ggplot** - static visualizations
- **leaflet**, **mapview**, **mapdeck** - interactive visualizations
- many more...
]

.pull-right[
<img src="figs/ecosystem.jpg" width="70%" style="display: block; margin: auto;" />
]

---
# Ecosystem(s)

.pull-left[
Spatial classes:

- **sf**
- **raster**
]

.pull-right[
Visualizations:

- **tmap**
]

---
# Ecosystem(s)

.pull-left[
Spatial classes:

- **sf**
- **raster**
- **sp**
- **stars**
- **terra**
- **spatial**
- ...
]

.pull-right[
Visualizations:

- **tmap**
- **ggplot** 
- **rasterVis**
- **leaflet**
- **mapview**
- **mapdeck**
- ...
]

---
# Prerequisites

.pull-left[
We can install the packages used in the workshop as follows:

```r
pkgs = c(
  "sf",                
  "raster",            
  "dplyr",            
  "tmap",              
  "landscapemetrics",
  "landscapetools",
  "usethis"
)
to_install = !pkgs %in% installed.packages()
if(any(to_install)) {
  install.packages(pkgs[to_install])
}
```
]

.pull-right[
The code below can be used to download all the scripts and data for the workshop:

```r
usethis::use_course("http://bit.ly/biogis19w")
```
]

---
# Exercises

1. Think about how R can help you? 
Why do you want to learn how to use it for spatial data analysis?
2. Have a look at [CRAN Task View: Analysis of Spatial Data](https://cran.r-project.org/web/views/Spatial.html) and check if there are any packages that may be useful to you in your daily work?

---
class: inverse, left, bottom
# Spatial data representation

---
# Spatial data representation

**sf**:
- The **sf** package is the successor of the **sp** package based on the OGC standard Simple Features.
- Combines the functionality of three previous packages: **sp**, **rgeos** and **rgdal**.
- Most of the functions in this package start with a prefix `st_`.
- This package handles additional vector data types (e.g. polygon and multipolygon are two separate classes), allows for easier data processing, and support sfor spatial databases such as PostGIS.
- https://r-spatial.github.io/sf/ and https://github.com/rstudio/cheatsheets/blob/master/sf.pdf

**raster**:
- The **raster** package contains classes and methods representing raster objects and operations.
- It allows raster data to be loaded and saved.
- It allows raster algebra and raster processing.
- It includes a number of additional functions, e.g. for analysis of terrain characteristics.
- It allows you to work on large sets of data.
- ?`raster-package` and http://www.rpubs.com/etiennebr/visualraster

---
# Spatial data reading - vector

.pull-left[

```r
library(sf)
polska = st_read("data/polska.gpkg")
```

```
## Reading layer `polska' from data source `/home/jn/Science/conferences/BioGIS_19/workshop/data/polska.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 171842.5 ymin: 132323.8 xmax: 861798.4 ymax: 775305.9
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
```
- The `st_write()` function can be used for spatial vector data writing.
]

.pull-right[

```r
plot(polska)
```

]

---
# Spatial data structure - vector

```r
polska
```

```
## Simple feature collection with 1 feature and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 171842.5 ymin: 132323.8 xmax: 861798.4 ymax: 775305.9
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
##     name  pop_est gdp_md_est                           geom
## 1 Poland 38476269    1052000 POLYGON ((487928.3 182549, ...
```

---
# Spatial data reading - raster

.pull-left[

```r
library(raster)
dem = raster("data/polska_srtm.tif")
```

- The `writeRaster()` function can be used for spatial raster data writing.
]

.pull-right[

```r
plot(dem)
```

<img src="index_files/figure-html/unnamed-chunk-10-1.png" style="display: block; margin: auto;" />
]

---
# Spatial data structure - raster

```r
dem
```

```
## class      : RasterLayer 
## dimensions : 702, 1202, 843804  (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : 14.125, 24.14167, 48.99167, 54.84167  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /home/jn/Science/conferences/BioGIS_19/workshop/data/polska_srtm.tif 
## names      : polska_srtm 
## values     : -9, 2220  (min, max)
```

---
# Exercises

1. Read the `wlkp_pn.gpkg` file containing the borders of the Wielkopolski National Park into R.
Display this object and view its structure.
What can you say about the contents of this file? 
What type of data does it store? 
What is the coordinate system used?
How many attributes does it contain?
What is its geometry?
2. Read the `wlkp_pn_dem.tif` file containing the digital elevation model for the Wielkopolska National Park area into R.
Display this object and view its structure.
What can you say about the contents of this file? 
What type of data does it store? 
What is the coordinate system used?
How many attributes does it contain?
How many dimensions does it have?
What is the data resolution?

---
class: inverse, left, bottom
# Spatial visualization

---
# Spatial visualization

```r
library(tmap)
```

- Thematic mapping package
- It allows you to create static maps, animated maps and interactive maps.
- It works by adding subsequent layers to visualize and then modifying them.
- https://github.com/mtennekes/tmap - list of additional materials about the **tmap** package

---
# Spatial visualization

.left-code[

```r
tm_shape(polska) + 
  tm_polygons()
```
]

.right-plot[
![](index_files/figure-html/tm1-1.png)
]

---
# Spatial visualization

.left-code[

```r
tm_shape(polska) + 
  tm_polygons() + 
  tm_scale_bar(position = c("left", 
                            "bottom")) 
```
]

.right-plot[
![](index_files/figure-html/tm2-1.png)
]

---
# Spatial visualization

.left-code[

```r
tm_shape(polska) + 
  tm_polygons() + 
  tm_scale_bar(position = c("left", 
                            "bottom")) + 
  tm_compass()
```
]

.right-plot[
![](index_files/figure-html/tm3-1.png)
]

---
# Spatial visualization

.left-code[

```r
tm_shape(polska) + 
  tm_polygons() + 
  tm_scale_bar(position = c("left", 
                            "bottom")) + 
  tm_compass() +
  tm_layout(title = "Poland")
```
]

.right-plot[
![](index_files/figure-html/tm4-1.png)
]

---
# Spatial visualization

```r
tmap_mode("view")
```

```r
tm_shape(polska) + 
  tm_polygons() + 
  tm_layout(title = "Poland")
```

```r
tmap_mode("plot")
```

---
# Spatial visualization

.left-code[

```r
tm_shape(dem) + 
  tm_raster()
```
]

.right-plot[
![](index_files/figure-html/tm5-1.png)
]

---
# Spatial visualization

.left-code[

```r
tm_shape(dem) + 
  tm_raster(style = "cont",
            midpoint = NA)
```
]

.right-plot[
![](index_files/figure-html/tm6-1.png)
]

---
# Spatial visualization

.left-code[

```r
tm_shape(dem) + 
  tm_raster(style = "cont", 
            midpoint = NA,
            palette = "-RdYlGn")
```

```r
tmaptools::palette_explorer()
```
]

.right-plot[
![](index_files/figure-html/tm7-1.png)
]

---
# Spatial visualization

.left-code[

```r
tm_shape(dem) + 
  tm_raster(style = "fixed",
            breaks = c(-99, 0, 300,
                       600, 9999), 
            midpoint = NA,
            palette = "-RdYlGn",
            title = "") + 
  tm_layout(legend.position = c("LEFT", 
                                "BOTTOM"))
```
]

.right-plot[
![](index_files/figure-html/tm8-1.png)
]

---
# Spatial visualization

.left-code[

```r
tm_shape(dem) + 
  tm_raster(style = "fixed",
            breaks = c(-99, 0, 300,
                       600, 9999), 
            labels = c("Depressions", 
                       "Plains", 
                       "Hills",
                       "Mountains"),
            midpoint = NA,
            palette = "-RdYlGn",
            title = "") + 
  tm_layout(legend.position = c("LEFT", 
                                "BOTTOM"))
```
]

.right-plot[
![](index_files/figure-html/tm9-1.png)
]

---
# Spatial visualization

.left-code[

```r
tm_shape(dem) + 
  tm_raster(style = "fixed",
            breaks = c(-99, 0, 300,
                       600, 9999), 
            labels = c("Depressions", 
                       "Plains", 
                       "Hills",
                       "Mountains"),
            midpoint = NA,
            palette = c("#5E8B73",
                        "#DAE97A", 
                        "#EADC70", 
                        "#AF8D5C"),
            title = "") + 
  tm_layout(legend.position = c("LEFT",
                                "BOTTOM"))
```
]

.right-plot[
![](index_files/figure-html/tm10-1.png)
]

---
# Spatial visualization

```r
map1 = tm_shape(dem) + 
  tm_raster(style = "fixed",
            breaks = c(-99, 0, 300, 600, 9999),
            labels = c("Depressions", "Plains", "Hills", "Mountains"),
            midpoint = NA, 
            palette = c("#5E8B73", "#DAE97A", "#EADC70", "#AF8D5C"),
            title = "") +
  tm_layout(legend.position = c("LEFT", "BOTTOM"))
```

```r
tmap_save(map1, "my_first_map.png")
```

---
# Exercises

1. Improve the `map1` object by e.g. adding a scale, north arrow, or title.
Also try to add the Polish border to this map.
2. Read the files `wlkp_pn.gpkg` and `wlkp_pn_gpkg_dem.tif` into R.
Create a map showing the terrain model for the Wielkopolski National Park area.
Save the obtained map to the file "WLKP_YOURNAME.png".
3. Additionally, try repeating the steps in the article [Geocomputation with R: maps extended](https://geocompr.github.io/geocompkg/articles/maps.html) in order to make a hillshade map of Poland.

---
class: inverse, left, bottom
# Data manipulation

---
# Data manipulation

```r
library(dplyr)
```

- The **dplyr** package allows for easy processing of tabular data.
- Its capabilities include selecting specific columns or rows, adding new columns, or performing summaries.
- More information can be found at https://dplyr.tidyverse.org/ and https://r4ds.had.co.nz/transform.html

---
# Attributes manipulation

```r
meteo_data = read.csv("data/polska_meteo_2017.csv", encoding = "UTF-8")
head(meteo_data)
```

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> Kod.stacji </th>
   <th style="text-align:left;"> Nazwa.stacji </th>
   <th style="text-align:right;"> Rok </th>
   <th style="text-align:right;"> Miesiac </th>
   <th style="text-align:right;"> Dzien </th>
   <th style="text-align:right;"> tavg </th>
   <th style="text-align:right;"> precip </th>
   <th style="text-align:right;"> vapor_pres </th>
   <th style="text-align:right;"> humidity </th>
   <th style="text-align:right;"> pressure </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 354150100 </td>
   <td style="text-align:left;"> KOŁOBRZEG </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 3.6 </td>
   <td style="text-align:right;"> 4.6 </td>
   <td style="text-align:right;"> 7.4 </td>
   <td style="text-align:right;"> 92.5 </td>
   <td style="text-align:right;"> 1013.5 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354150100 </td>
   <td style="text-align:left;"> KOŁOBRZEG </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:right;"> 3.5 </td>
   <td style="text-align:right;"> 4.7 </td>
   <td style="text-align:right;"> 6.0 </td>
   <td style="text-align:right;"> 76.9 </td>
   <td style="text-align:right;"> 1011.1 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354150100 </td>
   <td style="text-align:left;"> KOŁOBRZEG </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:right;"> 1.9 </td>
   <td style="text-align:right;"> 10.0 </td>
   <td style="text-align:right;"> 6.7 </td>
   <td style="text-align:right;"> 92.3 </td>
   <td style="text-align:right;"> 1007.1 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354150100 </td>
   <td style="text-align:left;"> KOŁOBRZEG </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:right;"> 2.7 </td>
   <td style="text-align:right;"> 2.3 </td>
   <td style="text-align:right;"> 5.8 </td>
   <td style="text-align:right;"> 76.0 </td>
   <td style="text-align:right;"> 993.9 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354150100 </td>
   <td style="text-align:left;"> KOŁOBRZEG </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> -1.7 </td>
   <td style="text-align:right;"> 6.2 </td>
   <td style="text-align:right;"> 3.9 </td>
   <td style="text-align:right;"> 74.5 </td>
   <td style="text-align:right;"> 1024.3 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354150100 </td>
   <td style="text-align:left;"> KOŁOBRZEG </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 6 </td>
   <td style="text-align:right;"> -4.5 </td>
   <td style="text-align:right;"> 4.0 </td>
   <td style="text-align:right;"> 3.5 </td>
   <td style="text-align:right;"> 77.3 </td>
   <td style="text-align:right;"> 1038.7 </td>
  </tr>
</tbody>
</table>

---
# Attributes manipulation

```r
meteo_stations = st_read("data/polska_stacje.gpkg")
```

```
## Reading layer `polska_stacje' from data source `/home/jn/Science/conferences/BioGIS_19/workshop/data/polska_stacje.gpkg' using driver `GPKG'
## Simple feature collection with 1998 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 172691.6 ymin: 139488.3 xmax: 854361.7 ymax: 846760.3
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +units=m +no_defs
```

```r
plot(st_geometry(meteo_stations))
```

---
# Attributes manipulation

```r
meteo_data_sel = meteo_data %>% 
  filter(Nazwa.stacji %in% unique(meteo_stations$NAZWA_ST))
```

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> Kod.stacji </th>
   <th style="text-align:left;"> Nazwa.stacji </th>
   <th style="text-align:right;"> Rok </th>
   <th style="text-align:right;"> Miesiac </th>
   <th style="text-align:right;"> Dzien </th>
   <th style="text-align:right;"> tavg </th>
   <th style="text-align:right;"> precip </th>
   <th style="text-align:right;"> vapor_pres </th>
   <th style="text-align:right;"> humidity </th>
   <th style="text-align:right;"> pressure </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 3.4 </td>
   <td style="text-align:right;"> 7.0 </td>
   <td style="text-align:right;"> 7.5 </td>
   <td style="text-align:right;"> 95.4 </td>
   <td style="text-align:right;"> 1013.5 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:right;"> 2.3 </td>
   <td style="text-align:right;"> 0.9 </td>
   <td style="text-align:right;"> 5.8 </td>
   <td style="text-align:right;"> 80.6 </td>
   <td style="text-align:right;"> 1010.7 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:right;"> 1.4 </td>
   <td style="text-align:right;"> 8.0 </td>
   <td style="text-align:right;"> 6.4 </td>
   <td style="text-align:right;"> 93.0 </td>
   <td style="text-align:right;"> 1007.0 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:right;"> 2.3 </td>
   <td style="text-align:right;"> 6.9 </td>
   <td style="text-align:right;"> 5.6 </td>
   <td style="text-align:right;"> 77.0 </td>
   <td style="text-align:right;"> 992.7 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> -2.9 </td>
   <td style="text-align:right;"> 1.5 </td>
   <td style="text-align:right;"> 2.8 </td>
   <td style="text-align:right;"> 55.6 </td>
   <td style="text-align:right;"> 1023.6 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 6 </td>
   <td style="text-align:right;"> -5.7 </td>
   <td style="text-align:right;"> 0.9 </td>
   <td style="text-align:right;"> 3.2 </td>
   <td style="text-align:right;"> 79.5 </td>
   <td style="text-align:right;"> 1038.8 </td>
  </tr>
</tbody>
</table>

---
# Attributes manipulation

```r
meteo_data_sel = meteo_data %>% 
  filter(Nazwa.stacji %in% unique(meteo_stations$NAZWA_ST)) %>% 
  mutate(data = as.Date(paste0(Rok, "-", Miesiac, "-", Dzien)))
```

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> Kod.stacji </th>
   <th style="text-align:left;"> Nazwa.stacji </th>
   <th style="text-align:right;"> Rok </th>
   <th style="text-align:right;"> Miesiac </th>
   <th style="text-align:right;"> Dzien </th>
   <th style="text-align:right;"> tavg </th>
   <th style="text-align:right;"> precip </th>
   <th style="text-align:right;"> vapor_pres </th>
   <th style="text-align:right;"> humidity </th>
   <th style="text-align:right;"> pressure </th>
   <th style="text-align:left;"> data </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 3.4 </td>
   <td style="text-align:right;"> 7.0 </td>
   <td style="text-align:right;"> 7.5 </td>
   <td style="text-align:right;"> 95.4 </td>
   <td style="text-align:right;"> 1013.5 </td>
   <td style="text-align:left;"> 2017-01-01 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:right;"> 2.3 </td>
   <td style="text-align:right;"> 0.9 </td>
   <td style="text-align:right;"> 5.8 </td>
   <td style="text-align:right;"> 80.6 </td>
   <td style="text-align:right;"> 1010.7 </td>
   <td style="text-align:left;"> 2017-01-02 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:right;"> 1.4 </td>
   <td style="text-align:right;"> 8.0 </td>
   <td style="text-align:right;"> 6.4 </td>
   <td style="text-align:right;"> 93.0 </td>
   <td style="text-align:right;"> 1007.0 </td>
   <td style="text-align:left;"> 2017-01-03 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:right;"> 2.3 </td>
   <td style="text-align:right;"> 6.9 </td>
   <td style="text-align:right;"> 5.6 </td>
   <td style="text-align:right;"> 77.0 </td>
   <td style="text-align:right;"> 992.7 </td>
   <td style="text-align:left;"> 2017-01-04 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> -2.9 </td>
   <td style="text-align:right;"> 1.5 </td>
   <td style="text-align:right;"> 2.8 </td>
   <td style="text-align:right;"> 55.6 </td>
   <td style="text-align:right;"> 1023.6 </td>
   <td style="text-align:left;"> 2017-01-05 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2017 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 6 </td>
   <td style="text-align:right;"> -5.7 </td>
   <td style="text-align:right;"> 0.9 </td>
   <td style="text-align:right;"> 3.2 </td>
   <td style="text-align:right;"> 79.5 </td>
   <td style="text-align:right;"> 1038.8 </td>
   <td style="text-align:left;"> 2017-01-06 </td>
  </tr>
</tbody>
</table>

---
# Attributes manipulation

```r
meteo_data_sel = meteo_data %>% 
  filter(Nazwa.stacji %in% unique(meteo_stations$NAZWA_ST)) %>% 
  mutate(data = as.Date(paste0(Rok, "-", Miesiac, "-", Dzien))) %>% 
  dplyr::select(-Rok, -Miesiac, -Dzien)
```

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> Kod.stacji </th>
   <th style="text-align:left;"> Nazwa.stacji </th>
   <th style="text-align:right;"> tavg </th>
   <th style="text-align:right;"> precip </th>
   <th style="text-align:right;"> vapor_pres </th>
   <th style="text-align:right;"> humidity </th>
   <th style="text-align:right;"> pressure </th>
   <th style="text-align:left;"> data </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 3.4 </td>
   <td style="text-align:right;"> 7.0 </td>
   <td style="text-align:right;"> 7.5 </td>
   <td style="text-align:right;"> 95.4 </td>
   <td style="text-align:right;"> 1013.5 </td>
   <td style="text-align:left;"> 2017-01-01 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2.3 </td>
   <td style="text-align:right;"> 0.9 </td>
   <td style="text-align:right;"> 5.8 </td>
   <td style="text-align:right;"> 80.6 </td>
   <td style="text-align:right;"> 1010.7 </td>
   <td style="text-align:left;"> 2017-01-02 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 1.4 </td>
   <td style="text-align:right;"> 8.0 </td>
   <td style="text-align:right;"> 6.4 </td>
   <td style="text-align:right;"> 93.0 </td>
   <td style="text-align:right;"> 1007.0 </td>
   <td style="text-align:left;"> 2017-01-03 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> 2.3 </td>
   <td style="text-align:right;"> 6.9 </td>
   <td style="text-align:right;"> 5.6 </td>
   <td style="text-align:right;"> 77.0 </td>
   <td style="text-align:right;"> 992.7 </td>
   <td style="text-align:left;"> 2017-01-04 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> -2.9 </td>
   <td style="text-align:right;"> 1.5 </td>
   <td style="text-align:right;"> 2.8 </td>
   <td style="text-align:right;"> 55.6 </td>
   <td style="text-align:right;"> 1023.6 </td>
   <td style="text-align:left;"> 2017-01-05 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 354160105 </td>
   <td style="text-align:left;"> KOSZALIN </td>
   <td style="text-align:right;"> -5.7 </td>
   <td style="text-align:right;"> 0.9 </td>
   <td style="text-align:right;"> 3.2 </td>
   <td style="text-align:right;"> 79.5 </td>
   <td style="text-align:right;"> 1038.8 </td>
   <td style="text-align:left;"> 2017-01-06 </td>
  </tr>
</tbody>
</table>

---
# Attributes manipulation

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> tavg </th>
   <th style="text-align:right;"> pressure </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 8.685298 </td>
   <td style="text-align:right;"> 960.5246 </td>
  </tr>
</tbody>
</table>

---
# Attributes manipulation

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> Kod.stacji </th>
   <th style="text-align:right;"> tavg </th>
   <th style="text-align:right;"> pressure </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 349190600 </td>
   <td style="text-align:right;"> 9.2780822 </td>
   <td style="text-align:right;"> 1017.513 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 349190625 </td>
   <td style="text-align:right;"> 6.5290411 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 349190650 </td>
   <td style="text-align:right;"> 0.1978082 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 349200660 </td>
   <td style="text-align:right;"> 9.2180822 </td>
   <td style="text-align:right;"> 1014.784 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 349210670 </td>
   <td style="text-align:right;"> 8.7728767 </td>
   <td style="text-align:right;"> 1017.553 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 349220690 </td>
   <td style="text-align:right;"> 8.3227397 </td>
   <td style="text-align:right;"> 1017.737 </td>
  </tr>
</tbody>
</table>

---
# Attribute joining

```r
meteo = meteo_stations %>% 
  left_join(meteo_data_sel, by = c("KOD_SZS" = "Kod.stacji"))
```

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> KOD_SZS </th>
   <th style="text-align:left;"> NAZWA_ST </th>
   <th style="text-align:right;"> tavg </th>
   <th style="text-align:right;"> pressure </th>
   <th style="text-align:left;"> geom </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 249190020 </td>
   <td style="text-align:left;"> BRZEŹNICA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:left;"> c(545976.260510018, 233149.457458902) </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 249199999 </td>
   <td style="text-align:left;"> OSIEK </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:left;"> c(519109.405047427, 231411.908795745) </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 249190040 </td>
   <td style="text-align:left;"> GIERAŁTOWICE </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:left;"> c(527197.290270776, 230151.641154815) </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 249190050 </td>
   <td style="text-align:left;"> RADZISZÓW </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:left;"> c(558306.218144104, 230885.909302064) </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 249190060 </td>
   <td style="text-align:left;"> SIEPRAW </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:left;"> c(569259.6149073, 226821.435334128) </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 249190070 </td>
   <td style="text-align:left;"> KĘTY </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:left;"> c(516182.656254972, 223533.556193128) </td>
  </tr>
</tbody>
</table>

---
# Attribute joining

```r
meteo = meteo_stations %>% 
  inner_join(meteo_data_sel, by = c("KOD_SZS" = "Kod.stacji"))
```

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> KOD_SZS </th>
   <th style="text-align:left;"> NAZWA_ST </th>
   <th style="text-align:right;"> tavg </th>
   <th style="text-align:right;"> pressure </th>
   <th style="text-align:left;"> geom </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 349190625 </td>
   <td style="text-align:left;"> ZAKOPANE </td>
   <td style="text-align:right;"> 6.5290411 </td>
   <td style="text-align:right;"> 0.000 </td>
   <td style="text-align:left;"> c(569804.515915351, 158924.036586795) </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 349190650 </td>
   <td style="text-align:left;"> KASPROWY WIERCH </td>
   <td style="text-align:right;"> 0.1978082 </td>
   <td style="text-align:right;"> 0.000 </td>
   <td style="text-align:left;"> c(571455.865783155, 152132.170449405) </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 349200660 </td>
   <td style="text-align:left;"> NOWY SĄCZ </td>
   <td style="text-align:right;"> 9.2180822 </td>
   <td style="text-align:right;"> 1014.784 </td>
   <td style="text-align:left;"> c(621910.953224735, 196894.713701549) </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 349210670 </td>
   <td style="text-align:left;"> KROSNO </td>
   <td style="text-align:right;"> 8.7728767 </td>
   <td style="text-align:right;"> 1017.553 </td>
   <td style="text-align:left;"> c(699590.639267207, 208053.055658539) </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 349220690 </td>
   <td style="text-align:left;"> LESKO </td>
   <td style="text-align:right;"> 8.3227397 </td>
   <td style="text-align:right;"> 1017.737 </td>
   <td style="text-align:left;"> c(742035.646054704, 183036.071495296) </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 350200570 </td>
   <td style="text-align:left;"> KIELCE-SUKÓW </td>
   <td style="text-align:right;"> 8.4849315 </td>
   <td style="text-align:right;"> 1016.895 </td>
   <td style="text-align:left;"> c(619185.547042676, 328424.723335625) </td>
  </tr>
</tbody>
</table>

---
# Reprojections

.pull-left[

```r
plot(dem)
```

<img src="index_files/figure-html/unnamed-chunk-37-1.png" style="display: block; margin: auto;" />
]

.pull-right[

```r
plot(st_geometry(meteo), axes = TRUE)
```

<img src="index_files/figure-html/unnamed-chunk-38-1.png" style="display: block; margin: auto;" />
]

---
# Reprojections

```r
meteo = meteo_stations %>% 
  inner_join(meteo_data_sel, by = c("KOD_SZS" = "Kod.stacji")) %>% 
  st_transform(4326)
```

---
# Reprojections

```r
tm_shape(dem) +
  tm_raster() +
  tm_shape(meteo) +
  tm_symbols(col = "tavg", palette = "RdBu")
```

---
# Vector-raster interactions

```r
meteo$elev = extract(dem, meteo)
head(meteo)
```

```
## Simple feature collection with 6 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 19.96032 ymin: 49.23253 xmax: 22.3417 ymax: 50.81048
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##     KOD_SZS        NAZWA_ST      tavg pressure                      geom
## 1 349190625        ZAKOPANE 6.5290411    0.000 POINT (19.96032 49.29382)
## 2 349190650 KASPROWY WIERCH 0.1978082    0.000 POINT (19.98182 49.23253)
## 3 349200660       NOWY SĄCZ 9.2180822 1014.784  POINT (20.6886 49.62714)
## 4 349210670          KROSNO 8.7728767 1017.553 POINT (21.76917 49.70674)
## 5 349220690           LESKO 8.3227397 1017.737  POINT (22.3417 49.46647)
## 6 350200570    KIELCE-SUKÓW 8.4849315 1016.895 POINT (20.69221 50.81048)
##   elev
## 1  844
## 2 1796
## 3  275
## 4  323
## 5  384
## 6  247
```

---
# Vector-raster interactions

.pull-left[

```r
lc = raster("data/polska_lc.tif")
plot(lc)
```

<img src="index_files/figure-html/unnamed-chunk-42-1.png" style="display: block; margin: auto;" />
]

.pull-right[

```r
wlkp_pn = st_read("data/wlkp_pn.gpkg", quiet = TRUE)
plot(st_geometry(wlkp_pn), axes = TRUE)
```

<img src="index_files/figure-html/unnamed-chunk-43-1.png" style="display: block; margin: auto;" />
]

---
# Vector-raster interactions

```r
wlkp_pn_lc = crop(lc, wlkp_pn)
plot(wlkp_pn_lc)
```

---
# Vector-raster interactions

```r
wlkp_pn_lc2 = mask(wlkp_pn_lc, wlkp_pn)
plot(wlkp_pn_lc2)
```

---
# Exercises

1. Read the 2017 meteorological data and the location of the meteorological stations.
Create a spatial object containing all measurements for three stations: `BIAŁYSTOK`, `KRAKÓW-BALICE`, and `POZNAŃ`. (Tip: use the `right_join()` function.
2. Calculate the average temperature in January and July for these three stations.
3. Determine the land cover of the three stations.
4. Additional: improve the average temperature map presented in the "Reprojections" section.

---
class: inverse, left, bottom
# Example: landscape metrics

---
# Example: landscape metrics

```r
library(landscapemetrics)
library(landscapetools)
```

```r
show_landscape(wlkp_pn_lc, discrete = TRUE)
```

---
# Example: landscape metrics

```r
list_lsm()
```

```
## # A tibble: 132 x 5
##    metric name                         type             level function_name
##    <chr>  <chr>                        <chr>            <chr> <chr>        
##  1 area   patch area                   area and edge m… patch lsm_p_area   
##  2 cai    core area index              core area metric patch lsm_p_cai    
##  3 circle related circumscribing circ… shape metric     patch lsm_p_circle 
##  4 contig contiguity index             shape metric     patch lsm_p_contig 
##  5 core   core area                    core area metric patch lsm_p_core   
##  6 enn    euclidean nearest neighbor … aggregation met… patch lsm_p_enn    
##  7 frac   fractal dimension index      shape metric     patch lsm_p_frac   
##  8 gyrate radius of gyration           area and edge m… patch lsm_p_gyrate 
##  9 ncore  number of core areas         core area metric patch lsm_p_ncore  
## 10 para   perimeter-area ratio         shape metric     patch lsm_p_para   
## # … with 122 more rows
```

---
# Landscape metrics - landscape level

- AI - Aggregation index

```r
lsm_l_ai(wlkp_pn_lc)
```

```
## # A tibble: 1 x 6
##   layer level     class    id metric value
##   <int> <chr>     <int> <int> <chr>  <dbl>
## 1     1 landscape    NA    NA ai      70.3
```

---
# Landscape metrics - class level

- ED - Edge density

```r
lsm_c_ed(wlkp_pn_lc)
```

```
## # A tibble: 12 x 6
##    layer level class    id metric  value
##    <int> <chr> <int> <int> <chr>   <dbl>
##  1     1 class    10    NA ed      8.27 
##  2     1 class    11    NA ed     10.6  
##  3     1 class    30    NA ed      3.12 
##  4     1 class    40    NA ed      0.767
##  5     1 class    60    NA ed      2.35 
##  6     1 class    70    NA ed      9.61 
##  7     1 class    90    NA ed      6.59 
##  8     1 class   100    NA ed      3.65 
##  9     1 class   130    NA ed      1.38 
## 10     1 class   180    NA ed      0.295
## 11     1 class   190    NA ed      5.32 
## 12     1 class   210    NA ed      1.69
```

---
# Landscape metrics - patch level

- PARA - Perimeter-Area ratio

```r
lsm_p_para(wlkp_pn_lc)
```

```
## # A tibble: 234 x 6
##    layer level class    id metric   value
##    <int> <chr> <int> <int> <chr>    <dbl>
##  1     1 patch    10     1 para   0.0169 
##  2     1 patch    10     2 para   0.0169 
##  3     1 patch    10     3 para   0.0117 
##  4     1 patch    10     4 para   0.00620
##  5     1 patch    10     5 para   0.0113 
##  6     1 patch    10     6 para   0.00893
##  7     1 patch    10     7 para   0.0169 
##  8     1 patch    10     8 para   0.0169 
##  9     1 patch    10     9 para   0.0127 
## 10     1 patch    10    10 para   0.0126 
## # … with 224 more rows
```

---
# Landscape metrics - all level

```r
calculate_lsm(wlkp_pn_lc, type = "aggregation metric")
```

```
## # A tibble: 428 x 6
##    layer level class    id metric value
##    <int> <chr> <int> <int> <chr>  <dbl>
##  1     1 class    10    NA ai      43.4
##  2     1 class    11    NA ai      85.0
##  3     1 class    30    NA ai      25.7
##  4     1 class    40    NA ai      40.6
##  5     1 class    60    NA ai      63.6
##  6     1 class    70    NA ai      68.7
##  7     1 class    90    NA ai      70.9
##  8     1 class   100    NA ai      33.5
##  9     1 class   130    NA ai      66.7
## 10     1 class   180    NA ai      37.5
## # … with 418 more rows
```

---
# Landscape metrics - vizualization

```r
mapa_p_para = get_lsm(wlkp_pn_lc, what = "lsm_p_para")
```

```
## Warning: only using 'what' argument
```

```r
show_landscape(mapa_p_para[[1]][[1]])
```

---
# Landscape metrics - vizualization

```r
show_lsm(wlkp_pn_lc, what = "lsm_p_para")
```

---
# Landscape metrics - samples

```r
punkt = st_sf(st_sfc(geom = st_point(c(347500, 493000))))
plot(wlkp_pn_lc)
plot(st_geometry(punkt), add = TRUE, cex = 2)
```

---
# Landscape metrics - samples

```r
sample_lsm(wlkp_pn_lc, punkt, shape = "circle", size = 2000, what = "lsm_c_ed")
```

```
## # A tibble: 10 x 8
##    layer level class    id metric  value plot_id percentage_inside
##    <int> <chr> <int> <int> <chr>   <dbl>   <int>             <dbl>
##  1     1 class    10    NA ed     11.5         1              100.
##  2     1 class    11    NA ed     19.9         1              100.
##  3     1 class    30    NA ed      5.26        1              100.
##  4     1 class    60    NA ed      2.14        1              100.
##  5     1 class    70    NA ed      9.34        1              100.
##  6     1 class    90    NA ed      7.50        1              100.
##  7     1 class   100    NA ed      7.80        1              100.
##  8     1 class   180    NA ed      0.795       1              100.
##  9     1 class   190    NA ed      2.14        1              100.
## 10     1 class   210    NA ed      0.795       1              100.
```

---
# Exercises

1. Borders of the national parks in Poland are located in the `Poland_pn.gpkg` file.
Read this file into R and select `"Białowieski Park Narodowy"` only.
2. Crop the land cover map from the `Polska_lc.tif` file to the extend of `"Białowieski Park Narodowy"`.
3. What land cover categories can be found in this area (you can use the `unique()` function for this purpose)?
4. Calculate the area of `"Białowieski Park Narodowy"` and then the areas the subsequent classes (`TA` and `CA`).
What are the three land cover categories with the largest area?
5. Check how many patches are in this area (`NP`).

---
class: inverse, left, bottom
# Resources

---
# Resources

.pull-left[
- http://geocompr.github.io/ - Lovelace R., Nowosad J., Muenchow J. 2019, **Geocomputation with R**. CRC Press
- **sf** - https://r-spatial.github.io/sf/
- **raster** - https://rspatial.org/spatial/index.html
- **tmap** - https://github.com/mtennekes/tmap
- **dplyr** - http://r4ds.had.co.nz/ - Wickham H., Goremund G. 2016. **R for Data Science**. O'Reilly Media
- **landscapemetrics** - https://r-spatialecology.github.io/landscapemetrics
]

.pull-right[
- Blogs (e.g., https://www.r-bloggers.com/)
- Courses (e.g., https://www.datacamp.com/)
- Websites (e.g., https://gis.stackexchange.com/)
- Twitter `#rspatial`, `#geocompr`
- Meetups
]

---
class: center, middle

.pull-left[
## About me:

Twitter: <i class="ion  ion-social-twitter "></i> jakub_nowosad

Email: nowosad.jakub@gmail.com

https://nowosad.github.io

https://geocompr.github.io/
]