Datos_eBird
Victoria Escobar
2024-03-01
Instalación de paquetes requeridos
Se emplearán diversos paquetes de R para acceder a datos de eBird, trabajar con datos espaciales, realizar procesamiento y manipulación de datos, así como para llevar a cabo el entrenamiento de modelos.
if (!requireNamespace("remotes", quietly = TRUE)) {
install.packages("remotes")
}
remotes::install_github("ebird/ebird-best-practices")## Skipping install of 'ebirdbestpractices' from a github remote, the SHA1 (b562653d) has not changed since last install.
## Use `force = TRUE` to force installationDatos SIG
A lo largo de la guía, se generarán mapas de distribuciones de especies, utilizando datos SIG para las fronteras políticas. Natural Earth se destaca como una fuente integral de datos SIG rasterizados y vectoriales, facilitando la creación de mapas cartográficos profesionales.
##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE# file to save spatial data
gpkg_file <- "data/gis-data.gpkg"
dir.create(dirname(gpkg_file), showWarnings = FALSE, recursive = TRUE)
# political boundaries
# land border with lakes removed
ne_land <- ne_download(scale = 50, category = "cultural",
type = "admin_0_countries_lakes",
returnclass = "sf") |>
filter(CONTINENT %in% c("North America", "South America")) |>
st_set_precision(1e6) |>
st_union()## Reading layer `ne_50m_admin_0_countries_lakes' from data source
## `C:\Users\victo\AppData\Local\Temp\RtmpEhcxNr\ne_50m_admin_0_countries_lakes.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 242 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
## Geodetic CRS: WGS 84# country boundaries
ne_countries <- ne_download(scale = 50, category = "cultural",
type = "admin_0_countries_lakes",
returnclass = "sf") |>
select(country = ADMIN, country_code = ISO_A2)## Reading layer `ne_50m_admin_0_countries_lakes' from data source
## `C:\Users\victo\AppData\Local\Temp\RtmpEhcxNr\ne_50m_admin_0_countries_lakes.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 242 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
## Geodetic CRS: WGS 84# state boundaries for united states
ne_states <- ne_download(scale = 50, category = "cultural",
type = "admin_1_states_provinces",
returnclass = "sf") |>
filter(iso_a2 == "US") |>
select(state = name, state_code = iso_3166_2)## Reading layer `ne_50m_admin_1_states_provinces' from data source
## `C:\Users\victo\AppData\Local\Temp\RtmpEhcxNr\ne_50m_admin_1_states_provinces.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 294 features and 121 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -46.96289 xmax: 180 ymax: 83.11611
## Geodetic CRS: WGS 84# country lines
# downloaded globally then filtered to north america with st_intersect()
ne_country_lines <- ne_download(scale = 50, category = "cultural",
type = "admin_0_boundary_lines_land",
returnclass = "sf") |>
st_geometry()## Reading layer `ne_50m_admin_0_boundary_lines_land' from data source
## `C:\Users\victo\AppData\Local\Temp\RtmpEhcxNr\ne_50m_admin_0_boundary_lines_land.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 390 features and 39 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -141.0021 ymin: -55.114 xmax: 145.9402 ymax: 70.06482
## Geodetic CRS: WGS 84lines_on_land <- st_intersects(ne_country_lines, ne_land, sparse = FALSE) |>
as.logical()
ne_country_lines <- ne_country_lines[lines_on_land]
# states, north america
ne_state_lines <- ne_download(scale = 50, category = "cultural",
type = "admin_1_states_provinces_lines",
returnclass = "sf") |>
filter(ADM0_A3 %in% c("USA", "CAN")) |>
mutate(iso_a2 = recode(ADM0_A3, USA = "US", CAN = "CAN")) |>
select(country = ADM0_NAME, country_code = iso_a2)## Reading layer `ne_50m_admin_1_states_provinces_lines' from data source
## `C:\Users\victo\AppData\Local\Temp\RtmpEhcxNr\ne_50m_admin_1_states_provinces_lines.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 581 features and 43 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -139.0565 ymin: -38.0716 xmax: 174.4685 ymax: 78.68672
## Geodetic CRS: WGS 84# save all layers to a geopackage
unlink(gpkg_file)
write_sf(ne_land, gpkg_file, "ne_land")
write_sf(ne_countries, gpkg_file, "ne_countries")
write_sf(ne_states, gpkg_file, "ne_states")
write_sf(ne_country_lines, gpkg_file, "ne_country_lines")
write_sf(ne_state_lines, gpkg_file, "ne_state_lines")Se revisa la versión de R.
## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.3.3 (2024-02-29 ucrt)
## os Windows 11 x64 (build 22631)
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate Spanish_Colombia.utf8
## ctype Spanish_Colombia.utf8
## tz America/Bogota
## date 2024-03-02
## pandoc 3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## bslib 0.5.1 2023-08-11 [1] CRAN (R 4.3.1)
## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.1)
## callr 3.7.3 2022-11-02 [1] CRAN (R 4.3.1)
## class 7.3-22 2023-05-03 [2] CRAN (R 4.3.3)
## classInt 0.4-10 2023-09-05 [1] CRAN (R 4.3.2)
## cli 3.6.1 2023-03-23 [1] CRAN (R 4.3.1)
## codetools 0.2-19 2023-02-01 [2] CRAN (R 4.3.3)
## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.3.1)
## curl 5.0.1 2023-06-07 [1] CRAN (R 4.3.1)
## DBI 1.1.3 2022-06-18 [1] CRAN (R 4.3.1)
## devtools 2.4.5 2022-10-11 [1] CRAN (R 4.3.1)
## digest 0.6.31 2022-12-11 [1] CRAN (R 4.3.1)
## dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.3.2)
## e1071 1.7-13 2023-02-01 [1] CRAN (R 4.3.1)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.1)
## evaluate 0.23 2023-11-01 [1] CRAN (R 4.3.2)
## fansi 1.0.4 2023-01-22 [1] CRAN (R 4.3.1)
## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.1)
## fs 1.6.2 2023-04-25 [1] CRAN (R 4.3.1)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.1)
## glue 1.6.2 2022-02-24 [1] CRAN (R 4.3.1)
## htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.3.1)
## htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.3.1)
## httpuv 1.6.11 2023-05-11 [1] CRAN (R 4.3.1)
## httr 1.4.7 2023-08-15 [1] CRAN (R 4.3.1)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.1)
## jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.3.1)
## KernSmooth 2.23-22 2023-07-10 [2] CRAN (R 4.3.3)
## knitr 1.45 2023-10-30 [1] CRAN (R 4.3.2)
## later 1.3.1 2023-05-02 [1] CRAN (R 4.3.1)
## lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.2)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.1)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.1)
## mime 0.12 2021-09-28 [1] CRAN (R 4.3.0)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.3.1)
## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.1)
## pkgbuild 1.4.2 2023-06-26 [1] CRAN (R 4.3.1)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.1)
## pkgload 1.3.2 2022-11-16 [1] CRAN (R 4.3.1)
## prettyunits 1.2.0 2023-09-24 [1] CRAN (R 4.3.1)
## processx 3.8.1 2023-04-18 [1] CRAN (R 4.3.1)
## profvis 0.3.8 2023-05-02 [1] CRAN (R 4.3.1)
## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.3.1)
## proxy 0.4-27 2022-06-09 [1] CRAN (R 4.3.1)
## ps 1.7.5 2023-04-18 [1] CRAN (R 4.3.1)
## purrr 1.0.1 2023-01-10 [1] CRAN (R 4.3.1)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.1)
## Rcpp 1.0.10 2023-01-22 [1] CRAN (R 4.3.1)
## remotes 2.4.2 2021-11-30 [1] CRAN (R 4.3.1)
## rlang 1.1.1 2023-04-28 [1] CRAN (R 4.3.1)
## rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.3.1)
## rnaturalearth * 1.0.1 2023-12-15 [1] CRAN (R 4.3.2)
## rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.1)
## s2 1.1.6 2023-12-19 [1] CRAN (R 4.3.2)
## sass 0.4.7 2023-07-15 [1] CRAN (R 4.3.1)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.1)
## sf * 1.0-15 2023-12-18 [1] CRAN (R 4.3.2)
## shiny 1.7.4 2022-12-15 [1] CRAN (R 4.3.1)
## stringi 1.7.12 2023-01-11 [1] CRAN (R 4.3.0)
## stringr 1.5.1 2023-11-14 [1] CRAN (R 4.3.2)
## terra 1.7-71 2024-01-31 [1] CRAN (R 4.3.2)
## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.1)
## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.3.1)
## units 0.8-5 2023-11-28 [1] CRAN (R 4.3.2)
## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.3.1)
## usethis 2.2.1 2023-06-23 [1] CRAN (R 4.3.1)
## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.3.1)
## vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.2)
## withr 3.0.0 2024-01-16 [1] CRAN (R 4.3.2)
## wk 0.9.1 2023-11-29 [1] CRAN (R 4.3.2)
## xfun 0.39 2023-04-20 [1] CRAN (R 4.3.1)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.1)
## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.3.0)
##
## [1] C:/Users/victo/AppData/Local/R/win-library/4.3
## [2] C:/Program Files/R/R-4.3.3/library
##
## ──────────────────────────────────────────────────────────────────────────────Importar datos de eBird a R
Lista de verificación de Georgia (donde anotan todos los pájaros que ven).
## auk 0.7.0 is designed for EBD files downloaded after 2023-10-25.
## No EBD data directory set, see ?auk_set_ebd_path to set EBD_PATH
## eBird taxonomy version: 2023##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(readr)
library(sf)
f_sed <- "C:/Users/victo/OneDrive - PUJ Cali/Documentos/Maestría/2. Proyecto aplicado/Replica datos eBird/ebd_US-GA_woothr_smp_relOct-2023/ebd_US-GA_woothr_smp_relOct-2023_sampling.txt"
checklists <- read_sampling(f_sed)
glimpse(checklists)## Rows: 1,120,132
## Columns: 31
## $ checklist_id <chr> "S78945522", "S100039674", "S7892023", "S789…
## $ last_edited_date <chr> "2021-04-13 15:28:08.442812", "2022-01-03 14…
## $ country <chr> "United States", "United States", "United St…
## $ country_code <chr> "US", "US", "US", "US", "US", "US", "US", "U…
## $ state <chr> "Georgia", "Georgia", "Georgia", "Georgia", …
## $ state_code <chr> "US-GA", "US-GA", "US-GA", "US-GA", "US-GA",…
## $ county <chr> "Monroe", "Monroe", "Chatham", "Chatham", "C…
## $ county_code <chr> "US-GA-207", "US-GA-207", "US-GA-051", "US-G…
## $ iba_code <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ bcr_code <int> 29, 29, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ usfws_code <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ atlas_block <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ locality <chr> "Rum Creek M.A.R.S.H. Project", "Rum Creek M…
## $ locality_id <chr> "L876993", "L876993", "L877616", "L877614", …
## $ locality_type <chr> "H", "H", "H", "H", "H", "H", "H", "H", "H",…
## $ latitude <dbl> 33.06119, 33.06119, 31.38881, 31.71668, 31.7…
## $ longitude <dbl> -83.87000, -83.87000, -80.21118, -80.43091, …
## $ observation_date <date> 1997-12-31, 1997-02-23, 1983-09-13, 1984-02…
## $ time_observations_started <chr> NA, NA, NA, NA, NA, NA, NA, "07:20:00", "15:…
## $ observer_id <chr> "obs1675368", "obs198993", "obs243631", "obs…
## $ sampling_event_identifier <chr> "S78945522", "S100039674", "S7892023", "S789…
## $ protocol_type <chr> "Historical", "Historical", "Incidental", "I…
## $ protocol_code <chr> "P62", "P62", "P20", "P20", "P20", "P20", "P…
## $ project_code <chr> "EBIRD", "EBIRD", "EBIRD", "EBIRD", "EBIRD",…
## $ duration_minutes <int> NA, NA, NA, NA, NA, NA, NA, 400, 120, NA, NA…
## $ effort_distance_km <dbl> NA, NA, NA, NA, NA, NA, NA, 64.374, 40.234, …
## $ effort_area_ha <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ number_observers <int> 1, NA, NA, NA, NA, NA, NA, 5, 15, NA, NA, NA…
## $ all_species_reported <lgl> TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALS…
## $ group_identifier <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ trip_comments <chr> NA, NA, "Joel McNeal entering accepted recor…# Filtrar las listas de verificación con información de distancia
checklists_con_distancia <- checklists %>%
filter(!is.na(effort_distance_km))
# Crear un histograma
histograma_distancias <- ggplot(checklists_con_distancia, aes(x = effort_distance_km)) +
geom_histogram(binwidth = 10, fill = "blue", color = "black", alpha = 0.7) +
labs(title = "Distribución de Distancias Recorridas",
x = "Distancia Recorrida (km)",
y = "Frecuencia") +
theme_minimal()
# Mostrar el histograma
print(histograma_distancias)
Contiene todas las observaciones (donde anotan todo lo que ven, incluso
si no ven nada).
f_ebd <- "C:/Users/victo/OneDrive - PUJ Cali/Documentos/Maestría/2. Proyecto aplicado/Replica datos eBird/ebd_US-GA_woothr_smp_relOct-2023/ebd_US-GA_woothr_smp_relOct-2023.txt"
observations <- read_ebd(f_ebd)
glimpse(observations)## Rows: 41,432
## Columns: 48
## $ checklist_id <chr> "G10000432", "G10001259", "G10002209", "G100…
## $ global_unique_identifier <chr> "URN:CornellLabOfOrnithology:EBIRD:OBS168581…
## $ last_edited_date <chr> "2023-04-16 08:18:20.744008", "2023-05-08 08…
## $ taxonomic_order <dbl> 27880, 27880, 27880, 27880, 27880, 27880, 27…
## $ category <chr> "species", "species", "species", "species", …
## $ taxon_concept_id <chr> "avibase-8E1D9327", "avibase-8E1D9327", "avi…
## $ common_name <chr> "Wood Thrush", "Wood Thrush", "Wood Thrush",…
## $ scientific_name <chr> "Hylocichla mustelina", "Hylocichla mustelin…
## $ exotic_code <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ observation_count <chr> "1", "2", "2", "2", "4", "1", "1", "2", "1",…
## $ breeding_code <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ breeding_category <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ behavior_code <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ age_sex <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ country <chr> "United States", "United States", "United St…
## $ country_code <chr> "US", "US", "US", "US", "US", "US", "US", "U…
## $ state <chr> "Georgia", "Georgia", "Georgia", "Georgia", …
## $ state_code <chr> "US-GA", "US-GA", "US-GA", "US-GA", "US-GA",…
## $ county <chr> "Bibb", "Laurens", "Wilkes", "Bartow", "Jasp…
## $ county_code <chr> "US-GA-021", "US-GA-175", "US-GA-317", "US-G…
## $ iba_code <chr> NA, NA, NA, NA, "US-GA_3108", NA, NA, NA, NA…
## $ bcr_code <int> 27, 27, 29, 28, 29, 29, 29, 29, 29, 29, 29, …
## $ usfws_code <chr> NA, NA, NA, NA, "USFWS_315", NA, NA, NA, NA,…
## $ atlas_block <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ locality <chr> "Bondsview Rd., Macon", "Wendell Dixon Rd. D…
## $ locality_id <chr> "L1286522", "L19345770", "L23606711", "L1118…
## $ locality_type <chr> "H", "P", "P", "H", "H", "H", "H", "H", "H",…
## $ latitude <dbl> 32.78092, 32.65948, 33.83174, 34.16456, 33.1…
## $ longitude <dbl> -83.57094, -83.02279, -82.80653, -84.73369, …
## $ observation_date <date> 2023-04-15, 2023-04-15, 2023-04-15, 2023-04…
## $ time_observations_started <chr> "08:18:00", "10:00:00", "09:11:00", "08:00:0…
## $ observer_id <chr> "obsr139202,obsr617813,obsr198476,obsr267268…
## $ sampling_event_identifier <chr> "S133806681,S133806682,S133857533,S133858448…
## $ protocol_type <chr> "Traveling", "Traveling", "Traveling", "Trav…
## $ protocol_code <chr> "P22", "P22", "P22", "P22", "P22", "P22", "P…
## $ project_code <chr> "EBIRD", "EBIRD", "EBIRD", "EBIRD", "EBIRD",…
## $ duration_minutes <int> 57, 45, 102, 287, 170, 99, 85, 71, 32, 10, 4…
## $ effort_distance_km <dbl> 2.849, 9.656, 8.005, 5.404, 11.955, 5.501, 2…
## $ effort_area_ha <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ number_observers <int> 11, 2, 3, 12, 2, 3, 2, 2, 1, 1, 2, 2, 2, 2, …
## $ all_species_reported <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TR…
## $ group_identifier <chr> "G10000432", "G10001259", "G10002209", "G100…
## $ has_media <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ approved <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TR…
## $ reviewed <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ reason <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ trip_comments <chr> NA, "Section of a 9.35 mile run I did today.…
## $ species_comments <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …Listas de verificación compartidas
checklists_shared <- read_sampling(f_sed, unique = FALSE)
# identify shared checklists
checklists_shared |>
filter(!is.na(group_identifier)) |>
arrange(group_identifier) |>
select(sampling_event_identifier, group_identifier)## [1] 1249905## [1] 1120132## [1] "S78945522" "S100039674" "S7892023" "S7892189" "S7891643"
## [6] "S7891963"## [1] "G7637089" "G7637100" "G7637324" "G7637248" "G7637403" "G7640468"# importar uno de los conjuntos de datos de ejemplo de auk sin acumular taxonomía
obs_ex <- system.file("extdata/ebd-rollup-ex.txt", package = "auk") |>
read_ebd(rollup = FALSE)
# rollup taxonomy
obs_ex_rollup <- auk_rollup(obs_ex)
# identificar las categorías taxonómicas presentes en cada conjunto de datos
unique(obs_ex$category)## [1] "domestic" "form" "hybrid" "intergrade" "slash"
## [6] "spuh" "species" "issf"## [1] "species"#sin resumen, hay cuatro observaciones
obs_ex |>
filter(common_name == "Yellow-rumped Warbler") |>
select(checklist_id, category, common_name, subspecies_common_name,
observation_count)#con resumen, se han combinado
obs_ex_rollup |>
filter(common_name == "Yellow-rumped Warbler") |>
select(checklist_id, category, common_name, observation_count)Filtrado para esrudiar región y temporada
El código facilita la selección exclusiva de la información relevante sobre el Wood Thrush en Georgia, abarcando los últimos 10 años y limitándose a las listas de verificación completas durante el mes de junio. Este proceso se asemeja a la analogía de tener un libro extenso y extraer únicamente las páginas esenciales para el estudio en cuestión.
# filter the checklist data
checklists <- checklists |>
filter(all_species_reported,
between(year(observation_date), 2014, 2023),
month(observation_date) == 6)
# filter the observation data
observations <- observations |>
filter(all_species_reported,
between(year(observation_date), 2014, 2023),
month(observation_date) == 6)Se está asegurando de incluir únicamente las observaciones de aves en tierra, excluyendo aquellas en agua. Esto se logra mediante la aplicación de un límite geográfico especial que define la ubicación de Georgia y considerando solo las observaciones cercanas. El proceso se asemeja a aplicar un filtro en un mapa para retener únicamente la información de interés. Para llevar a cabo esta tarea, se utiliza un “polígono límite”, una forma específica en el mapa que delimita con precisión la ubicación de Georgia, y se establece un criterio de inclusión para observaciones dentro de un radio de 1 kilómetro de distancia de Georgia.
# convert checklist locations to points geometries
checklists_sf <- checklists |>
select(checklist_id, latitude, longitude) |>
st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
# boundary of study region, buffered by 1 km
study_region_buffered <- read_sf("data/gis-data.gpkg", layer = "ne_states") |>
filter(state_code == "US-GA") |>
st_transform(crs = st_crs(checklists_sf)) |>
st_buffer(dist = 1000)
# spatially subset the checklists to those in the study region
in_region <- checklists_sf[study_region_buffered, ]
# join to checklists and observations to remove checklists outside region
checklists <- semi_join(checklists, in_region, by = "checklist_id")
observations <- semi_join(observations, in_region, by = "checklist_id")Se esta siendo cuidadosos y solo se guarda la información que todos están de acuerdo en haber visto, para que todo sea claro y no haya confusión. ¡Es como asegurarse de que todos estan hablando de lo mismo!
# remove observations without matching checklists
observations <- semi_join(observations, checklists, by = "checklist_id")En resumen, aunque se observe que alguien ha registrado avistamientos de un pájaro en varias ocasiones, si no se tiene conocimiento de cuántas salidas se realizaron en total, no es posible afirmar con certeza si ese pájaro es común o raro. Se asemeja a la situación en la que tus amigos salen en numerosas ocasiones, pero solo registran la presencia de ese pájaro en algunas de esas salidas.
El procedimiento para realizar esto se denomina “relleno cero” de los datos. Se utiliza una herramienta denominada auk_zerofill() para combinar la lista de avistamientos de aves con la lista completa de observaciones, permitiendo así obtener datos que indican qué pájaros fueron avistados y cuáles no.
A veces, cuando tus amigos escriben las horas o las distancias que caminaron, pueden estar escritas de una manera que no es tan fácil de entender. Entonces, usamos otra herramienta para convertir esas horas en números más fáciles de manejar (de 0 a 24), poner la distancia en 0 si están quietos, y hasta calculamos la velocidad.
Y, oh, a veces tus amigos escriben una “X” en lugar de decir cuántos pájaros vieron. Así que, para que sea más fácil trabajar con los datos, convertimos esas “X” en un espacio en blanco especial llamado “NA”, que significa que no sabemos cuántos pájaros vieron.
# function to convert time observation to hours since midnight
time_to_decimal <- function(x) {
x <- hms(x, quiet = TRUE)
hour(x) + minute(x) / 60 + second(x) / 3600
}
# clean up variables
zf <- zf |>
mutate(
# convert count to integer and X to NA
# ignore the warning "NAs introduced by coercion"
observation_count = as.integer(observation_count),
# effort_distance_km to 0 for stationary counts
effort_distance_km = if_else(protocol_type == "Stationary",
0, effort_distance_km),
# convert duration to hours
effort_hours = duration_minutes / 60,
# speed km/h
effort_speed_kmph = effort_distance_km / effort_hours,
# convert time to decimal hours since midnight
hours_of_day = time_to_decimal(time_observations_started),
# split date into year and day of year
year = year(observation_date),
day_of_year = yday(observation_date)
)## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `observation_count = as.integer(observation_count)`.
## Caused by warning:
## ! NAs introducidos por coerciónCada persona cuenta los pájaros de manera diferente, y algunos caminan mucho más que otros. Entonces, para que todos podamos contar los pájaros de una manera más justa, decidimos poner algunas reglas.
Decimos que solo vamos a contar los pájaros cuando:
Caminamos menos de 6 horas. Caminamos menos de 10 kilómetros. No vamos muy rápido, menos de 100 km/h. Hay 10 personas o menos contando pájaros juntas. Estas reglas nos ayudan a asegurarnos de que todos estemos contando los pájaros de manera similar.
# additional filtering
zf_filtered <- zf |>
filter(protocol_type %in% c("Stationary", "Traveling"),
effort_hours <= 6,
effort_distance_km <= 10,
effort_speed_kmph <= 100,
number_observers <= 10)La gran mayoría de las listas de verificación están muy por debajo del límite de 6 horas, y más de la mitad tienen menos de una hora de duración.
ggplot(zf_filtered) +
aes(x = effort_hours) +
geom_histogram(binwidth = 0.5,
aes(y = after_stat(count / sum(count)))) +
scale_y_continuous(limits = c(0, NA), labels = scales::label_percent()) +
labs(x = "Duration [hours]",
y = "% of eBird checklists",
title = "Distribution of eBird checklist duration")
Conjunto de Entrenamiento (80%):
Es como el conjunto de preguntas que practicamos. Aquí, el modelo aprende cómo responder a diferentes situaciones. Usamos el 80% de nuestros datos para esto.
Conjunto de Pruebas (20%):
Este es nuestro conjunto de preguntas reales. Usamos el 20% restante de nuestros datos para ver qué tan bien puede hacer el modelo en preguntas que no ha visto antes.
zf_filtered$type <- if_else(runif(nrow(zf_filtered)) <= 0.8, "train", "test")
# confirm the proportion in each set is correct
table(zf_filtered$type) / nrow(zf_filtered)##
## test train
## 0.1997576 0.8002424Queremos hacer espacio y quedarnos solo con lo que realmente necesitamos.
checklists <- zf_filtered |>
select(checklist_id, observer_id, type,
observation_count, species_observed,
state_code, locality_id, latitude, longitude,
protocol_type, all_species_reported,
observation_date, year, day_of_year,
hours_of_day,
effort_hours, effort_distance_km, effort_speed_kmph,
number_observers)
write_csv(checklists, "data/checklists-zf_woothr_jun_us-ga.csv", na = "")Vale la pena explorar el conjunto de datos para ver con qué estamos trabajando. Este mapa utiliza datos SIG disponibles.
Los datos usando mapas y gráficos
Mapa del Sesgo Espacial:
Imagina que tienes un mapa gigante de Georgia y colocas un punto cada vez que alguien ve aves. Pero, ¡sorpresa!, hay muchos puntos alrededor de Atlanta, la ciudad más grande. Esto es lo que llamamos “sesgo espacial”. Significa que más personas en esa área informan aves que en otros lugares. Puede ser porque hay más observadores allí o porque es más fácil ver aves en la ciudad. Este mapa nos ayuda a entender cómo están distribuidos los datos en el espacio.
Histogramas y Gráficos de Frecuencia:
Ahora, pensemos en las personas que salen a observar aves. Hacemos gráficos que nos muestran cosas como:
¿A qué hora del día la mayoría de la gente sale? (histograma) ¿Cuánto tiempo pasan observando aves? (histograma) ¿Cómo cambia la probabilidad de ver una especie según cuánto tiempo y esfuerzo ponen? (gráfico de frecuencia de detección) Es como si estuviéramos estudiando cómo la gente busca tesoros. ¿Van más de día o de noche? ¿Cuánto tiempo buscan? ¿Hay lugares donde encuentran tesoros más fácilmente?
# load gis data
ne_land <- read_sf("data/gis-data.gpkg", "ne_land") |>
st_geometry()
ne_country_lines <- read_sf("data/gis-data.gpkg", "ne_country_lines") |>
st_geometry()
ne_state_lines <- read_sf("data/gis-data.gpkg", "ne_state_lines") |>
st_geometry()
study_region <- read_sf("data/gis-data.gpkg", "ne_states") |>
filter(state_code == "US-GA") |>
st_geometry()
# prepare ebird data for mapping
checklists_sf <- checklists |>
# convert to spatial points
st_as_sf(coords = c("longitude", "latitude"), crs = 4326) |>
select(species_observed)
# map
par(mar = c(0.25, 0.25, 4, 0.25))
# set up plot area
plot(st_geometry(checklists_sf),
main = "Wood Thrush eBird Observations\nJune 2014-2023",
col = NA, border = NA)
# contextual gis data
plot(ne_land, col = "#cfcfcf", border = "#888888", lwd = 0.5, add = TRUE)
plot(study_region, col = "#e6e6e6", border = NA, add = TRUE)
plot(ne_state_lines, col = "#ffffff", lwd = 0.75, add = TRUE)
plot(ne_country_lines, col = "#ffffff", lwd = 1.5, add = TRUE)
# ebird observations
# not observed
plot(filter(checklists_sf, !species_observed),
pch = 19, cex = 0.1, col = alpha("#555555", 0.25),
add = TRUE)
# observed
plot(filter(checklists_sf, species_observed),
pch = 19, cex = 0.3, col = alpha("#4daf4a", 1),
add = TRUE)
# legend
legend("bottomright", bty = "n",
col = c("#555555", "#4daf4a"),
legend = c("eBird checklist", "Wood Thrush sighting"),
pch = 19)
box()
¡Vamos a explorar la hora del día en la que la gente ve aves!
Picos de Detección:
Piensa en las aves como tesoros. Algunas aves son más fáciles de encontrar en ciertos momentos del día. Por ejemplo, temprano en la mañana, durante el coro del amanecer, es como el mejor momento para buscar tesoros.
Intervalos de 1 Hora:
Ahora, para entender cuándo la gente encuentra más aves, dividimos el día en pedacitos de 1 hora. Es como dividir el día en muchas partes pequeñas para ver en cuál de ellas es más probable encontrar tesoros.
Condiciones para Graficar:
No queremos hacer trampa, así que solo mostraremos las horas en las que al menos 100 personas han buscado tesoros. Esto asegura que nuestras estimaciones sean confiables.
# summarize data by hourly bins
breaks <- seq(0, 24)
labels <- breaks[-length(breaks)] + diff(breaks) / 2
checklists_time <- checklists |>
mutate(hour_bins = cut(hours_of_day,
breaks = breaks,
labels = labels,
include.lowest = TRUE),
hour_bins = as.numeric(as.character(hour_bins))) |>
group_by(hour_bins) |>
summarise(n_checklists = n(),
n_detected = sum(species_observed),
det_freq = mean(species_observed))
# histogram
g_tod_hist <- ggplot(checklists_time) +
aes(x = hour_bins, y = n_checklists) +
geom_segment(aes(xend = hour_bins, y = 0, yend = n_checklists),
color = "grey50") +
geom_point() +
scale_x_continuous(breaks = seq(0, 24, by = 3), limits = c(0, 24)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Hours since midnight",
y = "# checklists",
title = "Distribution of observation start times")
# frequency of detection
g_tod_freq <- ggplot(checklists_time |> filter(n_checklists > 100)) +
aes(x = hour_bins, y = det_freq) +
geom_line() +
geom_point() +
scale_x_continuous(breaks = seq(0, 24, by = 3), limits = c(0, 24)) +
scale_y_continuous(labels = scales::percent) +
labs(x = "Hours since midnight",
y = "% checklists with detections",
title = "Detection frequency")
# combine
grid.arrange(g_tod_hist, g_tod_freq)
Como era de esperar, la detectabilidad del zorzal común es mayor
temprano en la mañana y disminuye rápidamente a medida que avanza el
día. La mayoría de los envíos de listas de verificación también se
realizan por la mañana; sin embargo, hay un número razonable de listas
de verificación entre las 6 a. m. y las 9 p. m. Es en esta región donde
las estimaciones de nuestro modelo serán más confiables.
Ahora, es como si estuviéramos estudiando a diferentes buscadores de tesoros. Algunos buscan durante solo un rato, mientras que otros buscan durante mucho más tiempo. Queremos ver cómo varía esto entre las personas que observan aves.
# summarize data by hour long bins
breaks <- seq(0, 6)
labels <- breaks[-length(breaks)] + diff(breaks) / 2
checklists_duration <- checklists |>
mutate(duration_bins = cut(effort_hours,
breaks = breaks,
labels = labels,
include.lowest = TRUE),
duration_bins = as.numeric(as.character(duration_bins))) |>
group_by(duration_bins) |>
summarise(n_checklists = n(),
n_detected = sum(species_observed),
det_freq = mean(species_observed))
# histogram
g_duration_hist <- ggplot(checklists_duration) +
aes(x = duration_bins, y = n_checklists) +
geom_segment(aes(xend = duration_bins, y = 0, yend = n_checklists),
color = "grey50") +
geom_point() +
scale_x_continuous(breaks = breaks) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Checklist duration [hours]",
y = "# checklists",
title = "Distribution of checklist durations")
# frequency of detection
g_duration_freq <- ggplot(checklists_duration |> filter(n_checklists > 100)) +
aes(x = duration_bins, y = det_freq) +
geom_line() +
geom_point() +
scale_x_continuous(breaks = breaks) +
scale_y_continuous(labels = scales::percent) +
labs(x = "Checklist duration [hours]",
y = "% checklists with detections",
title = "Detection frequency")
# combine
grid.arrange(g_duration_hist, g_duration_freq)
La mayoría de las listas de verificación duran una hora o menos y hay una rápida disminución en la frecuencia de las listas de verificación a medida que aumenta la duración. Pero aquí hay algo interesante: a medida que las búsquedas son más largas, hay una mayor probabilidad de encontrar al zorzal. Es como si las personas que buscan durante más tiempo tienen más posibilidades de encontrar a nuestro amigo el zorzal. Aunque generalmente, más tiempo significa más oportunidades de detectar aves, aquí notamos algo especial. Después de alrededor de 3.5 horas, la probabilidad de detectar al zorzal disminuye.
A continuación, esperamos a priori que cuanto mayor sea la distancia que recorra alguien, mayor será la probabilidad de encontrar al menos un Zorzal Bosque.
# summarize data by 1 km bins
breaks <- seq(0, 10)
labels <- breaks[-length(breaks)] + diff(breaks) / 2
checklists_dist <- checklists |>
mutate(dist_bins = cut(effort_distance_km,
breaks = breaks,
labels = labels,
include.lowest = TRUE),
dist_bins = as.numeric(as.character(dist_bins))) |>
group_by(dist_bins) |>
summarise(n_checklists = n(),
n_detected = sum(species_observed),
det_freq = mean(species_observed))
# histogram
g_dist_hist <- ggplot(checklists_dist) +
aes(x = dist_bins, y = n_checklists) +
geom_segment(aes(xend = dist_bins, y = 0, yend = n_checklists),
color = "grey50") +
geom_point() +
scale_x_continuous(breaks = breaks) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Distance travelled [km]",
y = "# checklists",
title = "Distribution of distance travelled")
# frequency of detection
g_dist_freq <- ggplot(checklists_dist |> filter(n_checklists > 100)) +
aes(x = dist_bins, y = det_freq) +
geom_line() +
geom_point() +
scale_x_continuous(breaks = breaks) +
scale_y_continuous(labels = scales::percent) +
labs(x = "Distance travelled [km]",
y = "% checklists with detections",
title = "Detection frequency")
# combine
grid.arrange(g_dist_hist, g_dist_freq)
La mayoría de las observaciones provienen de listas de verificación
cortas, es decir, donde las personas no caminan mucho, ¡menos de medio
kilómetro!Y aquí hay algo afortunado: dado que la mayoría de las listas
de verificación están en áreas pequeñas, podemos asumir que la
variabilidad en el hábitat (el tipo de entorno natural) en esas áreas no
es tan grande. Esto facilita nuestro trabajo al analizar y entender el
hábitat.
No entendí el segundo gráfico
Hay un aumento en la frecuencia de detección cuando se recorren mas kilometros.No entiendo porque cuando recorren menos que es la mayoria de las listas de verifcación hay una frecuencia de detección tan baja.
Lista de Verificación:
¿Qué es?: Es un registro de todas las especies de aves que se observan en un lugar y momento específicos.
Lista de Detección:
¿Qué es?: Registra solo las especies de aves que fueron detectadas de manera segura y confirmada.
A continuación hablaremos sobre el número de observadores:
Al principio, pensábamos que más amigos significaría más informes, descubrimos que hay un límite. Establecimos la regla de no considerar grupos con más de 30 amigos para que nuestros datos sean más confiables. ¡Así evitamos que los grupos muy grandes afecten nuestras conclusiones sobre las aves que están buscando!
# summarize data
breaks <- seq(0, 10)
labels <- seq(1, 10)
checklists_obs <- checklists |>
mutate(obs_bins = cut(number_observers,
breaks = breaks,
label = labels,
include.lowest = TRUE),
obs_bins = as.numeric(as.character(obs_bins))) |>
group_by(obs_bins) |>
summarise(n_checklists = n(),
n_detected = sum(species_observed),
det_freq = mean(species_observed))
# histogram
g_obs_hist <- ggplot(checklists_obs) +
aes(x = obs_bins, y = n_checklists) +
geom_segment(aes(xend = obs_bins, y = 0, yend = n_checklists),
color = "grey50") +
geom_point() +
scale_x_continuous(breaks = breaks) +
scale_y_continuous(labels = scales::comma) +
labs(x = "# observers",
y = "# checklists",
title = "Distribution of the number of observers")
# frequency of detection
g_obs_freq <- ggplot(checklists_obs |> filter(n_checklists > 100)) +
aes(x = obs_bins, y = det_freq) +
geom_line() +
geom_point() +
scale_x_continuous(breaks = breaks) +
scale_y_continuous(labels = scales::percent) +
labs(x = "# observers",
y = "% checklists with detections",
title = "Detection frequency")
# combine
grid.arrange(g_obs_hist, g_obs_freq)
La mayoría de las listas de verificación tienen uno o dos observadores y
parece haber un aumento en la frecuencia de detección con más
observadores. Sin embargo, es difícil distinguir un patrón discernible
en el ruido aquí, probablemente porque hay muy pocas listas de
verificación con más de 3 observadores.
Variables ambientales
Los modelos de distribución de especies son como detectives que buscan pistas para entender dónde viven los animales y por qué. Usan información sobre el lugar donde se encuentran y cosas como la altura del terreno, la vegetación y el agua para prever dónde podrían estar en lugares que aún no han sido explorados.
Cobertura del suelo
Usaremos datos especiales llamados “cobertura terrestre” para entender qué hay en cada parte del mapa.Estos datos provienen de algo llamado “MODIS MCD12Q1 v006”, que es como un súper detector de cómo se ve la tierra desde arriba. Cubre todo el mundo y nos muestra detalles cada 500 metros.
En R, usaremos el terrapaquete para trabajar con datos ráster.
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
## (status 2 uses the sf package in place of rgdal)## terra 1.7.71##
## Attaching package: 'tidyr'## The following object is masked from 'package:terra':
##
## extract## udunits database from C:/Users/victo/AppData/Local/R/win-library/4.3/units/share/udunits/udunits2.xml## Loading required package: viridisLite# load and inspect the landcover data
landcover <- rast("C:/Users/victo/OneDrive - PUJ Cali/Documentos/Maestría/2. Proyecto aplicado/Replica datos eBird/ebird-best-practices-data/data-raw/landcover_mcd12q1_umd_us-ga_2014-2022.tif")
print(landcover)## class : SpatRaster
## dimensions : 1174, 1249, 9 (nrow, ncol, nlyr)
## resolution : 463.3127, 463.3127 (x, y)
## extent : -8132528, -7553851, 3362724, 3906653 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs
## source : landcover_mcd12q1_umd_us-ga_2014-2022.tif
## names : 2014, 2015, 2016, 2017, 2018, 2019, ...
## min values : 0, 0, 0, 0, 0, 0, ...
## max values : 15, 15, 15, 15, 15, 15, ...# map the data for 2022
plot(as.factor(landcover[["2022"]]),
main = "MODIS Landcover 2022",
axes = FALSE)
La siguiente tabla sirve para comprender el mapa de la cobertura del
suelo.
lc_classes <- read_csv("C:/Users/victo/OneDrive - PUJ Cali/Documentos/Maestría/2. Proyecto aplicado/Replica datos eBird/ebird-best-practices-data/data-raw/mcd12q1_umd_classes.csv")## Rows: 15 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): name, label, description
## dbl (1): class
##
## ℹ 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.| class | name | label | description |
|---|---|---|---|
| 0 | Water bodies | water | At least 60% of area is covered by permanent water bodies. |
| 1 | Evergreen Needleleaf Forests | evergreen_needleleaf | Dominated by evergreen conifer trees (canopy >2m). Tree cover >60%. |
| 2 | Evergreen Broadleaf Forests | evergreen_broadleaf | Dominated by evergreen broadleaf and palmate trees (canopy >2m). Tree cover >60%. |
| 3 | Deciduous Needleleaf Forests | deciduous_needleleaf | Dominated by deciduous needleleaf (e.g. larch) trees (canopy >2m). Tree cover >60%. |
| 4 | Deciduous Broadleaf Forests | deciduous_broadleaf | Dominated by deciduous broadleaf trees (canopy >2m). Tree cover >60%. |
| 5 | Mixed Forests | mixed_forest | Dominated by neither deciduous nor evergreen (40-60% of each) tree type (canopy >2m). Tree cover >60%. |
| 6 | Closed Shrublands | closed_shrubland | Dominated by woody perennials (1-2m height) >60% cover. |
| 7 | Open Shrublands | open_shrubland | Dominated by woody perennials (1-2m height) 10-60% cover. |
| 8 | Woody Savannas | woody_savanna | Tree cover 30-60% (canopy >2m). |
| 9 | Savannas | savanna | Tree cover 10-30% (canopy >2m). |
| 10 | Grasslands | grassland | Dominated by herbaceous annuals (<2m). |
| 12 | Croplands | cropland | At least 60% of area is cultivated cropland. |
| 13 | Urban and Built-up Lands | urban | At least 30% impervious surface area including building materials, asphalt, and vehicles. |
| 15 | Non-Vegetated Lands | nonvegetated | At least 60% of area is non-vegetated barren (sand, rock, soil) or permanent snow and ice with less than 10% vegetation. |
| 255 | Unclassified | unclassified | Has not received a map label because of missing inputs. |
Las ubicaciones donde viste los pájaros a veces no son exactas, por lo que es mejor usar la información sobre cómo es el hábitat en el área cercana, no solo en ese punto específico. Para hacer esto, vamos a imaginar un círculo de 3 kilómetros de diámetro alrededor de cada lugar donde viste un pájaro. Existen muchas maneras de decir cómo es un vecindario, y en este caso, usaremos dos medidas sencillas: el porcentaje de tierra cubierta (pland), que nos dice qué tan común es cada tipo de paisaje (como bosques, pastizales o agua), y la densidad de bordes (ed), que nos dice cuánto borde hay entre esos lugares diferentes.
# ebird checklist locations
checklists <- read_csv("data/checklists-zf_woothr_jun_us-ga.csv") |>
# landcover data not availble for the full period, so we use the closest year
mutate(year_lc = as.character(pmin(year, 2022)))## Rows: 47863 Columns: 19
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (6): checklist_id, observer_id, type, state_code, locality_id, protoco...
## dbl (10): observation_count, latitude, longitude, year, day_of_year, hours_...
## lgl (2): species_observed, all_species_reported
## date (1): observation_date
##
## ℹ 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.# generate circular neighborhoods for all checklists
checklists_sf <- checklists |>
# identify unique location/year combinations
distinct(locality_id, year_lc, latitude, longitude) |>
# generate a 3 km neighborhoods
st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
buffers <- st_buffer(checklists_sf, dist = set_units(1.5, "km"))Después de encontrar esos lugares especiales y poner círculos mágicos alrededor de ellos, ahora queremos saber cómo es la tierra en esos círculos.Así que, para cada lugar especial, cortamos y cubrimos la información de cómo es la tierra (cobertura terrestre) que corresponde al mismo año en que se hizo la lista de aves. uego, usamos una herramienta llamada calculate_lsm() del paquete landscapemetrics. Es como un superpoder que nos ayuda a calcular dos cosas: el porcentaje de tierra cubierta (pland) y la densidad de bordes (edm) dentro de cada vecindario circular.
lsm <- NULL
for (i in seq_len(nrow(buffers))) {
buffer_i <- st_transform(buffers[i, ], crs = crs(landcover))
year <- as.character(buffer_i$year_lc)
# crop and mask landcover raster so all values outside buffer are missing
lsm[[i]] <- crop(landcover[[year]], buffer_i) |>
mask(buffer_i) |>
# calcualte landscape metrics
calculate_lsm(level = "class", metric = c("pland", "ed")) |>
# add variables to uniquely identify each point
mutate(locality_id = buffer_i$locality_id,
year_lc = buffer_i$year_lc) |>
select(locality_id, year_lc, class, metric, value)
}
lsm <- bind_rows(lsm)Hay un pequeño detalle: algunas clases de cobertura terrestre podrían no estar en la zona mágica alrededor de cada lugar. Entonces, para todas esas clases que no vimos en el vecindario, les diremos que su porcentaje de tierra cubierta y su densidad de bordes son cero, ¡como si no estuvieran allí!
Además, en lugar de usar números extraños para hablar sobre las clases de cobertura terrestre, usaremos nombres mágicos que nos dirán exactamente qué tipo de terreno es. Y todo esto lo haremos gracias a un archivo especial que nos da nombres importantes para las clases de cobertura terrestre.
lsm_wide <- lsm |>
# fill missing classes with zeros
complete(nesting(locality_id, year_lc),
class = lc_classes$class,
metric = c("ed", "pland"),
fill = list(value = 0)) |>
# bring in more descriptive names
inner_join(select(lc_classes, class, label), by = "class") |>
# transform from long to wide format
pivot_wider(values_from = value,
names_from = c(class, label, metric),
names_glue = "{metric}_c{str_pad(class, 2, pad = '0')}_{label}",
names_sort = TRUE) |>
arrange(locality_id, year_lc)Elevación
Ahora, queremos entender cómo es la elevación del suelo, es decir, qué tan alto o bajo está el terreno. Esto es importante porque afecta dónde eligen vivir las aves.Calcularemos dos cosas mágicas para cada lugar especial: la media y la desviación estándar de la elevación en un vecindario circular de 3 km alrededor de cada lugar donde viste aves.
# elevation raster
elevation <- rast("C:/Users/victo/OneDrive - PUJ Cali/Documentos/Maestría/2. Proyecto aplicado/Replica datos eBird/ebird-best-practices-data/data-raw/elevation_gmted_1km_us-ga.tif")
# mean and standard deviation within each circular neighborhood
elev_buffer <- exact_extract(elevation, buffers, fun = c("mean", "stdev"),
progress = FALSE) |>
# add variables to uniquely identify each point
mutate(locality_id = buffers$locality_id, year_lc = buffers$year_lc) |>
select(locality_id, year_lc,
elevation_mean = mean,
elevation_sd = stdev)Ahora, combinemos las variables de cobertura terrestre y elevación, unámoslas nuevamente a los datos de la lista de verificación y guárdelas para usarlas como predictores de modelos en los próximos capítulos.
# combine elevation and landcover
env_variables <- inner_join(elev_buffer, lsm_wide,
by = c("locality_id", "year_lc"))
# attach and expand to checklists
env_variables <- checklists |>
select(checklist_id, locality_id, year_lc) |>
inner_join(env_variables, by = c("locality_id", "year_lc")) |>
select(-locality_id, -year_lc)
# save to csv, dropping any rows with missing variables
write_csv(drop_na(env_variables),
"data/environmental-variables_checklists_jun_us-ga.csv",
na = "")Cuadrícula de predicción
Ahora, queremos hacer predicciones para toda nuestra área de estudio, para hacer esto, necesitamos algo especial llamado “cuadrícula de predicción”. Imagina que estás dividiendo todo el territorio en pequeños cuadros, como un gigantesco rompecabezas. Cada cuadro nos dirá cómo es el terreno en ese pedazo.
# lambert's azimuthal equal area projection for georgia
laea_crs <- st_crs("+proj=laea +lat_0=33.2 +lon_0=-83.7")
# study region: georgia
study_region <- read_sf("data/gis-data.gpkg", layer = "ne_states") |>
filter(state_code == "US-GA") |>
st_transform(crs = laea_crs)
# create a raster template covering the region with 3 km resolution
r <- rast(study_region, res = c(3000, 3000))
# fill the raster with 1s inside the study region
r <- rasterize(study_region, r, values = 1) |>
setNames("study_region")
# save for later use
r <- writeRaster(r, "data/prediction-grid_us-ga.tif",
overwrite = TRUE,
gdal = "COMPRESS=DEFLATE")A continuación, extraemos las coordenadas de los centros de las celdas del ráster que acabamos de crear, las convertimos en sfentidades de puntos y las almacenamos en zonas de influencia para generar vecindades circulares de 3 km.
# generate neighborhoods for the prediction grid cell centers
buffers_pg <- as.data.frame(r, cells = TRUE, xy = TRUE) |>
select(cell_id = cell, x, y) |>
st_as_sf(coords = c("x", "y"), crs = laea_crs, remove = FALSE) |>
st_transform(crs = 4326) |>
st_buffer(set_units(3, "km"))Ahora podemos calcular las variables de cobertura terrestre y elevación exactamente como lo hicimos para las listas de eBird en las dos secciones anteriores.En este caso utilizamos el año más reciente de datos de cobertura terrestre (es decir, 2022).
# estimate landscape metrics for each cell in the prediction grid
lsm_pg <- NULL
for (i in seq_len(nrow(buffers_pg))) {
buffer_i <- st_transform(buffers_pg[i, ], crs = crs(landcover))
# crop and mask landcover raster so all values outside buffer are missing
lsm_pg[[i]] <- crop(landcover[["2022"]], buffer_i) |>
mask(buffer_i) |>
# calcualte landscape metrics
calculate_lsm(level = "class", metric = c("pland", "ed")) |>
# add variable to uniquely identify each point
mutate(cell_id = buffer_i$cell_id) |>
select(cell_id, class, metric, value)
}
lsm_pg <- bind_rows(lsm_pg)
# transform to wide format
lsm_wide_pg <- lsm_pg |>
# fill missing classes with zeros
complete(cell_id,
class = lc_classes$class,
metric = c("ed", "pland"),
fill = list(value = 0)) |>
# bring in more descriptive names
inner_join(select(lc_classes, class, label), by = "class") |>
# transform from long to wide format
pivot_wider(values_from = value,
names_from = c(class, label, metric),
names_glue = "{metric}_c{str_pad(class, 2, pad = '0')}_{label}",
names_sort = TRUE,
values_fill = 0) |>
arrange(cell_id)Y ahora la media y la desviación estándar de elevación.
elev_buffer_pg <- exact_extract(elevation, buffers_pg,
fun = c("mean", "stdev"),
progress = FALSE) |>
# add variables to uniquely identify each point
mutate(cell_id = buffers_pg$cell_id) |>
select(cell_id, elevation_mean = mean, elevation_sd = stdev)Finalmente, combinamos las variables de cobertura terrestre y elevación y las guardamos en CSV.
# combine landcover and elevation
env_variables_pg <- inner_join(elev_buffer_pg, lsm_wide_pg, by = "cell_id")
# attach the xy coordinates of the cell centers
env_variables_pg <- buffers_pg |>
st_drop_geometry() |>
select(cell_id, x, y) |>
inner_join(env_variables_pg, by = "cell_id")
# save as csv, dropping any rows with missing variables
write_csv(drop_na(env_variables_pg),
"data/environmental-variables_prediction-grid_us-ga.csv",
na = "")Siempre podemos usar la plantilla ráster para convertir estas variables ambientales a un formato espacial, por ejemplo, si queremos mapearlas. Veamos cómo funciona esto para el porcentaje de cobertura de bosque caducifolio latifoliado (clase 4).
forest_cover <- env_variables_pg |>
# convert to spatial features
st_as_sf(coords = c("x", "y"), crs = laea_crs) |>
# rasterize points
rasterize(r, field = "pland_c04_deciduous_broadleaf")
# make a map
par(mar = c(0.25, 0.25, 2, 0.25))
plot(forest_cover,
axes = FALSE, box = FALSE, col = viridis(10),
main = "Deciduous Broadleaf Forest (% cover)")