```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```

##Overview
In this mini session, we will learn how to read in spatial data using the `sf` package, learn about the proporties of spatial data, and conduct a mini analysis using NYS county data and Census data from the `tidycensus` package. We will also create some maps with `ggplot2` and `mapview`.

##Required packages

```{r}
#install.packages("tidyverse")
#install.packages("tigris")
#install.packages("tidycensus")
#install.packages("sf")
#install.packages("spData")
#install.packages("mapview")

```
##Reading in data
+ read in vector data with the `sf` package and **st_read()** function 
+ read in raster data using the `raster` package and either the **raster()** for single-band rasters or **brick()** for multi-band rasters 

Here we will use the data from [NYS GIS dataset](http://gis.ny.gov/gisdata/inventories/details.cfm?DSID=927
) to read in a NYS county shapefile.

In the `sf` package:

+ spatial objects are stored as dataframes with some special properties
+ sf dataframes include spatial metadata like the coordinate reference system
+ geometry is stored in a list column 
+ tidyverse tools can be used with sf dataframes

### using st_read
```{r}
library(sf)

counties <- st_read("X:/spatial/NYS_Civil_Boundaries_SHP/Counties.shp")

counties

```
The short report printed mentions that there are 62 features (records, represented as rows) and 14 fields (attributes, represented as columns). The defined CRS in this data is the Universal Transverse Mercator (UTM), which uses a 2-dimensional Cartesian coordinate system to give locations on the surface of the Earth.

### using tigris package
```{r, results = FALSE}
library(tigris)

nys <- counties("New York", class= "sf")

nys2 <- counties("New York") #Output as Large SpatialPolygonsDataFrame (more complex, harder to decompose)

str(nys2) #compare to sf objects 
```

Methods of `sf` objects are:

```{r}
methods(class = "sf")
```

## Simple Plot
```{r}
plot(counties[4], reset = FALSE) # reset = FALSE: we want to add to a plot with a legend
plot(counties[1,1], col = 'grey', add = TRUE) #color first county as grey

```

## Using functions in sf package

+ **st_area()** returns the area of your features
+ **st_length()** returns the length of your features
+ **st_crs()** to assign a CRS to the data
+ **st_transform()** to transform CRS to match CRS for different layers

```{r}

library(tidyverse)


#compute areas of counties
areas <- st_area(counties)

summary(areas)
class(areas)

# Create a quick histogram of the areas using hist
hist(areas, xlim = c(0,  7.300e+09 ), breaks = 30)


#dplyr example
big_counties <- counties %>% filter(unclass(areas) > 3.233e+09)

plot(st_geometry(big_counties))

```

Note that the result of functions like **st_area()** and **st_length()** will not be a traditional vector. Instead the result has a class of units which means the vector result is accompanied by metadata describing the object's units. 



## Loading data from spData
```{r}
library("spData")

#get US-states data 
data(us_states)

class(us_states)

# Plot the us_states object using all defaults
plot(us_states)


# Plot just the total_pop_10 attribute of the us_states data
plot(us_states["total_pop_10"])

# Create a new object of just the us_states geometry
states_geo <- st_geometry(us_states)

# Plot the geometry of the us_states data
plot(states_geo)

```


## Using tidycensus
The United States Census Bureau has a number of different datasets that are available to anyone for free. `tidycensus` is an R package that allows users to interface with the US Census Bureau's decennial Census and five-year American Community APIs and return tidyverse-ready data frames, optionally with simple feature geometry included. 


API keys can be optained through (Census API website) [https://api.census.gov/data/key_signup.html].


```{r, include=FALSE}
library(tidycensus)

api_key <- "enter your api key here"
census_api_key(api_key)

Sys.getenv("CENSUS_API_KEY")
```



```{r}
#get data from 2010 decennial Census
state <- get_decennial(geography = "state", variables = "P001001")

state

#Get a dataset of median home values from the 2017 1-year ACS 
state_home <- get_acs(geography = "state",
                     variables = "B25077_001",
                     survey = "acs1",
                     geometry = TRUE,
                     shift_geo = TRUE)


```

### Variable search within tidycensus

To find Census variable IDs, use:

 + Online resources like Census Reporter
 + Built-in variable searching in tidycensus
 
```{r}

v17 <- load_variables(year = 2017,
           dataset = "acs1",
           cache = TRUE)

# Filter for table B19001
filter(v17, str_detect(name, "B19001"))

# Use median value to search for related variables
filter(v17, str_detect(label, fixed("median value", 
                                ignore_case = TRUE)))
```

### tidyverse example
```{r}
# Map through ACS1 estimates to see how they change through the years

ny_cities <- map_df(2012:2016, function(x) {
  get_acs(geography = "place", 
          variables = c(totalpop = "B01003_001"),  #total population
          state = "NY", 
          survey = "acs1", 
          year = x) %>%
    mutate(year = x)
})

ny_cities %>% arrange(NAME, year)
```


## Mapping data

```{r}
#simple choropleth plot of median home values
plot(state_home["estimate"])

#make a choropleth map with ggplot2
ggplot(state_home, aes(fill = estimate)) + geom_sf()

#use the viridis palettes for choropleth mapping
ggplot(state_home, aes(fill = estimate, col = estimate)) + 
  geom_sf() + 
  scale_fill_viridis_c() +  
  scale_color_viridis_c()
```

```{r}
ggplot(state_home, aes(fill = estimate, color = estimate)) + 
  geom_sf() + 
  scale_fill_viridis_c(labels = scales::dollar) +  
  scale_color_viridis_c(guide = FALSE) + 
  theme_minimal() + 
  coord_sf(crs = 26911, datum = NA) + 
  labs(title = "Median owner-occupied housing value by state", 
       subtitle = "USA", 
       caption = "Data source: 2017 ACS.\nData acquired with the R tidycensus package.", 
       fill = "ACS estimate")
```

### Graduated symbol maps

```{r}
#generate point centers for each US state
centers <- st_centroid(state_home)

#create graduated symbol map
ggplot() + 
  geom_sf(data = state_home, fill = "white") + 
  geom_sf(data = centers, aes(size = estimate), shape = 21, 
          fill = "lightblue", alpha = 0.7, show.legend = "point") + 
  scale_size_continuous(range = c(1, 20)) #the scale_size_continuous() function adjusts the range of sizes on the map.

```

### Interactive maps with mapview
```{r}
library(mapview)

m <- mapview(state_home, 
         zcol = "estimate")
m@map

```


References:

+ Pebesma, E., 2018. Simple Features for R: Standardized Support for Spatial Vector
  Data. The R Journal, https://journal.r-project.org/archive/2018/RJ-2018-009/
+ Kyle Walker (2019). tidycensus: Load US Census Boundary and Attribute Data as
  'tidyverse' and 'sf'-Ready Data Frames. R package version 0.9.
  https://CRAN.R-project.org/package=tidycensus
