Spatial Visualization

class: center, middle, inverse, title-slide

.title[
# Spatial Visualization
]
.author[
### Thiyanga S. Talagala
]

---

### R packages for spatial data analysis

.pull-left[

- rgdal

- sp

- rgeos

- raster

- sf

- tmap

- leaflet

- ggmap

]

.pull-right[

- maptools

- gstat

- spatstat

- stars

- geosphere

- RgoogleMaps

- rasterVis

]

---

### Geospatial vector data structures

**Point:** Individual longitude and latitude of locations

Eg: Locations, Buildings

**Lines:** Two or more vertices or points that are connected

Eg: Roads, rivers

**Polygons:** Three or more vertices and closed

Eg: Area of a country, state, district

---

.pull-left[

**What do we want to plot?** Globe - 3D

![](https://media.tenor.com/4A4UHbVtBIIAAAAC/globe-world.gif)

]
--

.pull-right[

**Where do we plot?** 
Computer screen or paper - 2D

![](sp_viz_files/figure-html/unnamed-chunk-1-1.png)

]

---
class: center, middle, inverse

# Challenge

## How to transform this three dimensional angular system to a two dimensional cartesian system?

--
> Solution : Spatial projection

---

## Types of map projections

- Cylindrical projections

- Conical projections

- Planar projections

> Visit [here](https://johnedevans.files.wordpress.com/2019/01/map-projections-explained.jpg) to see the visual illustration.

---
## **C**oordinate **R**eference **S**ystem (**CRS**)

Different types of CRSs:

- Mercator

- UTM

- Robinson

- Lambert

- Sinusoidal

- Alber

---

![](sp_viz_files/figure-html/unnamed-chunk-4-1.png)

---
![](sp_viz_files/figure-html/unnamed-chunk-5-1.png)
---

![](sp_viz_files/figure-html/unnamed-chunk-6-1.png)

---

![](sp_viz_files/figure-html/unnamed-chunk-7-1.png)

---

### Why some areas are stretched and some areas are compressed?

Projecting data from a 3D round surface onto a 2D flat surface results in visual modifications to the data.

--
### Why multiple CRS?

To optimize to best represent the:

shape

distance

area

direction

---

## Shapefiles

- Spatial data in vector format are often stored in a shapefile format

- The structure of points, lines, and polygons are different, each individual shapefile can only contain one vector type (all points, all lines or all polygons).

- You will not find a mixture of point, line and polygon objects in a single shapefile.

- Shapefiles contain geographic vector data that is expressed in a particular CRS

---

- A shapefile is made up of three required files with the same prefix ('spatial-data' in this case) but different extensions:

**spatial-data.shp:** main file that stores records of each shape geometries

**spatial-data.shx:** index of how the geometries in the main file relate to one-another

**spatial-data.dbf:** attributes of each record
---
## CRS

- geodetic datum (e.g, WGS84)

- the type of map projection (e.g., Mercator) 
- the parameters of the projection (e.g., location of the origin)

---

## Two types of maps

1. Simple feature map

2. Polygon maps

---
class: center, middle

# Simple feature map

---
## Example

```r
# install.packages("devtools")
devtools::install_github("thiyangt/ceylon")
library(ceylon)
library(tidyverse)
library(sp)
library(viridis)
```

```r
data(sf_sl_0)
class(sf_sl_0)
```

```
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"
```

```r
# work with spatial data; sp package will load with rgdal.
library(rgdal)
```

```
## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
```

```r
# for metadata/attributes- vectors or rasters
library(raster)
```

```
## 
## Attaching package: 'raster'
```

```
## The following object is masked from 'package:dplyr':
## 
##     select
```

```r
crs(sf_sl_0)
```

```
## [1] "PROJCRS[\"SLD99 / Sri Lanka Grid 1999\",\n    BASEGEOGCRS[\"SLD99\",\n        DATUM[\"Sri Lanka Datum 1999\",\n            ELLIPSOID[\"Everest 1830 (1937 Adjustment)\",6377276.345,300.8017,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",5233]],\n    CONVERSION[\"Sri Lanka Grid 1999\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",7.00047152777778,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",80.7717130833333,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9999238418,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"Sri Lanka - onshore\"],\n        BBOX[5.86,79.64,9.88,81.95]],\n    ID[\"EPSG\",5235]]"
```

---

```r
sf_sl_0
```

```
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 362203.3 ymin: 380301.9 xmax: 621918.1 ymax: 813560.9
## Projected CRS: SLD99 / Sri Lanka Grid 1999
## # A tibble: 1 × 2
##                                                                 geometry COUNTRY
##                                                       <MULTIPOLYGON [m]> <chr>  
## 1 (((481925.5 381353.7, 481922.9 381350.3, 481919 381348.2, 481914.6 38… SRI LA…
```

---

## 2. Visualizing using shapefiles: map of Sri Lanka

```r
ggplot(sf_sl_0) +
  geom_sf()
```

![](sp_viz_files/figure-html/unnamed-chunk-11-1.png)

---

```r
ggplot(sf_sl_0) +
  geom_sf(fill='beige') +
theme_minimal()
```

![](sp_viz_files/figure-html/unnamed-chunk-12-1.png)

---

```r
library(knitr)
sf_sl_0 %>% kable()
```

|geometry                       |COUNTRY   |
|:------------------------------|:---------|
|MULTIPOLYGON (((481925.5 38... |SRI LANKA |

---

## Your turn: Draw district and province maps

data:

```r
province
district
```

---

.pull-left[

Province

```r
ggplot(data=province) +
  geom_sf()
```

![](sp_viz_files/figure-html/unnamed-chunk-14-1.png)

]

.pull-right[

District

```r
ggplot(data=district) +
  geom_sf()
```

![](sp_viz_files/figure-html/unnamed-chunk-15-1.png)

]

---

```r
province
```

```
## Simple feature collection with 9 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 362203.3 ymin: 380301.9 xmax: 621918.1 ymax: 813560.9
## Projected CRS: SLD99 / Sri Lanka Grid 1999
## # A tibble: 9 × 4
##                                              geometry PROVINCE Status population
## *                                  <MULTIPOLYGON [m]> <chr>    <chr>       <dbl>
## 1 (((498211.2 611042.6, 498401.7 610897.1, 498415.9 … CENTRAL  Provi…    2781000
## 2 (((609877 559315.9, 609857.5 559313, 609845.2 5593… EASTERN  Provi…    1746000
## 3 (((501928.5 712099.6, 501970.4 712098.8, 502020 71… NORTH C… Provi…    1386000
## 4 (((393087.5 629959, 393098.1 629956.7, 393103.2 62… NORTH W… Provi…    2563000
## 5 (((405858.7 700083.7, 405855.3 700082, 405851.5 70… NORTHERN Provi…    1152000
## 6 (((432117.3 521875.5, 432083.5 521852.9, 432041.7 … SABARAG… Provi…    2070000
## 7 (((481925.5 381353.7, 481922.9 381350.3, 481919 38… SOUTHERN Provi…    2669000
## 8 (((522825 568220.9, 522872.6 568217.7, 522926.1 56… UVA      Provi…    1387000
## 9 (((411888.2 438189.4, 411886.4 438182.3, 411881.7 … WESTERN  Provi…    6165000
```

---

```r
district
```

```
## Simple feature collection with 26 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 362203.3 ymin: 380301.9 xmax: 621918.1 ymax: 813560.9
## Projected CRS: SLD99 / Sri Lanka Grid 1999
## # A tibble: 26 × 4
##                                              geometry DISTRICT Status population
##  *                                 <MULTIPOLYGON [m]> <chr>    <chr>       <dbl>
##  1 (((392952.6 594064.2, 392985.8 593971.5, 392986.2… [UNKNOW… <NA>           NA
##  2 (((551311.8 581603.6, 551324.5 581599.7, 551334.5… AMPARA   Distr…     736000
##  3 (((501928.5 712099.6, 501970.4 712098.8, 502020 7… ANURADH… Distr…     943000
##  4 (((522825 568220.9, 522872.6 568217.7, 522926.1 5… BADULLA  Distr…     886000
##  5 (((609877 559315.9, 609857.5 559313, 609845.2 559… BATTICA… Distr…     579000
##  6 (((400911.5 497888.4, 401032.5 497869.3, 401035.2… COLOMBO  Distr…    2455000
##  7 (((425038.8 404072.5, 425030.7 404072.5, 425022 4… GALLE    Distr…    1135000
##  8 (((430775.7 534473.9, 430789.2 534441.8, 430790.2… GAMPAHA  Distr…    2423000
##  9 (((592923.3 452679, 592949.2 452679, 592989.8 452… HAMBANT… Distr…     668000
## 10 (((362968.6 763698.5, 362883.6 763686.9, 362787.9… JAFFNA   Distr…     621000
## # ℹ 16 more rows
```

---

## Guess the output.

```r
ggplot() +
  geom_sf(data=province, lwd=2, col="black") + 
  geom_sf(data=district, linetype=21, col="red")
```

---

## Your turn:

Colour provinces according to the population.

---

## Annotations

.pull-left[

```r
ggplot(sf_sl_0) +
  geom_sf(fill='beige') +
ggspatial::annotation_north_arrow(location = "br")+
  ggspatial::annotation_scale(location = "bl")
```

]

.pull-right[

![](sp_viz_files/figure-html/unnamed-chunk-20-1.png)

]

---
## More work

```r
library(spData)
```

```
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
```

```
## 
## Attaching package: 'spData'
```

```
## The following object is masked _by_ '.GlobalEnv':
## 
##     world
```

---

```r
data("us_states") # from spData package
us_states
```

```
## Simple feature collection with 49 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.7042 ymin: 24.55868 xmax: -66.9824 ymax: 49.38436
## Geodetic CRS:  NAD83
## First 10 features:
##    GEOID        NAME   REGION             AREA total_pop_10 total_pop_15
## 1     01     Alabama    South 133709.27 [km^2]      4712651      4830620
## 2     04     Arizona     West 295281.25 [km^2]      6246816      6641928
## 3     08    Colorado     West 269573.06 [km^2]      4887061      5278906
## 4     09 Connecticut Norteast  12976.59 [km^2]      3545837      3593222
## 5     12     Florida    South 151052.01 [km^2]     18511620     19645772
## 6     13     Georgia    South 152725.21 [km^2]      9468815     10006693
## 7     16       Idaho     West 216512.66 [km^2]      1526797      1616547
## 8     18     Indiana  Midwest  93648.40 [km^2]      6417398      6568645
## 9     20      Kansas  Midwest 213037.08 [km^2]      2809329      2892987
## 10    22   Louisiana    South 122345.76 [km^2]      4429940      4625253
##                          geometry
## 1  MULTIPOLYGON (((-88.20006 3...
## 2  MULTIPOLYGON (((-114.7196 3...
## 3  MULTIPOLYGON (((-109.0501 4...
## 4  MULTIPOLYGON (((-73.48731 4...
## 5  MULTIPOLYGON (((-81.81169 2...
## 6  MULTIPOLYGON (((-85.60516 3...
## 7  MULTIPOLYGON (((-116.916 45...
## 8  MULTIPOLYGON (((-87.52404 4...
## 9  MULTIPOLYGON (((-102.0517 4...
## 10 MULTIPOLYGON (((-92.01783 2...
```

---

```r
data("us_states") # from spData package
us_states
```

---

## Plotting roads and other information

Click [here](https://www.emilyburchfield.org/courses/eds/making_maps_in_r)

---

## Map of Australia

```r
library(ozmaps)
library(sf)
oz_states <- ozmaps::ozmap_states
oz_states
```

```
## Simple feature collection with 9 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 105.5507 ymin: -43.63203 xmax: 167.9969 ymax: -9.229287
## Geodetic CRS:  GDA94
## # A tibble: 9 × 2
##   NAME                                                                  geometry
## * <chr>                                                       <MULTIPOLYGON [°]>
## 1 New South Wales              (((150.7016 -35.12286, 150.6611 -35.11782, 150.6…
## 2 Victoria                     (((146.6196 -38.70196, 146.6721 -38.70259, 146.6…
## 3 Queensland                   (((148.8473 -20.3457, 148.8722 -20.37575, 148.85…
## 4 South Australia              (((137.3481 -34.48242, 137.3749 -34.46885, 137.3…
## 5 Western Australia            (((126.3868 -14.01168, 126.3625 -13.98264, 126.3…
## 6 Tasmania                     (((147.8397 -40.29844, 147.8902 -40.30258, 147.8…
## 7 Northern Territory           (((136.3669 -13.84237, 136.3339 -13.83922, 136.3…
## 8 Australian Capital Territory (((149.2317 -35.222, 149.2346 -35.24047, 149.271…
## 9 Other Territories            (((167.9333 -29.05421, 167.9188 -29.0344, 167.93…
```

---
## Obtain the following map

![](sp_viz_files/figure-html/unnamed-chunk-25-1.png)

---

```r
ggplot(oz_states) + 
  geom_sf() + 
  coord_sf()
```

![](sp_viz_files/figure-html/unnamed-chunk-26-1.png)

---

## `coord_sf()`

- ensures that every layer in the plot uses the same projection

---

## Mark following locations on the map

.pull-left[

```r
oz_capitals <- tibble::tribble( 
  ~city,           ~lat,     ~lon,
  "Sydney",    -33.8688, 151.2093,  
  "Melbourne", -37.8136, 144.9631, 
  "Brisbane",  -27.4698, 153.0251, 
  "Adelaide",  -34.9285, 138.6007, 
  "Perth",     -31.9505, 115.8605, 
  "Hobart",    -42.8821, 147.3272, 
  "Canberra",  -35.2809, 149.1300, 
  "Darwin",    -12.4634, 130.8456, 
)
```

]

.pull-right[

```r
oz_capitals
```

```
## # A tibble: 8 × 3
##   city        lat   lon
##   <chr>     <dbl> <dbl>
## 1 Sydney    -33.9  151.
## 2 Melbourne -37.8  145.
## 3 Brisbane  -27.5  153.
## 4 Adelaide  -34.9  139.
## 5 Perth     -32.0  116.
## 6 Hobart    -42.9  147.
## 7 Canberra  -35.3  149.
## 8 Darwin    -12.5  131.
```

]

---
## Add other geoms

.pull-left[

```r
ggplot() + 
  geom_sf(data = oz_states, colour = "black", fill = NA) + 
  geom_point(data = oz_capitals, mapping = aes(x = lon, y = lat), colour = "red") 
```
]

.pull-right[

![](sp_viz_files/figure-html/unnamed-chunk-30-1.png)

]

---

## Add other geoms

.pull-left[

```r
ggplot() + 
  geom_sf(data = oz_states, colour = "black", fill = NA) + 
  geom_point(data = oz_capitals, mapping = aes(x = lon, y = lat), colour = "red") + coord_sf()
```
]

.pull-right[

![](sp_viz_files/figure-html/unnamed-chunk-32-1.png)

]

---
Your turn: Draw the electoral boundaries

```r
oz_votes <- rmapshaper::ms_simplify(ozmaps::abs_ced)
```

---

```r
library(ggplot2)
library(sf)
nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
point <- data.frame(id = "hi", lat = 36, long = -80) # point inside NC
nc
```

```
## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27
## First 10 features:
##     AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
## 1  0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1
## 2  0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0
## 3  0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5
## 4  0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1
## 5  0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9
## 6  0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7
## 7  0.062     1.547  1834    1834      Camden 37029  37029       15   286     0
## 8  0.091     1.284  1835    1835       Gates 37073  37073       37   420     0
## 9  0.118     1.421  1836    1836      Warren 37185  37185       93   968     4
## 10 0.124     1.428  1837    1837      Stokes 37169  37169       85  1612     1
##    NWBIR74 BIR79 SID79 NWBIR79                       geometry
## 1       10  1364     0      19 MULTIPOLYGON (((-81.47276 3...
## 2       10   542     3      12 MULTIPOLYGON (((-81.23989 3...
## 3      208  3616     6     260 MULTIPOLYGON (((-80.45634 3...
## 4      123   830     2     145 MULTIPOLYGON (((-76.00897 3...
## 5     1066  1606     3    1197 MULTIPOLYGON (((-77.21767 3...
## 6      954  1838     5    1237 MULTIPOLYGON (((-76.74506 3...
## 7      115   350     2     139 MULTIPOLYGON (((-76.00897 3...
## 8      254   594     2     371 MULTIPOLYGON (((-76.56251 3...
## 9      748  1190     2     844 MULTIPOLYGON (((-78.30876 3...
## 10     160  2038     5     176 MULTIPOLYGON (((-80.02567 3...
```

---
.pull-left[

## Data

```r
ggplot() +
  geom_sf(data = nc) +
  coord_sf(crs = st_crs(3347))
```

![](sp_viz_files/figure-html/unnamed-chunk-35-1.png)

]

.pull-right[

## Points

```r
ggplot() +
  geom_sf(data = nc) +
  geom_point(data = point, aes(x = long, y = lat))
```

![](sp_viz_files/figure-html/unnamed-chunk-36-1.png)

]

---
## Data and Points

```r
ggplot() +
  geom_sf(data = nc) +
  geom_point(data = point, aes(x = long, y = lat)) +
  coord_sf(crs = st_crs(3347))
```

![](sp_viz_files/figure-html/unnamed-chunk-37-1.png)

---

EPSG:3347

NAD83 / Statistics Canada Lambert

**Coordinate Reference System Formats: Three common formats include.**

- proj.4

- EPSG

- Well-known Text (WKT) formats.
---

```r
point <- st_as_sf(point, coords = c("long", "lat"), crs = 4326)

ggplot() +
  geom_sf(data = nc) +
  geom_sf(data = point) +
  coord_sf(crs = st_crs(3347))
```

![](sp_viz_files/figure-html/unnamed-chunk-38-1.png)

`coord_sf` knows what to transform from

---
class: center, middle

`coord_sf` can convert layers to common projections but they need to be sf objects for it to know how to do so.

---

## Projecting your data

```r
library(sf)
library(ggplot2)
# devtools::install_github("hrbrmstr/albersusa")
library(albersusa)
```

```
## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
##  Checking rgeos availability: TRUE
```

```
## rgeos version: 0.6-4, (SVN revision 699)
##  GEOS runtime version: 3.10.2-CAPI-1.16.0 
##  Please note that rgeos will be retired during October 2023,
## plan transition to sf or terra functions using GEOS at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html for details.
##  GEOS using OverlayNG
##  Linking to sp version: 2.0-0 
##  Polygon checking: TRUE
```

---

```r
crs_use = "+proj=laea +lat_0=30 +lon_0=-95"

d_points = data.frame(long = c(-110, -103, -84), 
                      lat  = c(45, 40, 41))
```

---

```r
A <- ggplot(data = usa_sf()) +
  geom_sf() + 
  geom_point(data = d_points, 
             aes(x = long, y = lat), 
             color = "red", size = 5) + 
  theme_minimal() + 
  ggtitle("(A) right point position, wrong projection")
```

---

```r
A
```

```
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
```

![](sp_viz_files/figure-html/unnamed-chunk-42-1.png)

---

```r
B <- ggplot(data = usa_sf()) +
  geom_sf() + 
  geom_point(data = d_points, 
             aes(x = long, y = lat), 
             color = "red", size = 5) + 
  coord_sf(crs = crs_use) + 
  theme_minimal() + 
  ggtitle("(B) right projection, wrong points using geom_point()")
```

---

```r
B
```

```
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
```

![](sp_viz_files/figure-html/unnamed-chunk-44-1.png)

---

```r
C <- ggplot() + 
  geom_sf(data = usa_sf()) + 
  geom_sf(data = st_as_sf(d_points,
                          coords = c("long", "lat"), crs = crs_use), 
          color = "red", size = 5) + 
  coord_sf(crs = crs_use) + 
  theme_minimal() + 
  ggtitle("(C) right projection, wrong points using geom_sf() points")
```

---

```r
C
```

```
## old-style crs object detected; please recreate object with a recent sf::st_crs()
```

![](sp_viz_files/figure-html/unnamed-chunk-46-1.png)

---

```r
#remotes::install_github("ropensci/USAboundaries")
#install.packages("USAboundariesData", repos = "https://ropensci.r-universe.dev", type = "source")
library(USAboundaries)
library(tidyverse)
crs_use <- "+proj=laea +lat_0=30 +lon_0=-95"

usa_sf <- us_boundaries(type="state", resolution = "low") %>% 
  dplyr::filter(!state_abbr %in% c("PR", "AK", "HI")) %>% 
  st_transform(crs = crs_use)

d_points <- data.frame(long = c(-110, -103, -84), 
                      lat  = c(45, 40, 41)) %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>% 
  st_transform(crs = crs_use)

D <- ggplot(data = usa_sf) +
  geom_sf() + 
  geom_sf(data = d_points, 
             color = "red", size = 5) + 
  theme_minimal()
```

---

```r
D
```

```
## function (expr, name) 
## .External(C_doD, expr, name)
## <bytecode: 0x7fbf8041fc40>
## <environment: namespace:stats>
```

---
`st_as_sf()`: what coordinate system do the coordinates you just gave me use?

`st_transform()`: what coordinate system should I transform those coordinates into?

**Analogy:** if you were translating a book,

the crs in `st_as_sf()` is "what language is the book written in?"

`st_transform()` is "what language am I translating this to?"

---

## What does CRS 4326 mean?

- The Three Most Common Coordinate Reference Systems

- WGS84 and EPSG:4326 are geographic CRS.

- It means that they use latitude and longitude coordinates to specify a location on the surface of the earth.

---

## Projecting your data

**Additional reading**

Click [here](https://www.earthdatascience.org/courses/earth-analytics/spatial-data-r/intro-to-coordinate-reference-systems/)

---

## Interactive chart

```r
library(tidyverse)
library(sf)
library(viridis)
#devtools::install_github("tylermorganwall/rayshader")
library(rayshader)

sightings = st_read(dsn = "~/Dropbox/R Code Repository/Sasquatch/WA_OR_Bigfoot.shp")
sightings = mutate(sightings, Sightings = Join_Count)

outputFile <- file.path("~/Dropbox/R Code Repository/Sasquatch/", "sasquatch_sightings.mp4")
```

Resources: https://geospatialtraining.com/3d-mapping-of-sasquatch-sightings-in-r-with-rayshader/
---

```r
sightings
```

---

```r
ggSightings <- ggplot(data = sightings) + 
  geom_sf(aes(fill = Sightings)) + 
  scale_fill_viridis("Sightings") +
  ggtitle("Sasquatch Sightings - Washington and Oregon") +
  theme_bw()
```

---

```r
ggSightings
```

---

```r
plot_gg(ggSightings,multicore=TRUE,width=5,height=5,scale=250,windowsize=c(1400,866),
        zoom = 0.55, phi = 30)
render_snapshot()
render_movie(filename = outputFile)
```

---
class: center, middle

# 2. Polygon map

---

```r
# devtools::install_github("UrbanInstitute/urbnmapr")
library(tidyverse)
library(magrittr)
```

```
## 
## Attaching package: 'magrittr'
```

```
## The following object is masked from 'package:raster':
## 
##     extract
```

```
## The following object is masked from 'package:purrr':
## 
##     set_names
```

```
## The following object is masked from 'package:tidyr':
## 
##     extract
```

```r
library(urbnmapr)
library(usmap) # to obtain population data
library(ggthemes)
library(scales)
```

```
## 
## Attaching package: 'scales'
```

```
## The following object is masked from 'package:viridis':
## 
##     viridis_pal
```

```
## The following object is masked from 'package:purrr':
## 
##     discard
```

```
## The following object is masked from 'package:readr':
## 
##     col_factor
```

```r
library(geofacet)
```

---

```r
statepop1 <- statepop %>% rename(state_name = full)
statepop2 <- full_join(statepop1, states, by = "state_name")
head(statepop2)
```

```
## # A tibble: 6 × 12
##   fips  abbr  state_name pop_2015  long   lat order hole  piece group state_fips
##   <chr> <chr> <chr>         <dbl> <dbl> <dbl> <int> <lgl> <fct> <fct> <chr>     
## 1 01    AL    Alabama     4858979 -88.5  31.9     1 FALSE 1     01.1  01        
## 2 01    AL    Alabama     4858979 -88.5  31.9     2 FALSE 1     01.1  01        
## 3 01    AL    Alabama     4858979 -88.5  31.9     3 FALSE 1     01.1  01        
## 4 01    AL    Alabama     4858979 -88.5  32.0     4 FALSE 1     01.1  01        
## 5 01    AL    Alabama     4858979 -88.5  32.0     5 FALSE 1     01.1  01        
## 6 01    AL    Alabama     4858979 -88.5  32.1     6 FALSE 1     01.1  01        
## # ℹ 1 more variable: state_abbv <chr>
```

---

```r
ggplot() + 
  geom_point(data = statepop2, mapping = aes(x = long, y = lat, group = group), fill = "grey", color = "black") 
```

![](sp_viz_files/figure-html/unnamed-chunk-56-1.png)

---

```r
ggplot() + 
  geom_polygon(data = statepop2, mapping = aes(x = long, y = lat, group = group),fill = "grey", color = "black") 
```

![](sp_viz_files/figure-html/unnamed-chunk-57-1.png)

---

Visit here: https://tstdataviz.netlify.app/slides/spatiotemporal.html
---

# Moving forward

### Visualizing using webmaps

```r
library(leaflet)
```

### Visualization of raster data

```r
library(raster)
```

---