# Make and save the visual as 'g1'
g1 <- ggplot() +
  # Visualize a nice background outline
  geom_sf(data = precincts, color = "darkgrey", size = 5) +
  # Visualize tan precinct polygons
  geom_sf(data = precincts, color = "white", fill = "tan") +
  # Visualize polling places as squares (shape = 22)
  geom_sf(data = polls, fill = "black", size = 1, shape = 22) +
  # visualize T-lines
  geom_sf(data = tlines, mapping = aes(color = line), size = 0.5) +
  # Vsiaulize t-stops
  geom_sf(data = tstops, mapping = aes(color = line), size = 1.15) +
  # Assign a color scheme to the t-lines and stops
  scale_color_manual(
    # Order the levels as desired
    breaks = c("GREEN", "ORANGE", "RED", "BLUE", "SILVER", "MULTIPLE"),
    # Now add the labels you want
    labels = c("Green", "Orange", "Red", "Blue", "Silver", "Multiple"),
    # And give each a color
    values = c("seagreen", "orangered", "firebrick", 
               "dodgerblue", "grey","purple")) +
  # Add a clear, simple theme
  theme_void(base_size = 30) +
  # And a basic legend
  labs(color = "Public\nTransit\nLine",
       subtitle = "Polling Places vs. Public Transit\nin Boston Precincts")

# Display the visual
g1

polls <- read_rds("polling_places.rds") %>%
  # Transform to Equal Distance Projection
  st_transform(crs = aed) 

tstops <- read_rds("train_stops_boston.rds") %>%
  # Transform to Equal Distance Projection
  st_transform(crs = aed) %>%
  # Save a copy of these train stop's geometry
  # in a new column, called geometry_stop
  mutate(geometry_stop = geometry)

# Identify which stops are closest to each polling places
nearest <- polls %>%
  # Join together the stops nearest
  st_join(tstops, join = st_nearest_feature) %>%
  # Extract the coordinates of our polls and our new t-stops'
         # st_coordinates returns a two-column matrix (x and y),
         # so we grab the first column for the x-coordiante,
         # and the second column for the y-coordinate
  mutate(x_poll = st_coordinates(geometry)[,1],
         y_poll = st_coordinates(geometry)[,2],
         x_stop = st_coordinates(geometry_stop)[,1],
         y_stop = st_coordinates(geometry_stop)[,2]) %>%
  # Convert to data.frame
  as.data.frame() %>%
  # Get rid of our now unnecessary geometry fields
  select(-geometry, -geometry_stop)

mydist <- nearest %>%
  # Zoom into just precincts where their coordinates
  # and their nearest t-stop coordinates could be obtained
  filter(!is.na(x_poll), !is.na(y_poll),
         !is.na(x_stop), !is.na(y_stop)) %>%
  # For each polling station and its precinct,
  group_by(id, ward_precinct) %>%
  # Build a new geometry column, with a 4-cell matrix
  mutate(geometry = matrix(
    # listing our polling place x and our t-stop x coordinate
    c(x_poll,  x_stop, 
    # and listing our polling place y and our t-stop y coordinate
      y_poll,  y_stop),
    # Make the matrix have two-rows and two-columns
    ncol = 2, nrow = 2) %>%
      # convert it to a line-string geometry
      st_linestring(dim = "XY") %>%
      # And bind it as a spatial data.frame,
      # with an Equal Distance Conic projection
      st_sfc(crs = aed)) %>%
    # Turn the whole data.frame into a spatial data.frame
  # with an Equal Distance Conic projection
  st_as_sf(crs = aed) %>%
  # Finally, calculate the distance of our linestrings!
  # It will calculate in the original units of our projections,
  # so meters, in this case.
  mutate(distance = st_length(geometry) %>% as.numeric())

mydist %>%
  head(3) %>%
  knitr::kable()

g1 +
  # Let's layer on top our new lines!
  geom_sf(data = mydist %>%
            # Remember to transform it to equal area conic
            st_transform(crs = aea), color = "black")

mydist %>%
  as.data.frame() %>%
  ggplot(mapping = aes(x = distance)) +
  geom_histogram(mapping = aes(y = ..density..), 
                 fill = "dodgerblue", color = "white") +
  geom_density(fill = "dodgerblue", alpha = 0.5) +
  labs(x = "Distance between Polling Place\nand Nearest T-Stop (meters)",
       y = "(%) Frequency\n(n = 256 polling places)") +
  theme_bw(base_size = 24)

mydist %>%
  write_sf("mydistances.kml", driver = "kml", delete_layer = TRUE)

mydist %>%
  write_sf("mydistances.shp", driver = "ESRI Shapefile", delete_layer = TRUE)