Geospatial Analysis with R

library(tidyverse)
library(rvest)

URL <- "https://history.house.gov/Institution/Party-Divisions/Party-Divisions/"
webpage <- read_html(URL)
cong <- as_tibble(html_table(webpage)[[1]])
onms <- colnames(cong)
newnms <- c("congress", "seats", "D", "R", "other", "delres")
cong <- cong %>% 
  slice((which(grepl("^Republican", `Anti-Administration`))[1] + 1):nrow(.)) %>%
  filter(`Congress (Years)` != colnames(cong)[1]) %>% 
  mutate(year = gsub("(*.*-)|(*.*–)|)|)2", "", `Congress (Years)`)) %>%
  mutate(year = as.numeric(year) - 2) %>% 
  rename_at(vars(onms), ~newnms) %>% select(-other, -delres) %>%
  mutate(seats = substr(seats, 1, 3)) %>% 
  mutate_at(.vars = vars(seats, D, R), as.numeric) %>% 
  mutate(swing = (D - R) / (D + R)) 
cong %>%  
  ggplot() + geom_line(aes(year, swing)) + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) + 
  geom_hline(yintercept = 0, lty = 2, col = "red") + 
  xlab("") + ylab("(D - R) / (D + R)") + 
  ggtitle("Normalized Party Control Index for US House") + theme_linedraw() + 
  scale_x_continuous(breaks = seq(1859, 2019, 10), expand = c(0, 2))

library(sf)
library(dplyr)
sz <- 100
dat <- lapply(1:3, function(x) {
  set.seed(x)
  d <- data.frame(
    id = 1:sz,
    v1 = runif(sz, min = 2, max = 12),
    grp = paste0("g", x)
  ) %>% mutate(v2 = v1 + rnorm(sz, mean = 2, sd = 2)) %>% 
    select(id, v1, v2, grp)
}) %>% bind_rows()

plot(dat$v1, dat$v2)

dat %>% 
  group_by(grp) %>% 
  # summarize(v2 = mean(v2)) %>%
  # summarize(across(v1:v2, mean)) %>%
  summarize(across(v1:v2, list("mu" = mean, "sd" = sd)))

t(sapply(unique(dat$grp), function(x) {
  colMeans(dat[dat$grp == x, c("v1", "v2")])
}))

ggplot(dat) + 
  geom_point(aes(x = v1, y = v2, color = grp))

ggplot(dat) + 
  geom_point(aes(x = v1, y = v2)) + 
  facet_wrap(~grp)
  # facet_wrap(~grp, nrow = 2, ncol = 2)

par(mfrow = c(1, 3))
for(i in unique(dat$grp)) {
  plot(v2 ~ v1, data = dat[dat$grp == i, ])
}

reg <- lm(v2 ~ v1, data = dat)
summary(reg)

plot(dat$v1, dat$v2)
abline(reg, col = "red")

ggplot(dat) + 
  geom_point(aes(x = v1, y = v2)) + 
  geom_smooth(aes(x = v1, y = v2), method = "lm")

# ? thoughts
ggplot(dat) + 
  geom_point(aes(x = v1, y = v2, color = grp)) + 
  geom_smooth(aes(x = v1, y = v2), method = "lm") + 
  # geom_smooth(aes(x = v1, y = v2, color = grp), method = "lm") #+
  facet_wrap(~grp)

library(tidyverse)
library(sf)
## read in csv
farmers <- system.file("extdata/farmer_spatial.csv", package = "geospaar") %>% 
  read_csv() 

## convert tibble to 'sf'
farmers_sf <- st_as_sf(farmers, coords = c("x", "y"), crs = 4326)

districts <- read_sf(system.file("extdata/districts.geojson", 
                                 package = "geospaar"))
roads <- st_read(system.file("extdata/roads.geojson", package = "geospaar"))

#st_write(districts, "districts_out.shp")
st_write(districts, "districts2.geojson")

st_crs(farmers_sf) <- 4326
farmers_sf %>% st_crs()

# or
farmers_sf %>% st_set_crs(NA)

farmers_sf %>% st_set_crs(4326)

farmers <- read_csv(system.file("extdata/farmer_spatial.csv", 
                                package = "geospaar"))
farmers <- st_as_sf(farmers, coords = c("x", "y")) %>% 
  st_set_crs(st_crs(districts))

farmers %>% st_crs() ## print new crs

roads <- st_transform(x = roads, crs = st_crs(districts))
st_crs(roads)

pt <- st_point(x = c(28, -14))
pt
#> POINT (28 -14)
#
# #2
pts <- st_multipoint(x = cbind(x = c(27.5, 28, 28.5), y = c(-14.5, -15, -15.5)))
pts
#> MULTIPOINT ((27.5 -14.5), (28 -15), (28.5 -15.5))
#
# #3
sline <- st_linestring(cbind(x = c(27, 27.5, 28), y = c(-15, -15.5, -16)))
sline
#> LINESTRING (27 -15, 27.5 -15.5, 28 -16)
#
# #4
pol <- st_polygon(list(cbind(x = c(26.5, 27.5, 27, 26, 26.5), 
                             y = c(-15.5, -16.5, -17, -16, -15.5))))
pol

par(mar = c(0, 0, 0, 0))
plot(districts %>% st_geometry(), col = "grey")
plot(pt, add = TRUE, col = "red", pch = 20)
plot(pts, add = TRUE, col = "blue", pch = 20)
plot(sline, add = TRUE, col = "cyan", lwd = 2)
plot(pol, add = TRUE, col = "orange", lwd = 2)

ggplot(districts) + geom_sf() +
  geom_sf(data = roads, col = "red") + 
  geom_sf(data = farmers_sf, col = "blue")

districts <- st_read(system.file("extdata/districts.geojson", 
                                 package = "geospaar"))
farmers_sf <- read_csv(system.file("extdata/farmer_spatial.csv", 
                                   package = "geospaar"))
farmers_s <- st_as_sf(farmers, coords = c("x", "y")) %>% 
  st_set_crs(st_crs(districts))

roads <- st_read(system.file("extdata/roads.geojson", package = "geospaar")) %>% 
  st_transform(x = roads, crs = st_crs(districts))

pt <- st_point(x = c(28, -14))
pts <- st_multipoint(x = cbind(x = c(27.5, 28, 28.5), y = c(-14.5, -15, -15.5)))
sline <- st_linestring(cbind(x = c(27, 27.5, 28), y = c(-15, -15.5, -16)))
pol <- st_polygon(list(cbind(x = c(26.5, 27.5, 27, 26, 26.5), 
                             y = c(-15.5, -16.5, -17, -16, -15.5))))

ggplot(districts) + geom_sf() +
  geom_sf(data = roads, col = "red") + 
  geom_sf(data = farmers_sf, col = "blue")

dist_areas <- districts %>% st_area()
dist_areas
dist_areas %>% units::set_units("km^2")

districts <- districts %>% 
  mutate(area = st_area(.) %>% units::set_units("km^2") %>% as.numeric()) 
districts

ggplot(districts) + geom_sf(aes(fill = area)) + 
  scale_fill_continuous(type = "viridis")

districts <- districts %>% 
  st_transform("ESRI:102022") %>% 
  mutate(area2 = as.numeric(st_area(.) / 10000)) %>% 
  st_transform(4326)

ggplot(districts) + 
  geom_histogram(aes(x = area - area2))