Este documento se basa en el taller “UPR Ponce Data Workshop” ofrecido por Arlene Megill (10/16/2024) REFERENCIA
Algunas líneas de códigos fueron creadas utilizando GitHub Copilot.
REFERENCIAS de metodologías y metadata de NEON
Los siguientes paquetes son necesarios para el análisis de datos de aves en NEON. Debe instalarlos primero si no están en su biblioteca de paquetes.
library(neonUtilities)
library(neonOS)
library(ggplot2)
library(dplyr)
library(tidyr)
library(auk)
library(ggmap)
library(readr)Para obtener los datos de aves en NEON, se puede descargar
directamente desde el portal de NEON o cargar los datos directamente a R
usando la función loadByProduct en la biblioteca
neonUtilities.
El siguiente enlace describe los parámetros y funciones de
neonUtilities:
guan_bird <- loadByProduct(dpID="DP1.10003.001",
site=c("GUAN"), # acrónimo de Guánica
package="expanded",
include.provisional=T,
check.size=F) # evita interrupciónEn nuestro caso dpID es DP1.10003.001 que
contiene datos de aves en NEON, se obtuvo de la página de NEON
Explore Data
Products.
El archivo guan_bird contiene los datos de aves en
Guánica en una estructura de datos de lista. Para extraer las diferentes
‘data frames’ de metadata y datos, se puede utilizar la función
list2env para convertir la lista en ‘data frames’ separadas
en el ambiente de R.
list2env(guan_bird, envir = .GlobalEnv)La ‘data frame’ guan_bird_data contiene los datos de
aves en Guánica. Para examinar los datos, se puede utilizar la función
head para ver las primeras filas de la ‘data frame’.
head(brd_countdata)Los tipos de variables en la ‘data frame’ brd_countdata
se pueden obtener utilizando la función str.
str(brd_countdata)Para conocer las especies que contiene la ‘data frame’
brd_countdata, se puede utilizar la función
unique.
unique(brd_countdata$scientificName) # variable de especie## [1] NA "Junco hyemalis"
## [3] "Sitta canadensis" "Catharus guttatus"
## [5] "Cyanocitta stelleri" "Poecile gambeli"
## [7] "Setophaga coronata auduboni" "Regulus calendula"
## [9] "Spinus pinus" "Setophaga coronata"
## [11] "Certhia americana" "Picoides villosus"
## [13] "Turdus migratorius" "Corvus brachyrhynchos"
## [15] "Selasphorus platycercus" "Empidonax oberholseri"
## [17] "Oreothlypis celata" "Empidonax occidentalis"
## [19] "Melospiza lincolnii" "Loxia curvirostra"
## [21] "Sialia currucoides" "Dendragapus obscurus"
## [23] "Perisoreus canadensis" "Cardellina pusilla"
## [25] "Pinicola enucleator" "Zonotrichia leucophrys"
## [27] "Nucifraga columbiana" "Eremophila alpestris"
## [29] "Anthus rubescens"Para seleccionar una especie de interés, se puede utilizar la función
filter de la biblioteca dplyr.
# cargar especie de interés a un objeto
species_interest <- c("Antrostomus noctitherus")
# filtrar la 'data frame' para la especie de interés
df_guabairo <- brd_countdata %>%
filter(scientificName %in% species_interest)Vamos ahora a examinar algunas estadísticas descriptivas de la especie de interés.
# Núnmero de observaciones
n_observations <- nrow(df_guabairo)
n_observations
# Ámbito de fechas
date_range <- range(df_guabairo$startDate)
date_range
# Tamaño de los grupos
cluster_summary <- summary(df_guabairo$clusterSize)
cluster_summary
# Métodos de detección
method_summary <- table(df_guabairo$detectionMethod)
method_summaryVamos crear una gráfica del número de observaciones por fecha.
# convertir datos de fecha a formato Date
df_guabairo$startDate <- as.Date(df_guabairo$startDate)
# gráfica de número de observaciones por fecha
ggplot(df_guabairo, aes(x = startDate)) +
geom_histogram(binwidth = 30, fill = "skyblue", color = "black") +
labs(
x = "Fecha",
y = "Número de observaciones")Figura 1. Número de observaciones de Guabairo en Guánica por fecha de muestreo.
Ahora un gráfico de puntos que muestre el promedio de observaciones por plotID con su línea de desviación estándar y por año.
# Load necessary libraries
library(dplyr)
library(ggplot2)
library(lubridate)
# Load the dataset
df_guabairo <- read.csv("df_guabairo.csv")
# Extract year from the startDate column
df_guabairo <- df_guabairo %>%
mutate(year = year(as.Date(startDate)))
# Step 1: Count the number of observations for each unique pointID within each plotID and year
point_counts <- df_guabairo %>%
group_by(plotID, year, pointID) %>%
summarise(
count_per_pointID = n(), # Count rows (observations) for each unique pointID
.groups = "drop"
)
# Step 2: Calculate mean and standard deviation of counts for each plotID and year
summary_stats <- point_counts %>%
group_by(plotID, year) %>%
summarise(
mean_obs = mean(count_per_pointID, na.rm = TRUE), # Mean of counts of unique pointIDs
sd_obs = sd(count_per_pointID, na.rm = TRUE), # Standard deviation of counts
.groups = "drop"
)
# Create the point graph with jittered points and fixed error bars
ggplot(summary_stats, aes(x = factor(year), y = mean_obs, color = plotID)) +
geom_point(size = 3, position = position_jitter(width = 0.2)) + # jittered points for mean
# geom_errorbar(
# aes(ymin = mean_obs - sd_obs, ymax = mean_obs + sd_obs),
# width = 0.2,
# position = position_jitter(width = 0.1)
# ) + # jittered error bars
labs(
title = "Mean Count of Unique PointIDs per PlotID per Year",
x = "Year",
y = "Mean Count of Unique PointIDs",
color = "PlotID"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1)
)Figura 2. Promedio de observaciones de Guabairo por PlotID y Año en Guánica.
Para crear un mapa debemos primero tener el archivo de localización
de los plots y puntos de muestreo. En primer lugar vamos a crear
archivos csv de brd_perpoint y
df_guabairo.
# crear archivo csv de brd_perpoint
write.csv(brd_perpoint, "brd_perpoint.csv")
# crear archivo csv de df_guabairo
write.csv(df_guabairo, "df_guabairo.csv")Creación del mapa de distribución de Guabairo en
# Load necessary libraries
library(ggplot2)
library(dplyr)
library(ggspatial)
library(rnaturalearth)
library(rnaturalearthdata)
# Load the datasets
brd_perpoint <- read.csv("brd_perpoint.csv")
df_guabairo <- read.csv("df_guabairo.csv")
# Summarize total counts of guabairo per plotID and pointID
guabairo_counts <- df_guabairo %>%
group_by(plotID, pointID) %>%
summarise(total_count = n(), .groups = "drop")
# Merge summarized counts with the coordinates from the brd_perpoint dataset
guabairo_map_data <- brd_perpoint %>%
select(plotID, pointID, decimalLatitude, decimalLongitude) %>%
distinct() %>%
left_join(guabairo_counts, by = c("plotID", "pointID"))
# Remove rows with NA counts
guabairo_map_data <- guabairo_map_data %>%
filter(!is.na(total_count))
# Load physical map data for the region (Puerto Rico)
world <- ne_countries(scale = "medium", returnclass = "sf")
# Plot the map with a physical background
ggplot() +
geom_sf(data = world, fill = "lightgray", color = "white") + # Physical map background
geom_point(data = guabairo_map_data, aes(x = decimalLongitude, y = decimalLatitude, size = total_count),
color = "blue", alpha = 0.7) + # Points sized by total count
labs(
x = "Longitude",
y = "Latitude",
size = "Total Count"
) +
theme_minimal() +
coord_sf(xlim = c(-66.88, -66.82), ylim = c(17.96, 17.99), expand = FALSE) + # Restrict to Puerto Rico
annotation_scale(location = "tr", width_hint = 0.5) + # Add a scale bar
annotation_north_arrow(location = "tl", which_north = "true", style = north_arrow_fancy_orienteering())Para analizar la biodiversidad de las especies de aves en los datos, podemos calcular el índice de diversidad de Simpson y el número de especies únicas por plotID y año.
# Load necessary libraries
library(dplyr)
library(ggplot2)
# Load the data
data <- read.csv("brd_coundata.csv")
# Extract year from startDate
data <- data %>%
mutate(year = as.numeric(format(as.Date(startDate, format = "%Y-%m-%d"), "%Y")))
# Count unique species (taxonID) per plotID per year
species_count <- data %>%
group_by(plotID, year) %>%
summarise(unique_species = n_distinct(taxonID), .groups = "drop")
# Create the ggplot2 graph
ggplot(species_count, aes(x = factor(year), y = unique_species, fill = plotID)) +
geom_bar(stat = "identity", position = "dodge") +
labs(
title = "Number of Unique Species per Plot per Year",
x = "Year",
y = "Number of Unique Species",
fill = "Plot ID"
) +
theme_minimal()# Load necessary libraries
library(dplyr)
library(gt)
# Load the data
data <- read.csv("brd_coundata.csv")
# Extract year from startDate
data <- data %>%
mutate(year = as.numeric(format(as.Date(startDate, format = "%Y-%m-%d"), "%Y")))
# Calculate Simpson's diversity index and unique species count (S) for each plotID and year
simpson_index <- data %>%
group_by(plotID, year, taxonID) %>%
summarise(count = n(), .groups = "drop") %>% # Count occurrences of each species
group_by(plotID, year) %>%
summarise(
total_count = sum(count), # Total number of individuals
unique_species = n_distinct(taxonID), # Number of unique species (S)
simpson_index = 1 - sum((count / sum(count))^2), # Simpson's index formula
.groups = "drop"
)
# Create a gt table to display the results
simpson_table <- simpson_index %>%
gt() %>%
tab_header(
title = "Simpson's Biodiversity Index with Unique Species Count",
subtitle = "Calculated for Each Plot and Year"
) %>%
cols_label(
plotID = "Plot ID",
year = "Year",
total_count = "Total Count",
unique_species = "Unique Species (S)",
simpson_index = "Simpson Index"
) %>%
fmt_number(
columns = vars(simpson_index),
decimals = 3
)
# Print the gt table
simpson_table