#replace my path with yours to the Lab 5 folder 
# setwd("/Users/xfliang/University of Michigan Dropbox/Xiaofan Liang/UM_Teaching/URP535_Urban_Informatics/W25/Lab/Lab10/")
setwd("G:/My Drive/MISC/MASTER/2023-24/Umich/Winter 25/URP 535/Week 4/Lab4")

#install.packages('tidyverse')
library(tidyverse)
## Warning: package 'readr' was built under R version 4.4.3
## Warning: package 'dplyr' was built under R version 4.4.3
#install.packages('tmap') 
library(tmap)
#install.packages('sf') 
library(sf)
#install.packages('tmap') 
library(osmdata)
#install.packages('tmap') 
library(osmdata)
#install.packages('leaflet') 
library(leaflet) # for visualizing leaflet maps
#install.packages("shiny")
library(shiny) # for writing a shiny app
## Warning: package 'shiny' was built under R version 4.4.3
# install.packages("rsconnect")  
library(rsconnect)# for deploying a shiny app

library(dplyr)

Q2: Implement a Choropleth

Explore leaflet documentation - Choropleths page. Use MI_counties.shp in Lab 4 and replicate a similar Choropleth as in the tutorial. The following features should present, while the rest is up to you:

# import shapefile into R
setwd("G:/My Drive/MISC/MASTER/2023-24/Umich/Winter 25/URP 535/Week 4/Lab4")
data <- read.csv('sample.csv')
cnty_shp <- st_read('MI_counties.shp')
## Reading layer `MI_counties' from data source 
##   `G:\My Drive\MISC\MASTER\2023-24\Umich\Winter 25\URP 535\Week 4\Lab4\MI_counties.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 83 features and 17 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -90.41839 ymin: 41.69612 xmax: -82.12297 ymax: 48.30606
## Geodetic CRS:  NAD83
cnty_shp <- st_transform(cnty_shp, crs = 4326)
cnty_shp %>% mutate(proof_of_concept = "you can use tidyverse!") %>% head()
## Simple feature collection with 6 features and 18 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -88.11761 ymin: 43.81368 xmax: -82.40872 ymax: 47.22327
## Geodetic CRS:  WGS 84
##   STATEFP COUNTYFP COUNTYNS GEOID        NAME           NAMELSAD LSAD CLASSFP
## 1      26      109 01622997 26109   Menominee   Menominee County   06      H1
## 2      26      007 01622946 26007      Alpena      Alpena County   06      H1
## 3      26      035 01622960 26035       Clare       Clare County   06      H1
## 4      26      103 01622994 26103   Marquette   Marquette County   06      H1
## 5      26      089 01622987 26089    Leelanau    Leelanau County   06      H1
## 6      26      153 01623017 26153 Schoolcraft Schoolcraft County   06      H1
##   MTFCC CSAFP CBSAFP METDIVFP FUNCSTAT      ALAND     AWATER    INTPTLAT
## 1 G4020   361  31940     <NA>        A 2704053268  761798738 +45.5362068
## 2 G4020  <NA>  10980     <NA>        A 1481127630 2908710755 +44.8949540
## 3 G4020  <NA>   <NA>     <NA>        A 1461700232   28500461 +43.9911368
## 4 G4020  <NA>  32100     <NA>        A 4685406068 4185352679 +46.6565964
## 5 G4020  <NA>  45900     <NA>        A  899241895 5659105307 +45.1461816
## 6 G4020  <NA>   <NA>     <NA>        A 3035145674 1844210517 +46.0211217
##       INTPTLON                       geometry       proof_of_concept
## 1 -087.5018115 POLYGON ((-87.49251 45.9855... you can use tidyverse!
## 2 -083.4265739 POLYGON ((-83.75924 45.2047... you can use tidyverse!
## 3 -084.8383253 POLYGON ((-85.0876 44.07366... you can use tidyverse!
## 4 -087.5840278 POLYGON ((-88.1166 46.50652... you can use tidyverse!
## 5 -086.0515740 POLYGON ((-85.56175 44.9522... you can use tidyverse!
## 6 -086.1970290 POLYGON ((-86.24015 46.5052... you can use tidyverse!
# only show the first six rows 
cnty_shp <- cnty_shp %>% 
  # only keep two columns and drop other irrelevant columns 
  select(c(GEOID, geometry)) %>% 
  # convert GEOID column to integer to match with GEOID type in data
  mutate(GEOID = as.integer(GEOID)) %>% 
  # join data to cnty_shape. Note that due to tidyverse syntax, the first argument in left_join is actually cnty_shp rather than data. copy=FALSE removes duplicated columns on the unique identifier. 
  left_join(data, by=c('GEOID'), copy=FALSE)
pal <- colorNumeric(
  palette = "YlOrRd",  
  domain = cnty_shp$population  
)

labels <- sprintf(
  "<strong>%s</strong><br/>%g people / mi<sup>2</sup>",
  cnty_shp$county, cnty_shp$population
) %>% lapply(htmltools::HTML)


leaflet(cnty_shp) %>%
  setView(-86, 45, 6) %>%
  addProviderTiles("MapBox", options = providerTileOptions(
    id = "mapbox.light",
    accessToken = Sys.getenv('MAPBOX_ACCESS_TOKEN'))) %>%
  addPolygons(
    fillColor = ~pal(population),
    weight = 2,
    opacity = 1,
    color = "white",
    dashArray = "3",
    fillOpacity = 0.7,
    highlightOptions = highlightOptions(
      weight = 5,
      color = "#666",
      dashArray = "",
      fillOpacity = 0.7,
      bringToFront = TRUE),
    label = labels,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "15px",
      direction = "auto")) %>%
  addLegend(pal = pal, values = ~population, opacity = 0.7, title = "Michigan Population",
    position = "bottomright")