1.1 Spatial file of US Counties

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

get_conus= function(data, var){
  filter(data, !get(var) %in%
           c('Hawaii', 'Puerto Rico', 'Guam', 
             'District of Columbia', 'Alaska'))
}
 
counties = st_transform(us_counties(), 5070) %>% 
  select(name, geoid, state_name) %>% 
  get_conus('state_name')

1.2 Anchors

county_centroid = st_centroid(counties) %>% 
  st_combine() %>% 
  st_cast("MULTIPOINT")

1.3 Tessellations/Coverages

vor = st_voronoi(county_centroid) %>% 
  st_cast() %>% 
  st_as_sf() %>% 
  mutate(id = 1:n())

triangle = st_triangulate(county_centroid) %>% 
  st_cast() %>% 
  st_as_sf() %>% 
  mutate(id = 1:n())

square = st_make_grid(county_centroid, n = c(70, 50)) %>% 
  st_as_sf() %>% 
  mutate(id = 1:n())

hexagon = st_make_grid(county_centroid, n = c(70, 50), square = FALSE) %>% 
  st_as_sf() %>% 
  mutate(id = 1:n())
county_centroid_union = st_union(county_centroid)

counties_u = counties %>% 
  st_union()

vor = vor %>% 
  st_intersection(counties_u)

triangle = triangle %>% 
  st_intersection(counties_u)

square = square %>% 
  st_intersection(counties_u)

hexagon = hexagon %>% 
  st_intersection(counties_u)

1.4

plot_tess = function(data, title){
  ggplot() + 
    geom_sf(data = data, fill = "white", col = "red", size = .2) +   
    theme_void() +
    labs(title = title, caption = paste("# of Tesselations:", nrow(data), "tiles" )) +
    theme(plot.title = element_text(hjust = .5, color =  "red", face = "bold"))
}
counties_resolved = counties %>% 
  ms_simplify(keep = 0.05)

Plots

plot_tess(counties, "Counties") #original

plot_tess(vor, "Voronoi") #voronoi

plot_tess(triangle, "Triangulated") #Triangulated

plot_tess(hexagon, "Hexagonial") #Hexagonal

plot_tess(square, "Square") #Square

2.1

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)
}
tess_summary(counties, "Original, No Tesselation")
##                tesselation features mean_area std_area tot_area
## 1 Original, No Tesselation     -Inf  2522.499 3404.614  7837404
tess_summary(vor, "Voroni Tesselation")
##          tesselation features mean_area std_area tot_area
## 1 Voroni Tesselation     3107  2522.499 2887.307  7837404
tess_summary(triangle, 'Triangular Tesselation')
##              tesselation features mean_area std_area tot_area
## 1 Triangular Tesselation     6194  1252.183 1575.912  7756021
tess_summary(hexagon, 'Hexagonal Cover')
##       tesselation features mean_area std_area tot_area
## 1 Hexagonal Cover     1639  3494.129 382.4741  5726877
tess_summary(square, "Grid Cover")
##   tesselation features mean_area std_area tot_area
## 1  Grid Cover     1648  3477.045 415.1009  5730171

2.3

tess_summary_bound = bind_rows(
  tess_summary(counties, "Original, No Tesselation"),
tess_summary(vor, "Voroni Tesselation"),
tess_summary(triangle, 'Triangular Tesselation'),
tess_summary(hexagon, 'Hexagonal Cover'),
tess_summary(square, "Grid Cover"),
  
)

2.4

#2.4
knitr::kable(tess_summary_bound,
caption = "Tesselated Surfaces",
col.names = c('Type', 'Features', 'Mean Area', 'SD of Area', 'Total Area'),
format.args = list(big.mark = ",")) %>%
kableExtra::kable_styling("basic", full_width = TRUE, font_size = 16)
Tesselated Surfaces
Type Features Mean Area SD of Area Total Area
Original, No Tesselation -Inf 2,522.499 3,404.6136 7,837,404
Voroni Tesselation 3,107 2,522.499 2,887.3075 7,837,404
Triangular Tesselation 6,194 1,252.183 1,575.9119 7,756,021
Hexagonal Cover 1,639 3,494.129 382.4741 5,726,877
Grid Cover 1,648 3,477.045 415.1009 5,730,171

3.1

library(readxl)
dam_data <- read_excel("/Users/lfinn443/github/geog-176A-labs/data/NID2019_U.xlsx") 
  

dam_2019 <- dam_data %>% 
  filter(!is.na(LONGITUDE), !is.na(LATITUDE)) %>% 
  st_as_sf(coords = c('LONGITUDE', 'LATITUDE'), crs = 4326) %>% 
  st_transform(5070)

3.2

point_in_polygon = 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()
}

3.3

county_dams = point_in_polygon(dam_2019, counties, 'geoid')
vor_dams = point_in_polygon(dam_2019, vor, 'id')
tri_dams = point_in_polygon(dam_2019, triangle, 'id')
sq_dams = point_in_polygon(dam_2019, square, 'id')
hex_dams = point_in_polygon(dam_2019, hexagon, 'id')

3.4

dam_plot= function(data=data, text)
  {ggplot() +
    geom_sf(data = data, aes(fill = log(n)), alpha = .8, size = .2, col = NA) +
    scale_fill_viridis_c() +
    theme_void() +
    theme(legend.position = 'none',
          plot.title = element_text(face = "bold", color = "black",  size = 24)) +
    labs(title = text,
         caption = paste0(sum(data$n), " Dams "))}

3.5

dam_plot(vor_dams, "US Dams: Voroni Tesselation")

dam_plot(tri_dams, "US Dams: Triangulation Tesselation")

dam_plot(sq_dams, "US Dams: Sqaure Cover")

dam_plot(hex_dams, "US Dams: Hexagonal Cover")

dam_plot(county_dams, "US Dams: County lines")

dam_plot_vor = dam_plot(vor_dams, "US Dams: Voroni Tesselation")

dam_plot_tri = dam_plot(tri_dams, "US Dams: Triangulation Tesselation")

dam_plot_sq = dam_plot(sq_dams, "US Dams: Grid Cover")

dam_plot_hex = dam_plot(hex_dams, "US Dams: Hexagonal Cover")

dam_plot_county = dam_plot(county_dams, "US Dams: County lines")