Working with spatial data in R

class: center, middle, inverse, title-slide

# Working with spatial data in R
## <a href="https://github.com/ecodatasci-tlv/spatial" class="uri">https://github.com/ecodatasci-tlv/spatial</a>
### 2019-05-22

---

# Topics we'll cover today:

* Types of spatial data
* Importing/Exporting spatial data
* Projections, extents, and units
* Comparison of different GIS tools out there
* Spatial data manipulation
* Intro to spatial analysis

---
class: inverse, middle, center

#Data types

.small[
https://www.gislounge.com/geodatabases-explored-vector-and-raster-data/
]  
---
.center[
.img-med[
![](pres_pictures/Screenshot 2019-05-22 13.08.45.png)
]]
---
class: middle, center

# **Vector data**

Equivalent to a spreadsheet with a geometry column that describes how to plot the data that represents that feature

---

### **Points**
Used to represent nonadjacent features 
.center[
.img-small[
![](pres_pictures/points-vector.png)
]]
---

### **Lines**
One dimensional and are only used to measure length
.center[
.img-small[
![](pres_pictures/line-vector.png)
]]

---
### **Polygons**
Two dimensional and are used for measurement of area and parameters of geographical feature. 
.center[
.img-small[
![](pres_pictures/polygon-vector.png)
]]

---

# File type

* Polygon and lines come in a shapefile format, point usually are stored in XY table (a spreadsheet with latitude and longitude columns)

* A shapefile is structured from 4 main files: 
 - **shp** - the geometry of each feature
 - **dbf** - a dBase file with the attribute data for each feature (can be opened in excel; stores the metadata of each feature)
 - **prj** - The projection file
 - **shx** - a spatial index that helps a GIS program find the features in the .shp file quickly

**IMPORTANT:** each shapefile can only store one type of geometry - point/ polygon/ line
---
class: middle, center

# **Raster data**

Cell-based data, represents surfaces. Each cell represents a value
---

## **Raster can represent two types of data**

### **Continuous** - temperature, elevation, etc.

### **Discrete** - population density, species richness, etc.

---

# **Rasters have three types of datasets:**

### **Thematic**
Used to represent discrete data, such as soil type or land use
.center[
.img-small[
![](pres_pictures/thematic_map.png)
]]
---
### **Spectral**
Aerial or setalite imagery used during digital map generation
.center[
.img-med[
![](pres_pictures/mx.jpg)
]]
---
### **Surface**
Represents continuous changes across a landscape (for example elevation)
.center[
.img-med[
![](pres_pictures/gmted2010.png)
]]
---

class: center

# Vector and raster comparison
<div class="tabwid"><table style='border-collapse:collapse;'><thead><tr><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Data_type</td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Pros</td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Cons</td></tr></thead><tbody><tr><td rowspan="3" style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Vector</td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">High geographical accuracy</td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Doesn’t work well with continuous data</td></tr><tr><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Follows topology rules which increases integrity</td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Processing intensive due to topology checks</td></tr><tr><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Used for measurements of proximity</td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"></td></tr><tr><td rowspan="3" style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Raster</td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Very good for representing remote sensing data</td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Less aesthetic pleasing</td></tr><tr><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Performs well with map algebra</td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Lacks attribute data flexibility </td></tr><tr><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"></td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;">Data size increases with the increase in resolution</td></tr></tbody></table></div>

---
class: center
# Importing Spatial Data
.center[
.img-med[
![](pres_pictures/dims.jpeg)
]]

---

Now that we've learned about different types of spatial data, let's load some. We'll use the `rgdal` package for this:

```r
pacman::p_load(rgdal)
pkmng_countries <- readOGR(dsn = "Data",
 layer = "pkmng_polygons")
```

```
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/marianovosolov/Dropbox/Studies/ecodatasci_tlv/spatial_analysis/Data", layer: "pkmng_polygons"
## with 246 features
## It has 12 fields
```

---

```r
head(pkmng_countries)
```

```
## class : SpatialPolygonsDataFrame 
## features : 6 
## extent : -61.88722, 50.37499, -18.01639, 42.61805 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables : 12
## # A tibble: 6 x 12
## FIPS ISO2 ISO3 UN NAME AREA POP2005 REGION SUBREGI LON
## <fct> <fct> <fct> <int> <fct> <int> <dbl> <int> <int> <dbl>
## 1 AC AG ATG 28 Anti~ 44 8.30e4 19 29 -61.8 
## 2 AG DZ DZA 12 Alge~ 238174 3.29e7 2 15 2.63
## 3 AJ AZ AZE 31 Azer~ 8260 8.35e6 142 145 47.4 
## 4 AL AL ALB 8 Alba~ 2740 3.15e6 150 39 20.1 
## 5 AM AM ARM 51 Arme~ 2820 3.02e6 142 145 44.6 
## 6 AO AO AGO 24 Ango~ 124670 1.61e7 2 17 17.5 
## # ... with 2 more variables: LAT <dbl>, pkmn_rc <int>
```

---
As usual, there is a `tidyverse` alternative to view spatial objects and edit the associated data
.small[

```r
pacman::p_load(spdplyr)
pkmng_countries
```

```
## class : SpatialPolygonsDataFrame 
## features : 246 
## extent : -180, 180, -90, 83.57027 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables : 12
## # A tibble: 246 x 12
## FIPS ISO2 ISO3 UN NAME AREA POP2005 REGION SUBREGI LON
## <fct> <fct> <fct> <int> <fct> <int> <dbl> <int> <int> <dbl>
## 1 AC AG ATG 28 Anti… 44 8.30e4 19 29 -61.8 
## 2 AG DZ DZA 12 Alge… 238174 3.29e7 2 15 2.63
## 3 AJ AZ AZE 31 Azer… 8260 8.35e6 142 145 47.4 
## 4 AL AL ALB 8 Alba… 2740 3.15e6 150 39 20.1 
## 5 AM AM ARM 51 Arme… 2820 3.02e6 142 145 44.6 
## 6 AO AO AGO 24 Ango… 124670 1.61e7 2 17 17.5 
## 7 AQ AS ASM 16 Amer… 20 6.41e4 9 61 -171. 
## 8 AR AR ARG 32 Arge… 273669 3.87e7 19 5 -65.2 
## 9 AS AU AUS 36 Aust… 768230 2.03e7 9 53 136. 
## 10 BA BH BHR 48 Bahr… 71 7.25e5 142 145 50.6 
## # … with 236 more rows, and 2 more variables: LAT <dbl>,
## # pkmn_rc <int>
```
]
---
With this package, you can use `tidyverse` functions to manipulate your spatial data:
.small[

```r
#filter out all countries without any Pokemon sightings, and let's rename the column while we're at it:
pkmng_countries %>%
  filter(pkmn_rc > 0) %>%
  rename(Country = NAME)
```

```
## class : SpatialPolygonsDataFrame 
## features : 72 
## extent : -179.999, 179.999, -55.9175, 83.10942 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables : 12
## # A tibble: 72 x 12
## FIPS ISO2 ISO3 UN Country AREA POP2005 REGION SUBREGI
## <fct> <fct> <fct> <int> <fct> <int> <dbl> <int> <int>
## 1 AO AO AGO 24 Angola 124670 1.61e7 2 17
## 2 AR AR ARG 32 Argent~ 273669 3.87e7 19 5
## 3 AS AU AUS 36 Austra~ 768230 2.03e7 9 53
## 4 BK BA BIH 70 Bosnia~ 5120 3.92e6 150 39
## 5 BM MM MMR 104 Burma 65755 4.80e7 142 35
## 6 BR BR BRA 76 Brazil 845942 1.87e8 19 5
## 7 BX BN BRN 96 Brunei~ 527 3.74e5 142 35
## 8 CA CA CAN 124 Canada 909351 3.23e7 19 21
## 9 CB KH KHM 116 Cambod~ 17652 1.40e7 142 35
## 10 CI CL CHL 152 Chile 74880 1.63e7 19 5
## # ... with 62 more rows, and 3 more variables: LON <dbl>, LAT <dbl>,
## # pkmn_rc <int>
```
]
---
Now let's load our point data:
.small[

```r
pkmng_points_dat <- read_csv("Data/pkmng_points.csv")
pkmng_points_dat
```

```
## # A tibble: 29,602 x 207
## pokemonId latitude longitude appearedLocalTime cellId_90m
## <dbl> <dbl> <dbl> <dttm> <dbl>
## 1 16 53.7 -0.443 2016-09-04 09:31:44 5.22e18
## 2 60 46.0 8.96 2016-09-06 07:57:06 5.15e18
## 3 19 46.2 15.3 2016-09-04 08:33:55 5.14e18
## 4 16 50.7 14.5 2016-09-04 09:23:24 5.12e18
## 5 23 44.5 11.3 2016-09-04 16:02:16 5.15e18
## 6 120 50.2 14.7 2016-09-03 15:51:17 5.12e18
## 7 46 33.6 -84.3 2016-09-07 01:05:11 9.87e18
## 8 10 46.0 8.95 2016-09-06 06:22:14 5.15e18
## 9 27 33.9 -118. 2016-09-06 05:13:44 9.28e18
## 10 23 34.3 -119. 2016-09-06 14:49:44 9.29e18
## # ... with 29,592 more rows, and 202 more variables:
## # cellId_180m <dbl>, cellId_370m <dbl>, cellId_730m <dbl>,
## # cellId_1460m <dbl>, cellId_2920m <dbl>, cellId_5850m <dbl>,
## # appearedTimeOfDay <chr>, appearedHour <dbl>,
## # appearedMinute <dbl>, appearedDayOfWeek <chr>, appearedDay <dbl>,
## # appearedMonth <dbl>, appearedYear <dbl>, terrainType <dbl>,
## # closeToWater <lgl>, city <chr>, continent <chr>, weather <chr>,
## # temperature <dbl>, windSpeed <dbl>, windBearing <dbl>,
## # pressure <dbl>, weatherIcon <chr>, sunriseMinutesMidnight <dbl>,
## # sunriseHour <dbl>, sunriseMinute <dbl>,
## # sunriseMinutesSince <dbl>, sunsetMinutesMidnight <dbl>,
## # sunsetHour <dbl>, sunsetMinute <dbl>, sunsetMinutesBefore <dbl>,
## # population_density <dbl>, urban <lgl>, suburban <lgl>,
## # midurban <lgl>, rural <lgl>, gymDistanceKm <dbl>,
## # gymIn100m <lgl>, gymIn250m <lgl>, gymIn500m <lgl>,
## # gymIn1000m <lgl>, gymIn2500m <lgl>, gymIn5000m <lgl>,
## # pokestopDistanceKm <dbl>, pokestopIn100m <lgl>,
## # pokestopIn250m <lgl>, pokestopIn500m <lgl>,
## # pokestopIn1000m <lgl>, pokestopIn2500m <lgl>,
## # pokestopIn5000m <lgl>, cooc_1 <lgl>, cooc_2 <lgl>, cooc_3 <lgl>,
## # cooc_4 <lgl>, cooc_5 <lgl>, cooc_6 <lgl>, cooc_7 <lgl>,
## # cooc_8 <lgl>, cooc_9 <lgl>, cooc_10 <lgl>, cooc_11 <lgl>,
## # cooc_12 <lgl>, cooc_13 <lgl>, cooc_14 <lgl>, cooc_15 <lgl>,
## # cooc_16 <lgl>, cooc_17 <lgl>, cooc_18 <lgl>, cooc_19 <lgl>,
## # cooc_20 <lgl>, cooc_21 <lgl>, cooc_22 <lgl>, cooc_23 <lgl>,
## # cooc_24 <lgl>, cooc_25 <lgl>, cooc_26 <lgl>, cooc_27 <lgl>,
## # cooc_28 <lgl>, cooc_29 <lgl>, cooc_30 <lgl>, cooc_31 <lgl>,
## # cooc_32 <lgl>, cooc_33 <lgl>, cooc_34 <lgl>, cooc_35 <lgl>,
## # cooc_36 <lgl>, cooc_37 <lgl>, cooc_38 <lgl>, cooc_39 <lgl>,
## # cooc_40 <lgl>, cooc_41 <lgl>, cooc_42 <lgl>, cooc_43 <lgl>,
## # cooc_44 <lgl>, cooc_45 <lgl>, cooc_46 <lgl>, cooc_47 <lgl>,
## # cooc_48 <lgl>, cooc_49 <lgl>, cooc_50 <lgl>, ...
```
]
---
We are going to create a shapefile from this table.

```r
pacman::p_load(sp)
xy <- tibble(longitude = pkmng_points_dat$longitude,
 latitude = pkmng_points_dat$latitude)
pkmng_points <- SpatialPointsDataFrame(coords = xy,
 data = pkmng_points_dat)
pkmng_points
```

```
## class       : SpatialPointsDataFrame 
## features    : 29602 
## extent      : -158.0198, 175.6162, -42.9845, 68.43698  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## variables   : 207
## names       : pokemonId,   latitude,   longitude,   appearedLocalTime,   cellId_90m,  cellId_180m,  cellId_370m,  cellId_730m, cellId_1460m, cellId_2920m, cellId_5850m, appearedTimeOfDay, appearedHour, appearedMinute, appearedDayOfWeek, ... 
## min values  :         1,   1.282526,    0.000046, 2016-09-02 21:50:11, 1.008879e+19, 1.008879e+19, 1.008879e+19, 1.008879e+19, 1.008879e+19, 1.008879e+19, 1.008879e+19,         afternoon,            0,              0,         dummy_day, ... 
## max values  :       149, -42.984504, -158.019758, 2016-09-08 03:51:27, 9.937698e+18, 9.937698e+18, 9.937698e+18, 9.937698e+18, 9.937698e+18, 9.937698e+18, 9.937698e+18,             night,           23,             59,         Wednesday, ...
```

---
# GeoJSON

So far we've been working with ESRI shapefiles. This is a format to encode geographic data, used widely in GIS systems.

Alternatively, there is the java-based GeoJSON format for encoding geographic data, which can support the features we've discussed (points, polygons).

There are advantages to GeoJSON, including reduced file sizes. Compare:

---

```r
pacman::p_load(geojsonio)
pacman::p_load(rmapshaper)

#we now convert our ESRI shapefile (a SpatialPoylgonsDataFrame object) to GeoJSON format:
pkmng_countries_json <- geojson_json(pkmng_countries) #this can be VERY time-consuming
```

Compare the file sizes: 6.4 Mb for the ESRI format, compared to 1.1 Mb for the GeoJSON format!

Read more about GeoJSON here: https://macwright.org/2015/03/23/geojson-second-bite.html

---
# Rasters
Finally, importing rasters is easily done with the `raster` package:

```r
pacman::p_load(raster)
pkmng_raster <- raster("Data/pkmng_raster.tif")

pkmng_raster
```

```
## class       : RasterLayer 
## dimensions  : 87, 180, 15660  (nrow, ncol, ncell)
## resolution  : 2, 2  (x, y)
## extent      : -180, 180, -90.42973, 83.57027  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/marianovosolov/Dropbox/Studies/ecodatasci_tlv/spatial_analysis/Data/pkmng_raster.tif 
## names       : pkmng_raster 
## values      : 0, 9.535896  (min, max)
```

---
# Projections
One of the most important things to remember when dealing with spatial data is the projection.

A projection is a mathematical transformation of coordinates from the surface of a sphere on a 2D plane.

.center[![](img/m-c.jpg)]

---
# Types of projections
## Projections by surface
Classification based on the type of surface onto which the globe is projected:
* Cylindrical
* Conic
* Azimuthal

etc.

---
## Projections by preservation of a metric property
Classification based on which property is maintained in the projection:
* Conformal (preserves shape, usually distorts area)
* Equal-area (preserves area, usually distorts shape)
* Equidistant (preserves distance from some point or line)

etc.

---
You can view the Coordinate Reference System of your spatial objects:

```r
crs(pkmng_countries)
```

```
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
```

You can see different values here, including the datum, units, etc.

---
Notice our points shapefile has no projection, because we didn't set it when creating the object:

```r
crs(pkmng_points)
```

```
## CRS arguments: NA
```

---
We can set the projection for this object, either by writing a string of our desired projection, or extracting the projection from an existing object:

```r
proj1 <- crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
proj1
```

```
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
```

```r
proj2 <- crs(pkmng_countries)
proj2
```

```
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
```

---

```r
proj4string(pkmng_points) <- proj1
pkmng_points
```

```
## class       : SpatialPointsDataFrame 
## features    : 29602 
## extent      : -158.0198, 175.6162, -42.9845, 68.43698  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 207
## names       : pokemonId,   latitude,   longitude, appearedLocalTime,           cellId_90m,          cellId_180m,          cellId_370m,         cellId_730m,         cellId_1460m,         cellId_2920m,         cellId_5850m, appearedTimeOfDay, appearedHour, appearedMinute, appearedDayOfWeek, ... 
## min values  :         1, -42.984504, -158.019758,        1472853011,    43172735374852096,    43172735576178688,    43172735307743232,   43172732086517760,    43172719201615872,    43172667662008320,    43172598942531584,         afternoon,            0,              0,         dummy_day, ... 
## max values  :       149,   68.43698,  175.616223,        1473306687, 12280882455157669888, 12280882455224778752, 12280882455493214208, 1.2280882456567e+19, 12280882443682054144, 12280882495221661696, 12280882701380091904,             night,           23,             59,         Wednesday, ...
```

---
You can also transform already existing objects to have different projections:

```r
pkmng_points_p <- spTransform(pkmng_points,
 crs("+proj=laea +lat_0=0 +lon_0=0"))
pkmng_points_p
```

```
## class       : SpatialPointsDataFrame 
## features    : 29602 
## extent      : -9285067, 11178209, -11993372, 9520128  (xmin, xmax, ymin, ymax)
## crs         : +proj=laea +lat_0=0 +lon_0=0 +ellps=WGS84 
## variables   : 207
## names       : pokemonId,   latitude,   longitude, appearedLocalTime,           cellId_90m,          cellId_180m,          cellId_370m,         cellId_730m,         cellId_1460m,         cellId_2920m,         cellId_5850m, appearedTimeOfDay, appearedHour, appearedMinute, appearedDayOfWeek, ... 
## min values  :         1, -42.984504, -158.019758,        1472853011,    43172735374852096,    43172735576178688,    43172735307743232,   43172732086517760,    43172719201615872,    43172667662008320,    43172598942531584,         afternoon,            0,              0,         dummy_day, ... 
## max values  :       149,   68.43698,  175.616223,        1473306687, 12280882455157669888, 12280882455224778752, 12280882455493214208, 1.2280882456567e+19, 12280882443682054144, 12280882495221661696, 12280882701380091904,             night,           23,             59,         Wednesday, ...
```

---
Why are projections so important?

Here's our points mapped on top of a shapefile of countries:

```r
plot(countries, axes = T)
plot(pkmng_points, add = T, pch = 1, col = "red")
```

---
Now, when the projections do not match:

```r
countries_p <- spTransform(countries,
 crs("+proj=laea +lat_0=0 +lon_0=0"))
plot(countries_p, axes = T)
plot(pkmng_points, add = T, pch = 1, col = "red")
```

---
Whereas matching projections give us all the coordinates in the right place:

```r
plot(countries_p, axes = T)
plot(pkmng_points_p, add = T, pch = 1, col = "red")
```

---
# Extent

Spatial objects also have extents - these are the minimum and maximum X and Y coordinates of the object. They don't have to have data in them. The coordinates of course depend on your projection and units:

```r
countries@bbox
```

```
##    min       max
## x -180 180.00000
## y  -90  83.57027
```

```r
countries_p@bbox
```

```
##         min      max
## x -12726751 12733710
## y -12641498 12299760
```

---
For raster objects, you find the extent with this syntax:

```r
pkmng_raster@extent
```

```
## class      : Extent 
## xmin       : -180 
## xmax       : 180 
## ymin       : -86.42973 
## ymax       : 83.57027
```

---
You can change the extent of objects using the `crop()` or `extend()` functions, depending on whether you want to decrease or increase the extent:

```r
pkmng_raster_c <- crop(pkmng_raster,
 extent(-50,
 100,
 0,
 80))
```

---

```r
plot(pkmng_raster)
plot(pkmng_raster_c, add = T, col = "red")
```

All the red cells are in the new, limited extent!

---
## Why is this important?

If you try to do spatial arithmetics on objects with mismatching extents and projections, you'll get errors. For example:

```r
raster::intersect(countries_p, pkmng_raster_c)
```

```
## Warning in intersect(x, y): non identical CRS
```

```
## Warning in intersect(x, y): polygons do not intersect
```

```
## NULL
```

---
class: center, middle
# Data manipulation

---
class: center, middle
# Geoprocessing
a GIS operation used to manipulate spatial data. A typical geoprocessing operation takes an input dataset, performs an operation on that dataset, and returns the result of the operation as an output dataset

---
### **Buffer points** - creating a polygon based on proximity
**For example:** buffer some of the pokemon points to 100km radius
We'll first subset the data to have less points

```r
sub_pkmng_p<- subset(pkmng_points,pokemonId %in%c(8:10))
plot(countries)
plot(sub_pkmng_p,add = T, col = "red")
```

<img src="index_files/figure-html/unnamed-chunk-2-1.png" style="display: block; margin: auto;" />
---
Now we'll buffer the points in a 100km radius to make them into polygons

```r
buffer_pkmng_points<- raster::buffer(sub_pkmng_p,width = 100000)
plot(countries)
plot(buffer_pkmng_points,add = T, col = "red")
```

---
### **Clip** - an overlay function that cuts out an input layer with the extent of a defined feature boundary. The result of this tool is a new clipped output layer.

There are a few ways of clipping polygons in R, we will look at using the rgdal package. 
---
Lets clip the buffered points to to have the ones that fall in Italy and to fit the boundaries of the Italian coastline.

First, well filter the countries shapefile to have a shapefile of Italy

```r
Italy <- countries %>% 
 filter(NAME == "Italy")
plot(Italy)
plot(buffer_pkmng_points,add = T, col = "red")
```

<img src="index_files/figure-html/unnamed-chunk-4-1.png" style="display: block; margin: auto;" />
---
# `rgeos` package example
Now, lets clip the the buffered points based on the Italian shapefile we created

```r
clip_pkmng_rgeos<- rgeos::gIntersection(Italy, buffer_pkmng_points, byid = TRUE, drop_lower_td = TRUE)
plot(Italy)
plot(clip_pkmng_rgeos,add = T, col = "red")
```

---

### **Intersect** - similar to the clip tool because the extents of input features defines the output. The only exception is that it preserves attributes from all the data sets that overlap each other in the output

Lets intersect the points with the country shapefile to find which points belong to which country

```r
intersect_pkmg<- raster::intersect(pkmng_points,countries) 
```

---
### **Merge** - combines data sets that are the same data type

```r
#subset the countries
sub_count<- countries %>% filter(NAME==c("France","Italy"))
sub_count_sf<- sf::st_as_sf(sub_count)
merge_count<- sf::st_combine(sub_count_sf)
plot(merge_count)
```

---
### **Dissolve** - The Dissolve Tool unifies boundaries based on common attribute values

```r
diss_sub_count<- rgeos::gUnaryUnion(sub_count)
plot(diss_sub_count)
```

---
# Zonal statistics
Vector-Raster analysis. Allows you to calculate statistics from raster data for each feature in the vector data

**For Example:** We want to calculate the temperature in each point that a pokemon was collected

---
Download and plot the mean annual temperature from the web using the getData function

```r
climate <- raster::getData('worldclim', var='bio', res=10)
annual_prec<- climate$bio12
plot(annual_prec)
```

---
Now lets calculate the mean temperature in each point

```r
poke_annual_prec<- raster::extract(annual_prec, pkmng_points,fun = mean,na.rm=T,df=T)
head(poke_annual_prec)
```

```
##   ID bio12
## 1  1   663
## 2  2  1326
## 3  3  1095
## 4  4   585
## 5  5   792
## 6  6   492
```

---
The data has no identification of which point it belongs to. To have that we need to merge it with the points identifiers

```r
poke_annual_prec_final<- cbind(pkmng_points@data$pokemonId,poke_annual_prec$bio12)
names(poke_annual_prec_final)<- c("pokemonId","annual_prec")
head(poke_annual_prec_final)
```

```
##      [,1] [,2]
## [1,]   16  663
## [2,]   60 1326
## [3,]   19 1095
## [4,]   16  585
## [5,]   23  792
## [6,]  120  492
```

---
# Spatial autocorrelation
Spatial data are not independent - violating assumptions of statistical tests.

Spatial autocorrelation decreases with distance - but the distance depends on your data.

So how do we test for spatial autocorrelation and what can we do about it?

---
## Testing for spatial autocorrelation

We'll use the `pgirmess` package:

```r
pacman::p_load(pgirmess)

#first we'll take a small sample of the points data, because this test is very time-consuming
pkmng_points_dat_s <- sample_frac(pkmng_points_dat, 0.05)
pgi.cor <- correlog(coords = pkmng_points_dat_s[,c("longitude","latitude")], #longitude and latitude, in that order
 z = pkmng_points_dat_s$temperature, #the variable we want to test for autocorrelation in
 method = "Moran")
```

---
The X axis is distance, and the Y axis is our measure of autocorrelation. Values higher than 0 mean there is spatial autocorrelation in that distance class.

```r
plot(pgi.cor)
```

---
And this is what it looks like when there's no spatial autocorrelation:

```r
pgi.cor2 <- correlog(coords = pkmng_points_dat_s[,c("longitude","latitude")], #longitude and latitude, in that order
 z = pkmng_points_dat_s$pokemonId, #the variable we want to test for autocorrelation in
 method = "Moran")
plot(pgi.cor2)
```

---
## Dealing with spatial autocorrelation
There are various methods to "eliminate" spatial autocorrelation for hypothesis testing and statistical inference. All are flawed, each in their own way.

We won't talk about them here, but you can read more in the literature:
https://onlinelibrary.wiley.com/doi/10.1111/j.2007.0906-7590.05171.x
https://onlinelibrary.wiley.com/doi/10.1111/jbi.12953

---
# Spatial non-stationarity
Often, relationships between variables can shift in space. So global coefficients of models can be very poor at describing the actual behavior of your data.

For example, let's look for a relationship between wind speed and temperature in our Pokemon Go sightings data:
.small[

```r
lm.model <- lm(windSpeed ~
 temperature,
 data = pkmng_points_dat_s)
summary(lm.model)
```

```
## 
## Call:
## lm(formula = windSpeed ~ temperature, data = pkmng_points_dat_s)
## 
## Residuals:
## Min 1Q Median 3Q Max 
## -7.071 -3.112 -0.737 2.241 49.873 
## 
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) 
## (Intercept) 5.06517 0.40037 12.651 < 2e-16 ***
## temperature 0.07405 0.01766 4.194 2.91e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.249 on 1478 degrees of freedom
## Multiple R-squared: 0.01176,	Adjusted R-squared: 0.01109 
## F-statistic: 17.59 on 1 and 1478 DF, p-value: 2.906e-05
```
]
---
Now, we'll map the residuals from this model:

```r
resids<-residuals(lm.model)
colours <- c("dark blue", "blue", "red", "dark red") 
map.resids <- SpatialPointsDataFrame(data = data.frame(resids),
 coords = cbind(pkmng_points_dat_s$longitude,
 pkmng_points_dat_s$latitude)
 ) 
```

---

```r
spplot(map.resids, cuts=quantile(resids), col.regions=colours, cex=1)
```

<img src="index_files/figure-html/gwr 3-1.png" style="display: block; margin: auto;" />
Notice how the red and blue dots are clustered together, and not randomly distributed - showing a spatial pattern in the residuals.

---
We can use a method called GWR (Globally Weighted Regression) to account for this - basically, you calculate a regression coefficient for each datum, weighing the contribution of all data based on their spatial distance from your focal datum.

.small[

```r
pacman::p_load(spgwr)
GWRbandwidth <- gwr.sel(windSpeed ~
 temperature,
 data = pkmng_points_dat_s,
 coords = cbind(pkmng_points_dat_s$longitude,
 pkmng_points_dat_s$latitude),
 adapt = T)
```

```
## Adaptive q: 0.381966 CV score: 25971.6 
## Adaptive q: 0.618034 CV score: 26470.16 
## Adaptive q: 0.236068 CV score: 25779.09 
## Adaptive q: 0.145898 CV score: 25429.72 
## Adaptive q: 0.09016994 CV score: 24672.23 
## Adaptive q: 0.05572809 CV score: 23923.28 
## Adaptive q: 0.03444185 CV score: 22690.01 
## Adaptive q: 0.02128624 CV score: 21834.46 
## Adaptive q: 0.01315562 CV score: 20909.26 
## Adaptive q: 0.008130619 CV score: 20374.53 
## Adaptive q: 0.005024999 CV score: 21990.94 
## Adaptive q: 0.01005 CV score: 20553.34 
## Adaptive q: 0.006944377 CV score: 21538.02 
## Adaptive q: 0.008863756 CV score: 20508.88 
## Adaptive q: 0.007677515 CV score: 20717.15 
## Adaptive q: 0.008410652 CV score: 20342.56 
## Adaptive q: 0.008357601 CV score: 20328.39 
## Adaptive q: 0.008304636 CV score: 20328.57 
## Adaptive q: 0.008357601 CV score: 20328.39
```
]
---

```r
gwr.model <- gwr(windSpeed ~
 temperature,
 data = pkmng_points_dat_s,
 coords = cbind(pkmng_points_dat_s$longitude,
 pkmng_points_dat_s$latitude),
 adapt = GWRbandwidth,
 hatmatrix = T,
 se.fit = T
 )
gwr.results <- as_tibble(gwr.model$SDF)
```

---
Now we can view the table of the results, and see the coefficient for each datum in our dataset - i.e. for each point locality. We can see they vary quite a bit.
.small[

```r
gwr.results
```

```
## # A tibble: 1,480 x 14
## sum.w X.Intercept. temperature X.Intercept._se temperature_se
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 30.7 31.2 -1.02 3.18 0.131 
## 2 13.8 6.22 -0.105 6.27 0.257 
## 3 14.4 11.4 -0.0643 7.80 0.390 
## 4 13.2 2.73 0.317 3.33 0.246 
## 5 63.8 5.35 0.0547 1.15 0.0454
## 6 45.0 6.98 0.0553 1.41 0.0763
## 7 20.7 10.4 -0.201 2.64 0.135 
## 8 13.9 1.26 0.365 3.11 0.182 
## 9 14.4 11.4 -0.0654 7.78 0.389 
## 10 13.1 -9.07 0.666 3.64 0.158 
## # ... with 1,470 more rows, and 9 more variables: gwr.e <dbl>,
## # pred <dbl>, pred.se <dbl>, localR2 <dbl>,
## # X.Intercept._se_EDF <dbl>, temperature_se_EDF <dbl>,
## # pred.se.1 <dbl>, coord.x <dbl>, coord.y <dbl>
```
]
---
Now let's plot this and see exactly what this non-stationarity looks like:
.small[

```r
pkmng_points_dat_s$coef <- gwr.results$temperature
countries_outline <- fortify(countries,
 region = "NAME")

mapWorld <- borders("world", colour=NA, fill="gray")
```

```
## Warning: package 'maps' was built under R version 3.5.2
```

```r
gwr.map <- ggplot(data = pkmng_points_dat_s,
 aes(x = longitude,
 y = latitude,
 colour = coef)
 ) + 
 mapWorld +
 geom_point() + 
 scale_colour_gradient2(low = "red",
 mid = "white",
 high = "blue",
 midpoint = 0,
 space = "rgb",
 na.value = "grey50",
 guide = "colourbar",
 guide_legend(title="Coefs")) +
 coord_equal()
```

```
## Warning: Non Lab interpolation is deprecated
```
]
---
And see how the regression coefficients change with space - blue for positive coefficients and red for negative ones:

```r
gwr.map
```

<img src="index_files/figure-html/gwr 8-1.png" style="display: block; margin: auto;" />
---
class: center, middle
# Different GIS tools out there

.center[
.img-med[
![](pres_pictures/Screenshot 2019-05-21 20.53.18.png)
]]
---
# TAKE HOME MESSAGES

* Be mindful of projections!
* Be mindful of spatial autocorrelation!
* Visualize spatial patterns!
* No one tool that fits them all