Grazers in a Protected Spanish Landscape

Author

Martin Cooper

Published

May 7, 2025

Context

Grasslands, covering large areas in multiple continents, are significant ecosystems that support biodiversity, regulate climate and provide many other ecosystem services (Lewińska et al. 2023). The long history of human settlement and societal development reinforces the importance of grasslands (Jouven et al. 2010). The Mediterranean biome in the South West of Europe is part of a global biodiversity hotspot with a large number of endemic plant species and habitats heavily influenced by historic grazing patterns (Porqueddu et al. 2016).

Not only are Mediterranean biomes influenced by grazing, but the dynamics of the habitats are influenced through fire disturbance. Both fire and grazing alter ecosystem dynamics and their combined influence creates distinct outcomes that differ to their individual effect (Noy‐Meir 1995). In the Mediterranean, the interaction between grazers and the habitat is well studied and changes to grazing has been shown to lead to changes to soil carbon capture (Peco et al. 2017), photosynthetic productivity and forage quality (Castillo-Garcia et al. 2022) and defense against invasion (Díaz Cando et al. 2025). These impacts go beyond the grasslands as the grasslands exisst within a wider landscape matrix, what happens in the grasslands will impact the distribution of both habitat types,# and the movement of organisms between them.

As with many habitats globally, there are factors that place Mediterranean habitats at risk. These include climate change, especially drought and other hydrological changes with an associated fire risk; land use changes and degradation; and the positive and negative influences of grazers (Zhang et al. 2025; Buisson et al. 2020).

Not only does grazing activity help with current conservation, both directly and via the indirect interactions with other risks, but they can also play key roles in ecosystem restoration. Since the translation of Frans Vera’s Grazing Ecology and Forest History (Vera 2000), interest in how vertebrate herbivores shape ecosystems has led to numerous studies and practical applications. Practical projects provide significant insight into the movement and ecology of herbivores. Two well known examples include Oostvaardersplassen in the Netherlands (Staatsbosbeheer 2025) and Knepp Wildland in the UK (Estate 2025).

Whilst these two well known sites are using herbivores to restore a variety of habitats directly, the understanding of how grazers support plant community assembly continues as well. The figure below, reproduced from Buisson et al. (2020), shows that not only is herbivory an important disturbance feature, the impact of grazers goes beyond this to include aspects such as seed dispersal too.

Figure 1. Large vertebrates identified as an efficient techniques for restoration, but also that more research is needed. Reproduced from Buisson et al. (2020)

All this makes it clear that understanding how grazers move and behave is important, they are significant actors in the ever shifting matrix of habitats. To support this, more than 90 of individuals in sheep and goat herds have been tagged and tracked by researchers at the Service for the Evaluation, Restoration and protection of Mediterranean Agroecosystems in Granada, Spain. The herds involved in the tagging were located on expansive livestock farms within protected areas of Andalusia, with data collected from October 2019 until August 2023 (Pérez-Luque et al. 2024).

Here I present the results of an analysis of part of this significant data set to examine how 9 goats utilize the area in which they have been tracked. This includes looking at their use of resources in the landscape and an investigation to see if we can highlight their movemement and across the landscape.

Analytical Approach

Before any analysis can begin, data must be retrieved, explored and basic statistics prepared. This will allow us to prepare a clean, appropriately formatted dataset ready for analysis. The data was retrieved from Movebank (study ID 3088763011) and exploration was undertaken following the details outlined in Smith et al. (2020).

workflow A Load Packages B Retrieve Data A->B C Load CSV (if needed) B->C D Clean & Transform Data - Select key columns - Format timestamps - Handle missing values - Remove duplicates C->D F Compute Summary Stats - Tracking duration - Individuals per species - Unique locations D->F E Exploratory Data Analysis - Structure & summary - Check distributions - Detect outliers - Visual representation G Filter by Species E->G F->E H Save Filtered Data G->H

The packages used for this were:

This workflow identified and eliminated missing and duplicates GPS locations and time stamps for ease of further analysis.

The dataset broke down as follows:

Following this analysis, the dataset was split to only contain the data around the goats. The same steps were then repeated for the goat only data. As duplicates and NAs had been removed prior to the split, these steps were not repeated.

Structure and Summary

Lets take a closer look at the steps in the workflow shown above. First a visual inspection of the data.

Summary of Goat Tracking Data
Species ID Total Days Tracked Total Locations
Capra hircus AV276 69 1045
Capra hircus AV277 309 6208
Capra hircus AV337 271 4515
Capra hircus AV338 266 3523
Capra hircus AV580 76 17734
Capra hircus AV581 173 41357
Capra hircus AV582 52 13432
Capra hircus AV583 171 24608
Capra hircus AV613 163 28865

There is significant variation between the counf of times and locations of the tracked goats which will impact the quality and validity of the analysis. As we are interested in the areas they are using, rather than comparisons between individuals, this variation is not as much of an issue as it could be. However, where there is less data there is a risk of underestimating areas.

Outliers

From the exploration of the full dataset, the data appears to be in 2 distinct areas. It is not possible to see outliers from scale that shows all the data, so a closer look is required. To do this the data is split into the two location groups.

Add seasonal data

As there are seasonal changes to foraging behavior (Abraham et al. 2022), and changes to plant response depending on the season (Petit Bon et al. 2021), it is worth investigating how the goats space use, resource selection and behavior vary by season. To do this, we can extract the data from the tracking information and bin into seasons. For months 12, 1, 2 were grouped as Winter,3, 4, 5 as Spring, 6, 7, 8 as Summer and 9, 10, 11 as Autumn. Whilst climate data may be the more appropriate classification, this data was not collected. Changes to the phenology of the ecosystem could make this classification less suitable in future studeis as well (Gordo and Sanz 2010).

Space Use

In this study, we are wanting to know what areas the goats are not only visiting, but the areas that they are using. After all, it is only by being there that the animals are going to be having an impact on the habitat. We are fortunate that the data we have is highly detailed GPS data, resulting in fine-scale location and time data. The study by Shakeri et al. (2021) found that goats in a different, more seasonal habitats, were consistent in their home range and utalisation distribution across years, but that there was variation by season. Both the importance of the area used and the differences found by Shakeri et al. (2021) suggest that a utalisation distribution would be more helpful that a standard home range. This leaves us with the following questions that this part of our analysis can answer:

  • How large an area do the goats use?

This will tell us about the scale of influence the animals have on the habitats.

  • Are there differences in the use by season?

This will help inform if certain habitats need more protection at different times of year. This is especially true for our study area as the level of stress (drought, heat etc.) will likely influence the impact of grazing. Further, as we know animals are key seed distributers, knowing where they may be dispersing plants seeds over winter and spring could help with conservation of habitats.

AKDE Utalisation Distributions

Due to the granularity in our data resulting from GPS location data, we can use an Autocorrelated Kernal Density Estimator (AKDE). Unlike KDEs, the AKDE incorporates the temporal nature of granular data, allowing for a much more accurate estimate of the space used (Fleming et al. 2015).Whist this could be computationally expensive, by focussing only on the goats in this analysis, this is achievable.

A significant advantage of using AKDEs is that the location points of the tracked animal are not treated as independent which is an assumption of KDE modelling. Since each animal is not independent of itself, or of the location it was previously in, it is clear that each location is not independent, violating this assumption. Another assumption of KDEs is that the data represents a random point and extrapolates the estimate of space use from this, again this is not true. With our more granular data, and sufficient computational resources to handle the extra data, an AKDE is the right method to use (Fleming et al. 2015).

To support our questions, we can do two AKDE analyses 1) with all data together to examine the goats space, and 2) dividing the data by season to see if there are differences in the space used by season. The workflow for this analysis is shown below.

AKDEworkflow A Load Packages - ctmm, dplyr, sf, tmap B Compute Full-Year AKDE UDs - Convert telemetry - Fit models - Generate UD - Save shapefiles A->B C Complete Seasonal AKDE UDs - Change data filters to incorperate seasons - Convert telemetry - Fit models - Generate UD - Save shapefiles B->C D Extract & Save Summary - Store UD areas - Compare seasonal shifts - Save CSV C->D

The parameters were set as follows for the AKDE analysis: Model fitting was completed with an inital guess, followed by a maximum likelihood estimation.

The packages used for this were:

Resource Use

Whilst AKDEs provide insights into where the goats are, they tell us nothing about the type of habitats the goats are in, which they prefer to stay in and move through. This leaves us with the question: - Do goats show any habitat type preference? As noted previously, goats, and other herbivores, can act as major influences on community structure.Our question ins:

  • Do goats show any habitat preference?

This is an important area to investigate, as the impact of grazers on Mediterranean habitats is well known and is crucial for the conservation of these habitats (Balata et al. 2022). Understanding, and identifying, habitats of heavy and low impact by grazers is important assessing whether certain habitats are disproportionately used and therefore require conservation interventions. In both the conservation of current target answer the question, we can evaluate the number of times the goats are found in each habitat class compared to the abundance of these classes in the wider landscape. For this, we can use a Resource Selection Function (RSF) (Boyce et al. 2002).

For this analysis, we compare the observed goat locations with their associated habitat class. This data has been collected via remote sensing, specifically the Copernicus Sentinel-2 satellite array, accessible via the Copernicus Data Space Ecosystem (Copernicus Land Monitoring Service 2021; European Space Agency (ESA) 2015).

Resource_use_workflow A Load Packages - dplyr, sf, sp, ctmm, ggplot2 B Download Habitat Raster - Access Sentinel-2 data via Copernicus - Retrieve 10m resolution land cover A->B C Prepare Habitat Data - Load raster - Reproject to EPSG:32630 - Crop to study area - Standardize classification B->C D Recall previously generated AKDE UDs - Extract individual space use polygons C->D E Run dBBMM Analysis - Estimate movement variance - Define probability surfaces of corridor use D->E F Extract Habitat Classes - Mask raster by AKDE outputs - Pair locations with available habitat data E->F G Create Available Locations - Random points at 3× ratio - Ensure spatial integrity F->G H Implement RSF Analysis - Logistic regression - Include mixed-effects for individual variability G->H I Model Validation - Extract coefficients - Generate ROC curve and AUC metrics H->I

For this workflow the following packages were used:

Movement Patterns and Behaviours

To expand on the understanding of how the goats are using the habitat patches within the landscape we can look at their behaviors and movement patterns. As a result of our comprehensive dataset, there are a variety of different approaches we could take. As we are interested in the impact of the goats on the habitat patches, we can now form a question:

  • What are the goats doing in the different habitat patches?

Knowing if the goats are grazing or just moving in the patches they have been found in will add to our previous analysis. The movement of animals in a landscape is not random and will be impacted by factors, such as the habitat type, and the internal state of the animal.

There are a wide variety of approaches to answer this question. Here we use Behavioral Change Point Analysis (BCPA). This method was developed by Gurarie et al. (2009) and allows us to elucidate where an animal’s behavior changes, which can inform us of areas that the animal uses that warrant further investigation. For the conservation of plant communities, this is important to planning conservation action. BCPA has advantages over other methods (such as HMMs) as it is much more effective at fine scale changes. This could indicate patterns of disturbance tied to these locations or other important ecological insights.

In addition to the BCPA, we can investigate the areas that the animals frequently return to. These areas could highlight important resources, such as water, high quality grazing habitat and regularly used movement corridors. This information can be further used to inform conservation action plans targeting plant communities by allowing planners to preserve movement corridors to facilitate grazing, or identify areas where barriers to dispersal could have a disproportionate impact on grazer behavior. Recursive movement patterns show us which areas are frequently revisited and the applications of this include highlighting areas essential for maintaining movement and areas where herbivore exclusion could be required to minimize habitat damage (Berger-Tal and Bar-David 2015).

Goat_Movement_Workflow A Load Packages - move2, bcpa, recurse, sf, dplyr, ggplot2, tmap B Load Dataset - Read CSV - Convert timestamps - Create Move2 object A->B C Data Validation - Check missing values - Validate track IDs B->C D Compute Movement Metrics - Speed, displacement, turning angles C->D E Filter Outliers - Remove time gaps > 1 day D->E F BCPA Analysis - Identify behavioral shifts - Detect movement transitions E->F G Movement Recursion Analysis - Detect revisited locations - Identify site fidelity F->G H Seasonal Segmentation - Compare movement across seasons G->H I Data Visualization - Map movement trajectories - Plot recursion hotspots H->I J Export Shapefile - Convert data to sf format - Save movement paths I->J K Habitat Preference Analysis - Overlay movement with habitat layers - Perform RSF modeling J->K

For this workflow the following packages were used:

Outcomes

Space Use

We have two questions around space use for the goats:

  • How large an area do the goats use?
  • Are there differences in the use by season?

To examine the size of the area used, the utilization distributions hold the answer. The full results of the area analysis can be expanded below. This is helpful information and tells us that, all together, the goats use 5,434.719 Ha of space.

     ID Patch     name  area_ha
1 AV337     1      Min 381.2336
2 AV337     1 Estimate 456.9145
3 AV337     1      Max 537.1660
4 AV338     1      Min 367.4255
5 AV338     1 Estimate 420.2961
6 AV338     1      Max 475.1128
      ID Patch     name   area_ha
1  AV276     2      Min 2135.3462
2  AV276     2 Estimate 3050.9645
3  AV276     2      Max 4121.8047
4  AV277     2      Min 1378.4726
5  AV277     2 Estimate 1584.8293
6  AV277     2      Max 1804.7078
7  AV580     2      Min 2161.5463
8  AV580     2 Estimate 3125.3552
9  AV580     2      Max 4264.6656
10 AV581     2      Min 1089.6836
11 AV581     2 Estimate 1143.7987
12 AV581     2      Max 1198.8787
13 AV582     2      Min 1056.8347
14 AV582     2 Estimate 1356.3599
15 AV582     2      Max 1692.6251
16 AV583     2      Min 2006.5272
17 AV583     2 Estimate 2448.4701
18 AV583     2      Max 2935.9072
19 AV613     2      Min  488.8662
20 AV613     2 Estimate  538.7624
21 AV613     2      Max  656.7628

The numbers are helpful, but it is easier to understand them on a map.

There is significant variation in the areas used by all animals.

Interestingly, the individuals in patch 1 (AV337 and AV338) both have smaller UDs compared to those in patch 2, although AV613 does not fit this trend. When looking at the time and amount of locations for these 3 individuals there is a large variation. From the data we have available, we can hypothesize that goats in patch 1 do not need as large an area to meet their needs, however, we only have two data points so further data collection would be required to validate this. It may be that there were less goats in this area and therefore this would also also explain the smaller utilization distribution.

The variations in size for goats in patch 2 suggest that the difference in area is not due to the quality of the habitat area, but that other factors may be at play. This data did not contain any further information the individuals, but research on other ungulates have found that male individuals have larger ranges (Vanp’e et al. 2008).

For differences in Seasons, we can plot the seasonal AKDEs to compare. Here are the results for partch 1

Summer Map

Summer_map <- tm_shape(Summer_masked) + tm_raster(palette = “Reds”, alpha = 0.7) + tm_text(text = “Summer”, size = 2, xmod = -0.4, ymod = 0.4) +
tm_basemap(server = “OpenStreetMap”) + tm_layout(legend.show = FALSE)

Autumn Map

Autumn_map <- tm_shape(Autumn_masked) + tm_raster(palette = “Oranges”, alpha = 0.7) + tm_text(text = “Autumn”, size = 2, xmod = -0.4, ymod = 0.4) +
tm_basemap(server = “OpenStreetMap”) + tm_layout(legend.show = FALSE)

Winter Map

Winter_map <- tm_shape(Winter_masked) + tm_raster(palette = “Blues”, alpha = 0.7) + tm_text(text = “Winter”, size = 2, xmod = -0.4, ymod = 0.4) +
tm_basemap(server = “OpenStreetMap”) + tm_layout(legend.show = FALSE)

View maps separately

tmap_arrange(Spring_map, Summer_map, Autumn_map, Winter_map)



## Resource Use

## Movement Patterns and Behaviours

#### 

## Present numerical and visual results of your analyses.

# Interpretation of Main Findings

There are several important results from our analysis, to summarise:

-   There are significant differences in the space different goats use
-   There are significant differences in how the goats use the habitat by season
-   The population studies has a preference for x habitat type
-   The goats behavior was more likely to change when in x habitat and stayed consistent when on x habitat
-   There were several areas that the goats regular returned to

This all represents a large amount of information about how the goats use the habitat, and therefore the impact on habitat patches by grazers. \## Discuss the significance of your findings.

# Conservation Perspectives

## Explore the implications of your findings for conservation/management efforts.

::: {.cell}

```{.r .cell-code}
# only run once!
# Define a vector of DOIs
doi_list <- c("10.1890/14-2010.1", "10.1186/s40462-020-00229-3", 
              "10.1111/j.1365-2656.2012.01955.x", "10.3390/rs14102322", 
              "10.1016/S0304-3800(02)00200-4")

# Retrieve BibTeX entries for all DOIs
bib_entries <- GetBibEntryWithDOI(doi_list)

# Print summary or save as a BibTeX file
print(bib_entries)
WriteBib(bib_entries, file = "bibliography.bib")

:::

Code 
library(RefManageR)

# Fetch BibTeX entry using DOI
bib_entry <- GetBibEntryWithDOI("10.1111/j.1365-2486.2009.02084.x")


# Convert to BibTeX format
bib_text <- toBibtex(bib_entry)

# Append to .bib file without reading
cat("\n", bib_text, "\n", file = "ref.bib", append = TRUE)

print(bib_text)

References

Abraham, J.O. et al., 2022. Seasonal strategies differ between tropical and extratropical herbivores [online]. Journal of Animal Ecology, 91 (3), 681–692. Available at: http://dx.doi.org/10.1111/1365-2656.13651.
Balata, D., Gama, I. and Domingos, T. et al., 2022. Using satellite NDVI time-series to monitor grazing effects on vegetation productivity and phenology in heterogeneous mediterranean forests [online]. Remote Sensing, 14 (10), 2322. Available at: http://dx.doi.org/10.3390/rs14102322.
Berger-Tal, O. and Bar-David, S., 2015. Recursive movement patterns: Review and synthesis across species [online]. Ecosphere, 6 (9), 1–12. Available at: http://dx.doi.org/10.1890/ES15-00106.1.
Boyce, M.S., Vernier, P.R. and Nielsen, S.E. et al., 2002. Evaluating resource selection functions [online]. Ecological Modelling, 157 (2–3), 281–300. Available at: http://dx.doi.org/10.1016/S0304-3800(02)00200-4.
Buisson, E. et al., 2020. Key issues in northwestern mediterranean dry grassland restoration [online]. Restoration Ecology, 29 (S1). Available at: http://dx.doi.org/10.1111/rec.13258.
Castillo-Garcia, M. et al., 2022. Understanding herbivore-plant-soil feedbacks to improve grazing management on mediterranean mountain grasslands [online]. Agriculture, Ecosystems &amp; Environment, 327, 107833. Available at: http://dx.doi.org/10.1016/j.agee.2021.107833.
Copernicus Land Monitoring Service, 2021. CLCplus backbone 2021 (raster 10 m), europe, 3-yearly. Available at: https://land.copernicus.eu/en/products/clc-backbone/clc-backbone-2021.
Díaz Cando, P.E. et al., 2025. Enemy behind the gates? Predicted climate change and land‐use intensification likely speed up <scp>C4</scp> grass invasions in europe [online]. Journal of Vegetation Science, 36 (2). Available at: http://dx.doi.org/10.1111/jvs.70023.
Dowle, M. and Srinivasan, A., 2023. Data.table: Fast aggregation of large data sets [eBook]. Available at: https://CRAN.R-project.org/package=data.table.
Dunnington, D., 2023. Ggspatial: Spatial data visualizations in ggplot2 [eBook]. Available at: https://CRAN.R-project.org/package=ggspatial.
Estate, K., 2025. Knepp wildland project. Available at: https://knepp.co.uk/.
European Space Agency (ESA), 2015. Sentinel-2: High-resolution optical imagery for land monitoring. Available at: https://dataspace.copernicus.eu/explore-data/data-collections/sentinel-data/sentinel-2.
Fleming, C.H. and Calabrese, J.M., 2024. Ctmm: Continuous-time movement modeling [eBook]. Available at: https://CRAN.R-project.org/package=ctmm.
Fleming, C.H., Fagan, W.F. and Mueller, T. et al., 2015. Rigorous home range estimation with movement data: A new autocorrelated kernel density estimator [online]. Ecology, 96 (5), 1182–1188. Available at: http://dx.doi.org/10.1890/14-2010.1.
Gordo, O. and Sanz, J.J., 2010. Impact of climate change on plant phenology in mediterranean ecosystems [online]. Global Change Biology, 16 (3), 1082–1106. Available at: http://dx.doi.org/10.1111/j.1365-2486.2009.02084.x.
Grolemund, G. and Wickham, H., 2024. Lubridate: Make dealing with dates easier [eBook]. Available at: https://CRAN.R-project.org/package=lubridate.
Gurarie, E., 2014. Bcpa: Behavioral change point analysis of animal movement [eBook]. Available at: https://cran.r-project.org/package=bcpa.
Gurarie, E., Andrews, R.D. and Laidre, K.L., 2009. A novel method for identifying behavioural changes in animal movement data [online]. Ecology Letters, 12 (5), 395–408. Available at: http://dx.doi.org/10.1111/j.1461-0248.2009.01293.x.
Jouven, M. et al., 2010. Rangeland utilization in mediterranean farming systems [online]. Animal, 4 (10), 1746–1757. Available at: http://dx.doi.org/10.1017/S1751731110000996.
Kranstauber, B. and Wehrens, R., 2023. Move: Visualizing and analyzing animal movement data [eBook]. Available at: https://CRAN.R-project.org/package=move.
Lewińska, K.E. et al., 2023. Beyond “greening” and “browning”: Trends in grassland ground cover fractions across eurasia that account for spatial and temporal autocorrelation [online]. Global Change Biology, 29 (16), 4620–4637. Available at: http://dx.doi.org/10.1111/gcb.16800.
Noy‐Meir, I., 1995. Interactive effects of fire and grazing on structure and diversity of mediterranean grasslands [online]. Journal of Vegetation Science, 6 (5), 701–710. Available at: http://dx.doi.org/10.2307/3236441.
Pebesma, E., 2024. Sf: Simple features for r [eBook]. Available at: https://CRAN.R-project.org/package=sf.
Peco, B. et al., 2017. Effects of grazing abandonment on soil multifunctionality: The role of plant functional traits [online]. Agriculture, Ecosystems &amp; Environment, 249, 215–225. Available at: http://dx.doi.org/10.1016/j.agee.2017.08.013.
Pérez-Luque, A.J. et al., 2024. vreGI Virtual Research Environment Grazing Intensity: Monitoring Grazing Patterns in Andalusia Natural Protected Areas [online]. Available at: https://github.com/serpam/vreGI_db.
Petit Bon, M. et al., 2021. Forage quality in tundra grasslands under herbivory: Silicon‐based defences, nutrients and their ratios in grasses [online]. Journal of Ecology, 110 (1), 129–143. Available at: http://dx.doi.org/10.1111/1365-2745.13790.
Picardi, S., 2020. Recurse: Analysis of movement recursions [eBook]. Available at: https://cran.r-project.org/package=recurse.
Porqueddu, C. et al., 2016. Grasslands in “old world” and “new world” mediterranean‐climate zones: Past trends, current status and future research priorities [online]. Grass and Forage Science, 71 (1), 1–35. Available at: http://dx.doi.org/10.1111/gfs.12212.
R., E. and A., M., 2023. move2: An improved framework for animal movement analysis [eBook]. Available at: https://CRAN.R-project.org/package=move2.
Roger Bivand, V.G.-R., Edzer Pebesma, 2024. Sp: Classes and methods for spatial data [eBook]. Available at: https://CRAN.R-project.org/package=sp.
Shakeri, Y.N., White, K.S. and Waite, J.N., 2021. Staying close to home: Ecological constraints on space use and range fidelity in a mountain ungulate [online]. Ecology and Evolution, 11 (16), 11051–11064. Available at: http://dx.doi.org/10.1002/ece3.7893.
Smith, C., Uzal, A. and Warren, M., 2020. Statistics in r for biodiversity conservation. S.I., s.n.
Staatsbosbeheer, 2025. Oostvaardersplassen nature reserve. Available at: https://www.staatsbosbeheer.nl/uit-in-de-natuur/locaties/oostvaardersplassen.
Tennekes, M., 2023. Tmap: Thematic maps in r [eBook]. Available at: https://CRAN.R-project.org/package=tmap.
Vanp’e, C. et al., 2008. Access to mates in a territorial ungulate is determined by the size of a male’s territory, but not by its habitat quality [online]. Journal of Animal Ecology, 78 (1), 42–51. Available at: http://dx.doi.org/10.1111/j.1365-2656.2008.01467.x.
Vera, F., 2000. Grazing ecology and forest history. Wallingford, UK: CABI Publishing.
Wickham, H., 2024a. ggplot2: Elegant data visualization using the grammar of graphics [eBook]. Available at: https://CRAN.R-project.org/package=ggplot2.
Wickham, H., 2024b. Tidyverse: Tidy data tools for r [eBook]. Available at: https://CRAN.R-project.org/package=tidyverse.
Wickham, H. and François, R., 2024. Dplyr: A grammar of data manipulation [eBook]. Available at: https://CRAN.R-project.org/package=dplyr.
Zhang, Z. et al., 2025. Editorial: Climate and environmental changes in circum-mediterranean regions [online]. Frontiers in Earth Science, 13. Available at: http://dx.doi.org/10.3389/feart.2025.1594165.
Zhu, H., 2023. kableExtra: Advanced tables in r [eBook]. Available at: https://CRAN.R-project.org/package=kableExtra.