In nhoteling/spaceheater: Convenience functions for spatial analysis

knitr::opts_chunk$set(echo = TRUE)

# For data handling
library(dplyr)
library(stringr)
library(readr)

# For spatial things
library(sf)
library(spatstat)             # for point pattern analysis
# library(lwgeom)             # needed to run st_distance
library(rnaturalearth)        # to get geo shapes
library(rnaturalearthhires)   # same as above
library(tidycensus)
library(tigris)               # to download counties data
library(osmdata)
library(maptools)

library(alphahull)            # detailed pgons from point data
library(dbscan)               # spatial clustering
library(solitude)  # for isolation forest stuff

# Plotting
library(ggplot2)

source("../R/spaceheater.R")
source("../R/spaceheater_data.R")

streets <- c("motorway", "trunk","primary","secondary","tertiary")
bb <- getbb("Portland, OR")
pdx_streets <- opq(bb) %>%
  add_osm_feature(key="highway", value=streets) %>%
  osmdata_sf()
pdx_waters <- opq(bb) %>%
  add_osm_feature(key="natural", value="water") %>%
  osmdata_sf()
saveRDS(pdx_streets, file="../data/pdx_streets.rds")
saveRDS(pdx_waters, file="../data/pdx_waters.rds")

pdx_streets <- readRDS("../data/pdx_streets.rds")
pdx_waters <- readRDS("../data/pdx_waters.rds")


epsg_code <- spaceheater_epsg(mean(st_coordinates(df.pdx$geometry)[,"X"]),
                              mean(st_coordinates(df.pdx$geometry)[,"Y"]))

bb1 <- parse_number(as.vector(str_split(pdx_streets$bbox, pattern=",",simplify=TRUE)))
bbox1 <- st_bbox(c("xmin"=bb1[2],"xmax"=bb1[4],
                  "ymin"=bb1[1],"ymax"=bb1[3]))
pgon1 <- st_as_sfc(bbox1) %>% 
  st_set_crs(st_crs(df.pdx$geometry)$input) %>% 
  st_transform(epsg_code)

df.streets1 <- pdx_streets$osm_lines %>% 
  st_transform(epsg_code) %>%
  st_crop(pgon1)

df.waters1 <- pdx_waters$osm_multipolygons %>%
  st_transform(epsg_code) %>%
  st_buffer(0) %>%
  st_crop(pgon1)

df.bwys <- df.pdx %>% 
  mutate(geometry = st_transform(geometry,epsg_code)) %>%
  filter(st_intersects(geometry,pgon1,sparse=FALSE))

df.untp <- readRDS("../data/dfuntp.rds") %>% 
  dplyr::select(byname,nratings, rating) %>%
  distinct() %>% 
  filter( (byname %in% df.bwys$byname) & nratings>100) %>%
  arrange(desc(rating)) %>%
  slice_head(n=10)

df.bwys1 <- df.bwys %>% 
  filter(byname %in% df.untp$byname) %>% 
  dplyr::select(byname,geometry) %>%
  left_join(df.untp,by=c("byname"))

hwy_vals <- c("motorway"=1.00, "trunk"=0.75, "primary"=0.50, "secondary"=0.3,"tertiary"=0.2)
p.pdx <- ggplot() + 
  geom_sf(data=st_buffer(pgon1,200), aes(geometry=geometry), fill="cornsilk") +
  geom_sf(data=df.waters1, aes(geometry=geometry), fill="aliceblue",size=0.2,color="grey60") +
  geom_sf(data=df.streets1, aes(geometry=geometry, size=highway), color="grey60") +
  geom_sf(data=df.bwys, aes(geometry=geometry), color="#de2d26", alpha=0.6, size=1) +
  #geom_sf(data=df.bwys1, aes(geometry=geometry), 
  #        shape=21, color="grey50",fill="#de2d26", alpha=0.8, size=2.5) +
  #lims(x=c(bbox1["xmin"],bbox["xmax"]),y=c(bbox["ymin"],bbox["ymax"])) +
  scale_size_manual(values=hwy_vals, guide="none") +
  theme_void()

p.pdx

#bwys2 <- st_as_sf(snapPointsToLines(as_Spatial(df.bwys$geometry),
#                                    as_Spatial(df.streets1$geometry)))

bwys2 <- snap_to_line(df.bwys$geometry, df.streets1$geometry)

p.pdx <- ggplot() + 
  geom_sf(data=pgon1, aes(geometry=geometry), fill="cornsilk") +
  geom_sf(data=df.waters1, aes(geometry=geometry), fill="aliceblue",size=0.2,color="grey60") +
  geom_sf(data=df.streets1, aes(geometry=geometry, size=highway), color="grey60") +
  geom_sf(data=bwys2, aes(geometry=geometry), color="#de2d26", alpha=0.6, size=1) +
  geom_text(data=df.bwys1, aes(x=st_coordinates(geometry)[,"X"],
                            y=st_coordinates(geometry)[,"Y"],
                            label=byname)) +
  #lims(x=c(bbox1["xmin"],bbox["xmax"]),y=c(bbox["ymin"],bbox["ymax"])) +
  scale_size_manual(values=hwy_vals, guide="none") +
  theme_void()

base_date <- "2021-06-19"
df.visits <- data.frame(byname=c("Cascade Brewing", "Great Notion Brewing",
                                   "Ecliptic Brewing", "Breakside Brewery",
                                   "Trap Door Brewing", "Fortside Brewing Company",
                                   "Level Beer", "Labyrinth Forge Brewing Company",
                                   "Threshold Brewing & Blending", "Wayfinder Beer"),
                        tm=c("12:00", "13:00", "14:00", "15:00", "16:00", "17:00",
                             "18:00", "19:00", "20:00", "21:00")) %>%
  mutate(timestamp=as.POSIXct(paste(base_date,tm))) %>%
  left_join(df.bwys %>% dplyr::select(byname,geometry), by="byname") %>%
  arrange(timestamp)

vvjust <- c(0.5, 0.5, 0.5, -0.2, 0.5, 1.2, 1.3, 1.4, 0.5, 0.5)
hhjust <- c(1.05, 1.05, 1.05, -0.02, 1.05, 0.2, 0.08, 0.1, 1.05, 1.05)
p.visits <- ggplot() + 
  geom_sf(data=st_buffer(pgon1,200), aes(geometry=geometry), fill="cornsilk") +
  geom_sf(data=df.waters1, aes(geometry=geometry), fill="aliceblue",size=0.2,color="grey60") +
  geom_sf(data=df.streets1, aes(geometry=geometry, size=highway), color="grey60") +
  geom_sf(data=df.bwys, aes(geometry=geometry), color="#de2d26", alpha=0.3, size=1) +
  geom_segment(data=df.visits, aes(x=st_coordinates(geometry)[,"X"],
                                   xend=lead(st_coordinates(geometry)[,"X"]),
                                   y=st_coordinates(geometry)[,"Y"],
                                   yend=lead(st_coordinates(geometry)[,"Y"])),
               size=1.5, alpha=0.6, color="#008ae6") +
  geom_sf(data=df.bwys1, aes(geometry=geometry), 
          shape=21, color="#008ae6",fill="#008ae6", alpha=0.8, size=2.5) +
  #geom_text(data=df.bwys1, aes(x=st_coordinates(geometry)[,"X"],
  #                          y=st_coordinates(geometry)[,"Y"],
  #                          label=byname), 
  #          hjust=hhjust, vjust=vvjust, color="#008ae6",size=4) +
  scale_size_manual(values=hwy_vals, guide="none") +
  theme_void()


#df %>% 
#  add_static_interval() %>%
#  interpolate_points() %>%
#  mutate(geometry=snap_to_road(geometry,roads)) %>%
#  segmentize()