#install.packages("plotly")
#install.packages("raster")
library(sf)
library(tmap)
## Warning: package 'tmap' was built under R version 4.4.3
library(readr)
library(dplyr)
library(purrr)
library(scales)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
library(plotly)
## Warning: package 'plotly' was built under R version 4.4.3
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.4.3
library(mapview)
library(stringr)
library(ggplot2)
library(tidyverse)
library(lubridate)
library(gganimate)
## Warning: package 'gganimate' was built under R version 4.4.3
# Load CSV
hurricane_data <- read_csv("./F_Project/HURRICANE/use_Hurricanes.csv") %>%
  mutate(
    Name = str_trim(Name),
    Year = as.integer(Year),
    STORM_ID = str_trim(toupper(STORM_ID))
  )
## Rows: 65 Columns: 15
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (8): Name, STORM_ID, Start_Date, Peak_Date, End_Date, Ocean, Areas_affec...
## dbl (7): Year, Duration_h, Max_Wind_Speed, Max_Wind_Speed_(mph), Highest_Cat...
## 
## ℹ 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.
meta_path <- "./F_Project/HURRICANE/L_shapefile_meta_Hurricanes.csv"

# Load metadata
meta_data <- read_csv("./F_Project/HURRICANE/L_shapefile_meta_Hurricanes.csv") %>%
  mutate(
    STORM_ID = str_trim(toupper(STORM_ID))  # Ensure uniform formatting
  )
## Rows: 65 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (8): Name, STORM_ID, Folder, Basin, PtsFile, LinFile, RadiiFile, windswa...
## dbl (1): Year
## 
## ℹ 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.
# Merge metadata and hurricane data
combined_meta <- meta_data %>%
  left_join(hurricane_data, by = "STORM_ID")
# Define base path and generate shapefile paths
base_path <- "./F_Project/HURRICANE/hurricanes"

combined_meta <- combined_meta %>%
  mutate(
    LinPath = file.path(base_path, Folder, LinFile),
    PtsPath = file.path(base_path, Folder, PtsFile),
    WindsPath = file.path(base_path, Folder, windswathFile)
  )

# Function to safely read shapefiles
safe_read_sf <- function(path) {
  if (!is.na(path) && file.exists(path)) {
    tryCatch(read_sf(path), error = function(e) NULL)
  } else {
    NULL
  }
}
# WindsPath

# Add shapefiles as spatial data columns
New_Hurricane_base <- combined_meta %>%
  mutate(
    PointsSF = map(PtsPath, safe_read_sf),
    LinesSF = map(LinPath, safe_read_sf),
    WindwathSF = map(WindsPath, safe_read_sf)
  )
# Filter for Hurricane Debby
Debby <- New_Hurricane_base %>% filter(STORM_ID == "DEBBY2024") # Ensure correct case

# Check if any rows were found for Debby

if (nrow(Debby) > 0) {
  # Check if LinesSF has at least one element and if it's not NULL
  
  if (length(Debby$LinesSF) > 0 && !is.null(Debby$LinesSF[[1]])) {
    # Check if PointsSF has at least one element and if it's not NULL
    
    if (length(Debby$PointsSF) > 0 && !is.null(Debby$PointsSF[[1]])) {
      
      if (length(Debby$WindwathSF) > 0 && !is.null(Debby$WindwathSF[[1]])) {
      # Check if WindwathSF has at least one element and if it's not NULL
        
      mapview(Debby$LinesSF[[1]], color = "blue") +
        mapview(Debby$PointsSF[[1]], col.region = "red") +
          mapview(Debby$WindwathSF[[1]], col.region = "yellow")
        } else {
          message("No valid point spatial data available for DEBBY2024.")
          }
      } else {
        message("No valid line spatial data available for DEBBY2024.")
        }
    } else {
      message("No valid Windwath spatial data available for DEBBY2024.")
      }
  } else {
    message("No data found for hurricane DEBBY2024.")
}
#WindwathSF = map(WindsPath, safe_read_sf)



# Add shapefiles as spatial data columns
New_Hurricane_base <- combined_meta %>%
  mutate(
    PointsSF = map(PtsPath, safe_read_sf),
    LinesSF = map(LinPath, safe_read_sf),
    WindwathSF = map(WindsPath, safe_read_sf)
  )

# Filter for Hurricane Debby
Debby <- New_Hurricane_base %>% filter(STORM_ID == "DEBBY2024") # Ensure correct case

# Check if any rows were found for Debby
if (nrow(Debby) > 0) {
  # Initialize an empty mapview object
  debby_map <- mapview()

  # Add line data if it exists
  if (length(Debby$LinesSF) > 0 && !is.null(Debby$LinesSF[[1]])) {
    debby_map <- debby_map + mapview(Debby$LinesSF[[1]], color = "blue", layer.name = "Debby Track")
  } else {
    message("No valid line spatial data available for DEBBY2024.")
  }

  # Add point data if it exists
  if (length(Debby$PointsSF) > 0 && !is.null(Debby$PointsSF[[1]])) {
    debby_map <- debby_map + mapview(Debby$PointsSF[[1]], col.region = "red", layer.name = "Debby Points")
  } else {
    message("No valid point spatial data available for DEBBY2024.")
  }

  # Add Windwath data if it exists
  if (length(Debby$WindwathSF) > 0 && !is.null(Debby$WindwathSF[[1]])) {
    debby_map <- debby_map + mapview(Debby$WindwathSF[[1]], color = "yellow", layer.name = "Debby Windwath")
  } else {
    message("No valid Windwath spatial data available for DEBBY2024.")
  }

  # Display the map if any layers were added
  if (length(debby_map@map) > 0) {
    debby_map
  } else {
    message("No spatial data to display for DEBBY2024.")
  }

} else {
  message("No data found for hurricane DEBBY2024.")
}
New_Hurricane_last <- New_Hurricane_base %>%
  select(-Name.x, -Year.x) %>%
  rename(
    Name = Name.y,
    Year = Year.y
  )
selected_year <- 2020

hurricanes_in_year <- New_Hurricane_last %>%
  filter(Year == selected_year)

print(hurricanes_in_year %>% select(STORM_ID, Name, Year))
## # A tibble: 9 × 3
##   STORM_ID    Name     Year
##   <chr>       <chr>   <int>
## 1 DELTA2020   Delta    2020
## 2 DOUGLAS2020 Douglas  2020
## 3 ETA2020     Eta      2020
## 4 HANNA2020   Hanna    2020
## 5 ISAIAS2020  Isaias   2020
## 6 LAURA2020   Laura    2020
## 7 MARCO2020   Marco    2020
## 8 SALLY2020   Sally    2020
## 9 ZETA2020    Zeta     2020
mapview::mapviewOptions(fgb = FALSE)  # Helps if performance is slow or viewer has issues
map_year <- mapview()
for (i in seq_len(nrow(hurricanes_in_year))) {
  hurricane_name <- hurricanes_in_year$Name[i] # Get the name of the current hurricane

  if (!is.null(hurricanes_in_year$LinesSF[[i]])) {
    map_year <- map_year +
      mapview(hurricanes_in_year$LinesSF[[i]], color = "blue",
              layer.name = hurricane_name) # Set layer name to hurricane name
  }

  if (!is.null(hurricanes_in_year$PointsSF[[i]])) {
    map_year <- map_year +
      mapview(hurricanes_in_year$PointsSF[[i]], col.region = "red",
              layer.name = paste(hurricane_name, "Points")) # Add "Points" for clarity
  }

  if (!is.null(hurricanes_in_year$WindwathSF[[i]])) {
    map_year <- map_year +
      mapview(hurricanes_in_year$WindwathSF[[i]], col.region = "yellow",
              layer.name = paste(hurricane_name, "Windwath")) # Add "Windwath" for clarity
  }
}
map_year
mapview::mapviewOptions(fgb = FALSE)  # Helps if performance is slow or viewer has issues
map_year <- mapview()
for (i in seq_len(nrow(hurricanes_in_year))) {
  hurricane_name <- hurricanes_in_year$Name[i] # Get the name of the current hurricane

  if (!is.null(hurricanes_in_year$LinesSF[[i]])) {
    map_year <- map_year +
      mapview(hurricanes_in_year$LinesSF[[i]], color = "blue",
              layer.name = hurricane_name) # Set layer name to hurricane name
  }
}
map_year
mapview::mapviewOptions(fgb = FALSE)  # Helps if performance is slow or viewer has issues
map_year <- mapview()
for (i in seq_len(nrow(hurricanes_in_year))) {
  hurricane_name <- hurricanes_in_year$Name[i] # Get the name of the current hurricane

  if (!is.null(hurricanes_in_year$WindwathSF[[i]])) {
    map_year <- map_year +
      mapview(hurricanes_in_year$WindwathSF[[i]], col.region = "yellow",
              layer.name = paste(hurricane_name)) # Add "Windwath" for clarity
  }
}
map_year
strongest_hurricanes <- New_Hurricane_last %>%
  filter(Max_Wind_Speed == max(Max_Wind_Speed, na.rm = TRUE))

print(strongest_hurricanes %>% select(Name, Year, Max_Wind_Speed, Min_Pressure, Damage, Deaths))
## # A tibble: 2 × 6
##   Name    Year Max_Wind_Speed Min_Pressure Damage         Deaths
##   <chr>  <int>          <dbl>        <dbl> <chr>           <dbl>
## 1 Dorian  2019            160          910 5,000,000,000      78
## 2 Wilma   2005            160          882 21,007,000,000     23
longest_duration <- New_Hurricane_last %>%
  filter(Duration_h == max(Duration_h, na.rm = TRUE))

print(longest_duration %>% select(Name, Year, Duration_h))
## # A tibble: 1 × 3
##   Name   Year Duration_h
##   <chr> <int>      <dbl>
## 1 Ivan   2004        522
New_Hurricane_stats <- New_Hurricane_last %>%
  mutate(
    Max_Wind_Speed_mph = Max_Wind_Speed * 1.15078, # Convert knots to mph if needed
    Category = cut(Max_Wind_Speed_mph,
                   breaks = c(-Inf, 38, 73, 96, 113, 136, Inf), # Tropical Depression/Storm, Cat 1-5
                   labels = c("TD/TS", "Cat 1", "Cat 2", "Cat 3", "Cat 4", "Cat 5"),
                   right = FALSE, # Intervals are [lower, upper) except the last
                   include.lowest = TRUE,
                   ordered_results = TRUE)
  )
# Interactive Plot: Number of Hurricanes per Year
ggplotly(ggplot(New_Hurricane_stats, aes(x = Year)) +
           geom_histogram(binwidth = 1, fill = "red", color = "white") +
           labs(title = "Number of US-Affected Hurricanes per Year",
                x = "Year",
                y = "Number of Hurricanes") +
           theme_minimal())
# Static Plot: Number of Hurricanes per Year
ggplot(New_Hurricane_stats, aes(x = Year)) +
  geom_histogram(binwidth = 1, fill = "steelblue", color = "white") +
  labs(title = "Number of US-Affected Hurricanes per Year",
       x = "Year",
       y = "Number of Hurricanes") +
  theme_minimal()

# Interactive Plot: Number of Hurricanes per Year
ggplotly(ggplot(New_Hurricane_stats, aes(x = Year)) +
           geom_histogram(binwidth = 1, fill = "green", color = "white") +
           labs(title = "Number of US-Affected Hurricanes per Year",
                x = "Year",
                y = "Number of Hurricanes") +
           theme_minimal())
# Static Plot: Number of Hurricanes by Saffir-Simpson Category
ggplot(New_Hurricane_stats %>% filter(!is.na(Category)), aes(x = Category, fill = Category)) +
  geom_bar() +
  labs(title = "Number of US-Affected Hurricanes by Saffir-Simpson Category",
       x = "Category",
       y = "Number of Hurricanes") +
  theme_minimal() +
  theme(legend.position = "none") # Hide legend if mapping fill to x-axis variable

# Group by year and list hurricane names
hurricanes_by_year <- New_Hurricane_stats %>%
  group_by(Year) %>%
  summarize(Hurricanes = paste(Name, collapse = ", "))

# Interactive Plot: Number of Hurricanes per Year with Names on Hover
interactive_year_plot_names <- ggplot(New_Hurricane_stats, aes(x = Year)) +
  geom_histogram(binwidth = 1, fill = "green", color = "white",
                 aes(text = paste("Year:", Year, "<br>Hurricanes:",
                                   hurricanes_by_year$Hurricanes[match(Year, hurricanes_by_year$Year)]))) +
  labs(title = "Number of US-Affected Hurricanes per Year",
       x = "Year",
       y = "Number of Hurricanes") +
  theme_minimal()
## Warning in geom_histogram(binwidth = 1, fill = "green", color = "white", :
## Ignoring unknown aesthetics: text
ggplotly(interactive_year_plot_names, tooltip = "text")
# Static Plot: Average Maximum Wind Speed per Year
avg_wind_speed_per_year <- New_Hurricane_stats %>%
  group_by(Year) %>%
  summarize(Avg_Max_Wind_Speed_mph = mean(Max_Wind_Speed_mph, na.rm = TRUE))

ggplot(avg_wind_speed_per_year, aes(x = Year, y = Avg_Max_Wind_Speed_mph)) +
  geom_line(color = "firebrick") +
  geom_point(color = "firebrick") +
  labs(title = "Average Maximum Wind Speed of US-Affected Hurricanes per Year",
       x = "Year",
       y = "Average Maximum Wind Speed (MPH)") +
  theme_minimal()

# Group hurricane names by wind speed bins
wind_speed_bins <- seq(floor(min(New_Hurricane_stats$Max_Wind_Speed_mph, na.rm = TRUE) / 10) * 10,
                       ceiling(max(New_Hurricane_stats$Max_Wind_Speed_mph, na.rm = TRUE) / 10) * 10,
                       by = 10)

hurricane_names_by_wind_speed <- New_Hurricane_stats %>%
  mutate(Wind_Speed_Bin = cut(Max_Wind_Speed_mph, breaks = wind_speed_bins, include.lowest = TRUE)) %>%
  group_by(Wind_Speed_Bin) %>%
  summarize(Hurricanes = paste(Name, collapse = ", "))

# Interactive Plot: Distribution of Maximum Wind Speed (MPH) with Names on Hover
interactive_wind_speed_hist <- ggplot(New_Hurricane_stats, aes(x = Max_Wind_Speed_mph)) +
  geom_histogram(binwidth = 10, fill = "darkorange", color = "white",
                 aes(text = paste("Wind Speed (MPH): [", floor(..x..), ", ", ceiling(..x..), ")<br>",
                                   "Count:", ..count.., "<br>",
                                   "Hurricanes:",
                                   hurricane_names_by_wind_speed$Hurricanes[match(as.character(cut(..x.., breaks = wind_speed_bins, include.lowest = TRUE)),
                                                                                  as.character(hurricane_names_by_wind_speed$Wind_Speed_Bin))]))) +
  labs(title = "Distribution of Maximum Wind Speeds (MPH)",
       x = "Maximum Wind Speed (MPH)",
       y = "Count") +
  theme_minimal()
## Warning in geom_histogram(binwidth = 10, fill = "darkorange", color = "white",
## : Ignoring unknown aesthetics: text
ggplotly(interactive_wind_speed_hist, tooltip = "text")
## Warning: The dot-dot notation (`..x..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(x)` instead.
## ℹ The deprecated feature was likely used in the base package.
##   Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Group hurricane names by duration bins
duration_bins <- seq(0, ceiling(max(New_Hurricane_stats$Duration_h, na.rm = TRUE) / 24) * 24, by = 24) # Bin by days

hurricane_names_by_duration <- New_Hurricane_stats %>%
  mutate(Duration_Bin = cut(Duration_h, breaks = duration_bins, include.lowest = TRUE)) %>%
  group_by(Duration_Bin) %>%
  summarize(Hurricanes = paste(Name, collapse = ", "))

# Interactive Plot: Distribution of Duration (Hours) with Names on Hover
interactive_duration_hist <- ggplot(New_Hurricane_stats, aes(x = Duration_h)) +
  geom_histogram(binwidth = 24, fill = "forestgreen", color = "white",
                 aes(text = paste("Duration (Hours): [", floor(..x..), ", ", ceiling(..x..), ")<br>",
                                   "Count:", ..count.., "<br>",
                                   "Hurricanes:",
                                   hurricane_names_by_duration$Hurricanes[match(as.character(cut(..x.., breaks = duration_bins, include.lowest = TRUE)),
                                                                                  as.character(hurricane_names_by_duration$Duration_Bin))]))) +
  labs(title = "Distribution of Hurricane Duration",
       x = "Duration (Hours)",
       y = "Count") +
  theme_minimal()
## Warning in geom_histogram(binwidth = 24, fill = "forestgreen", color = "white",
## : Ignoring unknown aesthetics: text
ggplotly(interactive_duration_hist, tooltip = "text")
# Static Plot: Total Duration of Hurricanes per Year
total_duration_per_year <- New_Hurricane_stats %>%
  group_by(Year) %>%
  summarize(Total_Duration_h = sum(Duration_h, na.rm = TRUE))

ggplot(total_duration_per_year, aes(x = Year, y = Total_Duration_h)) +
  geom_line(color = "purple") +
  geom_point(color = "purple") +
  labs(title = "Total Duration of US-Affected Hurricanes per Year",
       x = "Year",
       y = "Total Duration (Hours)") +
  theme_minimal()

# Scatter Plot: Max Wind Speed vs. Damage
# Convert 'Damage' to numeric, handling 'None' values
hurricane_data <- hurricane_data %>%
  mutate(Damage = ifelse(Damage == "None", NA, as.numeric(str_replace_all(Damage, ",", ""))))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Damage = ifelse(...)`.
## Caused by warning in `ifelse()`:
## ! NAs introduced by coercion
ggplot(hurricane_data, aes(x = `Max_Wind_Speed_(mph)`, y = Damage)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "red") + # Add a linear trend line
  labs(title = "Hurricane Damage vs. Maximum Wind Speed",
       x = "Maximum Wind Speed (mph)",
       y = "Damage (USD)") +
  theme_minimal() +
  scale_y_continuous(labels = scales::comma, na.value = NA) # Ensure y-axis is treated as numeric
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 14 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Scatter Plot: Year vs. Deaths
ggplot(hurricane_data, aes(x = Year, y = Deaths)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE, color = "blue") + # Loess for non-linear trends
  labs(title = "Hurricane-Related Deaths Over Time",
       x = "Year",
       y = "Number of Deaths") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

# Scatter Plot: Max Wind Speed vs. Duration
ggplot(hurricane_data, aes(x = `Max_Wind_Speed_(mph)`, y = Duration_h)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "purple") +
  labs(title = "Hurricane Duration vs. Maximum Wind Speed",
       x = "Maximum Wind Speed (mph)",
       y = "Duration (Hours)") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

# Function to visualize a hurricane's track
visualize_hurricane_track <- function(storm_id, data) {
  # Filter the data for the specified STORM_ID (case-insensitive comparison)
  hurricane <- data %>% filter(toupper(STORM_ID) == toupper(storm_id))

  # Check if any rows were found for the hurricane
  if (nrow(hurricane) > 0) {
    # Initialize an empty mapview object
    map_hurricane <- mapview()

    # Loop through the hurricane rows (should usually be 1 row, but robust)
    for (i in seq_len(nrow(hurricane))) {
      
      # Check for and add Windwath data
      if (!is.null(hurricane$WindwathSF[[i]]) && inherits(hurricane$WindwathSF[[i]], 'sf')) {
        map_hurricane <- map_hurricane +
          mapview(hurricane$WindwathSF[[i]], col.region = "yellow",
                  layer.name = paste(hurricane$Name[i], "Windspand"))
      } else {
        message(paste("No valid Windwath spatial data for", hurricane$Name[i],
                      "(", hurricane$STORM_ID[i], ")"))
      }
      
    }

    # Display the map
    print(map_hurricane) # Use print() to ensure it displays correctly in all environments

  } else {
    message(paste("No data found for hurricane with STORM_ID:", storm_id))
  }
}
visualize_hurricane_track("Laura2020", New_Hurricane_last)
# Function to visualize a hurricane's track
visualize_hurricane_track <- function(storm_id, data) {
  # Filter the data for the specified STORM_ID (case-insensitive comparison)
  hurricane <- data %>% filter(toupper(STORM_ID) == toupper(storm_id))

  # Check if any rows were found for the hurricane
  if (nrow(hurricane) > 0) {
    # Initialize an empty mapview object
    map_hurricane <- mapview()

    # Loop through the hurricane rows (should usually be 1 row, but robust)
    for (i in seq_len(nrow(hurricane))) {
      
      
      
      # Check for and add line data
      if (!is.null(hurricane$LinesSF[[i]]) && inherits(hurricane$LinesSF[[i]], 'sf')) {
        map_hurricane <- map_hurricane +
          mapview(hurricane$LinesSF[[i]], color = "blue",
                  layer.name = paste(hurricane$Name[i], "Track"))
      } else {
        message(paste("No valid line spatial data for", hurricane$Name[i],
                      "(", hurricane$STORM_ID[i], ")"))
      }
      
      # Check for and add point data
      if (!is.null(hurricane$PointsSF[[i]]) && inherits(hurricane$PointsSF[[i]], 'sf')) {
        map_hurricane <- map_hurricane +
          mapview(hurricane$PointsSF[[i]], col.region = "red",
                  layer.name = paste(hurricane$Name[i], "Points"))
      } else {
        message(paste("No valid point spatial data for", hurricane$Name[i],
                      "(", hurricane$STORM_ID[i], ")"))
      }
      
      # Check for and add Windwath data
      if (!is.null(hurricane$WindwathSF[[i]]) && inherits(hurricane$WindwathSF[[i]], 'sf')) {
        map_hurricane <- map_hurricane +
          mapview(hurricane$WindwathSF[[i]], col.region = "yellow",
                  layer.name = paste(hurricane$Name[i], "Points"))
      } else {
        message(paste("No valid Windwath spatial data for", hurricane$Name[i],
                      "(", hurricane$STORM_ID[i], ")"))
      }
      
    }

    # Display the map
    print(map_hurricane) # Use print() to ensure it displays correctly in all environments

  } else {
    message(paste("No data found for hurricane with STORM_ID:", storm_id))
  }
}
visualize_hurricane_track("Laura2020", New_Hurricane_last)
library(tidyverse)
library(sf)
library(mapview)
library(stringr)

# Function to visualize hurricane tracks by name and year
visualize_hurricane_track_by_name_and_year <- function(hurricane_name, data) {
  # Filter the data for hurricanes matching the name (case-insensitive)
  hurricanes <- data %>%
    filter(str_to_lower(Name) == str_to_lower(hurricane_name))

  # Check if any hurricanes were found
  if (nrow(hurricanes) > 0) {
    if (nrow(hurricanes) > 1) {
      # Prompt user to select year
      cat("Found multiple hurricanes with the name '", hurricane_name, "'. Please select the year:\n")
      for (i in seq_len(nrow(hurricanes))) {
        cat(i, ": ", hurricanes$Year[i], "\n", sep = "")
      }
      year_choice <- as.integer(readline(prompt = "Enter the number corresponding to the year: "))

      # Validate user input
      if (year_choice %in% seq_len(nrow(hurricanes))) {
        selected_hurricane <- hurricanes[year_choice, ]
        message(paste("Visualizing track for", selected_hurricane$Name, "(", selected_hurricane$Year, ")"))
        map_hurricane <- mapview()

        # Check and add line data
        if (!is.null(selected_hurricane$LinesSF[[1]]) && inherits(selected_hurricane$LinesSF[[1]], 'sf')) {
          map_hurricane <- map_hurricane +
            mapview(selected_hurricane$LinesSF[[1]], color = "blue",
                    layer.name = paste(selected_hurricane$Name, selected_hurricane$Year, "Track"))
        } else {
          message(paste("No valid line spatial data for", selected_hurricane$Name,
                        "(", selected_hurricane$STORM_ID, ")"))
        }

        # Check and add point data
        if (!is.null(selected_hurricane$PointsSF[[1]]) && inherits(selected_hurricane$PointsSF[[1]], 'sf')) {
          map_hurricane <- map_hurricane +
            mapview(selected_hurricane$PointsSF[[1]], col.region = "red",
                    layer.name = paste(selected_hurricane$Name, selected_hurricane$Year, "Points"))
        } else {
          message(paste("No valid point spatial data for", selected_hurricane$Name,
                        "(", selected_hurricane$STORM_ID, ")"))
        }

        # Check and add Windwath data
        if (!is.null(selected_hurricane$WindwathSF[[1]]) && inherits(selected_hurricane$WindwathSF[[1]], 'sf')) {
          map_hurricane <- map_hurricane +
            mapview(selected_hurricane$WindwathSF[[1]], col.region = "yellow",
                    layer.name = paste(selected_hurricane$Name, selected_hurricane$Year, "Windwath"))
        } else {
          message(paste("No valid Windwath spatial data for", selected_hurricane$Name,
                        "(", selected_hurricane$STORM_ID, ")"))
        }
        print(map_hurricane) #Print Map
      } else {
        message("Invalid year choice.")
      }
    } else {
      # If only one hurricane is found, proceed as before
      selected_hurricane <- hurricanes[1,]
      message(paste("Visualizing track for", selected_hurricane$Name, "(", selected_hurricane$Year, ")"))
      map_hurricane <- mapview()
      # Check and add line data
        if (!is.null(selected_hurricane$LinesSF[[1]]) && inherits(selected_hurricane$LinesSF[[1]], 'sf')) {
          map_hurricane <- map_hurricane +
            mapview(selected_hurricane$LinesSF[[1]], color = "blue",
                    layer.name = paste(selected_hurricane$Name, selected_hurricane$Year, "Track"))
        } else {
          message(paste("No valid line spatial data for", selected_hurricane$Name,
                        "(", selected_hurricane$STORM_ID, ")"))
        }

        # Check and add point data
        if (!is.null(selected_hurricane$PointsSF[[1]]) && inherits(selected_hurricane$PointsSF[[1]], 'sf')) {
          map_hurricane <- map_hurricane +
            mapview(selected_hurricane$PointsSF[[1]], col.region = "red",
                    layer.name = paste(selected_hurricane$Name, selected_hurricane$Year, "Points"))
        } else {
          message(paste("No valid point spatial data for", selected_hurricane$Name,
                        "(", selected_hurricane$STORM_ID, ")"))
        }

        # Check and add Windwath data
        if (!is.null(selected_hurricane$WindwathSF[[1]]) && inherits(selected_hurricane$WindwathSF[[1]], 'sf')) {
          map_hurricane <- map_hurricane +
            mapview(selected_hurricane$WindwathSF[[1]], col.region = "yellow",
                    layer.name = paste(selected_hurricane$Name, selected_hurricane$Year, "Windwath"))
        } else {
          message(paste("No valid Windwath spatial data for", selected_hurricane$Name,
                        "(", selected_hurricane$STORM_ID, ")"))
        }
      print(map_hurricane)
    }
  } else {
    message(paste("No hurricanes found with the name '", hurricane_name, "'."))
  }
}

# --- Interactive Part ---

# Assuming 'New_Hurricane_last' is your dataframe
user_input <- readline(prompt = "Enter the name of a hurricane to visualize: ")
## Enter the name of a hurricane to visualize:
# Check if the user entered something
if (str_trim(user_input) != "") {
  visualize_hurricane_track_by_name_and_year(user_input, New_Hurricane_last)
} else {
  message("No hurricane name entered.")
}
## No hurricane name entered.

EXTRA

library(tidyverse) library(sf) library(mapview) library(stringr)

Load hurricane data and metadata (assuming these are already loaded as in your script)

Load CSV

Trial_hurricane_data <- read_csv(“./F_Project/HURRICANE/use_Hurricanes.csv”) %>% mutate( Name = str_trim(Name), Year = as.integer(Year), STORM_ID = str_trim(toupper(STORM_ID)) )

meta_path <- “./F_Project/HURRICANE/L_shapefile_meta_Hurricanes.csv”

Load metadata

meta_data <- read_csv(“./F_Project/HURRICANE/L_shapefile_meta_Hurricanes.csv”) %>% mutate( STORM_ID = str_trim(toupper(STORM_ID)) # Ensure uniform formatting )

Merge metadata and hurricane data

Trial_combined_meta <- meta_data %>% left_join(Trial_hurricane_data, by = “STORM_ID”)

Define base path and generate shapefile paths

base_path <- “./F_Project/HURRICANE/hurricanes”

Trial_combined_meta <- Trial_combined_meta %>% mutate( LinPath = file.path(base_path, Folder, LinFile), PtsPath = file.path(base_path, Folder, PtsFile), WindsPath = file.path(base_path, Folder, windswathFile) )

Function to safely read shapefiles

safe_read_sf <- function(path) { if (!is.na(path) && file.exists(path)) { tryCatch(read_sf(path), error = function(e) NULL) } else { NULL } }

Add shapefiles as spatial data columns

Trial_New_Hurricane_base <- Trial_combined_meta %>% mutate( PointsSF = map(PtsPath, safe_read_sf), LinesSF = map(LinPath, safe_read_sf), WindwathSF = map(WindsPath, safe_read_sf) )

Load US states shapefile

us_states <- read_sf(“C:/Users/Alijah Anyagwosi/Downloads/SPRING 2025/RStudio/F_Project/HURRICANE/Contiguos_US/Contiguos_US.shp”)

Function to visualize hurricane track with intersecting states

visualize_hurricane_track_with_states <- function(storm_id, data, states_sf) { # Added states_sf # Filter hurricane data hurricane <- data %>% filter(toupper(STORM_ID) == toupper(storm_id))

# Check if data exists if (nrow(hurricane) > 0) { map_hurricane <- mapview()

for (i in seq_len(nrow(hurricane))) {
  # Add hurricane track
  if (!is.null(hurricane$LinesSF[[i]]) && inherits(hurricane$LinesSF[[i]], 'sf')) {
    # Ensure CRS consistency
    if (st_crs(hurricane$LinesSF[[i]]) != st_crs(states_sf)) {
      states_sf <- st_transform(states_sf, crs = st_crs(hurricane$LinesSF[[i]]))
    }
    map_hurricane <- map_hurricane +
      mapview(hurricane$LinesSF[[i]], color = "blue", layer.name = paste(hurricane$Name[i], "Track"))
  }

  #check for wind swath and add
  if (!is.null(hurricane$WindwathSF[[i]]) && inherits(hurricane$WindwathSF[[i]], 'sf')) {
     # Ensure CRS consistency
    if (st_crs(hurricane$WindwathSF[[i]]) != st_crs(states_sf)) {
      states_sf <- st_transform(states_sf, crs = st_crs(hurricane$WindwathSF[[i]]))
    }
    map_hurricane <- map_hurricane +
      mapview(hurricane$WindwathSF[[i]], col.region = "yellow",
              layer.name = paste(hurricane$Name[i], "Windswath"))
  }

  # Find intersecting states.  Use the LinesSF for intersection
  if (!is.null(hurricane$LinesSF[[i]]) && inherits(hurricane$LinesSF[[i]], 'sf')) {
    # Ensure CRS consistency
    if (st_crs(hurricane$LinesSF[[i]]) != st_crs(states_sf)) {
      states_sf <- st_transform(states_sf, crs = st_crs(hurricane$LinesSF[[i]]))
    }
    intersecting_states <- st_intersection(states_sf, hurricane$LinesSF[[i]])
    if (nrow(intersecting_states) > 0) {
      # Add intersecting states with a highlight color
      map_hurricane <- map_hurricane +
        mapview(intersecting_states, col.region = "red", layer.name = "Affected States")

      # Extract and print state names
      state_names <- intersecting_states$NAME # Assuming NAME is the column with state names
      message(paste("States intersected by", hurricane$Name[i], ":", paste(state_names, collapse = ", ")))
    } else {
      message(paste("No states intersected by", hurricane$Name[i]))
    }
  }
}
print(map_hurricane)

} else { message(paste(“No data found for hurricane with STORM_ID:”, storm_id)) } }

Example usage:

visualize_hurricane_track_with_states(“DEBBY2024”, Trial_New_Hurricane_base, us_states) #Added us_states

#colnames(combined_meta)

map_strongest <- mapview()

for (i in seq_len(nrow(strongest_hurricanes))) { if (!is.null(strongest_hurricanes\(LinesSF[[i]])) { map_strongest <- map_strongest + mapview(strongest_hurricanes\)LinesSF[[i]], color = “darkorange”) } if (!is.null(strongest_hurricanes\(PointsSF[[i]])) { map_strongest <- map_strongest + mapview(strongest_hurricanes\)PointsSF[[i]], col.region = “black”) } }

map_strongest

Load shapefiles

#counties_sf <- sf::st_read(“./F_Project/HURRICANE/USA/USA.shp”)

View them separately

#mapview(counties_sf, col.region = “white”, alpha.region = 0.5, color = “black”)

Load shapefiles

#usa_sf <- sf::st_read(“./F_Project/HURRICANE/USA/cb_2018_us_county_500k/cb_2018_us_county_500k.shp”)

View them separately

#mapview(usa_sf, col.region = “lightgray”, alpha.region = 0.3)

FIX

Histogram of Deaths

ggplot(hurricane_data, aes(x = Deaths)) + geom_histogram(binwidth = 5, fill = “darkgray”, color = “white”) + labs(title = “Distribution of Hurricane-Related Deaths”, x = “Number of Deaths”, y = “Frequency”) + theme_minimal()

Box Plots: Deaths by Highest Category

ggplot(hurricane_data %>% filter(!is.na(Highest_Category)), aes(x = Highest_Category, y = Deaths, fill = Highest_Category)) + geom_boxplot() + labs(title = “Hurricane-Related Deaths by Highest Category”, x = “Highest Category”, y = “Number of Deaths”) + theme_minimal() + theme(legend.position = “none”)

Convert ‘Damage’ to numeric (remove commas and handle “None”)

hurricane_data_clean <- hurricane_data %>% mutate( Damage_Numeric = as.numeric(gsub(“,”, ““, ifelse(Damage ==”None”, “0”, Damage))) )

Box Plots: Damage by Highest Category

ggplot(hurricane_data_clean %>% filter(!is.na(Highest_Category)), aes(x = Highest_Category, y = Damage_Numeric, fill = Highest_Category)) + geom_boxplot() + labs(title = “Hurricane Damage by Highest Category”, x = “Highest Category”, y = “Damage (USD)”) + theme_minimal() + scale_y_continuous(labels = scales::comma) + theme(legend.position = “none”)

hurricane_data_dates <- hurricane_data %>% mutate( Peak_Month = month(mdy(Peak_Date), label = TRUE) # Convert to date and extract month )

ggplot(hurricane_data_dates, aes(x = Peak_Month, fill = Peak_Month)) + geom_bar() + labs(title = “Number of Hurricanes by Month of Peak Intensity”, x = “Month of Peak Intensity”, y = “Number of Hurricanes”) + theme_minimal() + theme(legend.position = “none”)

— WARNING: This can be memory intensive for many hurricanes —

Attempt to combine all line shapefiles into a single sf object

Filter out NULL entries first

Need a base map of the US coastline/states

Get US state boundaries using maps or sf built-in data if available

library(maps) # if using base maps

us_states <- st_as_sf(map(“state”, fill=TRUE, plot = FALSE))

us_states <- st_transform(us_states, st_crs(all_hurricane_paths)) # Match CRS

Or use sf data if available and appropriate CRS

Example using built-in data, may need installation/download

us_states <- sf::st_as_sf(maps::map(“state”, fill=TRUE, plot = FALSE)) # Common way

Example: Plotting paths colored by Saffir-Simpson Category (if available in path data, which it might not be)

If category isn’t per-point in the shapefile, you’d need to join or color the whole line by max category

This requires more complex data wrangling to carry storm metadata into the combined sf object

Note on spatial plotting: Plotting all paths can be messy.

Consider plotting subsets (e.g., all Category 4+ storms, all storms in a specific decade).

This requires filtering New_Hurricane_stats before extracting and combining LinesSF.

valid_lines_sf <- New_Hurricane_stats\(LinesSF[!sapply(New_Hurricane_stats\)LinesSF, is.null)]

Check if there are valid spatial objects to combine

if(length(valid_lines_sf) > 0) { # Use do.call(rbind, …) or reduce(rbind, …) to combine sf objects # reduce is often safer with potential CRS issues all_hurricane_paths <- do.call(rbind, valid_lines_sf)

A simpler approach might be to get a low-res world map and filter to US region

library(rnaturalearth) library(rnaturalearthdata) world <- ne_countries(scale = “medium”, returnclass = “sf”) us_boundary <- world %>% filter(sovereignt == “United States of America”)

Ensure CRS matches - this is CRITICAL for plotting layers together

Assuming hurricane paths are likely in WGS84 (EPSG:4326)

target_crs <- st_crs(all_hurricane_paths) us_boundary <- st_transform(us_boundary, target_crs)

Static Plot: All Hurricane Paths

ggplot() + geom_sf(data = us_boundary, fill = “lightgray”, color = “white”) + geom_sf(data = all_hurricane_paths, aes(color = factor(Year)), linewidth = 0.5, alpha = 0.5) + # Color by year, adjust alpha/linewidth labs(title = “All US-Affected Hurricane Paths”, color = “Year”) + theme_minimal() + coord_sf(xlim = c(-130, -60), ylim = c(15, 50)) # Adjust limits to focus on relevant area

} else { message(“No valid hurricane line shapefiles found to plot an overview map.”) }

map_longest <- mapview()

for (i in seq_len(nrow(longest_duration))) { if (!is.null(longest_duration\(LinesSF[[i]])) { map_longest <- map_longest + mapview(longest_duration\)LinesSF[[i]], color = “purple”) } if (!is.null(longest_duration\(PointsSF[[i]])) { map_longest <- map_longest + mapview(longest_duration\)PointsSF[[i]], col.region = “green”) } }

map_longest

Assuming ‘New_Hurricane_last’ dataframe exists from your previous steps

Find the hurricane(s) with the longest duration

longest_duration <- New_Hurricane_last %>% filter(Duration_h == max(Duration_h, na.rm = TRUE))

Initialize an empty mapview object

map_longest <- mapview()

Loop through the hurricane(s) with the longest duration

for (i in seq_len(nrow(longest_duration))) { hurricane_name <- longest_duration$Name[i] # Get the name for layer naming

# Check if spatial data (LinesSF) exists for this hurricane and add to the map if (!is.null(longest_duration\(LinesSF[[i]]) && inherits(longest_duration\)LinesSF[[i]], ‘sf’)) { map_longest <- map_longest + mapview(longest_duration$LinesSF[[i]], color = “purple”, layer.name = hurricane_name) } else { message(paste(“No valid line spatial data for”, hurricane_name)) }

# Check if point spatial data (PointsSF) exists and add to the map if (!is.null(longest_duration\(PointsSF[[i]]) && inherits(longest_duration\)PointsSF[[i]], ‘sf’)) { map_longest <- map_longest + mapview(longest_duration$PointsSF[[i]], col.region = “green”, layer.name = paste(hurricane_name, “Points”)) } else { message(paste(“No valid point spatial data for”, hurricane_name)) }

# Check if spatial data (WindwathSF) exists for this hurricane and add to the map if (!is.null(longest_duration\(WindwathSF[[i]]) && inherits(longest_duration\)WindwathSF[[i]], ‘sf’)) { map_longest <- map_longest + mapview(longest_duration$WindwathSF[[i]], color = “orange”, layer.name = paste(hurricane_name, “Windwath”)) } else { message(paste(“No valid Windwath spatial data for”, hurricane_name)) } }

Display the map

map_longest

==========================================================