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 (Frans 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 (Knepp 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.
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 movement and across the landscape.
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).
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.
Lets take a closer look at the steps in the workflow shown above. First a visual inspection of the data.
| Species | ID | Total Days Tracked | Total Locations |
|---|---|---|---|
| Capra hircus | AV276 | 1085 | 1045 |
| Capra hircus | AV277 | 8476 | 6208 |
| Capra hircus | AV337 | 7817 | 4515 |
| Capra hircus | AV338 | 5152 | 3523 |
| Capra hircus | AV580 | 18245 | 17734 |
| Capra hircus | AV581 | 43715 | 41357 |
| Capra hircus | AV582 | 14251 | 13432 |
| Capra hircus | AV583 | 25509 | 24608 |
| Capra hircus | AV613 | 31038 | 28865 |
Table 2. Summary of data by individual
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.
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.
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 studies as well (Gordo and Sanz 2010).
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 utilisation 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 utilisation 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:
This will tell us about the scale of influence the animals have on the habitats.
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.
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.
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:
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:
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, D. 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).
For this workflow the following packages were used:
A 100m buffer was added to the UDs to ensure that variance in the GPS tags was accounted for and a 3x ratio was used as these was not a huge area to look at and more than this would result in an over saturation of random points leading to an inaccurate analysis.
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:
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. Step Selection Functions are a well used tool in conservation and ecology, however other methods better suit our question. Here we aim to use Behavioral Change Point Analysis (BCPA). This method was developed by Eliezer 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 (Oded Berger-Tal and Shirli Bar-David 2015). For this analysis the radius of the recursion was set at 1000m, with the accuracy of the GPS data being 10m, this means that patches of high return have a good chance of being identified.
For this workflow the following packages were used:
The full results of the utilisation distribution 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 |
|---|---|---|---|
| AV337 | 1 | Min | 381.234 |
| AV337 | 1 | Estimate | 456.915 |
| AV337 | 1 | Max | 537.166 |
| AV338 | 1 | Min | 367.425 |
| AV338 | 1 | Estimate | 420.296 |
| AV338 | 1 | Max | 475.113 |
| AV276 | 2 | Min | 2135.346 |
| AV276 | 2 | Estimate | 3050.965 |
| AV276 | 2 | Max | 4121.805 |
| AV277 | 2 | Min | 1378.473 |
| AV277 | 2 | Estimate | 1584.829 |
| AV277 | 2 | Max | 1804.708 |
| AV580 | 2 | Min | 2161.546 |
| AV580 | 2 | Estimate | 3125.355 |
| AV580 | 2 | Max | 4264.666 |
| AV581 | 2 | Min | 1089.684 |
| AV581 | 2 | Estimate | 1143.799 |
| AV581 | 2 | Max | 1198.879 |
| AV582 | 2 | Min | 1056.835 |
| AV582 | 2 | Estimate | 1356.360 |
| AV582 | 2 | Max | 1692.625 |
| AV583 | 2 | Min | 2006.527 |
| AV583 | 2 | Estimate | 2448.470 |
| AV583 | 2 | Max | 2935.907 |
| AV613 | 2 | Min | 488.866 |
| AV613 | 2 | Estimate | 538.762 |
| AV613 | 2 | Max | 656.763 |
Table 3. Utilisation Distribution areas by season and individual
The numbers are helpful, but it is easier to understand them on a map.
Figure 8. 95% UD Estimation for Patch 1
Figure 9. 95% UD Estimation for Patch 1
For differences in Seasons, we can plot the seasonal AKDEs to compare. Here are the results for patch 1, the areas of the different UDs are shown in the bottom corner of each plot.
Figure 10. Patch 1 95% UD estimates by season
And for patch 2:
Figure 11. Patch 1 95% UD estimates by season
Following data wrangling, 3 RSF models were fitted to assess goat habitat cover preference.
Initially, the model was fitted with all data.
The area under the curve for this model was 0.7967.
| effect | group | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | -6.932 | 0.051 | -136.761 | 0 |
| fixed | NA | patch_area | -0.082 | 0.007 | -10.947 | 0 |
| fixed | NA | isolation | 0.215 | 0.003 | 68.333 | 0 |
| fixed | NA | sealed | 5.435 | 0.033 | 166.843 | 0 |
| fixed | NA | needle_woody | 0.542 | 0.058 | 9.392 | 0 |
| fixed | NA | broadleaf_deciduous | -54.387 | 0.101 | -538.299 | 0 |
| fixed | NA | broadleaf_evergreen | 3.217 | 0.043 | 74.743 | 0 |
| fixed | NA | low_grow_woody | 3.695 | 0.029 | 127.163 | 0 |
| fixed | NA | herbaceous_perm | 5.409 | 0.027 | 197.065 | 0 |
| fixed | NA | herbaceous_period | 6.889 | 0.027 | 255.048 | 0 |
| fixed | NA | non_and_sparce_veg | 6.407 | 0.027 | 241.598 | 0 |
| ran_pars | id | sd__(Intercept) | 1.037 | NA | NA | NA |
Table 4. Summary Statistics for all habitat types
The estimate is the log-odds coefficient from the model gives us the relationship to the predictor, whilst the statistic value gives us the strength of this relationship.
As we are interested in the habitat types, another model was run removing sealed habitat and broadleaf deciduous as factors. The results of this model variation are below.
| effect | group | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | -1.512 | 0.056 | -26.930 | 0.000 |
| fixed | NA | patch_area | -0.086 | 0.008 | -10.949 | 0.000 |
| fixed | NA | isolation | 0.001 | 0.003 | 0.183 | 0.855 |
| fixed | NA | needle_woody | -4.886 | 0.050 | -97.072 | 0.000 |
| fixed | NA | broadleaf_evergreen | -2.217 | 0.034 | -66.136 | 0.000 |
| fixed | NA | low_grow_woody | -1.724 | 0.024 | -71.880 | 0.000 |
| fixed | NA | herbaceous_perm | -0.011 | 0.025 | -0.459 | 0.646 |
| fixed | NA | herbaceous_period | 1.463 | 0.022 | 67.660 | 0.000 |
| fixed | NA | non_and_sparce_veg | 0.986 | 0.021 | 46.377 | 0.000 |
| fixed | NA | water | -5.377 | 0.061 | -88.730 | 0.000 |
| ran_pars | id | sd__(Intercept) | 1.037 | NA | NA | NA |
Table 5. Summary Statistics for the reduced habitat type groups
A final evaluation was completed removing highly collinear predictors
The BCPA and Step Selection Functions were unable to be completed as the dataset had a wide range of values for both timelags between fixes and distance values which did not fit within realistic boundaries (one goat was found to have traveled at 4,326km/hr!). Despite this, some movement metrics were possible and they are displayed below:
| Statistic | Value |
|---|---|
| Min. | 3 |
| 1st Qu. | 306 |
| Median | 312 |
| Mean | 1337 |
| 3rd Qu. | 366 |
| Max. | 24548514 |
| NA’s | 9 |
Table 7. Spread of timeLags
As shown in this table, the timelags are very varied and filtering did not leave sufficient data for analysis.
| Statistic | Value |
|---|---|
| Min. | 0.000000 |
| 1st Qu. | 0.000085 |
| Median | 0.000186 |
| Mean | 0.000847 |
| 3rd Qu. | 0.000505 |
| Max. | 0.048539 |
| NA’s | 9.000000 |
Table 8. Spread of distance
As with the timelag data, the distance data was not reliable enough for analysis.
It was possible to assess how the goats were moving in their environment by examining the relative angles of movement.
Figure 15. Turning Angle Distribution
There is a clear pattern in this with the goats making decisions. The main types of behaviour being larger scale movements (peaks at 180 and -180 degrees) and feeding and resting (peak around 0).
It was not possible to perform the analysis in one go, processing was split into chunks and then recombined. The areas with highest revisits are shown on the maps below for each area. They clearly highlight some areas where the goats are found more often from the data.
The large amount of data for individuals in patch 2 means that along with the key areas, it is also possible to make out some key travel routes.
By extracting the habitat data from these highly used areas, we get the following information:
| Habitat | Patch.1 | Patch.2 |
|---|---|---|
| Broadleaf Deciduous | NA | 0.03 |
| Broadleaf Evergreen | 4.77 | 0.24 |
| Herbaceous Period | 3.78 | 12.11 |
| Herbaceous Permanent | 27.45 | 74.56 |
| Low-Grow Woody | 54.03 | 0.84 |
| Needle Woody | 1.36 | 0.06 |
| Non & Sparse Vegetation | 8.61 | 6.39 |
| Other | NA | 0.01 |
| Sealed | NA | 5.75 |
| Water | NA | 0.01 |
Table 9. Habitat cover percentages for highly revisited sites for Patch 1 & 2
It was unfortunate that not all of the planned analysis was possible. However, there are some results here that are worth exploring further.
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 also a large variation. From the data we have available, we can hypothesize that goats in patch 1 may not need as large an area to meet their needs, however, we only have two data sets so further data collection would be required to validate this.It maybe that as there a smaller number of animals in this area, explaining the smaller UD. No definite conclusion is possible for the reasons why the areas are smaller in patch 1 compared to 2.
Within patch 2, there is large variation in th size of individual animal UDs, this suggests that the difference in area is not due to the quality of the habitat area and 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).
There are also differences in the seasonal patterns of grazer areas, although as with above caution must be taken when looking at the patch 1 data. For patch 2, there is a clear reduction in the size of the area used with the summer area used (Figure 11). This change in seasonal grazing has been investigated for other species. Iris Schoenbaum et al. (2017) investigated how cattle showed differences in their habitat preference by season. Unlike the cattle in Iris Schoenbaum et al. (2017) which showed a preference for woody habitats, the goats studied here showed a preference for herbaceous habitat types. This difference highlights the importance of mixed grazing when planning conservation activities.
The cover of the habitat types in the recursive movement analysis reinforces what was found through the RSF, with the goats showing a preference for herbaceous habitat. It is interesting that patch 1 animals do not show this preference as strongly, however due to the smaller sample size, this should not be given much weight.
This answers most of the questions posed, although due to the incomplete analysis of movement data relating to behvaiour means that we have not been able to identity what the goats are doing in each area. Despite this, the combination of the turning angle distribution and the recursive movement does provide some insight. The turning angles suggest 3 main angles of movement, 2 at 180 and -180 representing directional reverses and a peak around 0 suggesting movements along linear features. When compared to the recursive sites, this is reinforced as there are clear broader patches where animals are repeatedly moving around an area and, whilst not as clear in patch 1, clear areas of directional movement. The combination of these results suggest that the goats are moving between frequently visited foraging and/or resting habitats with travel in between. These results should be taken with caution though as studies on wild goat populations in other parts of the world have found significant GPS bias for small scale movement analysis (Wells et al. 2011). Interestingly, this same study found that whilst fine scale movement data was unreliable, the habitat model did show potential. Whilst we cannot compare this directly due to the incomplete movement analysis here, it suggests that the results we have been able to gather represent important insights.
As this study, and the study by Wells et al. (2011), show, movement analysis has its value and its drawbacks as well as potential to help answer ecological questions. In this study, we have successfully identified the habitat types preferred by goats on large farms in a protected area in Andalucia, Spain. This was to to identify areas with high patterns of animal use, potentially leading to grazing and/or physical degradation of the habitat, as well as finding areas where grazing pressure and and insufficient animal faciliated seed dispersal could lead to nevgative outcomes for the habitats. insufficient graze activity.
Once areas where grazers need to be excluded, or areas where they need to be encouraged to frequent are identified, influencing the behaviour of the animals is necessary to impact on conservation outcomes. In large areas fencing and gates are not suitable and would have a negative impact on the landscape, they may also impact other ecological processes leading to undesirable outcomes. Virtual fencing has been identified as a possible solution to this challenge in other parts of the world (Harland et al. 2025). This offers potential in the area under investigation here as well as virtual fencing can be altered remotely allowing grazer managers to induce movement, and therefore habitat impact, changes. Given that the animals in this study show no impact on their movement preferences relating to patch isolation or size, the potential to encourage the animals to graze currently under grazed areas is stronger than if connectivity was vital to their movement. Despite this promise, there are challenges relating to managing goats with virtual fencing and results are currently inconclusive with significant variation between individual animals Wilms et al. (2025).
Exclusion of grazers is known to have a significant impact on grassland plant biodiversity (Xu and Guo 2015) and the impact of domestic grazers also influences wild grazer habitat use (DeGabriel et al. 2011), clearly demonstrating that grazing has many direct and indirect impacts on the biodiversity of plant habitats. Beyond plants, the presence, and species, of grazer has also been found to have an impact on bird diversity (Boyce et al. 2021). All this shows that the impact of grazers is varied and that management of grazers is important when planning conservation actions. From the results presented here, it is suggested that the status of the herbaceous habitats is assessed and that exclusion of grazers implemented where they are in poor condition. A comparative analysis of the status of those areas identified as under higher grazing pressure should be made with those where goats are not present. In a highly modified landscape like the Mediterranean, a delicate balance of activities are needed and they will largely depend on the aim of any activity. Overall, it is clear that grazers have an impact on the matrix of habitats in the landscape and their impact should be included in conservation planning.