THE POTENTIAL OF SENTINEL-1 DATA TO SUPPLEMENT HIGH RESOLUTION EARTH OBSERVATION DATA FOR MONITORING GREEN AREAS IN CITIES

Green areas play an important role within urban agglomerations due to their impact on local climate and their recreation function. For detailed monitoring, frameworks like the flora fauna habitat (FFH) classification scheme of the European Union’s Habitat Directive are broadly used. By date, FFH classifications are mostly expert-based. Within this study, a data-driven approach for FFH classification is tested. For two test areas in the municipality of Vienna, ALS point cloud data are used to derive predictor variables like terrain features, vegetation structure and potential insulation as well as reflection properties from full waveform analysis on a 1 m grid. In addition, Sentinel-1 C-Band time series data are used to increase the temporal resolution of the predicting features and to add phenological characteristics. For two 1.3 × 1.3 km test tiles, random forest classifiers are trained using different combinations (ALS, SAR, ALS+SAR) of input features. For all model test runs, the combination of ALS and SAR input features lead to best prediction accuracies when applied on test data.


INTRODUCTION
Green areas within urban agglomeration provide space for recreation and have considerable impact on local climate (Aram, et al., 2019, Tyrväinen, et al., 2005. Monitoring of green areas and biotope type mapping in the municipality of Vienna is required by nature conservation laws and national and international agreements on monitoring and reporting. The habitat types listed in Annex I of the Habitats Directive (European Union's Council Directive 92/43/EEC on the conservation of natural habitats of wild fauna and flora, classification, see European Commission, 2006) constitute a particularly relevant and prominent classification scheme. The scheme described therein comprises nine flora fauna habitat (FFH) groups and several sub-groups. To date, classification and mapping is mainly conducted by expert-based identification of general habitat types and plant species. It is labour-intensive, and therefore only conducted in a coarse temporal resolution. Datadriven approaches represent a less labour-intensive alternative that in addition may provide the possibility for a more detailed temporal monitoring. To implement such an approach, it is necessary to find features, which reflect the different FFH types. Morphological characteristics are known to be good predictors for different plant species (Hollaus, et al., 2009, Koenig and Höfle 2016, Puliti, et al., 2017. Additionally, phenological attributes and the temporal and spatial dynamics of the morphological parameters have been shown to be required to achieve a comprehensive classification (Martínková, et al., 2002, Dostálová, et al., 2018. The municipality of Vienna provides an extensive geodatabase and uses state-of-the-art methods for collecting geodata by surveying, airborne imaging and airborne laser scanning (ALS). Due to their limitations in terms of spatial resolution, satellite data have not been included in monitoring projects of the municipality of Vienna so far. * Corresponding author High-resolution earth observation (EO) data, especially point cloud data from airborne campaigns (e.g., ALS) provide valuable information for deriving structural parameters to describe morphological characteristics of vegetation (Hollaus, et al., 2006, Wagner, et al., 2008, Eysn, et al., 2012 and underlying terrain. Due to their low temporal resolution (for Vienna, there are currently only two comprehensive ALS data sets available -2007, 2015), the dynamics of structural parameters and the phenological characteristics cannot be extracted. The present study explores the possibility to fill this gap with satellite-based data. The two-satellite Sentinel-1 (S-1) constellation, operated as part of the European Copernicus programme, offers a combination of high spatial resolution in interferometric wide swath (IW) acquisition mode combined with a high temporal coverage of up to four acquisitions within a 12-days interval (two satellites and along ascending and descending orbits). As S-1 operates in the microwave portion of the electromagnetic spectrum, the data acquisition is not hampered by cloud cover. The current study aims to gain insight of the potential of EO data for FFH type mapping. Therefore, the replication of the FFH type-based classification of chosen study sites of the most recent biotope types mapping of the municipality of Vienna (Stadt Wien https://data.wien.gv.at, 2020) using data-driven approaches is targeted. For this purpose, different features are derived from both ALS-and synthetic aperture radar (SAR) data. This data is then analyzed in an explorative manner with reference to already existing FFH type-based classification of the selected areas of interest. In the next step, random forest (RF) models are trained using different combinations of predictors, including ALS-based predictors only, SAR-based predictors only and a combination of both. The accuracy of this fully-automated, supervised classification, and especially the added explanatory power of the S-1 data is investigated. Based on the derived results, the potential usage of S-1 data for future biotope mapping campaigns is discussed.

Habitats Directive classification
For the classification within this study, the targeted groups are the FFH types according Annex I of the Habitat Directive. The classification scheme differentiates a total of 229 FFH types: nine major groups with up to 31 subgroups per group (European Commission, 2006). Of the nine different major groups, five exist within the selected study areas in the municipality of Vienna. The research area comprises a total of 22 different FFH types. The mapping of the FFH classification for green areas within the municipality of Vienna is available as open data (Stadt Wienhttps://data.wien.gv.at, 2020). The most recent mapping was conducted in 2008. A comprehensive representation of the occurring types is shown in Table 1.

Study Sites
For this preliminary study, two major green areas within the municipality of Vienna are chosen for further investigation. The two test sites cover two main green area types in Vienna: (A) hilly, mainly forested areas in western Vienna and (B) riparian areas with river meadows and riparian forest vegetation along the Danube River. The five major groups in the areas are freshwater habitats (3), natural and semi-natural grassland formations (6), raised bogs and mires and fens (7), rocky habitats and caves (8) and forests (9). The first area (study site A) is located in the southwest of Vienna in the Vienna Woods. It covers 28.3 km². The area is hilly and cut by three major valleys in the south, north and northeast. The altitude is between 214 m and 515 m above the Adriatic (AA) with a mean value of 337 m AA. Study site A is dominated by forests (77.6% of the area). The major occurring forest types are Galio-Carpinetum oak-hornbeam forests (48.3%) and Asperulo-Fagetum beech forests (24.1%). 9.3% of the area are covered by natural and semi-natural grassland formations with lowland hay meadows (5.9%) semi-natural dry grasslands and scrubland facies on calcareous substrates (2.6%) as dominant types. 13.9% of the area is not assignable to any FFH type. The second area of interest (study site B) is situated in the southeast of Vienna in the riparian forests along the Danube River within the Donau-Auen National Park. The area covers 22.9 km². The area is mainly flat with indications of a former dominant and now partly regulated braided river system and a . 9% of the area is covered by semi-natural grassland formations with semi-natural dry grasslands and scrubland facies on calcareous substrates (4.6%) and lowland hay meadows (3.5%) as the most common types. 6.8% of the area is classified as freshwater habitat, with 6.5% natural eutrophic lakes with Magnopotamion or Hydrocharition -type vegetation.
The shares of the different FFH types for both study sites are also shown in Figure 1.

ALS data
The ALS data used was collected between November 9 th and November 24 th 2015 under leaf-off conditions in eight flight campaigns and covers the entire area of the municipality of Vienna. For the data acquisition, two different sensors were used: a Riegl LMS-Q680i and a Riegl LMS-Q560 (RIEGL Laser Measurement Systems, Horn, Austria). Both scanners provide full waveform analysis (Riegl, 2010, Riegl 2012). The acquisition resulted in a point density of > 16 echoes/m² for 97% of the whole city area.

Sentinel-1 data
The SAR features used are derived from data from S-1 in C band at a pixel spacing of 10 m. In detail, the Level 1b high resolution ground range detected (GRDH) and Level 1a single look complex (SLC) scenes acquired in the Interferometric Wide (IW) swath mode were used. The data from the relative orbits 22 (descending) and 146 (ascending) was used separately.

ALS features
All ALS features were derived directly or indirectly from the ALS point cloud described in 2.3 to 1 m resolution raster files. Terrain parameters such as aspect, slope, topographic positioning index (TPI) and topographic wetness index (TWI) describe the absolute as well as the relative position of a pixel in the surrounding terrain. The parameters were calculated using a pre-derived digital terrain model (DTM) from the ALS point cloud and gdaldem (Perry, et al., 2021) and SAGA Wetness Index (Conrad and Böhner, 2001). Aspect and slope give information about exposition and the inclination of the terrain. The TPI measures the difference of a central pixel and the elevation of a defined neighborhood. Positive TPI values imply a pixel location higher than the neighborhood. The pixel is therefore situated on a local peak or hilltop, or upslope. Negative values indicate a pixel location lower in comparison to its neighborhood. It can be seen as a sink or hill bottom (De Reu, et al., 2013). The TWI combines the upstream contribution area and the local slope.
High TWI values indicate a high potential runoff (Sørensen, et al., 2006). The insulation parameters give information about potential duration and intensity of insulation and incoming solar radiation. In a first step, a DSM was derived from the ALS point cloud using OPALS (Pfeifer, et al., 2014), and in a following step, the insolation features were derived using SAGA Potential Incoming Solar Radiation (Conrad, 2010). In combination, the features described above represent location characteristics for different habitat types.
To get information about the vertical and horizontal structure of the vegetation, structure features were derived. To get the absolute vegetation height, the normalized DSM (nDSM) was calculated by subtraction of the DTM from the DSM. For a detailed vertical structure description voxel-based analysis was applied. The voxel size was set to 1x1x1 m³. The understory height was set to the height of the voxel column above ground level, in which each voxel contains continuously at least one point. For the horizontal structure, a stand density index was calculated. Therefore, all pixels with a nDSM height < 0.2 m were defined as gaps, and the proportion of gaps within a 10 m radius neighborhood was calculated. The features described above are all based on point coordinates and the relative location. In addition to the location, the sensors used for the ALS campaign record the complete waveform of the backscattered echo, i.e., full-waveform lidar. This enables range measurements and the acquisition of additional information on the physical properties of the recorded objects (Wagner, et al., 2008). In order to obtain surface reflectance, a ratio between the amplitude of the ground echo and the surface echo was calculated.
To extract roughness information for the terrain and the surface,  The structure features as well as the full waveform features were calculated from the point cloud using OPALS (Pfeifer, et al., 2014). An overview of the metrics derived from ALS is shown in Table 2.

SAR features
From the SAR data, input features sensitive to forest canopy structure and phenological characteristics were derived. Therefore, seasonally averaged backscatter with VV (vertically polarized sent and vertically polarized received), VH (vertically polarized sent and horizontally polarized received) polarization and VH/VV cross ratio were calculated. Interferometric coherence was estimated between subsequent acquisitions with a temporal baseline of six days. Higher contributions of VH backscatter are typically seen as an indicator for volume scattering in vegetation canopies. The products were obtained separately for the relative orbits 22 (descending) and 146 (ascending). The sensitivity of the derived features to vegetation monitoring was already shown in several studies (e.g. Vreugdenhil, et al., 2018, Bruggisser, et al., 2021. Detailed information on the data processing is shown in Dostálová, et al. (2018). The data was harmonized by reprojecting and interpolated to a 1 m grid using gdalwarp (Warmerdam, et al., 2021) with bilinear resampling. Table 3 shows a detailed overview of the derived SAR features.

Exploratory statistics
For a first assessment of the usability of different predictor variables for supervised classification, the distribution of the values of the different SAR and ALS products was analyzed. Therefore, boxplots grouped by the present FFH classification were created for each of the two research areas separately. The potential of selected features as predictions for a FFH type-based classification is discussed.

RF training and prediction
Using multiple predictors, the supervised reproduction of the FFH classification was realized with a RF-model (Breiman, 2001). The model is trained and tested using the 'ranger' -Package (Wright, et al., 2020) implemented in R. The calculated features were imported, cropped and combined to raster stacks. For a first test, two tiles of 1.3 km × 1.3 km (one within each study site) were defined. The tiles were further spatially split in 75% training and 25% test area. Cases including missing values in at least one of the features were excluded.
The performance of the model is controlled by several parameters, including the number of trees (num.trees), the number of variables used for splitting at each node (mtry) and the maximal depth (max.depth) of the trees. In order to tune the model, the number of trees [250, 500] and the maximum depth [5,10,15,20] were varied. For each of the settings, three RF models were trained using different sets of predictors: the ALS features, the SAR features and both the ALS and SAR features. The trained RF were subsequently applied on the test data sets to predict the FFH types. The total prediction accuracy for each model was extracted by comparing the predicted values and the actual values for each pixel of the test data sets. For the chosen models, contingency tables were created to gain a more detailed overview of prediction accuracies on different classes. For a visual interpretation, the prediction and reference data was retransformed to a raster stack, exported as a GeoTIFF and plotted.

Exploratory statistics
For the FFH classification two research areas were investigated using descriptive statistics. The extent of the occurring FFH types varies widely among and within the research areas ( Figure 1). Figure 2 shows chosen results of the exploratory statistical analysis of the derived ALS and SAR products. The distribution of the feature values on all occurring FFH types is plotted, but further, only the dominant types covering more than 1% of the area are considered (Figure 1). Henceforth, the Natura 2000 Codes are used for addressing the different FFH types. The corresponding names can be looked up in Table 1. Within the study site A, the insulation differentiates between the groups of 6210, 6510 and 91F0 and the group of 9130 and 91E0. While the median potential annual solar insulation is between 43 kWh/m² and 47 kWh/m² for the first group, for the second group of areas, the potential energy input varies between 33 kWh/m² and 35 kWh/m². The stand density shows a perfect feature for distinguishing between forests and natural and seminatural grassland formations. While the median gap fraction (stand density) varies between 0.7 and 1.0 for grassland, the gap fraction of forest areas is close to zero.  In flat areas like the Donau-Auen National Park, the significance of the TWI to describe hydrological and topographic conditions is limited comparing to hilly areas, as the flow directions are rather undefined (Grabs, et al., 2009). The mean values of the forest classes are slightly lower than for the grasslands and the freshwater habitats. The clearest differences in SAR backscatter values are again found between forest and grassland areas, with forest areas showing higher values for the mean summer backscatter of VH orbit 146 as well as for VV orbit 022. The differences between the two forest types are negligible at the current scale of data analysis. The two grassland groups differ slightly, with 6120 showing higher reflectance values for both considered SAR features. The freshwater habitat 3150 shows values in between the two grassland classes for both SAR features.

Test tiles
For first tuning and testing of RF models, two tiles of the research areas were extracted. For an independent performance validation of the models, each tile was spatially split in 75% training area and 25% test area. The tiles were chosen by visual inspection of the research areas. The choice was based on heterogeneity of the mapped shapes in the tiles and a great variety of size and coverage of different FFH types. An overview of the chosen tiles is given in Figure 3.

RF-model on test tiles
After the definition of training and test areas, several test runs on RF tuning were performed using different parameter settings for the number of trees and the maximal depth of the trees. The models were trained with the input features cropped to the described training areas. With regard to the tiles, every single pixel was considered to be a discrete, unique entity. Thus, about 1.27 million pixels were used for training for each tile. The models were applied on 422 500 pixels of test data for each research area (the numbers might be slightly lower due to exclusion of missing values). The results of the model tests are summarized in Table . The prediction accuracy of models trained using a combination of both, ALS and SAR data was found to be slightly superior to  The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII- B3-2021XXIV ISPRS Congress (2021 those that are based on either ALS or SAR data alone. Considering the model parameters, for Tile A, 500 trees with a maximum depth of five lead to the largest amount of correctly predicted pixels when used on the test data. For Tile B, 500 trees with a maximum depth of 15 brought the highest correct classification output when used on the test data. These two results are chosen for more detailed investigations. The mapping of the model predictions is shown in Figure 3.
Plotting the data gives additional information on the spatial distribution of correctly and incorrectly classified pixels. The results of Tile A confirm the dominance of class 9130 and the misclassification of 9170. Basically all values of 9170, which are located in the northwestern corner of the test area, are classified as 9130, both with ALS and SAR predictors only as well as when a combination of both predictor sets is used. While the borders between the forest and the more fragmented grassland area in the southern part are visible, the different segments within the grassland were not properly classified using the presented approach. While the dominance of 6510 is still visible, the delineation between the different grassland types cannot be found   When taking a look at the training data, this can probably be traced back to an underrepresentation of this FFH type in the training data set. The homogeneity of the 91F0 area is mostly reproduced in the ALS data set, but in the SAR data set, patches of 91E0 are found in the 91F0 area, which is finally reflected in the output of the model using both ALS and SAR data. In general, using only ALS data for model training leads to speckled results. For SAR data only, the classification results show a rather patchier spatial distribution of the different FFH types. In combination, both of these characteristics are still visible, but smoothened.

REFLECTION AND OUTLOOK
A main challenge in this study is the temporal mismatch of the used data. While the biotope type mapping was performed in 2008, the ALS campaign was performed in 2015 and the used SAR data is from the year 2018. It is assumed that the major parts of the study sites did not change significantly within this timespan. For the two test tiles, 15 available orthophotos between 2010 and 2020 were examined and no major visually detectable changes in land cover were found.
The decision on which model to use for further investigation in this study was simply based on the prediction results. Although the amount of correctly classified pixels differs at a range of lower than 0.5 percentage points between 250 and 500 trees, the model using 500 trees was chosen. For future applications, especially on larger data sets, efficiency has to be considered as well, as the computation costs increase linearly with the number of trees when using RF-based algorithms (Scornet, 2018). Another challenge within this study occurs for the training and prediction of the type "no FFH type". There are no distinct classification characteristics known for this class and it is assumed to be very diverse. For future model training, the class will be excluded. The model performances values on the test data shown in this paper must be reflected critically. The quality of a RF model strongly depends on amount and selection of training data, both spatially and considering the number of predictor variables. The choice of test and training data set was not randomized within this study, and the choice of input predictors has not been optimized so far. Additionally, issues regarding the downsampling of the SAR data to a 1 m resolution are not considered here. What can be shown is, that combining ALS data with high spatial resolution and SAR data with a high temporal resolution can, at least for the two chosen test sites, increase the prediction accuracy of a RF classifier for FFH type classification. Further model tuning, subsampling of predictors and potential smoothing effects due to down-sampled SAR data and variations of the minimum mapping unit are subject of ongoing investigations.