The impact of hunting: Hunters’ inability to identify endangered and protected species

Authors: Dominik Marchowski, Robert Jurszo, Paweł Stańczak, Michał Jasiński, Sebastain Guentzel

Description

This document presents the code for the calculation part of the article devoted to waterfowl hunting on Lake Miedwie. The data comes from three hunting seasons: 2021/22, 2022/23 and 2023/24. To perform the analysis, raw data will be needed, which can be downloaded from:

https://marchowski.org/data/baza_V2.txt

hunt <- read.table(choose.files(), header = T) # Load the dataset - baza_V2.txt
head(hunt)
##   Pic   spec no HID seson prot
## 1  O1 ANACRE  1   1     O    0
## 2  O1 ANAPLA  3   1     O    0
## 3  O1 ANSALB 27   1     O    0
## 4  O1 ANSFAB  8   1     O    0
## 5  O1  ANSSP  8   1     O    0
## 6  O1 MARPEN 17   1     O    1

Packages

Load required packages

library(dplyr)
## 
## Dołączanie pakietu: 'dplyr'
## Następujące obiekty zostały zakryte z 'package:stats':
## 
##     filter, lag
## Następujące obiekty zostały zakryte z 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)

Group data by species and sum the ‘no’ column

summary <- hunt %>%
  group_by(spec) %>%
  summarise(total_count = sum(no))

# Display the results
print(summary)
## # A tibble: 16 × 2
##    spec    total_count
##    <chr>         <int>
##  1 ANAACU            9
##  2 ANACRE          102
##  3 ANAPLA          312
##  4 ANASPE            2
##  5 ANSALB          432
##  6 ANSANS          146
##  7 ANSFAB          153
##  8 ANSSP            24
##  9 AYTFER           59
## 10 AYTFUL            6
## 11 AYTMAR            1
## 12 FULATR            3
## 13 MARPEN           53
## 14 MARSTRE          26
## 15 MERALB            1
## 16 SPACLY            2

Add ‘code’ and ‘status’ columns to the data

summary <- summary %>%
  mutate(
    code = case_when(
      spec %in% c("MARPEN", "MARSTRE", "ANAACU", "SPACLY", "MERALB", "AYTMAR") ~ "P",
      spec == "AYTFER" ~ "VU",
      spec == "AYTFUL" ~ "NT",
      TRUE ~ ""
    ),
    status = case_when(
      code == "P" ~ "Protected",
      code == "VU" ~ "Vulnerable",
      code == "NT" ~ "Near Threatened",
      TRUE ~ "Other"
    )
  )

Create the plot

p <- ggplot(summary, aes(x = total_count, y = reorder(spec, total_count), fill = status)) +
  geom_bar(stat = "identity") +
  geom_text(
    data = subset(summary, code != ""), 
    aes(label = code),
    hjust = -0.5,
    size = 5,
    fontface = "bold"
  ) +
  labs(
    title = "",
    x = "Total number of birds hunted: seasons 2021/2022-2023/2024",
    y = "S   P   E   C   I   E   S"
  ) +
  theme_bw() +
  theme(
    legend.position = c(0.75, 0.5),     
    legend.background = element_rect(fill = "white", color = "black"),
    legend.text = element_text(size = 16, face = "bold"),
    legend.title = element_text(size = 18, face = "bold"),
    legend.key.size = unit(1.25, "lines"),
    legend.spacing = unit(1.25, "lines"),
    legend.margin = margin(15, 15, 15, 15),
    axis.title.x = element_text(size = 16),  
    axis.title.y = element_text(size = 16),  
    axis.text.x = element_text(size = 16),   
    axis.text.y = element_text(size = 14)
  ) +
  xlim(0, max(summary$total_count) * 1.2) +
  scale_fill_manual(
    values = c("Protected" = "red", "Vulnerable" = "orange", "Near Threatened" = "blue", "Other" = "gray44"),
    name = "Status"
  )
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(p)