R Notebook

library(terra)

## terra 1.7.71

library(leaflet)
library(tidyterra)

## 
## 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.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE

setwd("~/Desktop/Spatial R 370")

ny_rast <- rast("./NE1_50M_SR_W/NE1_50M_SR_W.tif")
# ny_rast <- project(ny_rast, "EPSG:32116", method = "bilinear")

ny_state <- st_read("~/Desktop/Spatial R 370/states/state_map.shp")

## Reading layer `state_map' from data source 
##   `/Users/kenlu/Desktop/Spatial R 370/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

plot(ny_rast)

us_rast <- crop(ny_rast, ny_state)
ggplot()+
  geom_spatraster(data = us_rast)+
  geom_sf(ny_state, mapping = aes(), fill = 'transparent', color = "grey")+
  theme_map()+
  scale_fill_whitebox_c(
    palette = "soft",
    na.value = "white")+
  labs(title = "United States",
     fill = "Fill")+
  theme_minimal()

## <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

# Education attainment in US dataset
# The data set provides the % of people in the states that actually attain a certain education degree
# We will be plotting the % that attained a bachelor's degree in leaflet
# I got 'st_as_sf' from https://r-spatial.github.io/sf/reference/st_as_sf.html 
# I couldnt get left_join to work so i just added a geometry column to my dataset and converted it to an sf object
library(tidyr)

## 
## Attaching package: 'tidyr'

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

setwd("~/Desktop/Spatial R 370")
education_data <-read.csv("education_data.csv")
geometry_data <- ny_state$geometry
state_data_full <- education_data %>% 
  mutate(geometry = geometry_data)
ny_state.sf <- st_as_sf(state_data_full)

ny_state.sf = st_transform(ny_state.sf, "EPSG:4326")
leaflet() %>%
  leaflet::addPolygons(data = ny_state.sf,
                       color = "green")

# Add pop up text so that when we click on a state it will display us the correct data 

p_popup <- paste0("<strong>Bachelor's Degree education attainment: </strong>", ny_state.sf$Bachelor.s.Degree)

pal_fun <- colorNumeric("YlOrRd", ny_state.sf$Bachelor.s.Degree)


leaflet() %>%
  leaflet::addPolygons(data = ny_state.sf, 
                       stroke = FALSE, 
                       fillColor = ~pal_fun(Bachelor.s.Degree), 
                       fillOpacity = 0.8, 
                       smoothFactor = 0.5, 
                       popup = p_popup) %>% 
  addLegend("bottomright",  
            pal=pal_fun,    
            values= ny_state.sf$Bachelor.s.Degree, 
            title = "Bachelor's Degree") %>% 
  addTiles()