#Libraries
library(tidyverse)
library(sf)
library(USAboundaries)
library(ggplot2)
library(gghighlight)
library(leaflet)
library(knitr)
library(readxl)
library(rmapshaper)
library(units)
library(dplyr)
library(leafpop)#Question 1
counties <- USAboundaries::us_counties() %>%
filter(!state_name %in% c("Hawaii", "Puerto Rico", "Alaska", "Guam")) %>%
st_transform(5070) %>%
st_as_sf() %>%
mutate(id = 1:n())
centroids <- counties %>%
st_centroid()
nrow(centroids)## [1] 3108cent_union <- centroids %>%
st_union()
cent_union## Geometry set for 1 feature
## Geometry type: MULTIPOINT
## Dimension: XY
## Bounding box: xmin: -2303763 ymin: 344857.6 xmax: 2200721 ymax: 3127770
## Projected CRS: NAD83 / Conus Albersboundary <- counties %>%
st_union() %>%
ms_simplify(keep = .025)
mapview::npts## function (x, by_feature = FALSE)
## {
## if (by_feature)
## nVerts(sf::st_geometry(x))
## else sum(nVerts(sf::st_geometry(x)))
## }
## <bytecode: 0x00000000458f1768>
## <environment: namespace:mapview>voronois <- st_voronoi(cent_union) %>%
st_cast() %>%
st_as_sf() %>%
mutate(id = 1:n()) %>%
st_intersection(boundary)
triangulate = st_triangulate(cent_union) %>%
st_cast() %>%
st_as_sf() %>%
mutate(id = 1:n()) %>%
st_intersection(boundary)
sqgrid = st_make_grid(cent_union, n = c(70)) %>%
st_cast() %>%
st_as_sf() %>%
mutate(id = 1:n()) %>%
st_intersection(boundary)
hex_grid = st_make_grid(counties, n = c(70), square = FALSE) %>%
st_as_sf() %>%
mutate(id = 1:n()) %>%
st_intersection(boundary)
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(voronois, "Voroni Coverage")
plot_tess(triangulate, "Triangulated Coverage")
plot_tess(sqgrid, "Square Coverage")
plot_tess(hex_grid, "Hexagonal Coverage")
plot_tess(counties, "Original Counties")
#Question 2
tess = function(data, description){
area = st_area(data)
area = set_units(area, "km^2")
area = as.numeric(area)
data.frame(attributes = c("Description", "Num of Features", "Mean Area Features", "Standard Deviation", "Tot Area"), values = c(description, mapview::npts(data), sd(area), mean(area), sum(area)))
}
tess(triangulate, "Original Centroids")## attributes values
## 1 Description Original Centroids
## 2 Num of Features 25324
## 3 Mean Area Features 1589.20171189929
## 4 Standard Deviation 1253.46031447994
## 5 Tot Area 7752652.04505845tess(voronois, "Voronoi Coverage")## attributes values
## 1 Description Voronoi Coverage
## 2 Num of Features 21498
## 3 Mean Area Features 2900.62811539088
## 4 Standard Deviation 2530.56543089913
## 5 Tot Area 7839691.70492551tess(sqgrid, "Square Coverage")## attributes values
## 1 Description Square Coverage
## 2 Num of Features 16334
## 3 Mean Area Features 511.262610251616
## 4 Standard Deviation 2406.3736032659
## 5 Tot Area 7818307.8370109tess(hex_grid, "Hexagonal Coverage")## attributes values
## 1 Description Hexagonal Coverage
## 2 Num of Features 15414
## 3 Mean Area Features 783.453809458224
## 4 Standard Deviation 3509.2621776748
## 5 Tot Area 7839691.70492551tess_summary = bind_rows(
tess(triangulate ,"triangulation"),
tess(voronois, "voroni"))
cov_summary = bind_rows(
tess(sqgrid, "square grid"),
tess(hex_grid, "hexagonal grid"))
count_summmary = bind_rows(
tess(counties, "counties"))
knitr::kable(tess_summary, caption = "Bound Rows")| attributes | values |
|---|---|
| Description | triangulation |
| Num of Features | 25324 |
| Mean Area Features | 1589.20171189929 |
| Standard Deviation | 1253.46031447994 |
| Tot Area | 7752652.04505845 |
| Description | voroni |
| Num of Features | 21498 |
| Mean Area Features | 2900.62811539088 |
| Standard Deviation | 2530.56543089913 |
| Tot Area | 7839691.70492551 |
knitr::kable(cov_summary, caption = "Bound Rows")| attributes | values |
|---|---|
| Description | square grid |
| Num of Features | 16334 |
| Mean Area Features | 511.262610251616 |
| Standard Deviation | 2406.3736032659 |
| Tot Area | 7818307.8370109 |
| Description | hexagonal grid |
| Num of Features | 15414 |
| Mean Area Features | 783.453809458224 |
| Standard Deviation | 3509.2621776748 |
| Tot Area | 7839691.70492551 |
knitr::kable(count_summmary, caption = "Bound Rows")| attributes | values |
|---|---|
| Description | counties |
| Num of Features | 51976 |
| Mean Area Features | 3404.32522035577 |
| Standard Deviation | 2521.74470782184 |
| Tot Area | 7837582.55191029 |
#Question 3
NID <- read_excel("../data/NID2019_U.xlsx") %>%
filter(!is.na(LONGITUDE)) %>%
filter(!is.na(LATITUDE)) %>%
st_as_sf(coords = c("LONGITUDE", "LATITUDE"), crs = 4326) %>%
st_transform(5070)
PIP <- function(points, polygons, id){
st_join(polygons, points) %>%
dplyr::count(.data[[id]])
}
NIDV = PIP(NID, voronois, "id")
NIDT = PIP(NID, triangulate, "id")
NIDS = PIP(NID, sqgrid, "id")
NIDH = PIP(NID, hex_grid, "id")
NIDC = PIP(NID, counties, "id")
plot_PIP <- function(data, title){
ggplot() +
geom_sf(data = data, aes(fill = n), size = 0.2, col = NA) +
scale_fill_viridis_c() +
theme_void() +
labs(title = title, caption = sum(data$n))
}
plot_PIP(NIDV, "Voronoi Tessellations")
plot_PIP(NIDT, "Triangulation Tessellations")
plot_PIP(NIDS, "Square Grid Tessellations")
plot_PIP(NIDH, "Hexagonal Tesssellations")
plot_PIP(NIDC, "Counties")
# Further on, I will be using NIDV (the voroni tessellations), because I like the shape of themQuestion 4
# I picked hydroelectric, irrigation, flood control, and water supply
NIDHYDRO = NID %>%
dplyr::filter(grepl("H", NID$PURPOSES))
hydro_pip = PIP(NIDHYDRO, voronois, "id")
NIDIR = NID %>%
dplyr::filter(grepl("I", NID$PURPOSES))
ir_pip = PIP(NIDIR, voronois, "id")
NIDFC = NID %>%
dplyr::filter(grepl("C", NID$PURPOSES))
fc <- PIP(NIDFC, voronois, "id")
NIDWS = NID %>%
dplyr::filter(grepl("S", NID$PURPOSES))
ws <- PIP(NIDWS, voronois, "id")
plot_PIP(hydro_pip, "Hydroeletric Voroni Tessellation")
plot_PIP(ir_pip, "Irrigation Voroni Tessellation")
plot_PIP(fc, "Flood Control Voroni Tessellation")
plot_PIP(ws, "Water Supply Voroni Tessellation")
# Extra Credit
mississippi <- USAboundaries::us_states() %>%
filter(state_name %in% "Mississippi")
read_sf("../data/rivers")## Simple feature collection with 98 features and 4 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -164.8874 ymin: -36.96945 xmax: 160.7636 ymax: 71.39249
## Geodetic CRS: WGS 84
## # A tibble: 98 x 5
## NAME SYSTEM MILES KILOMETERS geometry
## <chr> <chr> <dbl> <dbl> <MULTILINESTRING [°]>
## 1 Kolyma <NA> 2552. 4106. ((144.8419 61.75915, 144.8258 61.8036, 145.~
## 2 Parana Parana 1616. 2601. ((-51.0064 -20.07941, -51.02972 -20.2225, -~
## 3 San Fra~ <NA> 1494. 2404. ((-46.43639 -20.25807, -46.49835 -20.25306,~
## 4 Japura Amazon 1223. 1968. ((-76.71056 1.624166, -76.70029 1.688336, -~
## 5 Putumayo Amazon 890. 1432. ((-76.86806 1.300553, -76.86695 1.295, -76.~
## 6 Rio Mar~ Amazon 889. 1431. ((-73.5079 -4.459834, -73.79197 -4.621942, ~
## 7 Ucayali Amazon 1298. 2089. ((-73.5079 -4.459834, -73.51585 -4.506949, ~
## 8 Guapore Amazon 394. 634. ((-65.39585 -10.39333, -65.39578 -10.39358,~
## 9 Madre d~ Amazon 568. 914. ((-65.39585 -10.39333, -65.45279 -10.40111,~
## 10 Amazon Amazon 1890. 3042. ((-73.5079 -4.459834, -73.45141 -4.427785, ~
## # ... with 88 more rowsdams = read_excel("../data/NID2019_U.xlsx") %>%
filter(!is.na(LONGITUDE), !is.na(LATITUDE)) %>%
st_as_sf(coords = c("LONGITUDE", "LATITUDE"), crs = 4326)
plot(dams)
largest <- dams %>%
filter(!STATE %in% c("Hawaii", "Puerto Rico", " Alaska", "Guam")) %>%
slice_max(NID_STORAGE, n = 1) %>%
dplyr::select("DAM_NAME", "NID_STORAGE", "PURPOSES", "YEAR_COMPLETED")
leaflet() %>%
addProviderTiles(providers$CartoDB) %>%
addPolylines(data = mississippi) %>%
addCircleMarkers(data = largest,
radius = ~NID_STORAGE / 1500000,
color = "red",
fillOpacity = 1,
stroke = FALSE,
popup = leafpop::popupTable(st_drop_geometry(largest[1:4]),
feature.id = FALSE, row.numbers = FALSE))