A new approach to quantify grazing pressure under mediterranean pastoral systems using GIS and remote sensing

ABSTRACT Pastoral systems based on grazing itineraries are very common along the Mediterranean region and provide an opportunity to manage the fuel load and reduce fire risk in the ecosystem. Daily and seasonal movements of flocks bring on different grazing pressure (GP) (sheep ha−1) over the landscape. This study presents an approach to model sheep GP under a Mediterranean pastoral system in the Northeast of Portugal. The GP in a given location depends on the distance from the stable to the parish border, the distance to the stable, the heads of livestock and their preference for land use and cover (LUC). The geoprocessing we applied in this study integrated several spatial databases: the latest official Portuguese vector mapping of land use and cover (COS2015) and administrative boundaries (CAOP2018), the livestock stables location, and Sentinel-2 Multispectral Images. During the geoprocessing, the stocking density (SD) (sheep ha−1) were calculated and spatially interpolated. Homogeneous LUC units (permanent crops; annual crops; forest; shrubland; grassland; water bodies) were obtained by Random Forest supervised classification algorithm (kappa = 89.3%; global accuracy = 91.2%). Boolean overlapping of the LUC classes obtained by the supervised classifier with the mask created from COS2015 provides the potentially grazed LUC classes. Integrating LUC preferences with SD allows calculating and mapping the GP. The most common GP class is 0–0.25 sheep ha−1. Seeing the GP per LUC class, a value of 1.84 sheep ha−1 was found in permanent crops, 1.73 in annual crops, and 1.25 in grassland, 0.88 in grazed forests and 0.84 in shrublands. The GP modelling and mapping could assist in the implementation of herding programmes aimed at reducing fire hazards at a parish or at a regional scale.


Introduction
Since 1955 rural and livestock populations have shown a decreasing evolution in Portugal (Torres-Manso et al. 2014). This trend is also common to most countries along the northern Mediterranean rim (Hoggart and Paniagua 2001). The abandonment of rural areas particularly in mountain regions has led to an increase in shrublands and a decrease in pastures and cultivated areas (Ganteaume et al. 2013). As a result, there has a rise the fuel loads (Moreira, Francisco, and Ferreira 2001;Vilar et al. 2016), and in the dominance of simplified firevulnerable landscapes composed of extensive homogeneous communities (Lasanta et al. 2006;Vega-García et al. 2010;Azevedo, César, and Castro 2011). These landscape changes associated with an increase in extreme weather events have contributed to the dramatic rise in wildfires in most European countries along the northern Mediterranean basin (Ruiz-Mirazo, Martínez-Fernández, and Vega-García 2012).
Currently, Portugal has one of the highest forest fire risk rankings in Europe (Ruiz-Mirazo, Martínez-Fernández, and Vega-García 2012), and it has been the subject of numerous reviews by fire experts, who have made several recommendations for implementation of effective fuel management techniques. It seems that Portugal needs to improve various aspects in the fight against wildfires, among which experts point out the need to create a structural fire defence system of fuelbreaks and reduce the fuel load in critical areas (Beighley and Hyde 2009;Moreira et al. 2011;Marino et al. 2014). Therefore, it is essential to know how the livestock is distributed over the landscape at a parish scale (the smallest administrative unit in Portugal) so as to design grazing measures to promote fire prevention programmes. Apart from this, a new era of fire regimes demands a more comprehensive and balanced strategy, in which extensive livestock could be used as an instrument to reduce wildfire hazards (Fernandes et al. 2013). According to (Ruiz-Mirazo, Martínez-Fernández, and Vega-García 2012), the low number of livestock currently observed in many mountain areas requires a grazing intensification at a local scale, i.e. along fuelbreaks.
Pastoral farming systems have always adapted to the seasonal availability of forage resources and climate variability by moving animals. Through the ages, farmers and their animals have developed adaptive capacities to face spatial and temporal scarcity as well as the variability of pastoral resources (Schoenbaum et al. 2017). Animal mobility is one of the main strategies used in different forms: nomadism, transhumance or moving herds around outside the structural limits of the farm, a very common system in Northern Mediterranean areas, where small ruminant producers are landless farmers. In this way, in the North of Portugal, the livestock system is supported by grazing itineraries over unfenced or uncultivated parish territory (Baumont et al. 2000). The flocks are conducted by shepherds, crossing land with different uses considering their availability and feeding value, that is, their preference for land use and cover (LUC) (Castro 2004).
Electronic tracking devices like GNSS receivers ('global navigation satellite systems') have proven effective in monitoring the movement of flocks, allowing estimation of the time devoted to grazing in different land use and cover (Castro 2004;Bishop-Hurley et al. 2007;Putfarken et al. 2008;Tomkins et al. 2009;Brosh et al. 2010) and gathering detailed information related to animal movements and their spatial distribution (Perotto-Baldivieso et al. 2012). The combination of geographic information systems (GIS) with remotely sensed data allows the spatial distribution of grazing movements (Kawamura et al. 2005;Paz-Kagan et al. 2016). Sensors used in unmanned aerial vehicles (drones) are constantly evolving (e.g. Parrot Sequoia or Micasense). Although they do not yet match the spectral resolution of sensors mounted on orbiting platforms, such as Sentinel-2, Michez et al. (2019) successfully used multispectral information acquired by drones in grazing management, with great temporal, spatial and spectral advantages.
Digital remote sensing data has long been successfully used for LUC classification since the Landsat launch, which is frequently used on regional scales (Yuan et al. 2005;Weng 2006, 2007;De Fries et al. 1998;Wulder et al. 2018). Recently, the multispectral remote sensed data from Sentinel-2 has been providing accurate details about the landscape regarding spectral, temporal and spatial resolutions (Wulder et al. 2018;Steinhausen et al. 2018;Tavares et al. 2019), mainly by supervised digital classification. Maximum likelihood (ML) is one of the most common supervised classifiers of remotely sensed data. Over time, large classifier algorithms have been developed, such as the Random Forest (RF), and others (Spectral Angle Mapper, Machine Learning algorithms, such as the Support Vector Machines, and Neural Networks) (Kuching 2007;Son et al. 2018;Tavares et al. 2019).
The objective of this study is to design an approach for GP modelling, in pastoral systems in Mediterranean areas, using the north of Portugal as a case of study, by the integration of multispectral data from Sentinel-2 and geoprocessing. Our study sheds light on the potential of GIS and remote sensing to help construct a baseline for current GP in this region, which could be used in connection with policies or management strategies to address wildfires risk prevention.

Study area
The area under study is the Trás-os-Montes region, located in the northeast of Portugal (40°57ʹ N, 6°19ʹ W and 41°58ʹ N, 7°12ʹ W) (Figure 1), covering about 660,000 km 2 . This region is characterized by high climatic diversity. Due to being an interior region of the country, it has few maritime influences, but its diversified geomorphology results in a wide variation of rainfall, in the order of 400 mm in areas of the Upper Douro Valley and 1500 mm in mountain areas in the north (Castro and Fernández 2014). The region elevation ranges from 400 to 1486 m. The landscape is characterized by a mosaic of oak and pine forests, alternated with annual and permanent crops, and other non-forest land uses such as shrublands, grasslands and built-up areas/settlements (Azevedo, César, and Castro 2011).

Satellite data collection and processing
This study used data from Sentinel-2 Multispectral imagery of the European Space Agency. Sentinel-2 data Level 1 C, July 2018, were downloaded from the agency's Copernicus Sentinels Scientific Data Hub. We needed four Sentinel-2 scenes to cover the entire extent of the Bragança region. The cloud cover in all the scenes was very low ( Table 1). The monitoring LUC period coincided with the end of the cereal harvesting period (July) when the clean atmosphere is more frequent (cloud cover inferior to 2%) and it is possible to distinguish different LUC classes (e.g. permanent and annual crops).
The multispectral Sentinel-2 data in the visible (three bands in 10 m), near-infrared (five bands in 20 m) and shortwave infrared (two bands in 60 m) were used during the classification. The Level-1 C product was processed for radiometric and geometric corrections, including orthorectification and spatial registration on a global reference system (Sentinel-2Team 2015). The Sen2Cor processor was used to perform an atmospheric correction that transforms a Top-Of-Atmosphere (TOA) level-1 C into a Bottom-Of-Atmosphere (BOA) level-2A product (Müller-Wilm 2015). After that, all the bands for each scene were resampled to a 10 m × 10 m resolution and the four scenes were combined into a mosaic projected to UTM Zone 29 WGS84 and clipped by region of interest (ROI) extent. A flowchart of the data pre-processing and processing steps used in the study are summarized in Figure 2.
For the assessment of GP, eight land use and cover classes were considered: annual crops, permanent crops, grasslands, shrublands, grazed forests, ungrazed forests, urban and water bodies ( Table 2). The urban class was obtained directly from the latest official Portuguese land use and cover (DGT. 2015, Carta Administrativa Oficial de Portugal, e Carta de Uso e Ocupação Do Solo de Portugal Continental Para 2015). This thematic cartography was produced by national authorities based on the visual interpretation of high spatial resolution orthorectified aerial images. The urban area has stabilized regionally, according to periodic evaluations since 1995 (1995,2007,2010    In contrast, satellite image classification of annual crops has advantages compared to that obtained by COS2015 especially with images acquired at full crop maturity or at harvest time (e.g. cereals). Through satellite imagery, it is possible to classify the entire forest class, but not in grazed and ungrazed forests. In turn, the ungrazed forest class, which corresponds to the coniferous and eucalyptus forest, was easily identified by photo interpretation. It is therefore suggested to classify satellite images into global forest class and to remove ungrazed forest by a spatial overlap of COS2015. A homogeneous set of LUC class pixels were selected, and the algorithm was trained to classify the data based on this 'training data' (TD). Similarly, 'validation data' (VD) was created to validate the classification. The selection of these clusters for TD and VD creation was achieved by drawing polygons over homogeneous clusters based on the Sentinel-2 RGB colour composition (supplemental material), aided by the high-resolution images of Google Earth and COS2015 (disambiguation). A TD size of 49,946 pixels was obtained, distributed as follows: 8,455, 6,401, 8,218, 7,436, 11,581 and 7,855, for annual crops, permanent crops, shrublands, forest (grazed and ungrazed), and water bodies respectively.
Seeking for the best result of LUC to generate a GP map, RF supervised classifiers was used to discriminate the six classes initially proposed (Table 2). Vegetation and soil indices were calculated (Table 3) and integrated with spectral data (B 2 , B 3 , B 4 , B 5 , B 6 , B 7 ,B 8 , B 8a , B 11 , B 12 ). The RF algorithm was applied to spectral data considering two iterations: with (1st iteration) and without (2nd iteration) vegetation and soil indices. RF was first developed by Breiman (2001) as an ensemble machine learning algorithm based on regression trees. In both iterations, we used 500 trees as standard values and selected training data (TD) vectors ( Figure 2). Post-classification refinements were applied to the resulting image in a 3 × 3 kernel before the accuracy assessment ( Figure 2).
The discriminating quality of the classifications was evaluated using two coefficients: the overall accuracy, and the Cohan's Kappa coefficient (Cohen 1960;Foody 2002), which were derived from confusion matrices confronting the classification in lines with the ground truth reference data in columns (i, j). After getting the best classification, we used COS2015 to divide the forest into two classes (Table 2) and to add the urban class, which resulted in eight classes to be integrated into the GP modelling process.

Grazing pressure modelling
The sheep global grazing pressure (GP), understood as the number of sheep per ha, in the district of Bragança (northeast of Portugal), is modulated in the GIS environment with the multi-criteria evaluation to estimate the GP on different LUC classes Areas with shrubs and low grown trees as well as sparsely tree-covered areas. Land cover with coarse and irregular texture with a medium dark tone. Grazed Forests and Ungrazed Forests Forests are LUC classes with tree crown cover of more than 10 percent. Hardwoods exhibit higher overall reflectance than conifers. Forest areas, in general, have a coarse texture. In mixed forests, the texture is very coarse and irregular. Forest are divided in two classes, grazed and ungrazed forest. Urban Artificial surfaces intended for activities related to human societies. The appearance of urban areas in images may vary widely, depending on whether they are predominantly horizontal or vertical, continuous or discontinuous.

Water bodies
Rivers, reservoirs and lakes. Dark tone due to reduced overall reflectance.
( Table 2). The multi-criteria evaluation depends on: 1) the flock's daily movement distance (< 5 km); 2) the parish administrative limits restriction; and 3) the LUC preference by sheep given by the time spent in each land use in relation to the total time of the grazing itinerary (Castro 2004). Geoprocessing integrates optical data from Sentinel-2 (Table 3), Portugal LUC spatial database in 2015 (COS2015), the parishes administrative boundaries (CAOP2018), and a spatial database (from the 0108_OTSA_2_E project) which reports the stable location and livestock per stable ( Figure 1) (IPB 2010;DGT 2015). During the geoprocessing, the SD was initially calculated for the global area. Regarding the grazing preference (%) on each potentially grazed LUC classes obtained by (Castro 2004) for the same region, and the percentage of the ungrazed LUC classes (Water bodies and Urban areas, 13%), the SD was proportionally allocated to the grazed LUC classes (87%) ( Table 4). The Preference correction factor (F) was calculated by the product of preference value and the percentage of ungrazed classes. The LUC used by livestock (%) come out by the addition of Grazing preferences (%) and F (%).
The maximum distance potential travelled was considered at a radius of 5 km, defined from the stable and it not crossing the limits of the parish (Castro 2008;Castro and Fernández 2016). Considering the location of each stable and number of sheep per stable  (Pouget et al. 1990)  in the study region (Figure 1), the Euclidean distance to the nearest parish border was determined with a spatial resolution of 25 m (supplemental material). The values of each cell created in the raster surface were extracted to each stable. Similarly, the calculation of the distance around each stable was processed. As the herd moves away from its stable, its livestock is diluted over the gradually larger area. This effect was estimated through the multiple concentric rings (or buffers) of geometric increasing radius (100, 200, 400, 800, 1600, 3200, 6400 m), from each stable, considering the maximum distance a herd reaches on a daily itinerary before returning to the stable (Castro and Fernández 2016). The polygons derived gives the SD. We assumed that the herds have a variable itinerary and each location can be grazed by different herds of the same parish on the same day or on different days, consequently overlapping polygons will occur. The SD is accumulated where multiple polygons overlap (supplemental material).
Considering the average distance between neighbouring stables (850 m) and the way the SD of itinerant grazing smoothly varies across the landscape, we have created a network of systematic distribution points, with intervals of 1 km × 1 km, covering the entire study area, each point assuming the SD value of the intersected polygon (supplemental material). From this regularly distributed point network, representative of the regional variability of the SD, a spatial interpolation was performed for the entire area of interest, by ordinary kriging, an 'exact interpolator' where the interpolated value respects the original value at a known value position.
A GIS interpolation model should be a statistically good representation of spatial variability, with Mean Error (ME) (or mean prediction error) less than 1. The quality of predicting increased to near zero values of ME and Mean Standardized Error (MS), considering Standardized error (Stdd_Error) as the prediction errors (individual error) divided by their prediction standard errors (individual stderror). Beyond that, the prediction errors are normally distributed. The spatial resolution of the surface should be empirically tested considering the spatial variability (Houlding 2000;Schabenberger and Gotway 2017). By analysing prediction errors, the best spatial resolution for the spatially interpolated surface can be searched. It was generated a continuous surface of SD with a spatial resolution of 25 m. SD interpolation produces a continuous surface for the entire study area, including ungrazed classes (urban, water bodies, and ungrazed forest). The SD in ungrazed classes was later transferred to the grazed classes (Table 4) producing a new continuous surface of corrected stocking density (SD') (sheep ha −1 ). To generate the global GP dynamic through the Bragança region, we integrated SD' with the best result of the LUC classification using weighted spatial overlap geoprocessing techniques.

LUC classification
The RF supervised classifiers can be effectively applied to the Sentinel-2 images to segment the landscape in the desired LUC classes. The importance of the classifiers, as well as the training data (TD) and validation data (VD) should be highlighted (Vieira and Mather 2000;Foody 2002;Lillesand, Kiefer, and Chipman 2014). The interpretation of contingency matrices allows to highlight the confusion between the LUC classes and the overall accuracy of classification (Foody 2002;Congalton and Green 2008), helping to decide which algorithm is best and which bands and their transformations should be used in a particular landscape segmentation process.
The RF supervised classifiers demonstrated high discriminant capacity (Table 5), providing an overall accuracy of 90.8% and a kappa coefficient of 88.8%. The result of this classifier was further improved with the integration of vegetation and soil indices (RF 2st iteration) (Table 5), resulting in an overall accuracy of 91.2% and a kappa coefficient of 89.3%. The quality of the classification according to kappa was excellent (Cohen 1960;Rosenfield and Fitzpatrick 1986;Foody 2002;Congalton and Green 2008). In addition, the marginal accuracies were also excellent, from the point of view of both the producer and the user, with marginal kappa values always higher than 80%.
The results show that the choice made at the date of the images and the use of soil and vegetation indices were effective, and the quality of the classification was improved both overall and individually (Table 5). It is known that 'NDVI' enhances the contrast of light absorption by water from soil and vegetation, and has the advantage of being correlated with the level of chlorophyll activity in plant coverages (Rouse et al. 1974;Silleos et al. 2006;Saadat et al. 2011). It was also noted that the Brightness Index (BI) contrasts the bare soil (Escadafal 1989;Major, Baret, and Guyot 1990), especially in July, after annual crops harvesting, presenting this class a high reflectance in the visible and near-infrared region of electromagnetic radiation. In addition, the integration of these indices into RF 2st iteration improves the discrimination of the types of cover, particularly between the shrubland and the grassland.
The Annual crops (mainly cereals) are harvested in July (dry summer season), leaving bare brown soils behind, allowing their identification through the interpretation of remote sensing data. Knowledge of the dynamics of land use and cover in a region is an important factor to improve the accuracy of classification of images (Khatami and Mountrakis 2018), as well as a judicious choice of the date of acquisition. Also, Laborte, Maunahan, and Hijmans (2010) recommend the acquisition of images during the dry period.
A LUC map was created through the result of the RF 2st iteration classifier, also including the ungrazed classes (ungrazed forest and urban) that had previously been obtained through the COS2015 (Figure 3). The resulting LUC map shows that shrublands were the most abundant (188,315 ha, 28.5% of total area), with the remaining land uses being: permanent crops (161,724 ha, 24.5%); grassland (135,277 ha, 20.5%); ungrazed forest (69,781 ha, 10.6%); annual crops (49,814 ha, 7.5%); grazed forest (41,626 ha, 6.3%); urban (8,221 ha, 1.2%); and water bodies (5,242 ha, 0.8%) (Figures  3 and 4). According to Azevedo, César, and Castro (2011), the landscape of the Bragança region has gone through relevant modifications over the last decades, agriculture land has been replaced by shrublands and forests have decreased progressively from 22% of the landscape in 1958 to less than 5% in 2005. Also, annual crops have registered a gradual abandonment or reconversion to permanent crops, especially orchards (chestnut, almond, and olive trees).

Grazing pressure
Using the administrative boundaries (CAOP2018) and livestock stables location, was calculated the SD in vector format. Through kriging spatial interpolation, a smooth dispersion surface was produced of the SD. The spatial resolution of 25 m provided the smallest prediction errors and was adopted for all geoprocessing. The kriging model to interpolate SD is statistically a good representation of spatial variability, with prediction errors normally distributed, ME equal to 0.514 and MS equal to −0.018.
The integration of remote sensing data with advanced GIS tools allows us to model GP and to provide useful information on grazing intensity at a given point (Figure 4). The GP depends on the distance to the stable, the LUC class and the number of sheep per stable. Considering the whole territory, the most common GP is about 0-0.25 sheep ha −1 (25.5%). Seeing the GP per LUC, we found mean values of 1.84 sheep ha −1 in permanent crops, 1.73 in annual crops, 1.25 in grassland, 0.88 in the grazed forest and 0.84 in shrublands ( Figure 4).
Also, it should be stated that the agrarian LUC matrix is a GP that is compatible with a mosaic landscape, while the other LUC related to woodland matrixes showed a lower value of GP, which is not compatible with a fire-resistant landscape. In an attempt to mitigate the effects of the current low GP in remote Portuguese areas, the national forest authorities have introduced a fuel loads management programme by grazing which is focused on fuelbreaks and requires at least 1.4 sheep or goats per hectare (ICNF 2018). Similarly, in southern France, Etienne et al. (1993) reported a value between 0.6 and 1.4 sheep ha −1 in open areas, and up to 1.65 for forest ecosystems. In Spain, Evlagon et al. (2010) state 1.25-2.01 goat's ha −1 in woody vegetation, and 0.98-1.40 goat's ha −1 in grassland vegetation. Our results reflect the importance of understanding the spatial pattern of GP to reinforce this Portuguese programme, since grazing remains more effective to maintain biodiversity than other techniques such as prescribed burning (Pastro, Dickman, and Letnic 2011).
Land management practices (LMPs) improve the resilience of ecosystems to multiple disturbances such as wildfires (Jucker et al. 2018). The wildfire's dynamics in the Mediterranean basin is closely linked to shrub communities (Naveh 1994;Luis, Raventós, and González 2006;Moreira, Francisco, and Ferreira 2001;Vilar et al. 2016). In 2003 catastrophic fires in Portugal, 156,767 ha of shrublands and 132,407 ha of grasslands were burnt, emitting more than 20,000 tCO2 to the atmosphere, contributing to climate change (Ameray 2018). The increase of grazing intensity could be a useful practice to prevent wildfires in this type of Mediterranean pastoral system. This had already been claimed by Ruiz-Mirazo and Jabier (2012), in his studies to understand the linkages between pastoral wildfires and land cover patterns in similar managed Mediterranean landscapes.

Conclusions
Free access to current geographic databases and remote sensing images was particularly important for the success of SD and GP modelling. Classifying remote sensing imagery remains a challenge that depends on numerous factors such as the type and the acquisition date of the images, the classifier, and the complexity of the landscape in the study area. When we consider these factors, Sentinel-2 satellite images provide a LUC classification with high accuracy. RF showed to be a powerful machine learning classifier, particularly with the combination of the multispectral data with other classification criteria (e.g. vegetation and soil indices). From the digital processing of images Sentinel-2 resulted in a spatial resolution of 10 m, higher than the geoprocessing model used (25 m), opening up the possibility of applying the model to local studies in more detail.
The most common GP in the northeast landscapes of Portugal is low, about 0-0.25 sheep ha −1 . Seeing the GP per LUC, the mean values are 1.84 sheep ha −1 in permanent crops, 1.73 sheep ha −1 in annual crops, 1.25 sheep ha −1 in grassland, 0.88 sheep ha −1 in the grazed forest and 0.84 sheep ha −1 in shrublands. The agrarian LUC matrix is a GP reasonably compatible with a mosaic, while the LUC related to the woodland matrix has an extremely low GP, and therefore not compatible with a fire-resistant landscape. The GP modelling and mapping could assist in the implementation of herding programmes aimed at reducing fire hazards at both parish and regional scales.