*2.1. 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 e ffects 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◦4946.3"N, 2◦1737.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.

**Figure 1.** Upper images: the distribution of *Ziziphus lotus* priority habitat 5220\* in the Mediterranean area. Lower image: Dashed line shows the study area under urban and intensive agricultural pressure in 2016. UTM projection Zone 30N; WGS 1984 Datum. Map data: Google, Maxar.

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].

#### *2.2. 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.



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.

#### *2.3. 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].

## 2.3.1. Image Segmentation

To segmen<sup>t</sup> 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:

$$\text{ED2} = \sqrt{\text{(PSE)}^2 + \text{(NSR)}^2}.\tag{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:

$$\text{PSE} = \frac{\sum |\mathbf{s}\_i - r\_k|}{|r\_k|},\tag{2}$$

where rk is the area of the reference polygon and si is the overestimated area of the segmen<sup>t</sup> obtained during the segmentation;

$$\text{NSR} = \frac{\text{abs}(m-v)}{m},\tag{3}$$

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.*

#### 2.3.2. 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.

#### *2.4. 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.

#### *2.5. Shrub Area and Shape Dynamics*

Variations in the size and number of shrubs were determined by calculating the number of shrubs lost and di fferences 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 e ffect 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 di fferent 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 a ffected during this process. To evaluate whether the shrub population was a ffected by seawater intrusion between 1977 and 1984, the AMD and the ARD between the shrubs and the coastline were calculated as previously. Shrubs a ffected 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 e ffects 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.
