proy

— Libraries —

Lynx <- read.csv("Lynx.csv", sep="\t", quote = "")

Quebr <- read.csv("Quebrantahuesos.csv", sep="\t", quote = "")

4. Limpieza de datos

En esta sección se realiza la limpieza de los datos obtenidos de GBIF para ambas especies. El objetivo es construir una base homogénea y fiable para los análisis posteriores. Para ello, se aplican los mismos criterios de selección de variables, filtrado geográfico, revisión temporal y eliminación de duplicados en ambos conjuntos de datos.

4.1 Limpieza del lince ibérico

En primer lugar, se seleccionan únicamente las variables necesarias para el análisis: coordenadas geográficas, fecha del avistamiento, variables temporales y algunos campos descriptivos útiles para interpretar los registros.

lynx_clean <- Lynx[, c(
  "decimalLatitude",
  "decimalLongitude",
  "eventDate",
  "year",
  "month",
  "countryCode",
  "individualCount",
  "stateProvince",
  "locality"
)]

Tras esta primera reducción, se eliminan muchas columnas taxonómicas y administrativas que no aportan directamente al análisis principal, dejando una tabla más manejable y centrada en la información espacial y temporal.

A continuación, se filtran los registros localizados en España y se eliminan las observaciones sin latitud o longitud, ya que no pueden utilizarse ni en mapas ni en análisis espaciales. Además, se aplica un filtro espacial aproximado mediante un bounding box para conservar únicamente las observaciones situadas dentro del territorio español.

lynx_clean <- lynx_clean[!is.na(lynx_clean$countryCode) & lynx_clean$countryCode == "ES", ]
lynx_clean <- lynx_clean[
  !is.na(lynx_clean$decimalLatitude) &
  !is.na(lynx_clean$decimalLongitude),
]

lynx_clean <- lynx_clean[
  lynx_clean$decimalLatitude >= 35 &
  lynx_clean$decimalLatitude <= 45 &
  lynx_clean$decimalLongitude >= -10 &
  lynx_clean$decimalLongitude <= 5,
]

Después, se revisa la información temporal. Se eliminan los registros sin año y sin fecha válida, y también aquellos con años excesivamente antiguos, ya que podrían introducir ruido en el análisis comparativo reciente.

lynx_clean <- lynx_clean[
  !is.na(lynx_clean$year) & lynx_clean$year >= 1900,
]

Con ello se obtiene una base temporalmente más coherente, adecuada para representar la evolución de los avistamientos en las últimas décadas.

Finalmente, se revisan posibles duplicados, considerando como repetidos los registros con la misma fecha y las mismas coordenadas. En caso de coincidencia, se conserva una sola observación.

lynx_clean <- lynx_clean |> distinct(decimalLatitude, decimalLongitude, eventDate, .keep_all = TRUE)

4.2 Limpieza del quebrantahuesos

A continuación, se aplica el mismo procedimiento de limpieza al conjunto de datos del quebrantahuesos, utilizando los mismos criterios de selección de variables, filtrado espacial y depuración temporal.

quebr_clean <- Quebr[, c(
  "decimalLatitude",
  "decimalLongitude",
  "eventDate",
  "year",
  "month",
  "countryCode",
  "individualCount",
  "stateProvince",
  "locality"
)]

quebr_clean <- quebr_clean[!is.na(quebr_clean$countryCode) & quebr_clean$countryCode == "ES", ]
quebr_clean <- quebr_clean[
  !is.na(quebr_clean$decimalLatitude) &
  !is.na(quebr_clean$decimalLongitude),
]

quebr_clean <- quebr_clean[
  quebr_clean$decimalLatitude >= 35 &
  quebr_clean$decimalLatitude <= 45 &
  quebr_clean$decimalLongitude >= -10 &
  quebr_clean$decimalLongitude <= 5,
]

quebr_clean <- quebr_clean[
  !is.na(quebr_clean$year) & quebr_clean$year >= 1900,
]

quebr_clean <- quebr_clean |> distinct(decimalLatitude, decimalLongitude, eventDate, .keep_all = TRUE)

4.3 Resumen de la limpieza

En resumen, tras aplicar los filtros de selección, localización, calidad espacial, consistencia temporal y eliminación de duplicados, se obtienen dos conjuntos de datos más homogéneos y adecuados para el análisis. Esta depuración permite trabajar con observaciones comparables y reduce el efecto de registros incompletos o problemáticos.

5. Objetivo 1. Visualización de avistamientos

En este objetivo se analiza la evolución temporal y la distribución espacial de los avistamientos del lince ibérico y del quebrantahuesos en España durante los últimos 35 años.

5.1 Evolución temporal

En primer lugar, se analiza la evolución temporal del número de avistamientos en España durante el periodo 1990–2025. Se selecciona este intervalo porque se ajusta al horizonte temporal de 30–40 años planteado en el proyecto y permite evitar la influencia de registros demasiado antiguos y escasos.

lynx_obj1  <- lynx_clean[lynx_clean$year >= 1990, ]
quebr_obj1 <- quebr_clean[quebr_clean$year >= 1990, ]

lynx_years <- as.data.frame(table(lynx_obj1$year))
colnames(lynx_years) <- c("year", "n")
lynx_years$year <- as.numeric(as.character(lynx_years$year))

quebr_years <- as.data.frame(table(quebr_obj1$year))
colnames(quebr_years) <- c("year", "n")
quebr_years$year <- as.numeric(as.character(quebr_years$year))

ggplot(lynx_years, aes(x = year, y = n)) +
  geom_line() +
  geom_point() +
  labs(
    title = "Number of Iberian lynx sightings per year in Spain (1990–2025)",
    x = "Year",
    y = "Number of sightings"
  ) +
  theme_minimal()

En la serie temporal del lince ibérico se observa una tendencia general creciente en el número de avistamientos a lo largo del periodo analizado. Este aumento es especialmente visible a partir de la década de 2010.

Sin embargo, destaca un pico muy pronunciado en un año concreto, con un número de observaciones muy superior al resto. Este comportamiento no parece corresponder a una evolución natural de la población, sino más bien a un efecto relacionado con la recopilación de datos, como la incorporación masiva de registros o campañas de observación intensiva.

Por tanto, aunque existe una tendencia al alza, los resultados deben interpretarse con cautela debido a la posible influencia de este valor atípico.

ggplot(quebr_years, aes(x = year, y = n)) +
  geom_line() +
  geom_point() +
  labs(
    title = "Number of Bearded Vulture sightings per year in Spain (1990–2025)",
    x = "Year",
    y = "Number of sightings"
  ) +
  theme_minimal()

En el caso del quebrantahuesos, se observa una tendencia claramente creciente en el número de avistamientos desde principios de los años 2000. Este incremento es progresivo y se intensifica especialmente en la última década.

Este comportamiento puede reflejar tanto una mejora en el estado de conservación de la especie como un aumento en el esfuerzo de observación y en la disponibilidad de datos.

Sin embargo, en los últimos años se aprecia una caída brusca en el número de registros. Este descenso no debe interpretarse como una disminución real de la población, sino probablemente como consecuencia de que los datos más recientes aún no están completamente incorporados en la base de datos.

5.2 Distribución espacial

lynx_obj1$period <- cut(
  lynx_obj1$year,
  breaks = c(1990, 2000, 2010, 2020, 2026),
  labels = c("1990-1999", "2000-2009", "2010-2019", "2020-2025"),
  right = FALSE
)

lynx_obj1 <- lynx_obj1[!is.na(lynx_obj1$period), ]

quebr_obj1$period <- cut(
  quebr_obj1$year,
  breaks = c(1990, 2000, 2010, 2020, 2026),
  labels = c("1990-1999", "2000-2009", "2010-2019", "2020-2025"),
  right = FALSE
)

quebr_obj1 <- quebr_obj1[!is.na(quebr_obj1$period), ]

Tras analizar la evolución temporal, se estudia a continuación la distribución espacial de los avistamientos. Para facilitar la interpretación, los registros se agrupan en cuatro periodos temporales en lugar de utilizar el año exacto como escala continua.

La siguiente figura muestra la evolución temporal del número de avistamientos del lince ibérico en España.

spain_map <- map_data("world", region = "Spain")

ggplot() +
  geom_polygon(
    data = spain_map,
    aes(x = long, y = lat, group = group),
    fill = "gray95", color = "black"
  ) +
  geom_point(
    data = lynx_obj1,
    aes(x = decimalLongitude, y = decimalLatitude, color = period),
    alpha = 0.5, size = 1
  ) +
  labs(
    title = "Spatial distribution of Iberian lynx sightings in Spain by period",
    x = "Longitude",
    y = "Latitude",
    color = "Period"
  ) +
  coord_fixed(1.3) +
  theme_minimal()

A continuación se presenta la evolución temporal de los avistamientos del quebrantahuesos.

spain_map <- map_data("world", region = "Spain")

ggplot() +
  geom_polygon(
    data = spain_map,
    aes(x = long, y = lat, group = group),
    fill = "gray95", color = "black"
  ) +
  geom_point(
    data = quebr_obj1,
    aes(x = decimalLongitude, y = decimalLatitude, color = period),
    alpha = 0.5, size = 1
  ) +
  labs(
    title = "Spatial distribution of Bearded Vulture sightings in Spain by period",
    x = "Longitude",
    y = "Latitude",
    color = "Period"
  ) +
  coord_fixed(1.3) +
  theme_minimal()

En general, se observa que el lince ibérico presenta una distribución más concentrada en el sur de España, mientras que el quebrantahuesos se localiza principalmente en zonas montañosas del norte. Además, el número de observaciones aumenta en las últimas décadas, lo que puede reflejar tanto una mejora en la conservación como un aumento en la recopilación de datos.

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

OBJECTIVE 2: Analysis of the Impact of Protected Areas

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

# Disable S2 spherical geometry to avoid precision errors with OSM polygons
sf_use_s2(FALSE)

# ============================================================
# 6.1 Carga y Procesado (Optimizado con Caching .rds)
# ============================================================
processed_data_file <- "protected_areas_spain.rds"

if (file.exists(processed_data_file)) {
  # Si ya existe, lo cargamos
  protected_sf_spain <- readRDS(processed_data_file)
} else {
  # Si no existe, hacemos el trabajo pesado
  poly0 <- st_read("poly0/WDPA_Apr2026_Public_shp-polygons.shp", quiet = TRUE)
  poly1 <- st_read("poly1/WDPA_Apr2026_Public_shp-polygons.shp", quiet = TRUE)
  poly2 <- st_read("poly2/WDPA_Apr2026_Public_shp-polygons.shp", quiet = TRUE)
  
  protected_sf_spain <- rbind(poly0, poly1, poly2) %>%
    filter(ISO3 == "ESP")
  
  saveRDS(protected_sf_spain, processed_data_file)
  rm(poly0, poly1, poly2)
}

# ============================================================
# 6.2 Spatial Join & Map Preparation
# ============================================================
lynx_sf  <- st_as_sf(lynx_obj1, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326, remove = FALSE)
quebr_sf <- st_as_sf(quebr_obj1, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326, remove = FALSE)

lynx_sf$in_protected  <- lengths(st_intersects(lynx_sf,  protected_sf_spain)) > 0
quebr_sf$in_protected <- lengths(st_intersects(quebr_sf, protected_sf_spain)) > 0

lynx_obj2  <- st_drop_geometry(lynx_sf)
quebr_obj2 <- st_drop_geometry(quebr_sf)

# Simplificamos TODAS las zonas para el mapa (Bug visual arreglado)
protected_sf_map <- protected_sf_spain %>%
  st_simplify(dTolerance = 0.005)

# ============================================================
# 6.3 Spatial join: inside / outside protected area
# ============================================================

# Flag each sighting based on whether it falls within a protected area
lynx_sf$in_protected  <- lengths(st_intersects(lynx_sf,  protected_sf_spain)) > 0
quebr_sf$in_protected <- lengths(st_intersects(quebr_sf, protected_sf_spain)) > 0

# Drop geometry to work with standard data frames
lynx_obj2  <- st_drop_geometry(lynx_sf)
quebr_obj2 <- st_drop_geometry(quebr_sf)

# Quick summary check
cat("=== Iberian Lynx ===\n");   print(table(lynx_obj2$in_protected))

## === Iberian Lynx ===

## 
## FALSE  TRUE 
##   227   545

cat("=== Bearded Vulture ===\n"); print(table(quebr_obj2$in_protected))

## === Bearded Vulture ===

## 
## FALSE  TRUE 
##  4692  8900

# ============================================================
# 6.4 Plot: Annual evolution of sightings by zone
# ============================================================
plot_evolution <- function(df, species_name) {
  df %>%
    filter(year >= 1990) %>%
    mutate(Zone = ifelse(in_protected, "Protected", "Non-protected")) %>%
    group_by(year, Zone) %>%
    summarise(count = n(), .groups = "drop") %>%
    ggplot(aes(x = year, y = count, color = Zone)) +
    geom_line(linewidth = 1) +
    geom_point(size = 2) +
    scale_color_manual(values = c("Protected" = "#2E7D32", "Non-protected" = "#C62828")) +
    labs(
      title = paste("Evolution of", species_name, "sightings (1990–2025)"),
      x = "Year", y = "Number of sightings"
    ) +
    theme_minimal() +
    theme(legend.position = "bottom")
}

# Affichage côte à côte ou l'un sous l'autre
plot_evolution(lynx_obj2,  "Iberian Lynx")

plot_evolution(quebr_obj2, "Bearded Vulture")

# ============================================================
# 6.5 Heatmap: Period-to-period growth rate by zone
# ============================================================
target_years <- c(2007, 2010, 2013, 2018, 2020, 2021, 2022, 2023, 2024, 2025)

prepare_period_data <- function(df, species_name) {
  df %>%
    filter(year %in% target_years) %>%
    mutate(zone = ifelse(in_protected, "Protected", "Non-protected")) %>%
    group_by(year, zone) %>%
    summarise(total = n(), .groups = "drop") %>%
    complete(year = target_years, zone, fill = list(total = 0)) %>%
    mutate(species = species_name)
}

combined_periods <- bind_rows(
  prepare_period_data(lynx_obj2,  "Iberian Lynx"),
  prepare_period_data(quebr_obj2, "Bearded Vulture")
)

heatmap_data <- combined_periods %>%
  group_by(species, zone) %>%
  arrange(year) %>%
  mutate(
    prev_year  = lag(year),
    prev_total = lag(total),
    period     = ifelse(!is.na(prev_year), paste0(prev_year, "-", year), NA),
    delta_pct  = ifelse(!is.na(prev_total) & prev_total >= 5,
                        round((total - prev_total) / prev_total * 100, 1),
                        NA),
    label_text = case_when(
      is.na(prev_total)  ~ "N/A",
      prev_total < 5     ~ paste0(prev_total, " \u2192 ", total),
      delta_pct > 0      ~ paste0("+", delta_pct, "%"),
      TRUE               ~ paste0(delta_pct, "%")
    )
  ) %>%
  filter(!is.na(period)) %>%
  mutate(period = factor(period, levels = unique(period)))

ggplot(heatmap_data, aes(x = period, y = zone, fill = delta_pct)) +
  geom_tile(color = "white", linewidth = 0.8) +
  geom_text(
    aes(label = label_text, color = is.na(delta_pct) | abs(delta_pct) < 50),
    size = 3.5, fontface = "bold", show.legend = FALSE
  ) +
  scale_color_manual(values = c("TRUE" = "black", "FALSE" = "white")) +
  facet_wrap(~ species, ncol = 1) +
  scale_fill_gradient2(
    low = "#C62828", mid = "#FFF9C4", high = "#2E7D32",
    midpoint = 0, limits = c(-100, 150),
    oob = scales::squish, na.value = "grey90", name = "Growth (%)"
  ) +
  labs(
    title = "Period-to-period change in sightings by zone",
    subtitle = "Grey blocks: absolute counts shown when baseline < 5 sightings.",
    x = "Comparison period", y = NULL
  ) +
  theme_minimal(base_size = 12) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "right",
    strip.text = element_text(face = "bold", size = 13),
    panel.grid = element_blank()
  )

# ============================================================
# 6.6 % share: Protected vs Non-protected
# ============================================================
pct_data <- combined_periods %>%
  group_by(species, year) %>%
  mutate(pct = round(total / sum(total) * 100, 1)) %>%
  ungroup() %>%
  group_by(species, year) %>%
  filter(sum(total) > 0) %>%
  ungroup()

ggplot(pct_data, aes(x = factor(year), y = pct, fill = zone)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(
    aes(label = ifelse(pct >= 8, paste0(pct, "%"), "")),
    position = position_fill(vjust = 0.5),
    size = 3.2, fontface = "bold", color = "white"
  ) +
  scale_y_continuous(labels = percent_format()) +
  scale_fill_manual(
    values = c("Protected" = "#2E7D32", "Non-protected" = "#C62828"),
    name   = NULL
  ) +
  facet_wrap(~ species, ncol = 1) +
  labs(
    title = "Share of sightings inside vs outside protected areas",
    x = "Year", y = "Proportion of sightings"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "bottom",
    strip.text = element_text(face = "bold", size = 13),
    panel.grid.major.x = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

# ============================================================
# 6.7 Interactive map
# ============================================================
# Transformation pour Crosstalk : on enlève la géométrie sf complexe
map_data_coords <- bind_rows(
  lynx_sf  %>% filter(year %in% target_years) %>% mutate(species = "Lynx"),
  quebr_sf %>% filter(year %in% target_years) %>% mutate(species = "Vulture")
) %>%
  mutate(
    lng = st_coordinates(.)[,1],
    lat = st_coordinates(.)[,2],
    # Création d'une catégorie unique pour la couleur
    color_group = paste0(species, "_", in_protected) 
  ) %>%
  st_drop_geometry() %>%
  select(species, year, in_protected, lng, lat, color_group)

shared_map_data <- SharedData$new(map_data_coords)

# Palette Leaflet standard
pal <- colorFactor(
  palette = c("#006d2c", "#e6550d", "#2b8cbe", "#c51b8a"),
  domain = c("Lynx_TRUE", "Lynx_FALSE", "Vulture_TRUE", "Vulture_FALSE")
)

tagList(
  tags$div(
    style = "padding: 10px; background: #f8f9fa; border-bottom: 1px solid #ddd;",
    filter_slider(
      id         = "year_slider",
      label      = "Select observation year:",
      sharedData = shared_map_data,
      column     = ~year,
      step       = 1, sep = "", width = "100%"
    )
  ),
  leaflet(shared_map_data, width = "100%", height = 600) %>%
    addProviderTiles(providers$CartoDB.Positron) %>%
    addPolygons(
      data  = protected_sf_map,
      color = "#2ca25f", weight = 1, fillOpacity = 0.2,
      group = "Protected Areas",
      label = ~DESIG_ENG
    ) %>%
    addCircleMarkers(
      lng = ~lng, lat = ~lat,
      color = ~pal(color_group),
      radius = 4, stroke = TRUE, weight = 1, fillOpacity = 0.8,
      popup = ~paste0(
        "<b>", species, "</b> (", year, ")<br>",
        "Status: ", ifelse(in_protected, "Inside protected area", "Outside protected area")
      )
    ) %>%
    addLayersControl(
      overlayGroups = "Protected Areas",
      options = layersControlOptions(collapsed = FALSE)
    )
)

Select observation year: