Spatial data is special data

name:1
class: center middle main-title section-title-4

# Spatial data is special data

#### (Or - how does R keep track of all of this spatial stuff?)

**Session 3**

.light[HES597: Introduction to Spatial Data in R<br>
Boise State University Human-Environment Systems<br>
Fall 2021]

]

---

# Plan for today

Location, Location, Location

The Coordinate Reference System

Other elements that define location in R objects

Parting thoughts

---

# Questions we can answer with GIS

- Questions about geographic location
 
 - Questions about geographic distribution
 
 - Questions about geographic association
 
 - Questions about geographic interaction
 
 - Questions about geographic change

---

# Questions we can answer with GIS

.primary[- __Questions about geographic location__]
 
 .white[
 - Questions about geographic distribution
 
 - Questions about geographic association
 
 - Questions about geographic interaction
 
 - Questions about geographic change
 ]

---

layout: false
name: projections
class: center middle section-title section-title-2 animated fadeIn

# Location, Location, Location!!
### __(sure, but how do we actually describe where things are?)__
---

# Location
---
class: left middle

# Location

## __nominal: (potentially contested) place names__

## __absolute: the physical location on the earth's surface__

---
class: left middle

# Location

## __absolute: the physical location on the earth's surface__

---
# Coordinates: 2 or more measurements that specify location relative to a _reference system_

.pull-left[
- Cartesian coordinate system: grid formed by combining horizontal (*x*) and vertical (*y*) measurement scales

- _origin (O)_ = the point at which both measurement systems intersect

- By convention this is equal (0,0) and all measurements relative to the origin

- Adaptable to multiple dimensions (e.g. *z* for altitude)
]

.pull-right.center[
<figure>
  <img src="img/02/CartesianCoordinateSystem.png" alt="ZZZ" title="ZZZ" width="100%">
  <figcaption><a href="https://commons.wikimedia.org/wiki/File:CartesianCoordinateSystem.png">Svjo</a>, <a href="https://creativecommons.org/licenses/by-sa/4.0">CC BY-SA 4.0</a>, via Wikimedia Commons</figcaption>
</figure>
]
---

# Coordinates: 2 or more measurements that specify location relative to a _reference system_

.pull-left.center[
<figure>
  <img src="img/02/Latitude_and_Longitude.png" alt="ZZZ" title="ZZZ" width="100%">
  <figcaption><a href="https://commons.wikimedia.org/wiki/File:Latitude_and_Longitude_of_the_Earth_fr.svg">Djexplo</a>, CC0, via Wikimedia Commons</figcaption>
</figure>
]

- Geographic coordinate system uses longitude (*x*) and latitude (*y*)

- Horizontal location measured as longitude in degrees from prime meridian

- Vertical location measured as latitude in degrees from equator

- _Graticule_: the grid formed by the intersection of longitude and latitude

- The graticule is based on an ellipsoid model of earth's surface and contained in the _datum_

]
---

# Geographic Coordinate Systems: representing the earth's surface in lattitude and longitude

### __Geographic Coordinate Systems define *where* the data is located on the earth's surface__

- Generally round or ellipsoid so records location in angular units (e.g., degrees)
  - Multiple GCS because earth's surface is not perfectly round or smooth

### __The *datum* describes which ellipsoid to use and the precise relations between locations on earth's surface and Cartesian coordinates__

- Geodetic datums (e.g., `WGS84`): origin is based on the earth's center of gravity not optimized for any region.
  
  - Local datums (e.g., `NAD83`): allows for local variation in earth's surface to be modeled more accurately

---
# Projected Coordinate Systems (CRS)

### __Projected Coordinate Systems describe *how* the data should be translated to a flat surface__

- Necessary because maps, computer screens, and most R data structures are flat
  - Records locations in linear units (e.g., meters)
  - Three developable surfaces: planar (poles), conic (mid-latitudes), cylindrical (whole earth)
  
### __Projection necessarily induces some form of distortion (tearing, compression, or shearing__
   
   - Some projections minimize distortion of angle, area, or distance
   - Others attempt to avoid extreme distortion of any kind 
   - Distortions particularly challenging for raster data
   
### __How might you choose which projection you should use?__   
---
name: crs
# The Coordinate Reference System (CRS)

- The __CRS__ stores all of the necessary information for the projection calculations

- Includes: Datum, ellipsoid, units, and other information (e.g., False Easting, Central Meridian) to further map the projection to the GCS

- Not all projections have/require all of the parameters

- R stores these data in several formats ([EPSG](https://epsg.io/), [Proj](https://proj.org/), and [WKT](https://www.ogc.org/standards/wkt-crs]))

- Lots of projection info available at [spatialreference.org](https://spatialreference.org/)
---
# The Coordinate Reference System

### PROJ4 (deprecated)

`+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs`

### EPSG

`EPSG:4087`

### WKT

`PROJCS["WGS 84 / World Equidistant Cylindrical",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Equirectangular"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","4087"]]`

---
# The Coordinate Reference System (CRS)
Using WGS84 World Equidistant Cylindrical

```r
library(sf)
nc.sf <-  st_read(
  system.file("shapes/", package="maptools"), 
  "sids")
```

```
## Reading layer `sids' from data source 
##   `/Library/Frameworks/R.framework/Versions/4.0/Resources/library/maptools/shapes' 
##   using driver `ESRI Shapefile'
## 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
## CRS:           NA
```

```r
nc.sf <- nc.sf %>% 
  st_set_crs(4267) %>% 
  st_transform(4087)
#plot(st_geometry(nc.sf))
```
]

.pull-right[
<img src="02-slides_files/figure-html/plotproj1-1.png" width="504" style="display: block; margin: auto;" />

]
---
# The Coordinate Reference System (CRS)
Using NAD27 a local geodetic datum for North America

```r
library(sf)
nc.sf <-  st_read(
  system.file("shapes/", package="maptools"), 
  "sids")
```

```r
nc.sf <- nc.sf %>% 
  st_set_crs(4267) 
#plot(st_geometry(nc.sf))
```
]

.pull-right[
<img src="02-slides_files/figure-html/plotproj2-1.png" width="504" style="display: block; margin: auto;" />

]
---
# The Coordinate Reference System (CRS)
Using NAD83 adjusted for North Carolina

```r
library(sf)
nc.sf <-  st_read(
  system.file("shapes/", package="maptools"), 
  "sids")
```

```r
nc.sf <- nc.sf %>% 
  st_set_crs(4267) %>% 
  st_transform(3358)
#plot(st_geometry(nc.sf))
```
]

.pull-right[
<img src="02-slides_files/figure-html/plotproj3-1.png" width="504" style="display: block; margin: auto;" />

]
---
# The Coordinate Reference System (CRS)
Using UTM Zone 17N

```r
library(sf)
nc.sf <-  st_read(
  system.file("shapes/", package="maptools"), 
  "sids")
```

```r
nc.sf <- nc.sf %>% 
  st_set_crs(4267) %>% 
  st_transform(3747)
#plot(st_geometry(nc.sf))
```
]

.pull-right[
<img src="02-slides_files/figure-html/plotproj4-1.png" width="504" style="display: block; margin: auto;" />
]

---
.pull-left[
<figure>
  <img src="img/02/EPSG4087.png" alt="ZZZ" title="ZZZ" width="100%">
</figure>

<figure>
  <img src="img/02/EPSG4267.png" alt="ZZZ" title="ZZZ" width="100%">
</figure>
]

<figure>
  <img src="img/02/EPSG3747.png" alt="ZZZ" title="ZZZ" width="100%">
</figure>
]
---
# Accessing the CRS of Spatial* objects

The `sp` and `rgdal`packages are foundational to the `R` spatial ecosystem

```r
library(sp)
library(rgdal)
nc.sp <-  readOGR(
  system.file("shapes/", package="maptools"), 
  "sids", verbose = FALSE)
```
]

```r
proj4string(nc.sp)
```

```
## [1] NA
```

```r
proj4string(nc.sp) <- CRS("+init=epsg:4267")
proj4string(nc.sp)
```

```
## [1] "+proj=longlat +datum=NAD27 +no_defs"
```

```r
nc.sp.proj <- spTransform(nc.sp, CRS("+init=epsg:3747"))
proj4string(nc.sp.proj)
```

```
## [1] "+proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
```
]

The warning message you see is a function of the recent change from proj4 to proj6; details are [here](https://cran.r-project.org/web/packages/rgdal/vignettes/PROJ6_GDAL3.html)

---
name: elements
# Other important spatial attributes: Extent and Resolution

- **Extent** describes the area covered by the analysis

- **Resolution** describes the smallest unit being mapped

- Often combined (ambiguously) to describe the _scale_ of an analysis

---

# Extent in `R`

`R` has a very specific definition of extent: the rectangular region encompassed by the data

```r
nc.sf <-  st_read(
  system.file("shapes/", package="maptools"), 
  "sids", quiet = TRUE)

plot(st_geometry(nc.sf))
```

<img src="02-slides_files/figure-html/sfextent-1.png" width="50%" style="display: block; margin: auto;" />
]

```r
library(raster)
f <- system.file("external/test.grd", package="raster")
r <- raster(f)
plot(r)
```

<img src="02-slides_files/figure-html/rastext-1.png" width="50%" style="display: block; margin: auto;" />
]

---

# Extent in `R`

```r
nc.sf <-  st_read(
  system.file("shapes/", package="maptools"), 
  "sids", quiet = TRUE)
nc.sf %>% 
  st_set_crs(4267) %>% 
  st_transform(3358) %>% 
  st_bbox()
```

```
##      xmin      ymin      xmax      ymax 
## 123830.15  14739.33 930518.79 318254.74
```
]

```r
library(raster)
f <- system.file("external/test.grd", package="raster")
r <- raster(f)
extent(r)
```

```
## class      : Extent 
## xmin       : 178400 
## xmax       : 181600 
## ymin       : 329400 
## ymax       : 334000
```

```r
origin(r)
```

```
## [1] 0 0
```
]

---

# Resolution in `R`

- Thematically defined for `vector` datasets

- Based on pixel-size for `raster` datasets

```r
library(raster)
f <- system.file("external/test.grd", package="raster")
r <- raster(f)
crs(r)
```

```
## CRS arguments:
##  +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889
## +k=0.9999079 +x_0=155000 +y_0=463000 +datum=WGS84 +units=m +no_defs
```

```r
res(r)
```

```
## [1] 40 40
```

---

---

```r
nc.sp
```

```
## class       : SpatialPolygonsDataFrame 
## features    : 100 
## extent      : -84.32385, -75.45698, 33.88199, 36.58965  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=NAD27 +no_defs 
## variables   : 14
## names       :  AREA, PERIMETER, CNTY_, CNTY_ID,     NAME,  FIPS, FIPSNO, CRESS_ID, BIR74, SID74, NWBIR74, BIR79, SID79, NWBIR79 
## min values  : 0.042,     0.999,  1825,    1825, Alamance, 37001,  37001,        1,   248,     0,       1,   319,     0,       3 
## max values  : 0.241,      3.64,  2241,    2241,   Yancey, 37199,  37199,      100, 21588,    44,    8027, 30757,    57,   11631
```

---

```r
nc.sf
```

```
## 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
## CRS:           NA
## 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...
```

---

```r
r
```

```
## class      : RasterLayer 
## dimensions : 115, 80, 9200  (nrow, ncol, ncell)
## resolution : 40, 40  (x, y)
## extent     : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
## crs        : +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +datum=WGS84 +units=m +no_defs 
## source     : test.grd 
## names      : test 
## values     : 138.7071, 1736.058  (min, max)
```

---
# Final considerations

- `R` is **flexible**: you can manipulate most of these attributes

- The changes you make need to be **coherent** (often requires multiple changes)

- CRS choices for visualization are about accuracy, honesty, and clarity

- CRS choices for analysis based on consistency and minimizing locational error