Assignment 3; More Advanced Mapping

library(terra)

## Warning: package 'terra' was built under R version 4.3.3

## terra 1.7.71

library(tidyterra)

## Warning: package 'tidyterra' was built under R version 4.3.3

## 
## Attaching package: 'tidyterra'

## The following object is masked from 'package:stats':
## 
##     filter

library(sp)
library(ggplot2)
library(ggthemes)
library(urbnmapr)
library(sf)

## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE

#First we load in all the libraries needed for our project

setwd("C:/Users/lizzi/OneDrive/Desktop/DIDA 370/NE1_50M_SR_W/NE1_50M_SR_W")
my_rast<- rast("NE1_50M_SR_W.tif")
plot(my_rast)

#Changed directory and loaded in the tif file, did not re project raster object and plotted the earth map

setwd("C:/Users/lizzi/OneDrive/Desktop/DIDA 370/states/states")
states<- st_read("state_map.shp")

## Reading layer `state_map' from data source 
##   `C:\Users\lizzi\OneDrive\Desktop\DIDA 370\states\states\state_map.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 51 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -125.4117 ymin: 21.73154 xmax: -66.9499 ymax: 49.38436
## Geodetic CRS:  WGS 84

#Loaded in the states map file and was able to use the crop function since the projection matched.
us_rast<- crop(my_rast, states)
plot(us_rast)

ggplot()+
  geom_spatraster(data = us_rast)+
  geom_sf(states, mapping = aes())

## <SpatRaster> resampled to 500636 cells for plotting
## ! `tidyterra::geom_spatraster()`: Plotting 3 overlapping layers: NE1_50M_SR_W_1, NE1_50M_SR_W_2, and NE1_50M_SR_W_3. Either:
## 
##   Use `facet_wrap(~lyr)` for faceting or
## 
##   Use `aes(fill = <name_of_layer>)` for displaying single layers

#plotted both raster objects and the geom layer of the states on top.
#map doesn't look good so we need to add some customizations

us_map<- ggplot()+
  geom_spatraster(data = us_rast)+
  geom_sf(states, mapping = aes(), fill= "transparent", color = "black")+ #made the fill transparent to remove the white
  theme_map()+
  scale_fill_whitebox_c(
    palette = "deep",
    na.value = "white")+# tried out different palettes and decided on deep because it looked better and changed borders to black to be more apparent.
  labs(title = "United States of America",
       fill = "Fill")+
  theme(plot.title = element_text(hjust = 0.5, size = 16),
        legend.position="right")

## <SpatRaster> resampled to 500636 cells for plotting

us_map

## ! `tidyterra::geom_spatraster()`: Plotting 3 overlapping layers: NE1_50M_SR_W_1, NE1_50M_SR_W_2, and NE1_50M_SR_W_3. Either:
##   Use `facet_wrap(~lyr)` for faceting or
##   Use `aes(fill = <name_of_layer>)` for displaying single layers

counties <- get_urbn_map("counties", sf = TRUE) #Use counties map for this and extract only NYS
counties <- counties %>% 
  filter(state_abbv == "NY") %>% 
  st_transform("EPSG:32116")

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

## Warning in CPL_crs_from_input(x): GDAL Message 1: CRS EPSG:2163 is deprecated.
## Its non-deprecated replacement EPSG:9311 will be used instead. To use the
## original CRS, set the OSR_USE_NON_DEPRECATED configuration option to NO.

setwd("C:/Users/lizzi/OneDrive/Desktop/DIDA 370")
countydata1<- read.csv("2017_Unemployment_Rate.csv")
countydata2<- read.csv("Personal_Income_2017.csv")
#Had two data sets downloaded from the NYS data website, one showing unemployment rate in 2017 and another showing personal income. Cleaned both to ensure the merge can happen, changed personal income to numeric too.
library(tidyr)

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:terra':
## 
##     extract

counties <- counties %>% 
  separate(county_name, c("county_name", "county"), sep = " County") %>% 
  select(-county)
county_data_final<- counties %>%
  left_join(countydata1,by = c("county_name" = "County"))
county_data_final<- county_data_final %>% left_join(countydata2, by = c("county_name"= "County"))
plot(county_data_final)

## Warning in min(x): no non-missing arguments to min; returning Inf

## Warning in max(x): no non-missing arguments to max; returning -Inf

#Merged data sets and plotted them to see if the variables I wanted to plot were suitable.

unemp.sf<- st_transform(county_data_final, "EPSG:4326")
popup1 <- paste0("<strong>County: </strong>", unemp.sf$county_name," ",
                 "<strong> Unemployment_Rate </strong>", unemp.sf$Unemployment_Rate) #Added pop-up text to show both county name and unemployment rate. 
library(leaflet)

## Warning: package 'leaflet' was built under R version 4.3.3

pal_fun <- colorQuantile(palette= c("lightblue","yellow","orange","red"), domain= unemp.sf$Unemployment_Rate, n=4)
#Created my own color palette here with light blue being on the lowere end and red being on the higher end.
leaflet() %>%
  leaflet::addPolygons(data = unemp.sf, 
                       stroke = FALSE, 
                       fillColor = ~pal_fun(Unemployment_Rate), 
                       fillOpacity = 0.8, 
                       smoothFactor = 0.5, 
                       popup = popup1) %>% 
  addLegend("bottomright",
            pal=pal_fun,
            values=unemp.sf$Unemployment_Rate, 
            title = '2017 Unemployment Rate') %>% 
  addProviderTiles(providers$CartoDB.Positron) #Added the base map to make map look better. 

popup2 <- paste0("<strong>County: </strong>", unemp.sf$county_name," ",
                 "<strong> Personal_Income </strong>", unemp.sf$Personal_Income)

pal_fun1 <- colorNumeric(palette="Reds", domain= unemp.sf$Personal_Income)#created new palette function for the income since it wasn't in percentages and used reds. got this code from the github link shared in class. 
leaflet() %>%
  leaflet::addPolygons(data = unemp.sf, 
                       stroke = FALSE, 
                       fillColor = ~pal_fun1(Personal_Income), 
                       fillOpacity = 0.8, 
                       smoothFactor = 0.5, 
                       popup = popup2) %>% 
  addLegend("bottomright",
            pal=pal_fun1,
            values=unemp.sf$Personal_Income, 
            title = '2017 Personal_Income') %>% 
  addProviderTiles(providers$CartoDB.Positron)

leafletmap<- leaflet() %>%
  leaflet::addPolygons(data = unemp.sf, 
                       stroke = FALSE, 
                       fillColor = ~pal_fun1(Personal_Income), 
                       fillOpacity = 0.8, 
                       smoothFactor = 0.5, 
                       popup = popup2,
                       group= "2017 Personal_Income") %>% 
  leaflet::addPolygons(data = unemp.sf, 
                       stroke = FALSE, 
                       fillColor = ~pal_fun(Unemployment_Rate), 
                       fillOpacity = 0.8, 
                       smoothFactor = 0.5, 
                       popup = popup1,
                       group= "2017 Unemployment_Rate") %>% 
  addLegend("bottomright",
            pal=pal_fun1,
            values=unemp.sf$Personal_Income, 
            title = '2017 Personal_Income') %>%
   addLegend("bottomright",
            pal=pal_fun,
            values=unemp.sf$Unemployment_Rate, 
            title = '2017 Unemployment Rate') %>%
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addLayersControl(overlayGroups = c("2017 Personal_Income", "2017 Unemployment_Rate"),
                   options = layersControlOptions(collapsed = FALSE))
#Lastly, I combined both layers and added the Layers Control to select which layer to show on the leaflet map.
leafletmap