A Multi-Temporal Object-Based Image Analysis to Detect Long-Lived Shrub Cover Changes in Drylands

Climate change and human actions condition the spatial distribution and structure of vegetation, especially in drylands. In this context, object-based image analysis (OBIA) has been used to monitor changes in vegetation, but only a few studies have related them to anthropic pressure. In this study, we assessed changes in cover, number, and shape of Ziziphus lotus shrub individuals in a coastal groundwater-dependent ecosystem in SE Spain over a period of 60 years and related them to human actions in the area. In particular, we evaluated how sand mining, groundwater extraction, and the protection of the area affect shrubs. To do this, we developed an object-based methodology that allowed us to create accurate maps (overall accuracy up to 98%) of the vegetation patches and compare the cover changes in the individuals identified in them. These changes in shrub size and shape were related to soil loss, seawater intrusion, and legal protection of the area measured by average minimum distance (AMD) and average random distance (ARD) analysis. It was found that both sand mining and seawater intrusion had a negative effect on individuals; on the contrary, the protection of the area had a positive effect on the size of the individuals’ coverage. Our findings support the use of OBIA as a successful methodology for monitoring scattered vegetation patches in drylands, key to any monitoring program aimed at vegetation preservation.


Introduction
A fundamental fact of ecological observation is that most living organisms do not show random distributions. In fact, environmental controls and anthropogenic impacts are determinants of the spatial patterns of these organisms. This implies that it is possible to know the performance of ecosystems through the study of the spatial distribution patterns of the organisms that live in them [1,2]. This is particularly important in drylands where, as a result of water scarcity and edaphic limitations, vegetation appears to form isolated patches of one or more plant individuals [3,4]. In these ecosystems, it has been observed that the spatial pattern of these patches determines key aspects of ecosystem functioning such as primary production [5], water and nutrient cycles [6], and biotic interactions [7,8]. Tools to produce accurate vegetation maps at the appropriate spatial scale over time could be very useful for gaining knowledge about the health and dynamics of dryland ecosystems.
Remote sensing has proven to be the most useful tool for monitoring changes in vegetation, as it is cost-effective, allows repeated mapping, and produces information on a large scale [9][10][11]. Within this technique, pixel-based methods are the most commonly used. However, these methods show several limitations for producing accurate maps of vegetation patches or plant individuals in drylands. First, pixel-based methods do not consider the spatial context in which the pixels are framed, making it difficult to identify isolated image elements. Second, they often result in a final overlap of such elements from automatic classifications, when the analysis is based on high spatial resolution images [12]. In drylands, the land surface is characterized by scattered vegetation in a matrix of bare soil and scattered shrubs, so contextual information is very useful for image classification [13]. Both characteristics limit the possibility of identifying and classifying patches of vegetation and individual plant elements.
Several methods have emerged as an alternative to pixel-based methods for mapping individuals or vegetation patches. For example, in the case of forests, light detection and ranging (LiDAR) and very high frequency (VHF) synthetic aperture radar (SAR) images allow the characterization of various attributes of individual trees from their three-dimensional structure (e.g., [14][15][16]). However, this method is difficult to use when vegetation shows reduced aerial volume such as in drylands. In these cases, object-based image analysis (OBIA) can be a good solution for mapping patches of vegetation and individual plants [17], particularly because there is currently a wide variety of freely available high spatial resolution orthoimages. OBIA can provide a more accurate and realistic identification of scattered vegetation in drylands because of the combined spectral information of each pixel with the spatial context [18,19]. This method has yielded good results in the monitoring of spatial patterns, functioning, and structure of vegetation in these environments [20,21].
OBIA may be particularly useful for assessing the dynamics of populations of long-lived plants of conservation concern. In this case, it is difficult and costly to assess the environmental controls of population dynamics due to their high persistence and sometimes low rate of regeneration, which requires very long-term studies [22,23]. It has been proposed that the maintenance of long-lived plant populations is the result of a balance between regeneration (replacement of individuals by recruiting new recruits) and persistence (maintenance of individuals in space, physically and temporarily), or a combination of both strategies [24,25], depending on abiotic stress and biotic competition [26]. Monitoring populations of persistent individuals over time is complicated, as there are continuous disturbances in the environment that can alter their performance [24]. However, the availability of the analysis of historical aerial orthophotographs and high spatial resolution satellite images with OBIA provides a good opportunity to reconstruct the interannual dynamics of long-lived plant populations over long periods of time, thus enabling the evaluation of changes experienced by these shrub populations.
Ziziphus lotus (L.) Lam, a long-lived shrub from Mediterranean drylands [27], shows characteristics for a multi-temporal analysis of the spatial distribution in its populations with OBIA. This shrub species depends on groundwater [28], forms fertility islands, and is considered an engineering species [29] of an ecosystem of interest for conservation at the European level (Directive 92/43/CEE). The main European population of the shrub species is located in a flat coastal area surrounded by greenhouses in the Cabo de Gata-Níjar Natural Park, Spain. This population has been affected by several threats for many decades, including sand mining [30,31], reducing the amount of sand available to develop the Z. lotus fertility islands; urban pressure [32], which has reduced the potential distribution of Z. lotus; and the expansion of intensive agriculture [33,34], responsible for the decline in the level of the aquifer's water table, which may have caused the seawater intrusion [35,36]. Since 1944, several studies have evaluated this community of Z. lotus. Shrub patterns to identify groundwater dependence [28], the formations of shrubs in dunes [37], shrub spatial aggregation and consequences for reproductive Remote Sens. 2019, 11, 2649 3 of 17 success [29], and mutual positive effects between shrubs [38] have been researched. Yet, the monitoring of the shrub population dynamics has never been studied.
Despite most of the shrub population being located within the protected area, its temporal dynamics could be affected by several human-induced disturbances. However, due to the slow growth of shrubs and the inertia in the extinction of individuals, it is difficult to assess such dynamics without considering the population structure of the shrub species over the last several decades. This work proposes the use of remote sensing methods to map the spatial distributions of shrubs and to analyze their size and shape as a means of identifying anthropic disturbances. Our guiding hypothesis was that Z. lotus, phreatophytic shrubs, were affected by soil loss and seawater intrusion that decreased their cover area. On the contrary, after the legal protection of the area in 1987, the shrubs increased their cover area. Within this framework, the objectives of this work were as follows: (i) to make precision maps of scattered shrubs from historical remotely sensed images using OBIA and (ii) to extract information on changes in the shape, size, and spatial distribution of shrubs, and thus infer their relationships with human disturbances over a period of 60 years (1956-2016).

Study Case
We used a reliable and reproducible methodology to monitor structural changes in scattered vegetation of a dryland coastal zone using very high spatial resolution images and OBIA. The temporal dynamics of the Ziziphus lotus (L.) Lam population in a semi-arid coastal zone was evaluated to infer the effects of human disturbances on the shape and size of individuals over a period of 60 years. Two human disturbances were evaluated: (i) the extraction of coastal sands in the 1970s [30], which eliminated the aeolian sands found in the upper layer of the soil using heavy equipment and created roads and dirt tracks in the area, and (ii) the seawater intrusion in the mid-1980s caused by groundwater withdrawals for greenhouse irrigation. The withdrawals resulted in the water table of the main aquifer dropping by 30 m [35]. In addition, we evaluated the impacts of the protection of the area in 1987 in the shrub species.
The study area is located on a coastal aeolian plain in the Cabo de Gata-Níjar Natural Park, Spain (36 • 49 46.3"N, 2 • 17 37.1"W; Figure 1). This area is one of the driest in Europe, with a mean annual precipitation, temperature, and potential evapotranspiration (PET) of 200 mm, 19°C, and 1390 mm, respectively [39]. The area is a hydrogeological complex located between two wadis with numerous fractures [40,41] and shows the typical landscape of arid areas with bare soil and patches of shrubs dominated by Z. lotus shrubs [29,37]. This vegetation is supported by a shallow, unconfined coastal aquifer composed of gravel and sand deposits located in the discharge zone at the end of the two wadis. A major fault hydrologically separates this aquifer from the main regional one (Hornillo-Cabo de Gata, [32]). Consequently, inflows to the aquifer come mainly from scarce local rainfall. The scattered shrubs of Ziziphus lotus in SE Spain form the largest population of this shrub species in Europe. This population is protected by the Habitat Directive (5220* habitat, 92/43/CEE) and the Water Framework Directive (WFD) in Europe [31]. In this area, Z. lotus forms intricate structures of 1-3 m tall, accumulating sand under its cover called nebkas. This forms favorable microclimatic conditions under its cover with respect to the outside and increases the water availability due to hydraulic lift [38], carbon exchanges, and energy cycles [29], creating islands of fertility [42], which increases the diversity of animal and plant species. For this reason, Z. lotus is considered an ecosystem engineering species in this environment [31,33].

Datasets and Ground Truth
Eight orthoimages from two sources were used, namely, six orthoimages from the Andalusian Environmental Information Network (REDIAM) with a spatial resolution of 1 m/pixel from 1956 and 1977 (panchromatic images) and 0.5 m/pixel from 1984, 1997, 2004, and 2008 (multispectral images), and two Google Earth orthoimages from 2013 and 2016, with a resolution of 0.5 m/pixel (multispectral images). To work with the same spatial and spectral resolution, we homogenized the images to the lowest spatial resolution (i.e., 1 m/pixel) and transformed them into panchromatic images (1 band) from the three spectral band images with QGIS software v. 3.8 (manufacturer, city, state abbreviation, country). For sand mining mapping, we used airborne LiDAR data with 1 m point spacing obtained in 2011. A summary of the dataset is shown in Table 1. Table 1. Data sources, spatial resolution, band numbers and year of the data used in the object-based shrub mapping and sand extraction estimation from airborne laser scanning (light detection and ranging (LiDAR)). The scattered shrubs of Ziziphus lotus in SE Spain form the largest population of this shrub species in Europe. This population is protected by the Habitat Directive (5220* habitat, 92/43/CEE) and the Water Framework Directive (WFD) in Europe [31]. In this area, Z. lotus forms intricate structures of 1-3 m tall, accumulating sand under its cover called nebkas. This forms favorable microclimatic conditions under its cover with respect to the outside and increases the water availability due to hydraulic lift [38], carbon exchanges, and energy cycles [29], creating islands of fertility [42], which increases the diversity of animal and plant species. For this reason, Z. lotus is considered an ecosystem engineering species in this environment [31,33].

Datasets and Ground Truth
Eight orthoimages from two sources were used, namely, six orthoimages from the Andalusian Environmental Information Network (REDIAM) with a spatial resolution of 1 m/pixel from 1956 and 1977 (panchromatic images) and 0.5 m/pixel from 1984, 1997, 2004, and 2008 (multispectral images), and two Google Earth orthoimages from 2013 and 2016, with a resolution of 0.5 m/pixel (multispectral images). To work with the same spatial and spectral resolution, we homogenized the images to the lowest spatial resolution (i.e., 1 m/pixel) and transformed them into panchromatic images (1 band) from the three spectral band images with QGIS software v. 3.8 (manufacturer, city, state abbreviation, country). For sand mining mapping, we used airborne LiDAR data with 1 m point spacing obtained in 2011. A summary of the dataset is shown in Table 1. Table 1. Data sources, spatial resolution, band numbers and year of the data used in the object-based shrub mapping and sand extraction estimation from airborne laser scanning (light detection and ranging (LiDAR)).

Data Source Spatial Resolution Band Number Year
Andalusian Two hundred perimeters of Z. lotus and 200 points of bare soil with scarce vegetation were randomly taken as the ground truth. A submeter precision GPS (Leica GS20 Professional Data Mapper; Leica, Wetzlar, Germany) was used. To do this, 12 longitudinal transects along the coast with a separation of 150 m between them were followed. The perimeter was taken with a distance of 1 m between nodes and the bare soil points were taken with a separation of at least 2 m from the nearest shrub. In addition, 200 shrub perimeters were digitized in each historical image with a distance of 1 m between nodes coinciding with the pixel size of the orthoimages in QGIS software v. 3.8.

Object-Based Image Analysis
OBIA consists of two phases, namely, the segmentation of the image into almost homogeneous objects and its subsequent classification based on similarities of shape, spectral information, and contextual information [17]. In the segmentation phase, it is necessary to establish an appropriate scale level depending on the size of the object studied in the image [43]; for example, low values for small vegetation and high values for large constructions [44,45]. During the classification, the segmented objects are classified to obtain cartographies of the classes of interest using algorithms such as nearest neighbor [46]. The success of the classification depends on the accuracy of the previous segmentation [47].

Image Segmentation
To segment the images, we used the multiresolution segmentation algorithm implemented in eCognition v. 8.9 software (Definiens, Munich, Germany). This algorithm depends on three parameters: (i) Scale, which controls the amount of spatial variation within objects and therefore their output size; (ii) Shape, which considers the form and color of objects; if it is set to high values, the form will be considered and if it is close to 0, the color will be considered instead; (iii) Compactness, a weighting to represent the smoothness of objects formed during the segmentation; if it is set to high values, the compactness will be considered complex and if it is set to values close to 0, the smoothness will be considered as simple [48]. To obtain the optimal value for each segmentation parameter, we used a ruleset in eCognition v8.9 that segmented the image by systematically increasing the Scale parameter in steps of 5 and the Shape and Compactness parameters in steps of 0.1 [49]. The Scale ranged from 5 to 50, and the Shape and the Compactness ranged from 0.1 to 0.9. A total of 6480 shapefiles were generated with possible segmentations of Z. lotus shrubs in a computer with a Core i7-4790K, 4 GHz and 32G of RAM memory (Intel, Santa Clara, CA, USA).
To evaluate the accuracy of all segmentations, we developed an R script to calculate the Euclidean Distance v.2 (ED2; [50]; Equation (1)), measuring the arithmetic and geometric discrepancies between the 200 reference polygons of Z. lotus and the corresponding segmented objects: (1) ED2 optimizes geometric and the arithmetic discrepancies with the "Potential Segmentation Error" (PSE; Equation (2)) and the "Number-of-Segments Ratio" (NSR; Equation (3)), respectively. According to [50], values of ED2 close to 0 indicate good arithmetic and geometric coincidence, whereas high values indicate a mismatch between them: where r k is the area of the reference polygon and s i is the overestimated area of the segment obtained during the segmentation; where NSR is the arithmetic discrepancy between the polygons of the resulting segmentation and the reference polygons and abs is the absolute value of the difference of the number of reference polygons, m, and the number of segments obtained, v.

Classification and Validation of Segments
We used the k-nearest neighbors algorithm to classify the best segmentations (lowest ED2 values) in two classes, that is, (i) Ziziphus lotus shrub (Z) and (ii) Bare soil with sparse vegetation patches (S). In order to train the classification algorithm, 70% of the ground-truth samples (140 Z and 140 S) and the features of greatest separability (J) between them, obtained using the separability and threshold (SEaTH) algorithm, were used [51,52]. The remaining 30% of the ground-truth samples (60 Z and 60 S) were used as the validation set [53,54], and the accuracy of the classifications was evaluated using error and confusion matrices, extracting Cohen's kappa index of accuracy (KIA) [55] and the overall accuracy (OA) of them. Finally, errors in shrub segmentation were evaluated by estimation of the root-mean-square Error (RMSE) and the mean bias error (MBE) between reference polygons and segments classified as Z. lotus shrubs.

Sand Extraction Curvature Analysis
The evaluation of areas affected by sand extraction within the study area in the 1970s was performed using a geomorphometric analysis of the land surface [56]. The analyses were based on a LiDAR-derived digital elevation model (DEM) dataset, generated using an ArcGIS toolbox for multiscale DEM geomorphometric analysis. This toolbox allows the generation of a number of curvature-related land surface variables [57], including plane, profile, mean, minimum profile, maximum profile, tangential, non-sphericity, and total Gaussian curvature; positive and negative openness; and signed average relief. Several maps were derived for each curvature variable at different spatial scales. The sizes of the analysis window ranged from 3 × 3 m to 101 × 101 m with a 14 m interval. Univariate and bivariate statistics were calculated for variables related to curvature [56,58].
A window size of 61 m was selected for the geomorphometric analysis of the curvature of the surface, which provided a good compromise between the size of land surface depressions resulting from sand mining operations and spatial generalization. The sand extraction areas were located and digitized on a final map and validated with a field survey. In addition, an estimation of the volume of soil loss resulting from sand extractions was performed. To calculate the volume of soil loss, a new digital surface model was generated without the soil loss zones extracted with the previous curvature analysis. Then, the difference was applied to the initial surface model with the areas identified as soil loss and to the digital surface model without soil loss, obtaining the volume of the previously identified soil loss areas.

Shrub Area and Shape Dynamics
Variations in the size and number of shrubs were determined by calculating the number of shrubs lost and differences in shrub cover area between consecutive image pairs. To calculate losses and gains in the coverage of the individuals, we assumed that a resulting negative area meant a loss of surface coverage, whereas a positive area meant a gain of coverage. To determine the edge effect and the health indicator on shrubs [59], the round shape index was calculated as the ratio between the cover area and the perimeter of each shrub in different years [60]. In order to evaluate whether the shrub cover reduction that occurred between 1956 and 1977 was related to sand extraction, the average minimum distance (AMD) and average random distance (ARD) between shrubs and sand extraction areas were calculated using PASSaGE v.2 software (The Biodesign Institute, Arizona State University, Tempe, AZ, USA) [61] with 999 permutations. We assumed that the shrubs in the 1956 image rather than the 1977 image were removed during sand mining, and those that reduced their cover were affected during this process. To evaluate whether the shrub population was affected by seawater intrusion between 1977 and 1984, the AMD and the ARD between the shrubs and the coastline were calculated as previously. Shrubs affected by sand extractions (those appearing in the 1956 image but not in the 1977 one) were not included in this analysis. When calculating the AMD, the shrubs that showed a reduction in cover over the corresponding period were used, whereas for the calculation of the ARD simulated shrubs were used. In order to evaluate the effects of protecting the shrubs within the Natural Park in 1987, reduction in shrub cover and number of shrubs in the 1984-2016 period was determined.

Segmentation Accuracy
The average values of the Scale and the ED2 were 25 and 0.45, respectively ( Table 2)

Classification and Characteristics of Ziziphus lotus Shrubs
The analysis of class separability and threshold with the SEaTH algorithm showed that the best features for discriminating between classes (i.e., those with the highest separability) were mainly related to texture (i.e., the family of features related to the Gray-Level Co-Occurrence Matrix (GLCM)) and brightness of objects (Table 3). Table 3. Features used in the classifications and separability between them using the separability and threshold (SEaTH) algorithm. In bold, the two features for each year with the highest separability used to classify the images. All the classifications were highly accurate, with values of OA and KIA close to 1 ( Table 4)

Shrub Number, Area, and Shape Dynamics
During the 60-year period evaluated, the number of shrubs decreased by 742. The moment of highest shrub population was 1977, with 2625 shrubs. Conversely, the lowest number of shrubs was detected in 2016, with 1883 shrubs (Table 5). However, the total shrub area between 1956 and 2016 increased by 3692 m 2 . In addition, we observed an increase in the maximum cover area value of shrubs after 1997. Finally, the most circular shrubs appeared in 1956 (i.e., the lowest values of the round shape index) and the high values of the round shape index increased over the years (Table 5). In general, the cover area of shrubs between pairs of years showed an increase, with a trend of smaller individuals to lose more cover area than larger shrubs. In the period 1984-1997, 1423 shrubs reduced their cover area. In the period 1977-1984 (Table 6), 1650 shrubs increased their cover area. Table 6. Change of cover and frequency of the difference in Ziziphus lotus area in the studied years . The negative and positive areas are the result of the subtraction between the year and its predecessor. The highest lost area, the largest positive area, the balance between greater areas, the positive frequency and the negative frequency of shrubs are highlighted in bold type.

Sand Extraction Mapping and Curvature Analysis
The results of the analyses indicated that more than 187 m 3 of sand were extracted in the study area (4.2 km 2 ). According to [30], a visual analysis of resulting maps also suggested that sand extractions were distributed spatially following a connected network and following existing roads in the area (Figure 2).

Spatial Relationships of Shrubs with Sand Extractions, Coastline (Seawater Intrusion), and Protected Area
Between 1956-1977, 752 shrubs reduced their cover area in the sand mining event. The AMD between the shrubs and the zones of sand extractions presented an average minimum distance of 25.57 ± 37.49 m. The ARD analysis showed an average minimum distance of 127.48 m ± 23.68 m between the random simulated shrubs and the zones of the sand extractions ( Figure 2).   In the period 1984-2016, 551 individuals were lost (Figure 4), but in the area there was a total gain of more than 23 m 2 (Table 5), coinciding with the protection of the study zone under the Cabo de Gata-Níjar Natural Park. In the period 1984-2016, 551 individuals were lost (Figure 4), but in the area there was a total gain of more than 23 m 2 (Table 5), coinciding with the protection of the study zone under the Cabo de Gata-Níjar Natural Park. In the period 1984-2016, 551 individuals were lost (Figure 4), but in the area there was a total gain of more than 23 m 2 (Table 5), coinciding with the protection of the study zone under the Cabo de Gata-Níjar Natural Park.

Discussion
The first step in inferring human disturbances in vegetation was to generate an accurate object-based map of the study area. The high accuracies obtained in the segmentations, evaluated with the ED2 index, facilitated the classifications of the shrubs, obtaining similar accuracies to

Discussion
The first step in inferring human disturbances in vegetation was to generate an accurate object-based map of the study area. The high accuracies obtained in the segmentations, evaluated with the ED2 index, facilitated the classifications of the shrubs, obtaining similar accuracies to previous studies in which other species of shrubs with OBIA were detected [21,62]. In the classification step, the segments of Z. lotus and bare soil with sparse vegetation showed high separability, using features related to their brightness and texture as a consequence of clear spectral differences between vegetation and bare soil [63]. According to [64] and [52], the best spectral-related features for discriminating between vegetation and bare soil with sparse vegetation were Brightness, and the GLCM family of features, which present high separability values in well-differentiated classes, such as vegetation and soil [52,65,66]. The worst features to discriminate (i.e., lowest separability) were the geometry-related ones (e.g., Roundness, Area). This is contrary to previous studies, in which Roundness has been suggested as a potential feature to discriminate between rounded shrubs and bare soil [21,66,67]. This discrepancy with previous studies could be explained by the high heterogeneity in the vegetation form that can present after disturbance [5]. However, care must be taken interpreting these results, since shadows generated by large individuals may result in an overestimation of shrub cover area [68,69], and individuals growing together may be underestimated as a result of appearing in the images as one [62].
The spatial distribution of sand extraction areas was unrelated to the topography of the area but was related to the location of older roads and tracks. This suggests that sand extractions were preferentially located to previous sand extractions in an effort to minimize labor costs. This observation is consistent with a previous study on sand extraction in the region (e.g., [30]). According to [70], the negative effect on shrubs by sand mining was shown in the low values of AMD calculated in the period 1956-1977. This reduction in population cover area could be related to sand mining during the 1960s and 1970s [30], which might confirm the positive effect that sands have on the health of Z. lotus shrubs in the study area [31] and in other areas of North Africa [71,72].
The lower value of the AMD between shrubs with reduced cover in the period 1977-1984 to the coastline suggested by the ARD that this reduction could be related to a decrease in the freshwater table and the intrusion of seawater into the aquifer [35,73]. The smallest shrubs were the most affected, which can be related to difficulties of access to groundwater due to their smaller roots compared to larger and more developed individuals [28,74]. These results agree with previous studies evaluating the negative effect that seawater intrusion has on vegetated communities and groundwater-dependent ecosystems (e.g., [75,76]).
In addition, the results of this work could be affected by other natural conditions or affections, for example, shrubs could be affected by herbivory [5,77], climate change [78], or uncontrolled use of pesticides [79]. We argue that in order to better understand the results obtained in this work, it is necessary to complement remote sensing techniques with in situ works. For example, complementing the results obtained with the presence-absence of isotopes and relating them with the seawater intrusion [80,81] could provide a better understanding of how this phreatophytic community responds to anthropic perturbations over time.
Although 742 Z. lotus individuals were lost during the study period, their average size and the round shape of the shrubs were higher and bigger at the end than at the beginning of the period. However, the variability of these characteristics also increased over time, which means that a greater variety of shapes and sizes was observed in the population. This could be explained by the 1987 declaration regarding the Cabo de Gata-Níjar Natural Park, where the study area is located. This protection, in addition to a slow recovery of the aquifer after undergoing seawater intrusion between 1977 and 1984, could have contributed to a slow but continued development in time by adults, which might have better access to fresh water from the aquifer due to a more developed root system [82] up to a length of 60 m [83]. Furthermore, the fact that the largest shrubs were the most developed in time supports the longevity character of this species through longevity [27], which is an important strategy for its survival in the Mediterranean region [26]. This, together with the anthropic pressure on the system, may explain the development of adult individuals, but the lack of recruitment of juveniles, as observed by [27] not only in this area but also in other regions in SE Spain.

Conclusions
The combination of very high-resolution historical images and OBIA is a powerful tool for identifying and monitoring communities of sparse vegetation in drylands [62]. Our results suggest that monitoring changes in the number and the cover of a shrub community in a semi-arid ecosystem could help to infer anthropogenic disturbances that affect its health. The vegetation conditions showed that the loss of sandy substrate affected Z. lotus negatively, either by reducing its cover or by eliminating individuals by direct sand extraction processes. In addition, seawater intrusion into the aquifer influenced the cover and structure of the shrubs close to the coastline in a period of massive groundwater extraction [35], negatively affecting the smallest shrubs for the most part. However, the legal protection of the area had a positive effect on the health of the remaining individuals, which increased their coverage. The implementation of semi-automatic methods to infer the effects of human activities on shrub populations, such as the one evaluated in this study, could help improve the monitoring programs of existing protected areas. This could reduce the cost of these activities, not only in economic terms but also from a human perspective, which is key to the long-term preservation of any protected area.