Geospatial Analysis with R

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

roads <- st_read(
  system.file("extdata/roads.geojson", package = "geospaar")
) 
roads <- roads %>% 
  mutate(road_length = as.numeric(units::set_units(st_length(.), "km"))) 
ggplot() +
  geom_sf(data = districts, color = "black", fill = NA) +
  geom_sf(data = roads, aes(color = road_length), lwd = 1) + 
  scale_color_continuous(type = "viridis") + 
  theme_bw()

ggplot(districts) + 
  geom_sf(color = "black", fill = NA) +
  st_read(system.file("extdata/roads.geojson", package = "geospaar")) %>% 
  mutate(road_length = as.numeric(units::set_units(st_length(.), "km"))) %>% 
  geom_sf(data = ., aes(color = road_length), lwd = 1) + 
  scale_color_continuous(type = "viridis") +
  theme_bw()

p <- ggplot(districts) + geom_sf()
p + roads %>% 
  mutate(road_length = as.numeric(units::set_units(st_length(.), "km"))) %>% 
  filter(., road_length <= 500 & road_length > 49) %>% 
  geom_sf(data = ., color = 'red') + theme_minimal()

roads_ss <- roads %>% 
  mutate(road_length = as.numeric(units::set_units(st_length(.), "km"))) %>% 
  filter(., road_length <= 500 & road_length > 49)

p + geom_sf(data = roads_ss, color = 'red') + 
  theme_minimal()

coords <- cbind("x" = c(25, 25.1, 26.4, 26.4, 26.4, 25), 
                "y" = c(-15.5, -14.5, -14.6, -14.8, -15.6, -15.5))
pol <- st_polygon(x = list(coords)) %>% 
  st_sfc %>% 
  st_sf(ID = 1, crs = 4326)
roads_rpj <- st_transform(roads, crs = st_crs(districts))

par(mar = rep(0, 4))
plot(st_geometry(districts), col = "grey")
plot(pol, col = "purple", add = TRUE)
plot(roads_rpj$geometry, col = 'red', add = TRUE)

st_intersects(x = roads_rpj, y = districts)

highway <- roads_rpj[401, ]
dists_int <- districts %>% 
  slice(st_intersects(x = highway, y = .)[[1]])
# unlist(st_intersects(x = highway, y = districts))
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
dists_int
plot(st_geometry(dists_int), col = 'grey')
plot(highway, col = 'red', add = TRUE)

pol_int_dists <- st_intersection(x = pol, y = districts)

par(mfrow = c(1, 2), mar = rep(0, 4))
plot(st_geometry(districts), col = "grey")
plot(st_geometry(pol_int_dists), col = rainbow(n = nrow(pol_int_dists)), 
     add = TRUE)
plot(st_geometry(districts), col = "grey")
districts %>% 
  slice(unlist(st_intersects(x = pol, y = .))) %>% 
  plot(add = TRUE)

d1 <- st_difference(x = districts, y = pol)
d2 <- st_difference(x = pol, y = districts)

par(mfrow = c(1, 2), mar = c(0, 0, 1, 0))
d1 %>% 
  st_geometry %>% 
  plot(col = "grey")
d2 %>% 
  st_geometry %>%
  plot(col = "grey")

uni1 <- st_union(x = pol, y = districts)
uni1

uni2 <- st_union(x = districts, y = pol)
uni2

par(mfrow = c(1, 2), mar = c(0, 0, 1, 0))
plot(uni1[2], reset = FALSE, main = "Union of pol with districts")
plot(uni2[1], reset = FALSE, main = "Union of districts with pol")

## project everything to Albers
farmers_alb <- farmers_sf %>% 
  st_transform(st_crs(roads))
long_roads <- roads %>% 
  filter(as.numeric(st_length(.)) > 500000)
districts_alb <- districts %>% 
  st_transform(st_crs(roads))

roads_buff_20km <- st_buffer(long_roads, 20000)
plot(st_geometry(districts_alb), col = "grey")
plot(st_geometry(roads_buff_20km), col = "blue", add = TRUE)
plot(st_geometry(long_roads), col = "red", add = TRUE)

par(mar = rep(0, 4))
districts %>% 
  st_union %>% 
  plot(col = "grey")
set.seed(1)
districts %>% 
  st_union %>% 
  st_sample(size = 25, exact = TRUE) %>% 
  plot(pch = 20, add = TRUE)

par(mar = rep(0, 4))
districts %>% 
  st_geometry %>% 
  plot(col = "grey")
set.seed(1)
districts %>% 
  st_sample(size = rep(5, nrow(.)), by_polygon = TRUE) %>% 
  plot(pch = 20, add = TRUE)

set.seed(1)
pts <- st_union() %>% 
  st_sample() %>% 
  st_transform(st_crs())
dists_in_pts <- districts_alb %>% 
  slice(unlist(st_intersects(st_buffer(, dist = 50000), .)))
ggplot() + 
  geom_sf(data = ) + 
  geom_sf(data = , fill = "blue") + 
  geom_sf(data = )

farmers_alb %>% 
  slice(unlist(st_intersects(st_buffer(), .))) %>% 
  ?

# code to adapt
roads %>% mutate(km = as.numeric(st_length(.)) / 1000) %>% 
  filter(km > 200) %>% 
  st_buffer(dist = 20000) %>% 
  st_transform(crs = 4326)
st_difference(larger_object, smaller_object)
ggplot() + geom_sf

buffered_object %>% st_sample()

library(ggplot2)
mypoly <- st_polygon(list(cbind(x = c(26.5, 27.5, 27, 26, 26.5), 
                                y = c(-15.5, -16.5, -17, -16, -15.5)))) %>%
  st_sfc(crs = 4326) %>% st_intersection(., districts) %>% st_as_sf() %>% 
  mutate(ID = seq(1:nrow(.)))
#st_crs(mypoly) <- 4326
smallest_centroid <- districts %>% 
  filter(distName == "Lusaka") %>% 
  st_centroid()


p <- ggplot(districts) + 
  geom_sf() + 
  geom_sf(data = smallest_centroid, color = "red")

p1 <- p + geom_sf(data = mypoly, aes(color = as.factor(ID)),
                  fill = "transparent") + labs(color = "ID")
# small_pol2 <- pol %>% st_sfc() %>% st_sf(crs = 4326) %>% mutate(ID = 1)
p2 <- p + geom_sf(data = mypoly, aes(color = as.factor(ID)), 
                  fill = "transparent") + 
  labs(color = "ID")

gridExtra::grid.arrange(p1, p2, ncol = 2)

g1 <- cowplot::plot_grid(p1 + theme(legend.position = "none"), 
                         p2 + theme(legend.position = "none"))
cowplot::plot_grid(g1, cowplot::get_legend(p1), rel_widths = c(2, 0.3))

library(patchwork)

p1 + p2

p1 + p2 + plot_layout(guides = "collect")

p1 + p1 / p
(p1 + p2) / (p1 + p2) + plot_layout(guides = "collect")