#1.1 Spatial Conus file
conus = USAboundaries::us_counties() %>% 
  filter(!state_name %in% c("Hawaii", "Puerto Rico", "Alaska", "Guam", "District of Columbia")) %>% 
  st_transform(5070)

#1.2 Anchoring
get_conus = function(data, var){
  filter(data, !get(var) %in%
           c("Hawaii", "Puerto Rico", "Alaska", "Guam", "District of Columbia"))
}

counties = us_counties() %>% 
  get_conus("state_name") %>% 
  st_cast() %>% 
  st_as_sf() %>% 
  mutate(id = 1:n())  %>% 
  st_transform(5070)

conus_cent = st_centroid(counties) %>%
  st_union() 
#1.3 Tesselations
counties_1 = conus_cent %>% 
  st_cast() %>% 
  st_as_sf() %>% 
  mutate(id = 1:n()) 

conus_voroni = conus_cent %>% 
  st_voronoi() %>% 
  st_cast() %>% 
  st_as_sf() %>% 
  mutate(id = 1:n()) 
conus_tri = conus_cent %>% 
  st_triangulate() %>% 
  st_cast() %>% 
  st_as_sf() %>% 
  mutate(id = 1:n()) 
#1.3 Coverages
conus_grd = counties %>% 
  st_make_grid(n = c(70,70)) %>% 
  st_as_sf() %>% 
  mutate(id = 1:n()) 

conus_hex = counties %>% 
  st_make_grid(n = c(70,70), square = FALSE) %>% 
  st_as_sf() %>% 
  mutate(id = 1:n())
counties_union = counties %>% 
  st_union()

conus_voroni = conus_voroni %>% 
  st_intersection(counties_union)

conus_tri = conus_tri %>% 
  st_intersection(counties_union)

conus_grd = conus_grd %>% 
  st_intersection(counties_union)

conus_hex = conus_hex %>% 
  st_intersection(counties_union)
#1.4 Cutting Tesselations
resolved_counties = counties_union %>% 
  ms_simplify(keep = 0.05)

conus_voroni = conus_voroni %>% 
  st_intersection(resolved_counties)

conus_tri = conus_tri %>% 
   st_intersection(resolved_counties)
#1.6 Tessellation Function
plot_tess = function(arg1, title){
  ggplot()+
    geom_sf(data = arg1, col = "white", size = 0.2)+
    theme_void()+
    theme(legend.position = "none",
          plot.title = element_text(face = "bold", color = "darkblue", hjust = 0.5, size = 24))+
    labs(title = paste0(title),
         caption = paste0(mapview::npts(arg1), "features"))
}
#1.7 Plotting Tessellations
plot_tess(counties_union, "Original")

plot_tess(conus_voroni, "Voroni")

plot_tess(conus_tri, "Triangulation")

plot_tess(conus_hex, "Hexagonal")

#2.1 Summary Function
tess_summary = function(sf_object, descrip){
  object_area = st_area(sf_object) %>% 
    set_units("km^2") %>% 
    drop_units()
  area_df= data.frame(tesselation = descrip, features = max(sf_object$id), mean_area = mean(object_area), std_area = sd(object_area), tot_area = sum(object_area))
return(area_df)
}
#2.2 Summarizing Each Tesselation
Voroni_dataframe = tess_summary(conus_voroni, "Voroni Tessellation")

Triangulation_datafram = tess_summary(conus_tri, "Triangulation Tessellation")

Grid_dataframe = tess_summary(conus_grd, "Grid Cover")

Hexgrid_datafram = tess_summary(conus_hex, "Hexagonal Cover")

Original_dataframe = tess_summary(counties_union, "No Tesselation")
#2.3 Binding rows
tess_summary2 = bind_rows(
  tess_summary(conus_tri, 'Triangulation Tesselation'),
  tess_summary(conus_voroni, 'Voroni Tesselation'),
  tess_summary(conus_grd, 'Grid Cover'),
  tess_summary(conus_hex, 'Hexagonal Cover'),
  tess_summary(counties, 'County')
  
)
#2.4 Printing Table
knitr::kable(tess_summary2, caption = 'US Counties Tesselation', col.names = c('Tesselation', '# of Features', 'Mean Area', 'Standard Deviation of Area', 'Total Area'))
US Counties Tesselation
Tesselation # of Features Mean Area Standard Deviation of Area Total Area
Triangulation Tesselation 6194 1247.199 1574.5349 7723900
Voroni Tesselation 3107 2510.354 2883.6265 7794648
Grid Cover 3108 2521.687 606.5924 7837404
Hexagonal Cover 2271 3451.080 870.9742 7837404
County 3107 2522.499 3404.6136 7837404

#2.5 The triangulation tessellation has the greatest number of features, whereas the Hexagonal cover has the greatest average area per tesselation. Although the hexagonal cover has the greatest average area, the County tessellation has the greatest standard deviation. Lastly, the Voroni tesselation has the largest total area.

#3.1 loading dams dataset
dams <- read_excel("~/github/geog-176A-labs/data/NID2019_U.xlsx")

dams_sf = dams %>% 
  filter(!is.na(LONGITUDE), !is.na(LATITUDE)) %>% 
  st_as_sf(coords = c('LONGITUDE', 'LATITUDE'), crs = 4326) %>% 
  st_transform(5070)
#3.2 Creating point in polygon function
pip = function(points, polygon, bar){
  st_join(polygon, points) %>% 
    st_drop_geometry() %>% 
    count(get(bar)) %>% 
    setNames(c(bar, "n")) %>% 
    left_join(polygon, by = bar) %>% 
    st_as_sf()
}
plot(conus_voroni)

#3.3 applying point-in-polygon function to each tessellation
dams_voroni = pip(dams_sf, conus_voroni, 'id')
dams_tri = pip(dams_sf, conus_tri, 'id')
dams_grd = pip(dams_sf, conus_grd, 'id')
dams_hex = pip(dams_sf, conus_hex, 'id')
dams_county = pip(dams_sf, counties, 'id')
#3.4 extending pip function to a new plotting function
pip_plot = function(data = data, text){
  ggplot()+
    viridis::scale_fill_viridis()+
    theme_void()+
    theme(plot.title = element_text(face = 'bold', color = 'black', size = 24), plot.subtitle = element_text(size = 13),
          plot.caption = element_text(face = 'bold', size = 13), legend.title = element_text('bold'),
          legend.text = element_text(face = 'bold'))+
    labs(title = text,
         subtitle = 'National Inventory of Dams',
         fill = '# of Dams',
         caption = paste0(sum(data$n), "dams"))+
    theme(aspect.ratio = 5)
}
#3.5 applying new plotting function to each of the tessellations 

pip_plot(dams_voroni, "Voroni Tesselation of US Dams")

pip_plot(dams_tri, "Triangulation Tesselation of US Dams")

pip_plot(dams_grd, "Grid Tesselation of US Dams")

pip_plot(dams_hex, "Hexagonal Tesselation of US Dams")

pip_plot(dams_county, "County Lines of US Dams")

#3.6 Going forward we will only use the Voroni Tessellation because it provides us with the most information.

#4.1 Creating a 
nid_classifier = data.frame(abbr = c("I", "H", "C", "N", "S", "R", "P", "F", "D", "T", "G", "O"),
                            purpose = c("Irrigation", "Hydroelectric", "Flood Control", "Navigation", "Water Supply", "Recreation", "Fire Protection", "Fish and Wildlife", "Debris Control", "Tailings", "Grade Stabilization", "Other"))

Dam_Frq = strsplit(dams$PURPOSES, split = "") %>% 
  unlist() %>% 
  table() %>% 
  as.data.frame() %>% 
  setNames(c("abbr", "count")) %>% 
  left_join(nid_classifier) %>% 
  mutate(Lab = paste0(purpose, "\n(", abbr, ")"))

recdams = pip(dams_sf[grepl("R", dams_sf$PURPOSES),], dams_voroni, 'id')
flood_control = pip(dams_sf[grepl("C", dams_sf$PURPOSES),], dams_voroni, 'id')
fire_dams = pip(dams_sf[grepl("F", dams_sf$PURPOSES),], dams_voroni, 'id')
water_supply = pip(dams_sf[grepl("S", dams_sf$PURPOSES),], dams_voroni, 'id')
irrigation = pip(dams_sf[grepl("I", dams_sf$PURPOSES),], dams_voroni, 'id')
hydroelectric =pip(dams_sf[grepl("H", dams_sf$PURPOSES),], dams_voroni, 'id')
#4.1 Creating a new point-in-polygon function to count the number of dams per purpose.
pip_plot2 = function(data, text){
  ggplot()+
    geom_sf(data = data, aes(fill = n), alpha = 5, size = 2)+
    gghighlight(n > mean(n) + sd(n))+
    viridis::scale_fill_viridis()+
    theme_void()+
    theme(plot.title = element_text(face = "bold", color = "black", size = 24), plot.subtitle = element_text(size = 11), legend.text = element_text(face = "bold"))+
    labs(title = text,
         subtitle = 'National Inventory of Dams',
         fill = '# of Dams', caption = paste0(sum(data$n), "dams"))+
    theme(aspect.ratio = .5)
}
#4.2 using the plotting function from Q3 to map the counts of dams 

#recdams_plot = pip_plot2(recdams, "Recreational Dams")
#plot(recdams_plot)

#error: Problem with `mutate()` input `fill`. x Input `fill` must be a vector, not a function. ℹ Input `fill` is `tryCatch(n, error = function(e) NA)`.
#Extra Credit


library(leaflet)



rivers = read_sf("~/github/geog-176A-labs/data/majorrivers_0_0")

rivers = rivers %>% 
  filter(SYSTEM == "Mississippi")

biggest_dams = dams_sf %>% 
  filter(HAZARD == "H", grepl("C", PURPOSES)) %>% 
  group_by(STATE) %>% 
  slice_max("NID_STORAGE") %>% 
  select("DAM_NAME", "NID_STORAGE", "PURPOSES", "YEAR_COMPLETED")

#Leaflet

leaflet() %>% 
  addProviderTiles(providers$CartoDB) %>% 
  addPolylines(data = rivers) %>% 
  addCircleMarkers(data = st_transform(biggest_dams, 4326), fillOpacity =1, radius = ~NID_STORAGE/1500000, color = "red", stroke = FALSE, popup = leafpop::popupTable(st_drop_geometry(biggest_dams), feature.id = FALSE, row.numbers = FALSE))