Projections & Point Data

class: left, middle, inverse
background-image: url("https://live.staticflickr.com/65535/50362989122_a8ee154fea_k_d.jpg")
background-size: cover

# .orange[Spatial Data]

### .fancy[Projections & Points <svg aria-hidden="true" role="img" viewBox="0 0 384 512" style="height:1em;width:0.75em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:orange;overflow:visible;position:relative;"><path d="M172.268 501.67C26.97 291.031 0 269.413 0 192 0 85.961 85.961 0 192 0s192 85.961 192 192c0 77.413-26.97 99.031-172.268 309.67-9.535 13.774-29.93 13.773-39.464 0zM192 272c44.183 0 80-35.817 80-80s-35.817-80-80-80-80 35.817-80 80 35.817 80 80 80z"/></svg>]

---
class: bottom
background-image: url("https://live.staticflickr.com/65535/50437121667_ac9c3e7f84_c_d.jpg")
background-size: 55%

.footnote[ [Knight *et al.* (2005) *Nature*.  **437**, 880-883.](https://www.researchgate.net/deref/http%3A%2F%2Fdx.doi.org%2F10.1038%2Fnature03962)]

---

# Toblers First Law of Geography

.blue[
> Everything is related to everything else, but near things are more related to each other.
]

.pull-left[
### Physical Measurements
- Location

- Distance

- Network

- Neighborhoods & Regions
]

.pull-right[
### Emerging Properties

- Spatial Heterogeneity

- Spatial Dependence

- Objects & Fields

]

.footnote[[Tobler, W. R. 1970. *Economic Geography*, **46**, 234–240.](https://doi.org/10.1002/9781118786352.wbieg1011)]

---
class: sectionTitle, inverse

# .orange[Projections]

## .fancy[Representing Location On .redinline[*This*] Planet.]

---

# A Projection

A projection is a mathematical mapping of points measured on this surface of this earth that can be represented on things like computer screens.

.center[![Example](https://live.staticflickr.com/65535/50437002623_e8df5b3c9f_c_d.jpg) ]

---
background-image: url("https://live.staticflickr.com/65535/50437046573_20efc64bdb_o_d.png")
background-position: right
background-size: auto
background-color: #00013B

# .blueinline[Earth == Lumpy Bumpy Potato?]

&nbsp;

## .redinline[Yes!    
&nbsp; &nbsp; &nbsp; - J. Ciminelli]

---

# Ellipsoids

Models of the physical structure of the surface of the planet.

- NAD83/GRS80 - Satellite measurements of distance from core to surface of earth.

- WGS84 - Model built on global GPS system.

---

# Example Data - Maps

For these examples, I'm going to be using the `maps` library1. .footnote[1If you do not have it, install it by `install.packages("maps")`.]For

```r
library( ggplot2 )
library( maps )
states <- map_data( "state" )
head(states, n=3)
```

```
## long lat group order region subregion
## 1 -87.46201 30.38968 1 1 alabama <NA>
## 2 -87.48493 30.37249 1 2 alabama <NA>
## 3 -87.52503 30.37249 1 3 alabama <NA>
```

A basic map, notice the use of `group`.

```r
ggplot( states, aes(x=long, y=lat, group=group) ) + 
  geom_polygon( fill="lightgray", color="black", lwd=0.25) + 
  theme_void() -> p
```

---

# Azimuth Projections

.pull-left[

Projected onto a 2-dimensional plane that is tangential to a single point on the surface of the earth (commonly  a pole, though not a necessity).

![Azequidistant](https://live.staticflickr.com/65535/50437120363_d8e0686d38_w_d.jpg)
]
.pull-right[

```r
p + coord_map( "azequalarea")
```

]

---

# Cylindrical Projections

.pull-left[

```r
p + coord_map( "cylindrical" )
```

<img src="slides_files/figure-html/unnamed-chunk-4-1.png" width="504" style="display: block; margin: auto;" />
]
.pull-right[
Parallels are straight lines and horizontal up and down the plot from latitude = 0

.center[
![](https://live.staticflickr.com/65535/50437120498_8dd67df3f1_w_d.jpg)
]
]

---

# Conic Projections

.pull-left[

Symmetric around the Prime Meridian and parallels are segments of concentric circles.

![](https://live.staticflickr.com/65535/50437120428_6da48bed81_w_d.jpg)

]
.pull-right[

```r
p + coord_map( "conic", lat0 = 30)
```

]

---
class: sectionTitle, inverse

# .green[Datum]

---

# Coordinate Systems

The *Datum* are the coordinate system used on the ellipsoid.  Common types include:

- Longitude & Latitude - The East/West & North/South position on the surface of the earth.

- Prime Meridian (0° Longitude) passes thorugh the [Royal Observatory](https://en.wikipedia.org/wiki/Royal_Observatory,_Greenwich) in Greenwich England, with positive values of longitude to the east and negative to the west.
  
  - Equator (0° Latitude) and is defined as the point on the planet where both northern and southern hemisphers have equal amounts of day and night at the [equinox](https://en.wikipedia.org/wiki/Equinox) (Sept. 21 & March 21).
  
  - Richmond, Virginia: 37.533333 Latitude, -77.466667 Longitude
  
--

- Universal Trans Mercator - A division of the earth into 60 zones (~6°longitude each, labeled 01 - 60) and 20 bands each of which is ~8° latitude (labeled C-X excluding I & O with A & B dividing up Antartica).  See image [here](https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#/media/File:Universal_Transverse_Mercator_zones.svg).

- Coordinates include Zone & band designation as well as coordinates in Easting and Northing (planar coordinates within the zone) measured in meters.
  
  - Richmond, Virginia: 18S 282051 4156899

---

# Functional Consequences

&nbsp;

.center[ 
.red[DATA NEED TO BE IN THE SAME ELLIPSOID & DATUM

.red[FROM the BEGINNING]]
]

---
class: sectionTitle, inverse

# .blue[Points]

---

# Beetle Data

The Bark Beetle, *Araptus attenuatus*, is a Sonoran Desert endemic parasite that lives on within the plant *Euphorbia lomelii*.

.pull-left[

&nbsp;

![Araptus attenuatus](https://live.staticflickr.com/65535/50441339417_74e04216fa_w_d.jpg)
]

.pull-right[
.center[![Euphorbia lomelii](https://live.staticflickr.com/65535/50441175211_ba3b9df2ea_w_d.jpg)]
]

---

# Sampling Sites

```r
library( readr )
url <- "https://raw.githubusercontent.com/dyerlab/ENVS-Lectures/master/data/Araptus_Disperal_Bias.csv"
read_csv( url ) %>%
 select( Site, Longitude, Latitude, everything() ) %>%
 arrange( Latitude ) -> data 
summary( data )
```

```
##      Site             Longitude         Latitude         Males      
##  Length:31          Min.   :-114.3   Min.   :23.29   Min.   : 9.00  
##  Class :character   1st Qu.:-113.1   1st Qu.:24.95   1st Qu.:16.00  
##  Mode  :character   Median :-112.0   Median :26.64   Median :21.00  
##                     Mean   :-112.0   Mean   :26.44   Mean   :25.68  
##                     3rd Qu.:-110.8   3rd Qu.:27.78   3rd Qu.:31.50  
##                     Max.   :-109.3   Max.   :29.33   Max.   :64.00  
##     Females       Suitability        MFRatio        GenVarArapat   
##  Min.   : 5.00   Min.   :0.0563   Min.   :0.5938   Min.   :0.0500  
##  1st Qu.:15.50   1st Qu.:0.2732   1st Qu.:0.8778   1st Qu.:0.1392  
##  Median :21.00   Median :0.3975   Median :1.1200   Median :0.2002  
##  Mean   :23.52   Mean   :0.4276   Mean   :1.1598   Mean   :0.2006  
##  3rd Qu.:29.00   3rd Qu.:0.5442   3rd Qu.:1.3618   3rd Qu.:0.2592  
##  Max.   :63.00   Max.   :0.9019   Max.   :2.2000   Max.   :0.3379  
##   GenVarEuphli   
##  Min.   :0.0500  
##  1st Qu.:0.1777  
##  Median :0.2171  
##  Mean   :0.2203  
##  3rd Qu.:0.2517  
##  Max.   :0.5122
```

---

# A Simple Map - Leaflet

.pull-left[
The `leaflet` library allows you to make quick, interactive maps.

```r
library( leaflet )
data %>%
 mutate( Label = paste( "Site:", Site, 
 "<hr>\nFemales:", Females, 
 " Males: ", Males,
 " Suitability:", Suitability) ) %>%
 leaflet() %>%
 addMarkers( ~Longitude, ~Latitude, popup = ~Label ) %>%
 addProviderTiles( "OpenTopoMap" )
```

]

.pull-right[
<div id="htmlwidget-bc421d5954aed62db17e" style="width:504px;height:504px;" class="leaflet html-widget"></div>
<script type="application/json" data-for="htmlwidget-bc421d5954aed62db17e">{"x":{"options":{"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}}},"calls":[{"method":"addMarkers","args":[[23.2855,24.00789,24.0195,24.13389,24.2115,24.21441,24.58843,24.87611,25.0247,25.34819,25.55757,25.60521,25.91409,26.0155,26.20876,26.63783,26.94589,27.0367,27.18232,27.2028,27.3632,27.40498,27.52944,28.03661,28.22308,28.40846,28.66056,28.72796,28.96651,29.01457,29.32541],[-110.10429,-109.85071,-110.096,-110.46236,-110.951,-110.27252,-110.74599,-110.69175,-111.675,-111.60056,-111.21563,-111.32638,-112.08062,-111.35475,-111.37833,-109.327,-112.04613,-112.986,-112.66552,-112.408,-112.964,-112.52959,-113.31609,-113.39991,-113.18263,-112.86985,-113.99141,-113.48973,-113.66787,-113.94486,-114.29353],null,null,null,{"interactive":true,"draggable":false,"keyboard":true,"title":"","alt":"","zIndexOffset":0,"opacity":1,"riseOnHover":false,"riseOffset":250},["Site: Aqu <hr>\nFemales: 9 Males: 12 Suitability: 0.7217","Site: 73 <hr>\nFemales: 5 Males: 11 Suitability: 0.1455","Site: 157 <hr>\nFemales: 30 Males: 26 Suitability: 0.881","Site: 153 <hr>\nFemales: 41 Males: 35 Suitability: 0.7325","Site: 163 <hr>\nFemales: 21 Males: 21 Suitability: 0.4329","Site: 48 <hr>\nFemales: 27 Males: 18 Suitability: 0.6195","Site: 75 <hr>\nFemales: 18 Males: 16 Suitability: 0.8656","Site: 77 <hr>\nFemales: 16 Males: 19 Suitability: 0.4226","Site: const <hr>\nFemales: 11 Males: 18 Suitability: 0.1744","Site: 51 <hr>\nFemales: 11 Males: 9 Suitability: 0.4125","Site: 168 <hr>\nFemales: 25 Males: 28 Suitability: 0.4965","Site: 64 <hr>\nFemales: 19 Males: 16 Suitability: 0.2876","Site: 166 <hr>\nFemales: 26 Males: 19 Suitability: 0.2673","Site: 58 <hr>\nFemales: 9 Males: 11 Suitability: 0.9019","Site: 169 <hr>\nFemales: 28 Males: 47 Suitability: 0.3132","Site: 32 <hr>\nFemales: 27 Males: 40 Suitability: 0.0563","Site: 93 <hr>\nFemales: 21 Males: 25 Suitability: 0.1627","Site: 161 <hr>\nFemales: 63 Males: 64 Suitability: 0.2791","Site: 12 <hr>\nFemales: 21 Males: 24 Suitability: 0.3519","Site: 162 <hr>\nFemales: 41 Males: 57 Suitability: 0.6136","Site: sfran <hr>\nFemales: 16 Males: 26 Suitability: 0.3211","Site: 160 <hr>\nFemales: 36 Males: 48 Suitability: 0.3652","Site: 159 <hr>\nFemales: 15 Males: 22 Suitability: 0.188","Site: 89 <hr>\nFemales: 19 Males: 18 Suitability: 0.3975","Site: 171 <hr>\nFemales: 39 Males: 38 Suitability: 0.4614","Site: 173 <hr>\nFemales: 32 Males: 19 Suitability: 0.5243","Site: 177 <hr>\nFemales: 50 Males: 49 Suitability: 0.2618","Site: 175 <hr>\nFemales: 8 Males: 13 Suitability: 0.4214","Site: 84 <hr>\nFemales: 17 Males: 13 Suitability: 0.3953","Site: 9 <hr>\nFemales: 10 Males: 11 Suitability: 0.5641","Site: 88 <hr>\nFemales: 18 Males: 23 Suitability: 0.2186"],null,null,null,null,{"interactive":false,"permanent":false,"direction":"auto","opacity":1,"offset":[0,0],"textsize":"10px","textOnly":false,"className":"","sticky":true},null]},{"method":"addProviderTiles","args":["OpenTopoMap",null,null,{"errorTileUrl":"","noWrap":false,"detectRetina":false}]}],"limits":{"lat":[23.2855,29.32541],"lng":[-114.29353,-109.327]}},"evals":[],"jsHooks":[]}</script>
]

---

# The Simple Features (`sf`) library.

.orangeinline[Simple Features] are an open standard developed by the Open Geospatial Consortium ([OGC](https://ogc.org)).  They define the following basic types:

.pull-left[

### Basic Geometry Types

- POINT

- LINESTRING

- POLYGON  
]

.pull-right[
### Corresponding 'multi' versions

-  MULTIPOINT

- MULTILINESTRING

- MULTIPOLYGON

- GEOMETRYCOLLECTION

]

Each of these basic types can be represented within a column of a `data.frame`.

---

# Setting `sf` objects

To convert from a normal `data.frame` with columns of numerical values to an `sf` object, you .redinline[**must**] indicate:  
- Names of the `x` and `y` coordinates, and   
- A reference to the Coordinate Reference System

```r
data %>% 
  st_as_sf( coords=c("Longitude","Latitude"), crs=4326 ) -> data
head( data )
```

```
## Simple feature collection with 6 features and 7 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -110.951 ymin: 23.2855 xmax: -109.8507 ymax: 24.21441
## Geodetic CRS: WGS 84
## # A tibble: 6 × 8
## Site Males Females Suitability MFRatio GenVarArapat GenVarEuphli
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Aqu 12 9 0.722 1.33 0.120 0.0968
## 2 73 11 5 0.146 2.2 0.137 0.253 
## 3 157 26 30 0.881 0.867 0.150 0.191 
## 4 153 35 41 0.732 0.854 0.333 0.276 
## 5 163 21 21 0.433 1 0.298 0.338 
## 6 48 18 27 0.620 0.667 0.115 0.213 
## # … with 1 more variable: geometry <POINT [°]>
```

---
background-image: url("https://live.staticflickr.com/65535/50445538736_bc1c3456e5_k_d.jpg)
background-size: cover

---

# Reprojecting Points

We can reproject data (that **already** has a coordinate reference system) into any other projection.

```r
data %>%
  st_transform( 6372 ) %>%
  st_bbox()
```

```
##    xmin    ymin    xmax    ymax 
## 1307745 1274010 1773676 1968473
```

---
class: sectionTitle, inverse

# .red[Plotting `sf` Objects]

---

.pull-left[
# Plotting `sf` Objects

Eash of the data columns are associated with the `geometry` column we made.  If we plot the whole `data.frame`, it makes replicate plots of each data column displayed in the coordinate space of the `geometry`.

- Color will represent gradient in values for each plot.

- Normal `plot()` parameters can be used to customize.
]

.pull-right[

```r
plot( data , pch=16, cex=2)
```

<img src="slides_files/figure-html/unnamed-chunk-11-1.png" width="504" style="display: block; margin: auto;" />
]

---

# Linked data and `geometry`

.pull-left[

```r
plot( data$Suitability, pch=16, cex=2, bty="n" )
```

<img src="slides_files/figure-html/unnamed-chunk-12-1.png" width="432" style="display: block; margin: auto;" />
]

.pull-right[

```r
plot( data["Suitability"], pch=16, cex=2)
```

<img src="slides_files/figure-html/unnamed-chunk-13-1.png" width="432" style="display: block; margin: auto;" />
]

---

# Hello `ggplot` My Old Friend

```r
ggplot( data ) + geom_sf( aes(size=Suitability) ) 
```

---

.pull-left[
# GGPlot and Labels

```r
ggplot( data ) + 
  geom_sf_label( aes(label=Site)) + 
  theme_void() + 
  coord_map()
```
]

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

---
class: sectionTitle, inverse

# .blue[ Map Overlays 🗺]

## .fancy[Making it look better]

---

# Maps

As with other `ggplot` activities, we can mix the data being plot by changing `data` and `aes` mappings.

.pull-left[

Grab the Mexico map data like before.

```r
map_data("world") %>%
  filter( region == "Mexico") -> map

ggplot( ) + 
  geom_polygon( aes( x=long, 
                     y=lat, 
                     group=group ), 
                data=map, 
                fill="grey" ) + 
  geom_sf( data=data, 
           aes(color=Suitability), 
           size=2) +
  xlab("Longitude") + 
  ylab("Latitude") + 
  theme_bw( base_size = 12 ) +
  coord_sf( xlim = c(-115, -105),
            ylim = c(20, 30) )
```

]

.pull-right[ 
<img src="slides_files/figure-html/unnamed-chunk-18-1.png" width="432" style="display: block; margin: auto;" />
]

---

# Natural Earth

.pull-left[
<img src="slides_files/figure-html/unnamed-chunk-19-1.png" width="504" style="display: block; margin: auto;" />
]

.pull-right[

Loading `sf` object for background map from `naturalearth` libraries and then plotting.  *Notice:* `xlim` and `ylim` are configured withing `coord_sf()`.

```r
library( rnaturalearth )
library( rnaturalearthdata )

world <- ne_countries( scale = "medium", returnclass = "sf")

ggplot() + 
  geom_sf( data=world ) + 
  geom_sf( data=data, 
           aes(color=Suitability)) +
  coord_sf( xlim = c(-115, -105),
            ylim = c(20, 30) ) + 
  xlab("Longitude") + 
  ylab("Latitude")
```
]

---
class: sectionTitle, inverse

# .green[Spatial Operations]

## .fancy[What can we do with this now?]

---

# Geospatial Operations

### aka the `Spatial Toolbox`

Once the data are in `sf` features, there are a lot of different spatial operations that we can use, including:

.pull-left[
- Finding the bounding box

```r
data %>% 
  st_bbox()
```

```
##       xmin       ymin       xmax       ymax 
## -114.29353   23.28550 -109.32700   29.32541
```

- Finding the convex hull

```r
data %>% 
  st_union() %>% 
  st_convex_hull() -> hull
```

]

.pull-right[
- Estimating its centroid & and area.

```r
hull %>%
  st_centroid()
```

```
## Geometry set for 1 feature 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -111.3417 ymin: 26.37741 xmax: -111.3417 ymax: 26.37741
## Geodetic CRS:  WGS 84
```

```r
library( units )
hull %>%
  st_area() %>%
  set_units( km^2 )
```

```
## 122130.5 [km^2]
```

]

---

class: middle
background-image: url("images/contour.png")
background-position: right
background-size: auto

.center[

# 🙋🏻‍♀️ Questions?

![Peter Sellers](https://live.staticflickr.com/65535/50382906427_2845eb1861_o_d.gif+)
]

&nbsp;

.bottom[ If you have any questions for about the content presented herein, please feel free to [submit them to me](mailto://rjdyer@vcu.edu) and I'll get back to you as soon as possible.]