NFLD ports.R

library(tidyverse)
library(dplyr)
library(gganimate)
library(gifski)
library(wesanderson)


#home port data
hp <- read.csv('./data/home port.csv')
head(hp)
summary(hp)
structure(hp)
unique(hp$Species.Code)
unique(hp$Home.Port.Code)

# Drop Labrador ports
ports_to_drop <- c("Makkovik", "Black Tickle", "Square Islands", 
                   "Williams Harbour", "St. Lewis", "Red Bay", 
                   "Capstan Island", "L'Anse Au Loup", "Forteau")

# Filter out the rows for these ports
hp <- hp %>%
  filter(!(Home.Port %in% ports_to_drop)) %>%
  filter(Landing.Year < 2024)

#filter out ports that have no records
hp_no_zeros <- hp %>%
  filter(!(Home.Fishers == 0 & Home.Vessels == 0))

#Identify ports that have a nonzero sum of fishers and vessels
ports_with_any_activity <- hp %>%
  group_by(Home.Port) %>%
  summarise(
    total_fishers = sum(Home.Fishers, na.rm = TRUE),
    total_vessels = sum(Home.Vessels, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  # Keep only ports that have at least some fishers AND some vessels
  filter(total_fishers > 0 & total_vessels > 0)

#Filter the main dataset to keep only those ports
hp_filtered <- hp %>%
  filter(Home.Port %in% ports_with_any_activity$Home.Port) 


# summarize total fishers by year
fishers_by_year <- hp_filtered %>%
  group_by(Landing.Year) %>%
  summarise(total_fishers = sum(Home.Fishers, na.rm = TRUE)) %>%
  ungroup()

head(fishers_by_year)


# summarize total vessels by year
vessels_by_year <- hp_filtered %>%
  group_by(Landing.Year) %>%
  summarise(total_vessels = sum(Home.Vessels, na.rm = TRUE)) %>%
  ungroup()

#filter ports that don't have at least 10 observations

#Identify ports that have at least 10 observations
ports_at_least_10 <- hp %>%
  group_by(Home.Port) %>%
  summarise(n_obs = n(), .groups = "drop") %>%
  filter(n_obs >= 10)

#Filter the main dataset
hp_filtered10 <- hp %>%
  filter(Home.Port %in% ports_at_least_10$Home.Port)

# summarize total ports by year
ports_by_year <- hp_filtered %>%
  group_by(Landing.Year) %>%
  summarise(total_ports = n_distinct(Home.Port, na.rm = TRUE)) %>%
  ungroup()


# Count how many *distinct* species are being fished in each port per year
hp_port_year_species <- hp_filtered %>%
  group_by(Landing.Year, Home.Port) %>%
  filter(Landing.Year < 2024) %>%
  summarise(
    num_species = n_distinct(Species.Name),
    .groups = "drop"
  )

# Summarize by year: average number of species across ports
hp_avg_species_year <- hp_port_year_species %>%
  group_by(Landing.Year) %>%
  summarise(
    avg_num_species = mean(num_species),
    .groups = "drop"
  )


###  TOP 25 PORTS BY VESSELS  ###

#filter top 25 ports by vessel number over entire time series

# Compute the total number of vessels per port
ports_by_vessels <- hp %>%
  group_by(Home.Port) %>%
  summarise(total_vessels = sum(Home.Vessels, na.rm = TRUE), .groups = "drop")

# Sort the ports in descending order by total vessels
ports_by_vessels <- ports_by_vessels %>%
  arrange(desc(total_vessels))

# Select the top 25 ports
ports_top_25 <- ports_by_vessels %>%
  slice_head(n = 25)
#export as csv
write.csv(ports_top_25, "ports_top_25.csv", row.names = FALSE)

# Filter the main dataset to keep only rows from these top 25 ports
hp_filtered25 <- hp %>%
  filter(Home.Port %in% ports_top_25$Home.Port)

# Combine groundfish and seal groups
hp_filtered25 <- hp %>%
  filter(Home.Port %in% ports_top_25$Home.Port) %>%
  mutate(
    Species.Name = case_when(
      str_detect(Species.Name, regex("Groundfish", ignore_case = TRUE)) ~ "Groundfish",
      str_detect(Species.Name, regex("Seal",       ignore_case = TRUE)) ~ "Harp Seal",
      str_detect(Species.Name, regex("Shrimp",       ignore_case = TRUE)) ~ "Shrimp",
      TRUE ~ as.character(Species.Name)
    ),
    Species.Name = factor(Species.Name, levels = sort(unique(Species.Name)))  #alphabetical
  )

#rename port
hp_filtered25 <- hp_filtered25 %>%
  mutate(Home.Port = ifelse(Home.Port == "Little Harbour East, Placentia Bay",
                            "Little Harbour East",
                            Home.Port))

#if Home.Port is a factor, drop/refresh levels
hp_filtered25 <- hp_filtered25 %>%
  mutate(Home.Port = factor(Home.Port))

#alphabetical factor order for legend/colour mapping
species_levels <- sort(unique(hp_filtered25$Species.Name))
hp_filtered25 <- hp_filtered25 %>%
  mutate(Species.Name = factor(Species.Name, levels = species_levels))

#make n discrete colors from the zissou1 palette (interpolated)
n_species <- length(species_levels)
zissou_cols <- colorRampPalette(wesanderson::wes_palette("Zissou1", 5))(n_species)
names(zissou_cols) <- species_levels

#ensure Species.Name is alphabetical (after your Groundfish/Seal collapsing)
species_levels <- sort(unique(as.character(hp_filtered25$Species.Name)))
hp_filtered25 <- hp_filtered25 %>%
  mutate(Species.Name = factor(Species.Name, levels = species_levels))

n_species <- length(species_levels)
k_cols <- ceiling(n_species / 2)

#make k zissou colours, then reuse each colour twice
zissou_base <- colorRampPalette(wesanderson::wes_palette("Zissou1", 5))(k_cols)

species_cols <- rep(zissou_base, each = 2)[seq_len(n_species)]
names(species_cols) <- species_levels

#solid/dashed for alternating species
species_lty <- rep(c("solid", "dashed"), times = k_cols)[seq_len(n_species)]
names(species_lty) <- species_levels

p <- ggplot(hp_filtered25, aes(x = Landing.Year, y = Home.Vessels, group = Species.Name)) + 
  geom_line(aes(colour = Species.Name, linetype = Species.Name)) +
  scale_colour_manual(values = species_cols, breaks = species_levels, name = "Species") +
  scale_linetype_manual(values = species_lty, breaks = species_levels, name = "Species") +
  labs(x = "Year", y = "Number of Vessels") +
  theme_bw() +
  facet_wrap(vars(Home.Port)) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    strip.background = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA),
    axis.text.x = element_text(angle = 300)
  ) +
  guides(
    linetype = "none",
    colour = guide_legend(
      ncol = 1,
      override.aes = list(
        linetype = unname(species_lty),  #make dashed show in legend
        size = 1
      )
    )
  )
p

ggsave("plots/Figure 4.jpg", dpi=300, width = 25, height = 25, units = "cm")


#repeat for all ports

#use all ports (no top-25 filter)
hp_allports <- hp_filtered10 %>%
  mutate(
    Species.Name = case_when(
      stringr::str_detect(Species.Name, stringr::regex("Groundfish", ignore_case = TRUE)) ~ "Groundfish",
      stringr::str_detect(Species.Name, stringr::regex("Seal",       ignore_case = TRUE)) ~ "Harp Seal",
      TRUE ~ as.character(Species.Name)
    ),
    Home.Port = ifelse(Home.Port == "Little Harbour East, Placentia Bay",
                       "Little Harbour East",
                       as.character(Home.Port))
  )

hp_allports <- hp_allports %>%
  mutate(
    Home.Port = trimws(as.character(Home.Port)),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Goose Cove", Home.Port, ignore.case = TRUE), "Goose Cove", Home.Port),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("North Harbour", Home.Port, ignore.case = TRUE), "North Harbour", Home.Port),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Petty Harbour", Home.Port, ignore.case = TRUE), "Petty Harbour", Home.Port),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Sally's Cove", Home.Port, ignore.case = TRUE), "Sally's Cove", Home.Port),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Sallys Cove",  Home.Port, ignore.case = TRUE), "Sally's Cove", Home.Port), #covers missing apostrophe
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Seal Cove", Home.Port, ignore.case = TRUE), "Seal Cove", Home.Port),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Wild Cove", Home.Port, ignore.case = TRUE), "Wild Cove", Home.Port),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Little Harbour East", Home.Port, ignore.case = TRUE), "Little Harbour East", Home.Port),
    Home.Port = factor(Home.Port)  #refresh factor levels
  )

#alphabetical levels for per-species legend
species_levels <- sort(unique(hp_allports$Species.Name))
hp_allports <- hp_allports %>%
  mutate(Species.Name = factor(Species.Name, levels = species_levels))

#reuse each zissou colour twice + solid/dashed alternating
n_species <- length(species_levels)
k_cols <- ceiling(n_species / 2)

zissou_base <- colorRampPalette(wesanderson::wes_palette("Zissou1", 5))(k_cols)
species_cols <- rep(zissou_base, each = 2)[seq_len(n_species)]
names(species_cols) <- species_levels

species_lty <- rep(c("solid", "dashed"), times = k_cols)[seq_len(n_species)]
names(species_lty) <- species_levels

p_allports <- ggplot(hp_allports, aes(x = Landing.Year, y = Home.Vessels, group = Species.Name)) + 
  geom_line(aes(colour = Species.Name, linetype = Species.Name)) +
  scale_colour_manual(values = species_cols, breaks = species_levels, name = "Species") +
  scale_linetype_manual(values = species_lty, breaks = species_levels, name = "Species") +
  labs(x = "Year", y = "Number of Vessels") +
  theme_bw() +
  facet_wrap(vars(Home.Port)) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    strip.background = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA),
    axis.text.x = element_text(angle = 300)
  ) +
  guides(
    linetype = "none",
    colour = guide_legend(
      ncol = 1,
      override.aes = list(linetype = unname(species_lty), size = 1)
    )
  )
p_allports

ggsave("plots/Figure S2.jpg", dpi=300, width = 40, height = 40, units = "cm")


# Count how many *distinct* species are being fished in each port per year
hp_port_year_species25 <- hp_filtered25 %>%
  group_by(Landing.Year, Home.Port) %>%
  filter(Landing.Year < 2024) %>%
  summarise(
    num_species = n_distinct(Species.Name),
    .groups = "drop"
  )

# Summarize by year: average number of species across ports
hp_avg_species_year25 <- hp_port_year_species25 %>%
  group_by(Landing.Year) %>%
  summarise(
    avg_num_species = mean(num_species),
    .groups = "drop"
  )



###  NUMBER OF SPECIES CAUGHT BY PORT FOR TOP 25 PORTS  ###


# Count how many *distinct* species are being fished in each port per year
hp_port_year_species25 <- hp_filtered25 %>%
  group_by(Landing.Year, Home.Port) %>%
  filter(Landing.Year < 2024) %>%
  summarise(
    num_species = n_distinct(Species.Name),
    .groups = "drop"
  )


# Summarize by year: average number of species across ports
hp_avg_species_year25 <- hp_port_year_species25 %>%
  group_by(Landing.Year) %>%
  summarise(
    avg_num_species = mean(num_species),
    .groups = "drop"
  )


#calculate the number of home ports through time

# Summarize mean (or total) fishers & vessels by landing year
year_summary <- hp %>%
  group_by(Landing.Year) %>%
  summarise(
    mean_fishers = mean(Home.Fishers),
    sum_fishers  = sum(Home.Fishers),
    mean_vessels = mean(Home.Vessels),
    sum_vessels  = sum(Home.Vessels),
    .groups = "drop"
  )

head(year_summary)


port_summary <- hp %>%
  group_by(Home.Port) %>%
  summarise(
    mean_fishers = mean(Home.Fishers),
    sum_fishers  = sum(Home.Fishers),
    mean_vessels = mean(Home.Vessels),
    sum_vessels  = sum(Home.Vessels),
    .groups = "drop"
  )

head(port_summary)


species_summary <- hp %>%
  group_by(Species.Code, Species.Name) %>%
  summarise(
    mean_fishers = mean(Home.Fishers),
    sum_fishers  = sum(Home.Fishers),
    .groups = "drop"
  )

head(species_summary)


###  NUMBER OF FISHERS BY YEAR  ###


# 1. Summarize data by year and port
hp_port_year <- hp %>%
  group_by(Landing.Year, Home.Port) %>%
  summarise(
    TotalFishers = sum(Home.Fishers, na.rm = TRUE),
    .groups = "drop"
  )



###  TOTAL NUMBER OF FISHERS PER YEAR IN TOP 25 BIGGEST PORTS  ###

year_summary25 <- hp_filtered25 %>%
  group_by(Landing.Year) %>%
  summarise(
    mean_fishers = mean(Home.Fishers),
    sum_fishers  = sum(Home.Fishers),
    mean_vessels = mean(Home.Vessels),
    sum_vessels  = sum(Home.Vessels),
    .groups = "drop"
  )

head(year_summary25)


port_summary25 <- hp_filtered25 %>%
  group_by(Home.Port) %>%
  summarise(
    mean_fishers = mean(Home.Fishers),
    sum_fishers  = sum(Home.Fishers),
    mean_vessels = mean(Home.Vessels),
    sum_vessels  = sum(Home.Vessels),
    .groups = "drop"
  )

head(port_summary25)


species_summary25 <- hp_filtered25 %>%
  group_by(Species.Code, Species.Name) %>%
  summarise(
    mean_fishers = mean(Home.Fishers),
    sum_fishers  = sum(Home.Fishers),
    .groups = "drop"
  )

head(species_summary25)



###  NUMBER OF FISHERS BY YEAR BY PORT  ###


# Summarize data by year and port
hp_port_year <- hp %>%
  group_by(Landing.Year, Home.Port) %>%
  summarise(
    TotalFishers = sum(Home.Fishers, na.rm = TRUE),
    .groups = "drop"
  )



###  NUMBER OF PORTS FISHING THROUGH TIME  ###


# Sum Home.Fishers per port-year
hp_port_year <- hp %>%
  group_by(Landing.Year, Home.Port) %>%
  summarise(
    total_fishers = sum(Home.Fishers, na.rm = TRUE),
    .groups = "drop"
  )

# For each year, count how many ports have >= 1 fisher
port_count_by_year <- hp_port_year %>%
  filter(total_fishers >= 1) %>%
  filter(Landing.Year < 2024) %>%
  group_by(Landing.Year) %>%
  summarise(
    n_ports_with_fishers = n(),
    .groups = "drop"
  )


### NUMBER OF PORTS FISHING EACH SPECIES  ###


# Count the number of distinct ports fishing each species in each year
hp_species_year <- hp %>%
  group_by(Landing.Year, Species.Name) %>%
  summarise(
    num_ports = n_distinct(Home.Port),  # how many different ports fished this species
    .groups = "drop"
  )


# Count how many *distinct* species are being fished in each port per year
hp_port_year_species <- hp %>%
  group_by(Landing.Year, Home.Port) %>%
  filter(Landing.Year < 2024) %>%
  summarise(
    num_species = n_distinct(Species.Name),
    .groups = "drop"
  )


# Summarize by year: average number of species across ports
hp_avg_species_year <- hp_port_year_species %>%
  group_by(Landing.Year) %>%
  summarise(
    avg_num_species = mean(num_species),
    .groups = "drop"
  )


# Count how many species are being fished in each port per year
hp_port_year_species <- hp %>%
  group_by(Landing.Year, Home.Port) %>%
  filter(Landing.Year < 2024) %>%
  summarise(
    num_species = n(Species.Name),
    .groups = "drop"
  )


# Summarize by year: average number of species across ports
hp_avg_species_year <- hp_port_year_species %>%
  group_by(Landing.Year) %>%
  summarise(
    avg_num_species = mean(num_species),
    .groups = "drop"
  )


# Count the number of distinct species harvested for each port-year
hp_port_year_species <- hp %>%
  group_by(Landing.Year, Home.Port) %>%
  summarise(
    num_species = n_distinct(Species.Name),  # or Species.Code if that’s your identifier
    .groups = "drop"
  )

# For each year, find the maximum number of species harvested among all ports
hp_max_species_year <- hp_port_year_species %>%
  filter(Landing.Year < 2024) %>%
  group_by(Landing.Year) %>%
  summarise(
    max_species = max(num_species),
    .groups = "drop"
  )


#species list
species_levels <- hp %>%
  mutate(Species.Name = as.character(Species.Name)) %>%
  filter(!is.na(Species.Name)) %>%
  distinct(Species.Name) %>%
  arrange(Species.Name)

write.csv(species_levels, "data/Species list.csv", row.names = FALSE)



###  MAP PORTS  ###

lats <- read.csv('./data/NLCE Ports lat lons.csv')
summary(lats)
head(lats)

# Perform a left join, keeping all rows from 'hp' and matching lat/long from 'lats'
hp_with_coords <- hp %>%
  left_join(lats, by = c("Home.Port" = "Port.Desc"))

# Check the joined data
head(hp_with_coords)

# Convert coordinates to decimal degrees
hp_with_coords$Port.Latitude <- hp_with_coords$Port.Latitude /100
hp_with_coords$Port.Longitude <- hp_with_coords$Port.Longitude /-100
head(hp_with_coords)
structure(hp_with_coords)


# Function to convert a coordinate in degree minute "DD.MM" format to decimal degrees
dm_to_dd <- function(dm) {
  # work with the absolute value first
  dm_abs <- abs(dm)
  
  # extract degrees (integer part)
  deg <- floor(dm_abs)
  
  # extract minutes: the fractional part multiplied by 100
  minutes <- (dm_abs - deg) * 100
  
  # convert to decimal degrees
  dd <- deg + minutes / 60
  
  # reapply the sign of the original value
  dd * sign(dm)
}


#convert degree minute to decimal degrees
hp_with_coords <- hp_with_coords %>%
  mutate(
    Port.Latitude.dec = dm_to_dd(Port.Latitude),
    Port.Longitude.dec = dm_to_dd(Port.Longitude)
  )

head(hp_with_coords)


#standardize port locations
hp_with_coords_std <- hp_with_coords %>%
  group_by(Home.Port) %>% 
  mutate(
    # For each port, compute the mean latitude and longitude
    Port.Latitude.dec = mean(Port.Latitude.dec, na.rm = TRUE),
    Port.Longitude.dec = mean(Port.Longitude.dec, na.rm = TRUE)
  ) %>%
  ungroup()

# Check the results for one port:
hp_with_coords_std %>% 
  filter(Home.Port == "Daniel's Harbour") %>% 
  select(Home.Port, Port.Latitude.dec, Port.Longitude.dec) %>% 
  head() %>%
  tail()

head(hp_with_coords_std)


#sum number of vessels for each port and year
hp_port_year <- hp_with_coords_std %>%
  group_by(Home.Port, Landing.Year) %>%
  summarise(
    total_vessels = sum(Home.Vessels, na.rm = TRUE),
    # If all rows in the group have the same coordinates, using first() is sufficient.
    Port.Latitude.dec = first(Port.Latitude.dec),
    Port.Longitude.dec = first(Port.Longitude.dec),
    .groups = "drop"
  )

# View the first few rows of the summarized data
head(hp_port_year)


#sum number of vessels by port - no year
hp_port <- hp_with_coords_std %>%
  group_by(Home.Port) %>%
  summarise(
    total_vessels_all = sum(Home.Vessels, na.rm = TRUE),
    # If all rows in the group have the same coordinates, using first() is sufficient.
    Port.Latitude.dec = first(Port.Latitude.dec),
    Port.Longitude.dec = first(Port.Longitude.dec),
    .groups = "drop"
  )



###  CORRECTED VALUE OF NL ANNUAL LANDINGS  ###

value <- read.csv('./data/NL landings value.csv')
head(value)

# Plot the corrected NL landings value by year
e <- ggplot(value, aes(x = Year, y = Adjusted.value)) +
  geom_line() +
  labs(
    x = "Year",
    y = "Adjusted $ value in millions"
  ) +
  ylim(0, 1500) +
  theme_bw() + 
  theme(
    panel.grid.major = element_blank(),  # remove major gridlines
    panel.grid.minor = element_blank())   # remove minor gridlines
e


###  PLOT LANDINGS BY GROUP  ###

landings <- read.csv('./data/NL_provincial_quantities_subtotals_1998_2023_FULL.csv')
head(landings)

library(scales)
scale_y_continuous(labels = comma)


# read data
df <- read_csv("./data/NL_provincial_quantities_subtotals_1998_2023_FULL.csv")

# reshape to long format
df_long <- df %>%
  pivot_longer(
    cols = c(groundfish_t, pelagic_other_t, shellfish_t, others_t),
    names_to = "category",
    values_to = "quantity_t"
  )

# rename categories
df_long$category <- recode(df_long$category,
                           groundfish_t = "Groundfish",
                           pelagic_other_t = "Pelagic",
                           shellfish_t = "Shellfish",
                           others_t = "Others"
)

# make Zissou1 palette for 4 categories
pal <- wes_palette("Zissou1", 100, type = "continuous")
pal4 <- pal[seq(1, 100, length.out = 4)]

# stacked area plot without outlines, with comma formatting
f <- ggplot(df_long, aes(x = year, y = quantity_t, fill = category)) +
  geom_area(alpha = 0.95) +
  theme_bw() +
  labs(
    x = "Year",
    y = "Landings (tonnes)",
    fill = NULL   # removes legend title
  ) +
  scale_fill_manual(values = pal4) +
  scale_y_continuous(labels = comma) +  # <-- comma formatting
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

f



###  MAP  ###

# Install packages if not already installed:
# install.packages(c("rnaturalearth", "rnaturalearthdata", "sf", "ggplot2", "gganimate"))

library(sf)                   # for spatial data handling
library(rnaturalearth)        # for downloading natural earth data
library(rnaturalearthdata)    # supporting data
library(ggplot2)             # for plotting
library(gganimate)           # for animation
library(ggrepel)             # for labelling

# Load coastline data as an sf object (at the "medium" scale)
coastline <- ne_coastline(scale = "medium", returnclass = "sf")

# Calculate bounding box
lon_min <- min(hp_port_year$Port.Longitude.dec, na.rm = TRUE) - 0.2
lon_max <- max(hp_port_year$Port.Longitude.dec, na.rm = TRUE) + 0.2
lat_min <- min(hp_port_year$Port.Latitude.dec,  na.rm = TRUE) - 0.2
lat_max <- max(hp_port_year$Port.Latitude.dec,  na.rm = TRUE) + 0.2


wes_palette("Zissou1")
print("Zissou1")
pal <- wes_palette("Zissou1", 5)
print(pal)


# Basic static map

pal <- wes_palette("Zissou1", 100, type = "continuous")
scale_color_gradientn(colors = pal, name = "Vessels")


map <- ggplot() +
  
  # Plot the coastline
  geom_sf(data = coastline, color = "black", fill = NA) +
  
  # Plot standardized port points, sized and colored by total vessels
  geom_point(
    data = hp_port,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec, 
        size = total_vessels_all, colour = total_vessels_all),
    alpha = 0.7
  ) +
  
  # Add port labels using geom_text_repel
  geom_text_repel(
    data = hp_port,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec, label = Home.Port),
    size = 3,
    color = "black",
    nudge_y = 0.05
  ) +
  
  # Adjust point size range and remove size legend
  scale_size(range = c(1, 8), guide = "none") +
  
  # Add the custom continuous color palette
  scale_color_gradientn(colors = pal, name = "Vessels") +
  
  # Use coord_sf with bounding box (assuming lon_min, lon_max, lat_min, lat_max are defined)
  coord_sf(
    xlim = c(lon_min, lon_max),
    ylim = c(lat_min, lat_max),
    expand = FALSE
  ) +
  
  # Labels and theme
  labs(x = "Longitude", y = "Latitude") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),  # remove major gridlines
    panel.grid.minor = element_blank())   # remove minor gridlines

map



#map with only top 25 port labelled

#rename port
hp_port <- hp_port %>%
  mutate(
    Home.Port = trimws(as.character(Home.Port)),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Goose Cove", Home.Port, ignore.case = TRUE), "Goose Cove", Home.Port),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("North Harbour", Home.Port, ignore.case = TRUE), "North Harbour", Home.Port),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Petty Harbour", Home.Port, ignore.case = TRUE), "Petty Harbour", Home.Port),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Sally's Cove", Home.Port, ignore.case = TRUE), "Sally's Cove", Home.Port),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Sallys Cove",  Home.Port, ignore.case = TRUE), "Sally's Cove", Home.Port), #covers missing apostrophe
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Seal Cove", Home.Port, ignore.case = TRUE), "Seal Cove", Home.Port),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Wild Cove", Home.Port, ignore.case = TRUE), "Wild Cove", Home.Port),
    Home.Port = ifelse(!is.na(Home.Port) & grepl("Little Harbour East", Home.Port, ignore.case = TRUE), "Little Harbour East", Home.Port),
    Home.Port = factor(Home.Port)  #refresh factor levels
  )


# get top 25 ports by total_vessels_all
hp_port_labels <- hp_port %>%
  arrange(desc(total_vessels_all)) %>%
  slice_head(n = 25)

map <- ggplot() +
  
  # coastline
  geom_sf(data = coastline, color = "black", fill = NA) +
  
  # all ports as points
  geom_point(
    data = hp_port,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec, 
        size = total_vessels_all, colour = total_vessels_all),
    alpha = 0.7
  ) +
  
  # labels only for top 25
  geom_text_repel(
    data = hp_port_labels,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec, label = Home.Port),
    size = 3,
    color = "black",
    nudge_y = 0.05
  ) +
  
  # adjust point size range and remove size legend
  scale_size(range = c(1, 8), guide = "none") +
  
  # custom continuous color palette with legend label
  scale_color_gradientn(colors = pal, name = "Vessels") +
  
  # map extent
  coord_sf(
    xlim = c(lon_min, lon_max),
    ylim = c(lat_min, lat_max),
    expand = FALSE
  ) +
  
  # labels and theme
  labs(x = "Longitude", y = "Latitude") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

map


# Newfoundland label 
nf_lon <- -56.1
nf_lat <- 48.8

map <- ggplot() +
  geom_sf(data = coastline, color = "black", fill = NA) +
  geom_point(
    data = hp_port,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec,
        size = total_vessels_all, colour = total_vessels_all),
    alpha = 0.7
  ) +
  geom_text_repel(
    data = hp_port_labels,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec, label = Home.Port),
    size = 3,
    color = "black",
    nudge_y = 0, 
      segment.color = NA   #remove label leader lines
  ) +
  annotate(
    "text",
    x = nf_lon, y = nf_lat,
    label = "Newfoundland",
    size = 7,          #bigger than port labels
    color = "black"
  ) +
  scale_size(range = c(1, 8), guide = "none") +
  scale_color_gradientn(colors = pal, name = "Vessels") +
  coord_sf(
    xlim = c(lon_min, lon_max),
    ylim = c(lat_min, lat_max),
    expand = FALSE
  ) +
  labs(x = "Longitude", y = "Latitude") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

map

ggsave("plots/Fig 1.jpg", dpi=300, width = 20, height = 20, units = "cm")
ggsave("plots/Fig 1 med.jpg", dpi=300, width = 30, height = 30, units = "cm")
ggsave("plots/Fig 1 lrg.jpg", dpi=300, width = 40, height = 40, units = "cm")


# Basic static map with no labels




# plot four maps by time period, 1998-2003; 2004-2010; 2010-2016; 2017-2023

#Create a Time Period Variable

hp_with_coords_std <- hp_with_coords_std %>%
  mutate(time_period = case_when(
    Landing.Year >= 1998 & Landing.Year <= 2003 ~ "1998-2003",
    Landing.Year >= 2004 & Landing.Year <= 2010 ~ "2004-2010",
    Landing.Year >= 2011 & Landing.Year <= 2016 ~ "2011-2016",
    Landing.Year >= 2017 & Landing.Year <= 2023 ~ "2017-2023",
    TRUE ~ NA_character_
  ))


#Sum Vessels by Port and Time Period

hp_port_period <- hp_with_coords_std %>%
  group_by(Home.Port, time_period) %>%
  summarise(
    total_vessels_all = sum(Home.Vessels, na.rm = TRUE),
    Port.Latitude.dec = first(Port.Latitude.dec),
    Port.Longitude.dec = first(Port.Longitude.dec),
    .groups = "drop"
  )


#filter NAs
hp_port_period <- hp_port_period %>%
  filter(!is.na(time_period))

write.csv(hp_port_period, "ports_ranked_period.csv", row.names = FALSE, fileEncoding = "UTF-8")


#MAP 

pal <- wes_palette("Zissou1", 100, type = "continuous")

map_faceted <- ggplot() +
  
  # Plot the coastline
  geom_sf(data = coastline, color = "black", fill = NA) +
  
  # Plot the port points (sized and colored by total vessels)
  geom_point(
    data = hp_port_period,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec, 
        size = total_vessels_all, colour = total_vessels_all),
    alpha = 0.7
  ) +
  
  # Adjust point size range and remove the size legend
  scale_size(range = c(1, 8), guide = "none") +
  
  # Use the custom continuous color palette for vessel totals
  scale_color_gradientn(colors = pal, name = "Vessels") +
  
  # Set coordinate limits (ensure lon_min, lon_max, lat_min, lat_max are defined)
  coord_sf(
    xlim = c(lon_min, lon_max),
    ylim = c(lat_min, lat_max),
    expand = FALSE
  ) +
  
  # Axis labels and theme
  labs(x = "Longitude", y = "Latitude") +
  theme_bw() +
  
  # Facet by time period and remove the facet label boxes/grey background
  facet_wrap(~ time_period, drop = TRUE) +
  theme(strip.background = element_blank(),
        strip.text = element_text(color = "black", size = 10))
map_faceted


ggsave("plots/Map of ports by period.jpg", dpi=300, width = 20, height = 20, units = "cm")
ggsave("plots/Map of ports by period med.jpg", dpi=300, width = 30, height = 30, units = "cm")
ggsave("plots/Map of ports by period lrg.jpg", dpi=300, width = 40, height = 40, units = "cm")




#no minor gridlines

map_faceted <- ggplot() +
  
  # Plot the coastline
  geom_sf(data = coastline, color = "black", fill = NA) +
  
  # Plot the port points (sized and colored by total vessels)
  geom_point(
    data = hp_port_period,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec, 
        size = total_vessels_all, colour = total_vessels_all),
    alpha = 0.7
  ) +
  
  # Remove the size legend
  scale_size(range = c(1, 8), guide = "none") +
  
  # Apply the custom continuous color palette for vessel totals
  scale_color_gradientn(colors = pal, name = "Vessels") +
  
  # Set coordinate limits (make sure lon_min, lon_max, lat_min, lat_max are defined)
  coord_sf(
    xlim = c(lon_min, lon_max),
    ylim = c(lat_min, lat_max),
    expand = FALSE
  ) +
  
  # Axis labels
  labs(x = "Longitude", y = "Latitude") +
  
  # Use theme_bw() and remove minor gridlines
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),              # Remove minor gridlines
    strip.background = element_blank(),              # Remove the background box behind facet labels
    strip.text = element_text(color = "black", size = 10)  # Keep facet labels with custom styling
  ) +
  
  # Facet the map by time period
  facet_wrap(~ time_period, drop = TRUE)

map_faceted


ggsave("plots/Map of ports by period no minor.jpg", dpi=300, width = 20, height = 20, units = "cm")
ggsave("plots/Map of ports by period med no minor.jpg", dpi=300, width = 30, height = 30, units = "cm")
ggsave("plots/Map of ports by period lrg no minor.jpg", dpi=300, width = 40, height = 40, units = "cm")


#pick the top 5 ports per time period
top_ports <- hp_port_period %>%
  group_by(time_period) %>%
  slice_max(order_by = total_vessels_all, n = 5, with_ties = FALSE) %>%
  ungroup()

#map
map_faceted <- ggplot() +
  geom_sf(data = coastline, color = "black", fill = NA) +
  geom_point(
    data = hp_port_period,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec,
        size = total_vessels_all, colour = total_vessels_all),
    alpha = 0.7
  ) +
  ggrepel::geom_text_repel(
    data = top_ports,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec, label = Home.Port),
    inherit.aes = FALSE,
    size = 3,
    seed = 1,
    max.overlaps = Inf
  ) +
  scale_size(range = c(1, 8), guide = "none") +
  scale_color_gradientn(colors = pal, name = "Vessels") +
  coord_sf(
    xlim = c(lon_min, lon_max),
    ylim = c(lat_min, lat_max),
    expand = FALSE
  ) +
  labs(x = "Longitude", y = "Latitude") +
  theme_bw() +
  theme(
    panel.grid.minor = element_blank(),     #no minor gridlines
    strip.background = element_blank(),
    strip.text = element_text(color = "black", size = 10)
  ) +
  # Use theme_bw() and remove minor gridlines
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),              # Remove minor gridlines
    strip.background = element_blank(),              # Remove the background box behind facet labels
    strip.text = element_text(color = "black", size = 10)  # Keep facet labels with custom styling
  ) +
  facet_wrap(~ time_period, drop = TRUE)

map_faceted


#increase font size
#map
map_faceted <- ggplot() +
  geom_sf(data = coastline, color = "black", fill = NA) +
  geom_point(
    data = hp_port_period,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec,
        size = total_vessels_all, colour = total_vessels_all),
    alpha = 0.7
  ) +
  ggrepel::geom_text_repel(
    data = top_ports,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec, label = Home.Port),
    inherit.aes = FALSE,
    size = 4,          #label font size (increase as needed)
    seed = 1,
    max.overlaps = Inf
  ) +
  scale_size(range = c(1, 8), guide = "none") +
  scale_color_gradientn(colors = pal, name = "Vessels") +
  coord_sf(
    xlim = c(lon_min, lon_max),
    ylim = c(lat_min, lat_max),
    expand = FALSE
  ) +
  labs(x = "Longitude", y = "Latitude") +
  theme_bw(base_size = 16) +   #increases all text globally
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(color = "black", size = 16)  #facet label size
  ) +
  facet_wrap(~ time_period, drop = TRUE)

map_faceted

ggsave("plots/Map of ports by period labels no minor.jpg", dpi=300, width = 20, height = 20, units = "cm")
ggsave("plots/Map of ports by period labels med no minor.jpg", dpi=300, width = 30, height = 30, units = "cm")
ggsave("plots/Map of ports by period labels lrg no minor.jpg", dpi=300, width = 40, height = 40, units = "cm")


###  MULTIPLOT  ###


library(patchwork)

# arrange 3 rows x 3 columns
(a | b | c) /
  (d | e | f) +
  plot_annotation(tag_levels = 'A')  # adds A, B, C, D


ggsave("plots/Figure 3i.jpg", dpi=300, width = 25, height = 15, units = "cm")


# Create the animated plot
animated_map <- ggplot() +
  
  # Plot the coastline
  geom_sf(data = coastline, color = "black", fill = NA) +
  
  # Plot your points, sized by the number of vessels
  geom_point(
    data = hp_port_year,
    aes(x = Port.Longitude.dec, y = Port.Latitude.dec, size = total_vessels),
    color = "blue",      # adjust color to preference
    alpha = 0.7
  ) +
  
  # Adjust the point size range
  scale_size(range = c(1, 8), name = "Total Vessels") +
  
  # Use coord_sf with bounding box
  coord_sf(
    xlim = c(lon_min, lon_max),
    ylim = c(lat_min, lat_max),
    expand = FALSE
  ) +
  
  # Labels and minimal theme
  labs(
    title = "Number of Vessels by Port - Year: {frame_time}",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_bw()