*2.2. Methods*

Several methodologies for change detection and classification can be applied, and selecting the appropriate technique is related to the objective of the study [46]. Techniques such as single image differing or ratioing provide binary information (change/no-change), and if detailed information is required, such as the change direction, classification techniques are preferable. Another technique providing information on change type and direction is index differencing. These are mathematical transformations in multispectral mode and are produced separately; then, other change detection techniques (e.g., differencing or ratioing) can be applied [46]. Concerning the unit of analysis, on the one hand, in [46], it is explained that pixel-based change detection methods have been used traditionally, the main advantages of this unit of analysis being its suitability for large pixels; it does not generalize the data; and it is an effective methodology, especially when the relationship between pixel intensity and the land cover changes under investigation is strong [14]. On the other hand, the object-based approach allows the exploitation of the spatial context, reduces the noisy outputs of isolated changed pixels and allows direct object change detection (DOCD) by comparing spectral information [46]. One of the object-based units of analysis is the vector polygon, which is extracted from existing geodatabases; they group together pixels that are suitable for statistical analysis, the result of which may indicate changes within the corresponding polygons. On one hand, vector polygons provide a cartographically 'clean' basis for analysis [14], allow the exploitation of additional thematic information about the objects to obtain better results and enhance the interpretation of the image [47]. They also provide important information on the location of the objects to investigate for change detection. In [47], this type of object-based approach, combined with spectral indices, was used for the automatic change detection of buildings in an urban environment as it can handle the complexity of urban environments. On the other hand, vector polygons generalize the data, and the size and shape of objects cannot be compared [14]. As regards the current study, we opted for this latter methodology, where a set of features (multi-spectral indices and radar backscattering) are used to create what we could define as the temporal signatures of the RDSs. The methodology responds to the need to monitor the RDS polygons at a regional scale, and to have generalized information of changes detected for the three types of classes (vegetation, building and soil), thereby reducing the manual work [14].

The proposed overall methodology, whose main goal is to provide a shortlist of the sites that are likely to have changed and for which an on-field visit would be required, is shown in Figure 4. The first main block is the feature extraction, where the Sentinel-1 and Sentinel-2 images available in Terrascope and described in the previous section are processed to obtain the above-mentioned temporal profiles of the RDSs (each RDS has multiple temporal profiles—one per feature). The second main block is devoted to the characterization of the changes, which is carried out in two steps: (i) the change detection, which flags a site as changed (or not) and provides an estimate of the change date(s)—this is carried out once every two months; (ii) the change classification, which is divided into two separated processes. First, when a change date is detected, a rule-based classification is performed in order to provide additional information on the type of change: vegetation, building or soil. Second, the same methodology is applied once per year, considering summer average features in order to detect gradual changes.

**Figure 4.** Workflow of the automatic change detection and classification of the RDSs.

The final output of the service is a csv file that is automatically delivered to the operator. For each RDS, this report includes: (1) information on whether a change has occurred or not, (2) the type of change and (3) the estimated date of the change (if available).

### 2.2.1. Features Extraction and Temporal Profiles

For each Sentinel-1 acquisition (more specifically, the VH band, which was found to be the most suitable for our scope) that contained the site of interest within the desired time frame, the average backscatter (sigma0) for that site was computed and used to populate the corresponding temporal profile. Since a site can be typically seen from 3 to 4 different viewing angles (considering both ascending and descending orbits), separate profiles were created for each satellite pass and then averaged to obtain a unique "sigma0VH" feature.

Regarding Sentinel-2 data, all the L2A tiles over the area were analyzed. Only the tiles presenting less than 25% of clouds were selected, which greatly reduced the number of undetected cloud pixels. Then, each image was clipped based on the RDS vector polygons file. Image co-registration was ensured during this process. Then, the Scene Classification layer, a classification map generated via the Sen2Cor ESA processor that accompanies every L2A image and is directly available in Terrascope, was used to remove every single pixel classified as "No\_Data", "Cloud\_Shadows", "Cloud\_Medium\_Probability", "Cloud\_High\_Probability", "Thin\_Cirrus" and "Snow". This allowed us to remove, site per site, the dates for which no data, shadows, clouds, or snow pixels were present.

In the feature extraction step, six widely used spectral indices were calculated that were to be used in the next processes: (1) the Built-Up Areas Index (BAI) [47], (2) the Brightness Index (BI) [48], (3) the Second Brightness Index (BI2) [48], (4) the Normalized Vegetation Index (NDVI) [49], (5) the second Normalized Difference Water Index (NDWI2) [50] and (6) the soil brightness index (SBI) [47]. The selection of the spectral indices was motivated by their widespread application in the literature and by considering that most built-up indices require SWIR bands, which are available only in a coarse resolution for Sentinel-2. The BI2 index has been tested for built-up detection after applying NDVI and NDWI2 to mask vegetation and water [51]. BAI has proven to be useful to detect asphalt and concrete surfaces [47], and SBI has been successfully investigated by [47] and [52]. For each index, each RDS and each available image since 2015, the average per RDS was calculated and used to generate the Sentinel-2 time series:

$$\text{BAI} = (\text{(BO2} - \text{BOS)}) / ((\text{BO2} + \text{BOS})) \tag{1}$$

$$\text{BI} = \sqrt{((\text{B04\*B04}) + (\text{B03\*B03}))/2} \tag{2}$$

$$\text{BI2} = \sqrt{\left(\left(\text{(BO4\*B04)} + \left(\text{BOS\*B03}\right) + \left(\text{BOS\*B08}\right)\right)/3\right)}\tag{3}$$

$$\text{NDVI} = ((\text{BOS} - \text{BO4})) / ((\text{BOS} + \text{BO4})) \tag{4}$$

$$\text{NDWTL2} = \left( \left( \text{BO3} - \text{BO8} \right) \right) / \left( \left( \text{BO3} + \text{BO8} \right) \right) \tag{5}$$

$$\text{SBI} = \sqrt{((\text{BO4\*B04}) + (\text{BOS\*BOS}))}\tag{6}$$

where B0n corresponds to the n-th Sentinel-2 band used for the calculation, here B02, B03, B04 and B08, all with a 10 m resolution.

To create the final temporal profiles (each RDS has multiple profiles, one per feature), a linear interpolation to fill in the gaps (1 data point per day) in the data and a smoothing using a Gaussian kernel with a standard deviation of 61 were performed.
