# Load Libraries

library(leaflet)
library(sf)
## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(readr)
library(ggplot2)
library(urbnmapr)
sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
# Load and Clean EFAM Data
efam <- read_csv("~/Downloads/efam_data_geocoded2024.csv") %>%
  rename(
    State = starts_with("State"),
    City  = starts_with("City"),
    Zip   = starts_with("Zip")
  ) %>%
  mutate(
    Latitude  = as.numeric(Latitude),
    Longitude = as.numeric(Longitude)
  ) %>%
  filter(
    !is.na(Latitude),
    !is.na(Longitude),
    !(Latitude == 0 & Longitude == 0),
    !(Longitude > -20 & Longitude < 20 & Latitude > -20 & Latitude < 20)
  ) %>%
  filter(!(Type == "Summer Meal Sites" &
           grepl("BINGHAMTON", Name, ignore.case = TRUE)))
## New names:
## • `Street` -> `Street...3`
## • `City` -> `City...4`
## • `State` -> `State...5`
## • `Zip` -> `Zip...6`
## • `Street` -> `Street...20`
## • `City` -> `City...23`
## • `State` -> `State...24`
## • `Zip` -> `Zip...26`
## Rows: 179 Columns: 28
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (22): Type, Name, Street...3, City...4, State...5, Hours of Operation, A...
## dbl  (6): Zip...6, Latitude, Longitude, Accuracy Score, Unit Number, Zip...26
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Convert EFAM to SF (WGS84)
efam_sf <- st_as_sf(
  efam,
  coords = c("Longitude", "Latitude"),
  crs = 4326,
  remove = FALSE
)


# Load Flood Shapefiles
flood06 <- st_read("~/Downloads/Flood_2006/2006_Flood.shp")
## Reading layer `2006_Flood' from data source 
##   `/Users/malaykamudassar/Downloads/flood_2006/2006_Flood.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 935338.9 ymin: 728907.9 xmax: 1045917 ymax: 802227.6
## Projected CRS: NAD83 / New York Central (ftUS)
flood11 <- st_read("~/Downloads/Flood_2011/Flood_2011_09_08_High_Water.shp")
## Reading layer `Flood_2011_09_08_High_Water' from data source 
##   `/Users/malaykamudassar/Downloads/flood_2011/Flood_2011_09_08_High_Water.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 52 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 1265455 ymin: 15264730 xmax: 1425718 ymax: 15344070
## Projected CRS: NAD83 / UTM zone 18N
# Convert All Layers to Projected CRS (EPSG 26918)

flood06_26918 <- st_transform(flood06, 26918)
flood11_26918 <- st_transform(flood11, 26918)
efam_26918    <- st_transform(efam_sf, 26918)

# Create 500 m Buffers
flood06_buf_26918 <- st_buffer(flood06_26918, dist = 500)
flood11_buf_26918 <- st_buffer(flood11_26918, dist = 500)


# Spatial Intersections: Identify EFAM Locations Inside Buffers
efam_26918 <- efam_26918 %>%
  mutate(
    within2006 = lengths(st_intersects(., flood06_buf_26918)) > 0,
    within2011 = lengths(st_intersects(., flood11_buf_26918)) > 0
  )


# Transform Everything Back to WGS84 for Leaflet Mapping

efam_ll        <- st_transform(efam_26918, 4326)
flood06_ll     <- st_transform(flood06_26918, 4326)
flood11_ll     <- st_transform(flood11_26918, 4326)
flood06_buf_ll <- st_transform(flood06_buf_26918, 4326)
flood11_buf_ll <- st_transform(flood11_buf_26918, 4326)

# Create Color Palette for EFAM Types
type_vals <- sort(unique(efam_ll$Type))
pal <- colorFactor("Set2", domain = type_vals)

# Build Final Leaflet Map
leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%

  # 2006 Flood
  addPolygons(
    data = flood06_ll,
    color = "sienna",
    fillColor = "darkorange",
    fillOpacity = 0.5,
    weight = 1,
    group = "2006 Flood"
  ) %>%
  addPolygons(
    data = flood06_buf_ll,
    color = "orange",
    fillColor = "gold",
    fillOpacity = 0.25,
    weight = 1,
    group = "2006 500 m Buffer"
  ) %>%

  # 2011 Flood
  addPolygons(
    data = flood11_ll,
    color = "firebrick",
    fillColor = "maroon",
    fillOpacity = 0.5,
    weight = 1,
    group = "2011 Flood"
  ) %>%
  addPolygons(
    data = flood11_buf_ll,
    color = "red",
    fillColor = "pink",
    fillOpacity = 0.25,
    weight = 1,
    group = "2011 500 m Buffer"
  ) %>%

  # EFAM Points
  addCircleMarkers(
    data = efam_ll,
    radius = 6,
    stroke = TRUE, weight = 1,
    color = "#333333",
    fillColor = ifelse(efam_ll$within2006 | efam_ll$within2011,
                       "red",      # within 500 m
                       "grey60"),  # outside
    fillOpacity = 0.9,
    label = ~Name,
    group = "EFAM Sites"
  ) %>%

  # Layer Control
  addLayersControl(
    overlayGroups = c("EFAM Sites",
                      "2006 Flood", "2006 500 m Buffer",
                      "2011 Flood", "2011 500 m Buffer"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%

  # Legend
  addLegend(
    "bottomright",
    pal = pal,
    values = ~Type,
    data = efam_ll,
    title = "EFAM Type",
    opacity = 1
  )