Geography 176A

# SPDS
library(tidyverse)
library(sf)
library(units)

# Data
library(USAboundaries)
library(rnaturalearth)

# Visualization
library(gghighlight)
library(ggrepel)
library(knitr)
library(rmapshaper)
library(maps)
library(leaflet)
library(mapview)

Question 1

#1.1
conus = USAboundaries::us_counties() %>%
  filter(!state_name %in% c("Puerto Rico",
                            "Alaska",
                            "Hawaii")) %>% 
  st_transform(5070)
st_geometry(conus)

#1.2
mp_counties= conus %>% 
  group_by(state_name) %>% 
  summarise() 

c_counties= mp_counties %>% 
  st_centroid() %>% 
  st_union()

#1.3
v_grid = st_voronoi(mp_counties) %>% 
  st_cast() %>% 
  st_as_sf() %>% 
  mutate(id = 1:n())

t_grid = st_triangulate(mp_counties) %>% 
  st_cast() %>% 
  st_as_sf() %>% 
  mutate(id = 1:n())

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

hex_grid = st_make_grid(c_counties, n = 70, square = FALSE) %>% 
  st_as_sf() %>% 
  mutate(id = 1:n())

#1.4 & 1.5
boundary = conus %>% 
  st_union() %>% 
  st_cast("MULTILINESTRING") %>% 
  ms_simplify(keep = .0000000000001)
print(mapview::npts(boundary))
plot(boundary)

#1.6
v_grid = st_intersection(v_grid, boundary)

t_grid = st_intersection(t_grid, boundary)

#1.7
plot_tess = function(data, title){
  ggplot() + 
    geom_sf(data = data, fill = "white", col = "navy", size = .2) +   
    theme_void() +
    labs(title = title, caption = paste("This tesselation has:", nrow(data), "tiles" )) +
    theme(plot.title = element_text(hjust = .5, color =  "navy", face = "bold"))
}
plot_tess(hex_grid, "hex")
plot_tess(v_grid, "voronoi")
plot_tess(t_grid, triangulate)
plot_tess(sq_grid, "square")

## Error: <text>:33:1: unexpected symbol
## 32: 
## 33: hex_grid
##     ^

Question 2

#2.1
sum_tess = function(data, var) {
 data = data %>% 
    mutate(area = data %>% 
             st_area() %>% 
             set_units("km^2") %>% 
             drop_units())
  var <- data.frame("tesselation" = var, "features" = count(data$id), "mean_area" = mean(data$area), "std" = sd(data$area, na.rm = FALSE), "total_area" = sum(data$area))
  return(var)
  }

#2.2
sum_tess(hex_grid, "hex")

## Error in eval(lhs, parent, parent): object 'hex_grid' not found

sum_tess(v_grid, "voronoi")

## Error in eval(lhs, parent, parent): object 'v_grid' not found

sum_tess(t_grid, "triangulate")

## Error in eval(lhs, parent, parent): object 't_grid' not found

sum_tess(sq_grid, "square")

## Error in eval(lhs, parent, parent): object 'sq_grid' not found

#2.3
tess_summary = bind_rows(
  sum_tess(hex_grid, "hex"),
  sum_tess(v_grid, "voronoi"),
  sum_tess(t_grid, "triangulate"),
  sum_tess(sq_grid, "square"),
  sum_tess(mp_counties, "raw"))

## Error in eval(lhs, parent, parent): object 'hex_grid' not found

#2.4
knitr::kable(tess_summary, 
             caption = "Summary of Tesellations",
             col.names = c("Tesselation", "# of Features", "Mean Area (km^2", "Standard Deviation", "Total Area (km^2"),
             format.args = list(big.mark = ","))

## Error in knitr::kable(tess_summary, caption = "Summary of Tesellations", : object 'tess_summary' not found

#2.5
#The voronoi tessellation returns back an area much bigger than the raw data, potentially leading to more calculations. The Triangulate and raw data are quite similar so triangulation is likely a good way of representing data. Square and hex have the least area and so are likely easiest to calculate.

Question 3

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

sf_dams = dams[!is.na(dams$LONGITUDE)&!is.na(dams$LATITUDE),] 

sf_dams = st_as_sf(sf_dams,
           coords=c("LONGITUDE", "LATITUDE"),
           crs= 5070) 
st_geometry(sf_dams)

## Geometry set for 91223 features 
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -176.6873 ymin: 0 xmax: 144.7063 ymax: 71.26667
## projected CRS:  NAD83 / Conus Albers
## First 5 geometries:

#3.2
point_in_polygon3 = function(points, polygon, group){
  st_join(polygon, points) %>%
    st_drop_geometry() %>%
    count(group) %>%
    setNames(c(group, "n")) %>%
    left_join(polygon, by = group) %>%
    st_as_sf()
}


#3.3
hex_pip = point_in_polygon3(sf_dams, hex_grid, id)

## Error in st_join(polygon, points): object 'hex_grid' not found

v_pip = point_in_polygon3(sf_dams, v_grid, id)

## Error in st_join(polygon, points): object 'v_grid' not found

sq_pip = point_in_polygon3(sf_dams, sq_grid, id)

## Error in st_join(polygon, points): object 'sq_grid' not found

t_pip = point_in_polygon3(sf_dams, t_grid, id)

## Error in st_join(polygon, points): object 't_grid' not found

raw_pip = point_in_polygon3(sf_dams, mp_counties, RECORDID)

## Error in st_join(polygon, points): object 'mp_counties' not found

#3.4
plot_pip = function(data){
  ggplot() + 
    geom_sf(data = data, aes(fill = count(n)), alpha = .9, size = .2) + 
    scale_fill_gradient("viridis") + 
    theme_void() + 
    theme(legend.position = 'none',
          caption = paste0(sum(data$n), "total dams")) 
}

#3.5
plot_pip(hex_pip)

## Error in fortify(data): object 'hex_pip' not found

plot_pip(v_pip)

## Error in fortify(data): object 'v_pip' not found

plot_pip(t_pip)

## Error in fortify(data): object 't_pip' not found

plot_pip(sq_pip)

## Error in fortify(data): object 'sq_pip' not found

plot_pip(raw_pip)

## Error in fortify(data): object 'raw_pip' not found

#3.6
#I was unable to get my code to work, however I would likely choose the hex tesselation because it is aesthetically pleasing and relatively quick to calculate.

Question 4

#4.1
#I will choose recreation, water supply, fish and wildlife, and irrigation because these are all human needs that can potentially be in conflict with eachother.
dams_R = sf_dams %>% 
  filter(grepl("R", sf_dams$PURPOSES) == "TRUE")
dams_S = sf_dams %>% 
  filter(grepl("S", sf_dams$PURPOSES) == "TRUE")
dams_F = sf_dams %>% 
  filter(grepl("F", sf_dams$PURPOSES) == "TRUE")
dams_I = sf_dams %>% 
  filter(grepl("I", sf_dams$PURPOSES) == "TRUE")

#4.2
plot_pip2 = function(data){
  ggplot() + 
    geom_sf(data = data, aes(fill = count(n)), alpha = .9, size = .2) + 
    scale_fill_gradient("viridis") + 
    theme_void() + 
    theme(legend.position = 'none',
          caption = paste0(sum(data$n), "total dams")) +
    gghighlight(count(n)> (mean(data)+sd(data)))
}

plot_pip2(dams_R)

## Error: All calculations failed! Please provide a valid predicate.

plot_pip2(dams_S)

## Error: All calculations failed! Please provide a valid predicate.

plot_pip2(dams_F)

## Error: All calculations failed! Please provide a valid predicate.

plot_pip2(dams_I)

## Error: All calculations failed! Please provide a valid predicate.

#4.3
#The dams would likely be distributed around large water systems such as the mississippi river system. lab-05.rmd

Geography 176A

Lab 04: Tesselations, Point-in-Polygon

Kyler Plouffe

Question 1

Question 2

Question 3

Question 4