Using Landsat Vegetation Indices to Estimate Impervious Surface Fractions for European Cities

Abstract

Impervious surfaces (IS) are a key indicator of environmental quality, and mapping of urban IS is important for a wide range of applications including hydrological modelling, water management, urban and environmental planning and urban climate studies. This paper addresses the accuracy and applicability of vegetation indices (VI), from Landsat imagery, to estimate IS fractions for European cities. The accuracy of three different measures of vegetation cover is examined for eight urban areas at different locations in Europe. The Normalized Difference Vegetation Index (NDVI) and Soil Adjusted Vegetation Index (SAVI) are converted to IS fractions using a regression modelling approach. Also, NDVI is used to estimate fractional vegetation cover (FR), and consequently IS fractions. All three indices provide fairly accurate estimates (MAEs ≈ 10%, MBE’s < 2%) of sub-pixel imperviousness, and are found to be applicable for cities with dissimilar climatic and vegetative conditions. The VI/IS relationship across cities is examined by quantifying the MAEs and MBEs between all combinations of models and urban areas. Also, regional regression models are developed by compiling data from multiple cities to examine the potential for developing and applying a single regression model to estimate IS fractions for numerous urban areas without reducing the accuracy considerably. Our findings indicate that the models can be applied broadly for multiple urban areas, and that the accuracy is reduced only marginally by applying the regional models. SAVI is identified as a superior index for the development of regional quantification models. The findings of this study highlight that IS fractions, and spatiotemporal changes herein, can be mapped by use of simple regression models based on VIs from remote sensors, and that the method presented enables simple, accurate and resource efficient quantification of IS.

1. Introduction

Urban land-use is typically characterized by high levels of temporal and spatial dynamics (e.g., due to rising populations) and a large degree of heterogeneity as roads, buildings or other paved areas, different kinds of vegetation, water and areas of bare soil are located within small distances. The detailed composition of the urban environment to a large extent defines the hydrological and thermal processes in cities, and hence mapping accurately the urban structure and development at the proper spatial and temporal scales required, e.g., for future urban planning, climate adaptation and cloudburst management is essential [1,2,3]. For many such applications, remote sensing techniques including satellite imagery provides superior temporal and spatial coverage facilitating systematic, accurate and resource efficient (semi- or fully- automated) mapping of urban land cover at various scales as compared to, e.g., traditional field surveys with GPS or manual digitizing of hard copy maps [2,4]. In general, significant changes in urban land cover are found at annual to decadal timescales [5]. While state-of-the-art high resolution (HR) (<5 m) satellite imagery only dates back to the late 1990s, medium resolution (MR) (5–100 m) and low resolution (LR) data are available for the past 30–40 years, allowing for extended time series of urban change and development monitoring at a low cost of spatial resolution.

Impervious surfaces (IS) defined as man-made sealed surfaces through which water cannot infiltrate typically dominate urban environments and may be used as a proxy for many aspects of urbanization [4,5]. IS are a key indicator of environmental quality as they have important implications for many bio-physical processes at the regional and local scale [6]. For urban areas, these processes are primarily linked to the hydrological cycle (run-off volumes and velocity, transport of non-point pollutants, groundwater recharge), the surface energy budget (reflective properties of surfaces, changing heat fluxes, urban heat islands) and biological functions (flora and fauna distribution, biodiversity) [7,8]. Since the 1970s, a number of different pixel and feature based methods have been developed for estimating the quantity and location of urban IS from satellite imagery, including the use of vegetation fractions, regression modelling, object/feature-based methods, Spectral Mixture Analysis (SMA), Maximum Likelihood Classifiers (MLC) and Artificial Neural Networks (ANN), etc. (a complete review may be found in [4]). It is generally agreed by the scientific community that resolving the finer scales of the urban heterogeneity poses a challenge for remote sensing and that only analyses based on high-resolution imagery provides the accuracy needed to, e.g., overcome the problem of “mixed pixels”. On the other hand, in terms of routine applications the high complexity, acquisition and analysis expenses, limited spatial and temporal extent, and high costs of sequential coverage often limit the use of such methods [9,10]. Likewise, the wide variety of urban remote sensing techniques, image sources, and effectiveness and how they are evaluated restricts our ability to compare different studies, making it difficult if not impossible to systematically assess, e.g., cities’ resilience to flooding at scales ranging from regional to even continental scales. In this paper, we explore whether a simple and transferable methodology may bridge this gap by facilitating the inter-comparison of estimated IS fractions across cities in a robust and systematic way from easily accessible Landsat imagery.

VIs like the Normalized Difference Vegetation Index (NDVI) [11,12] and the Soil Adjusted Vegetation Index (SAVI) [13] have a long record of merit within the remote sensing of vegetation and vegetation indices (VI)s have been demonstrated by several authors to provide relatively accurate estimates of the quantity and distribution of IS and urban land cover [14,15,16,17]. The use of VIs as a tool to estimate urban IS fractions is based on the assumption of a strong inverse relationship between vegetation cover and IS, i.e., it is implicitly assumed that non-IS within urban areas are covered with green vegetation. Limitations in using VIs for mapping IS include the influence of bare soils, shadow effects from buildings and trees and tree crowns covering IS [17]. Arguably, the most critical of the three, bare soil has similar spectral characteristics as urban fabrics and hence is easily confused with IS [4,18]. Also, the presence of bare soil causes misclassifications to occur when using automatic/semi-automatic classification techniques to create IS data from HR imagery for validations purposes. For urban areas, Bauer et al., 2008 has proposed using growing season images only to quantify IS fractions to reduce the presence of bare soil. In general, however, using VIs for mapping IS should be done with caution in cities with large areas of bare soil, e.g., as found in many developing countries.

The Landsat archive covers the past 30–40 years and is freely available, making it highly applicable and well suited for monitoring changes to IS and other aspects of urban development. Specifically, we have compared the applicability of Landsat 8 NDVI and SAVI, and the fractional vegetation cover (FR) approach in quantifying sub-pixel IS fractions for eight cities in Europe and their transferability to each other. We apply a well-known regression-based approach as recommended by [2] and [14], which is found to be suitable for deriving time series of urban IS fractions from VIs. We then investigate whether the models based on VI’s that are derived for a specific urban area may be reused at different geographical locations while asserting the loss of accuracy [15]. Lastly, we examine different regional models aggregating data from different cities.

2. Materials and Methods

Three different IS datasets are developed based on information on vegetation cover (NDVI, SAVI and fractional vegetation cover) from Landsat 8 imagery. Eight cities are analyzed in the following, representing different vegetative and climatic conditions in Europe while allowing us to investigate the spatial transferability of city-specific regression models (see below), and the development of regional quantification models (Figure 1). The accuracies of the methods are evaluated at a 30 m spatial resolution.

2.1. Landsat Data and Processing

Radiometrically and geometrically corrected Landsat 8 imagery (Level 1T product) were downloaded from the USGS Global Visualization Viewer and converted to top of atmosphere reflectance using the rescaling coefficients provided in the included metadata files [20]. Only images with a clear sky view of the urban areas were included in the analyses, and atmospheric corrections were performed using dark object subtraction to estimate the reflectance at the surface. As Landsat 8 has only recently been launched (data are available from April 2013 onwards), a limited number of clear sky view scenes were available for analyses (Table 1).

Figure 1. Location of case cities within major terrestrial ecoregions in Europe [19].

Figure 1. Location of case cities within major terrestrial ecoregions in Europe [19].

Table 1. Landsat 8 image paths, rows and acquisition dates.

Country	Norway	Germany	France	Austria	UK	Denmark	France	Spain
City	Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona
Path/Row	197/18 197/18	195/23 195/23 196/23	195/26	189/27	204/25 204/25	195/21	194/30	197/31 197/31 198/31
			195/26	190/26		195/21	195/30
			196/26	190/26		195/21	195/30
			196/26	190/27		195/22	195/30
Date (mm/yyyy)	07/2013 09/2013	07/2013 08/2013 10/2013	06/2013	07/2013	06/2013 07/2013	05/2013	05/2013	05/2013 07/2013 08/2013
			07/2013	06/2013		06/2013	07/2013
			04/2013	08/2013		07/2013	08/2013
			07/2013	08/2013		07/2013	09/2013

Table 1. Landsat 8 image paths, rows and acquisition dates.

Country	Norway	Germany	France	Austria	UK	Denmark	France	Spain
City	Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona
Path/Row	197/18 197/18	195/23 195/23 196/23	195/26	189/27	204/25 204/25	195/21	194/30	197/31 197/31 198/31
			195/26	190/26		195/21	195/30
			196/26	190/26		195/21	195/30
			196/26	190/27		195/22	195/30
Date (mm/yyyy)	07/2013 09/2013	07/2013 08/2013 10/2013	06/2013	07/2013	06/2013 07/2013	05/2013	05/2013	05/2013 07/2013 08/2013
			07/2013	06/2013		06/2013	07/2013
			04/2013	08/2013		07/2013	08/2013
			07/2013	08/2013		07/2013	09/2013

2.2. NDVI, SAVI and Fractional Vegetation Cover

The Normalized Difference Vegetation Index (NDVI) is a measure of “greenness” that indicates the relative abundance and activity of vegetation (Equation (1)) [11,12,21]. Upon download of the Landsat 8 Level 1T product, the pixel based NDVI values were calculated based on top of atmosphere reflectance information available in band 4 (${\rho }_{red}$) and band 5(${\rho }_{nir}$).

$$NDVI=\frac{{\rho }_{nir}-{\rho }_{red}}{{\rho }_{nir}+{\rho }_{red}}$$

(1)

The Soil Adjusted Vegetation Index (SAVI) was developed in order to adjust for the influence of variations in soil background color on the ρ_nir/ρ_red relationship for areas with partial vegetation cover [13]. The additional quantity L in SAVI is a soil background adjustment factor, which varies from 0 to 1 for areas with different vegetation density and leaf area index (LAI) (Equation (2)). In the current study, L was set to 0.5 as this has been found to provide a robust index for a wide range of vegetative conditions [13].

$$SAVI=\frac{(1+L)\times {\rho }_{nir}-{\rho }_{red}}{{\rho }_{nir}+{\rho }_{red}+L}$$

(2)

Accurate measurements of sub-pixel fractional vegetation cover (FR) can be derived from remote sensing measurements of NDVI or by use of SMA, where vegetation is considered as one of the end members [15,22]. In the current paper, FR is based on the approach by [15] and is calculated from the relationship between values of actual NDVI (NDVI) and the values for dense vegetation (NDVI_S) and bare soil (NDVI₀) (Equation (3)). NDVI₀ and NDVI_S are not constant in space/time and different values must be applied when calculating FR for different cities. In general, NDVI₀ is determined by topography and the color of non-vegetated surfaces and the satellite sensor used [23]. The value of NDVI_S is mainly affected by vegetation type and should in general be set at the point where the surface is completely covered by vegetation [23,24,25].

$$FR={\left(\frac{NDVI-NDV{I}_0}{NDV{I}_S-NDV{I}_0}\right)}^2$$

(3)

2.3. Reference Impervious Surface Data

Reference IS data, used to validate the accuracy of the Landsat estimates, were created by manual digitalization of HR aerial photography. After digitalization, the data were resampled to 30 m spatial resolution to match the Landsat 8 pixel size. The digitalization was performed in ArcMap using the world imagery, which is available in the ArcGIS software [26]. The spatial resolution and acquisition dates of the HR imagery used for creating the reference IS data varies between 0.1 and 0.4 m from 2009–2012 for the different urban areas (Table 2). To reduce the influence of variability in sub-area characteristics and size, sample areas with comparable mean imperviousness and size were chosen for the analyses.

2.4. Using Vegetation Indices to Map Impervious Surface Fractions

Information on the extent and quantity of vegetation cover may be used to provide fairly accurate estimates of sub-pixel IS fractions for urban areas [14,16]. This is based on the assumption of a strong inverse relationship between vegetation cover and IS, i.e., it is implicitly assumed that non-IS within urban areas are covered with green vegetation. The theoretical relationship between sub-pixel imperviousness and vegetation cover, as measured by NDVI/SAVI and FR, is shown in Figure 2. As NDVI and SAVI have been found to saturate for surfaces characterized by dense vegetation cover, there is an inherent non-linear relationship between NDVI/SAVI and sub-pixel imperviousness for partially vegetated surfaces, including urban areas [27,28].

Table 2. Acquisition dates and spatial resolution of HR imagery used to create reference IS data, mean imperviousness and size of urban sub-areas as measured by the number of Landsat 30 m pixels.

Country	Norway	Germany	France	Austria	UK	Denmark	France	Spain
City	Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona
Resolution (m)	0.3	0.4	0.3	0.3	0.3	0.1	0.3	0.3
Date (mm/yyyy)	08/2011	04/2009	05/2012	08/2011	07/2011	06/2012	11/2012	09/2011
Mean imperviousness of sub-areas	60%	66%	74%	50%	68%	69%	65%	70%
Number of 30 m pixels included	729	777	752	729	722	744	745	727

Table 2. Acquisition dates and spatial resolution of HR imagery used to create reference IS data, mean imperviousness and size of urban sub-areas as measured by the number of Landsat 30 m pixels.

Country	Norway	Germany	France	Austria	UK	Denmark	France	Spain
City	Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona
Resolution (m)	0.3	0.4	0.3	0.3	0.3	0.1	0.3	0.3
Date (mm/yyyy)	08/2011	04/2009	05/2012	08/2011	07/2011	06/2012	11/2012	09/2011
Mean imperviousness of sub-areas	60%	66%	74%	50%	68%	69%	65%	70%
Number of 30 m pixels included	729	777	752	729	722	744	745	727

Figure 2. Conceptual relationship between vegetation indices from Landsat imagery and sub-pixel imperviousness (adapted from [17]).

Figure 2. Conceptual relationship between vegetation indices from Landsat imagery and sub-pixel imperviousness (adapted from [17]).

There are a number of limitations in using VIs for mapping of IS. As mentioned above, these include the presence of bare soils, shadow effects from buildings and trees and tree crowns covering IS [17]. Bare soil has similar spectral characteristics as urban fabrics and hence is easily confused with IS [4,18]. As a result, areas characterized by large shares of bare soil may tend to have low NDVI/SAVI values and will be perceived (falsely) as having a high level of imperviousness. This in turn implies that the presence of bare soils decrease the overall accuracy of VI based estimates. However, when using growing season images to quantify IS fractions for urban areas only, relatively little bare soil is present [17]. In the current analyses, the influence of bare soils is considered to be negligible as bare soils generally only cover a very small area within the central parts of major cities in Europe. This was, e.g., the case for the eight cities probed in this study. Shadows from buildings and tall trees covering vegetation reduce the VI signal from these areas and therefore influences the strength of the relationship between VIs and IS fractions. The NDVI/SAVI signal is inflated for pixels where tree crowns cover IS. The consequence hereof is likewise a reduced accuracy of the VI based estimates.

Assuming that urban areas are comprised only by a combination of vegetation cover and IS, FR can be used to estimate impervious surface fractions [29] (Equation (4)). FR is calculated for each Landsat pixel from measurements of NDVI, NDVI₀ and NDVI_s (Equation (3)) [15]. In the current study NDVI₀ and NDVI_s denote the average values for 100% and 0% imperviousness, respectively, and are estimated individually for each urban area through visual inspection of high resolution images.

$$IS=1-FR$$

(4)

2.5. Regression Modelling

To exclude major areas characterized by bare soil, we establish city boundaries prior to developing and training of the regression models. Urban masks were created by manual digitalization of the Landsat images for the individual cities. Also, the final regression models are only accurate for urban areas where the area characterized by bare soil can be considered to be negligible (an urban mask should be created prior to the use of VIs for mapping of IS fractions), and the approach therefore requires some pre-knowledge of the urban areas where the method is to be applied. As a consequence, hereof it may not be appropriate to use VIs for mapping of IS for many cities in developing countries.

In the initial preprocessing stage (Figure 3A), Landsat imagery is used to create Maximum Value Composite (MVC) NDVI and SAVI images using two to four clear sky view summer images for each urban area (Table 1). MVCs are used to reduce the influence of intra-annual variations in greenness hereby enabling an optimal comparison of the city-specific quantification models. SAVI was tested as an alternative to NDVI to investigate if the inclusion of a background adjustment factor would reduce the impact of intercity variations in construction materials and background color on the relationship between IS and VIs. Geo-rectification was conducted prior to the regression analyses to perfectly align the Landsat data and the high resolution images. Linear and second order polynomial regression analyses (Figure 3B) were performed for each urban sub-area using the Landsat VI composites (MVC NDVI and MVC SAVI) and the reference IS data as the independent and dependent variables, respectively. This yields a total of 16 regression models for each VI (one linear and one second order polynomial for each city). Inverse calibration following [30] was performed for the NDVI/SAVI based regression models to reduce the overestimation of pixels with low IS fractions and the underestimation of pixels with high IS fractions [17]. Inverse calibration is conducted by using the slope and intercept of the linear relation between observed (reference IS) and Landsat estimated IS fractions to calibrate the final quantification. Two sample areas, one from the central parts of the cities and one from the fringes, were selected for training and validation of the regression models for each city, corresponding to a total number of ≈700–800 × 30 m pixels (Table 2). Of these, 50% were used to train the regression models and 50% were used for validation. The location of the sample areas were carefully selected in order to represent most of the heterogeneity in surface features that exist within the urban areas. Also, the location of the sample areas were chosen to ensure a minimum number of 50 pixels for each level of imperviousness (0%–30%, 31%–60%, 61%–90%, 90%–100%) and thus that all levels of IS fractions were represented adequately, i.e., to realistically represent variations in urban land cover. Visual inspections of Landsat TM imagery were performed for all the sample areas at the acquisition time of the aerial images to ensure that major land cover changes had not occurred during the period between the acquisition times of the aerial images (reference data) and the Landsat 8 imagery as this would influence the relationship between vegetation cover and NDVI/SAVI.

Figure 3. Flowchart of image processing and analyses procedures. (A) Preprocess, (B) Regression modelling, (C) Impervious surface data, (D) Accuracy assessment, (E) Spatial transferability analysis.

Figure 3. Flowchart of image processing and analyses procedures. (A) Preprocess, (B) Regression modelling, (C) Impervious surface data, (D) Accuracy assessment, (E) Spatial transferability analysis.

The regression models were used to derive four IS datasets (Figure 3C) for the individual urban sub-areas (two SAVI based products and two NDVI based products). Since regression modeled estimates may be found to exceed 100% and to be less than 0%, IS values greater than 100% and less than 0% were set to 100% and 0%, respectively. The FR was converted to IS fractions to create an additional dataset for each urban sub-area (Equations (3) and (4)).

2.6. Accuracy Assessment and Model Comparison

The performance of the linear and polynomial regression models is evaluated by calculating the Mean Absolute Errors (MAE) and the Mean Bias Errors (MBE) of the VI based IS estimates as compared to the reference IS fractions (Figure 3D) (Equations (5) and (6)). To evaluate the spatial transferability of the regression models, we specifically estimate the loss of accuracy when using a transferred model (e.g., a model developed for another city) as compared to a local model. (Figure 3E). This is done by quantifying the MAEs and MBEs for all possible combinations of IS models and urban sub-areas. Also, regional regression models are developed by aggregating NDVI/SAVI and reference IS data from the cities located within the same terrestrial ecoregions (Figure 1). The accuracies of three different regional regression models are examined for both NDVI and SAVI, namely aEuropean model (all cities), a Central model (Odense, Hamburg, Exeter, Strasbourg and Vienna) and a Southern model (Nice and Barcelona).

The MAE is defined as the absolute difference between the observed (reference IS) and the predicted values (VI based estimates) (Equation (5)). $y_i$ and ${\widehat{y}}_i|$ is the predicted and observed pixel value respectively, and n is the total number of pixels.

$$MAE=\frac{1}{n}{\displaystyle \sum _{j=1}^n}|y_j-{\widehat{y}}_j|$$

(5)

Conversely, the MBE measures the average model “bias”, i.e., how much the models under- or overestimates the mean imperviousness for the urban sub-areas (Equation (6)). $\overline{P}$ is the predicted mean IS fraction for the urban sub-area while $\overline{O}$ is the observed mean IS fraction as estimated from the reference data.

$$MBE=\left(\overline{P}-\overline{O}\right)$$

(6)

The linear and polynomial regression models are compared using standard statistical tests, i.e., the F test and the Bayesian Information Criterion (BIC) to discern which of the two models are statistically superior. Assuming normally distributed variances and independency the F statistic is defined by

$$F=\frac{({RSS}_{\text{1}}-{RSS}_{\text{2}})/({df}_{\text{1}}-{df}_{\text{2}})}{{RSS}_{\text{2}}/{df}_{\text{2}}}$$

(7)

Similarly, BIC is defined by

BIC = n · ln(RSS/n) + k · ln(n)

(8)

where index 1 refers to the simple linear model, 2 to the second order polynomial, n is the total amount of data points, RSS the residual sum-of-squares, k the number of fitted model parameters (e.g., 3 for the second order polynomial) and df = n – k is the number of free parameters. Since the linear model is a special instance of the second order polynomial the F test can be used to evaluate whether the more complex model provides a significantly better fit to the data. Similarly, BIC yields a criterion for model selection among a finite set of models by penalizing the increase in the likelihood achieved by adding additional parameters to a model fit by the number of parameters in the fit.

3. Results

3.1. Accuracy of Impervious Surface Estimates

For all eight cities, a strong inverse relationship between NDVI/SAVI and reference IS fractions was observed (Figure 4). Linear and second order polynomial regression models were both found to provide a good description of the relationship between NDVI/SAVI and reference IS fractions. The degree of non-linearity was seen to vary dependent on the vegetative conditions of the individual cities. Higher maximum NDVI and SAVI values were found in the urban areas located in the northern and central parts of Europe (Oslo, Hamburg, Strasbourg, Vienna, Exeter and Odense) as compared to those in southern Europe (Nice and Barcelona). Likewise, we observed a strong linear relation between IS fractions and FR for all the urban sub-areas (Figure 4), which confirms our assumption that bare soils are only marginally present within the selected urban sub-areas.

Figure 4. Relationship between reference IS fractions and NDVI/SAVI/FR for the eight urban sub-areas.

Figure 4. Relationship between reference IS fractions and NDVI/SAVI/FR for the eight urban sub-areas.

In general, NDVI, SAVI and FR all perform well as indicators of IS fractions (average R² = 0.74–0.78, average RMSE = 13.2%–14.6%) (Table 3). The polynomial regression models are slightly better (higher R² and lower RMSE values) in describing the relationship between the two VI´s and the reference IS fractions as compared to the linear models (Table 3). In all cases, both F and BIC tests (results not shown) clearly indicate that the improvements due to the polynomial models are statistically significant with a probability of this result being due to random scatter virtually zero. That is the polynomial models are statistically found to outperform the linear models. Some inter-city variation was observed with RMSE’s ranging from 11.7%–15.7% for the polynomial models, and from 11.8%–16.2% for the linear models. Corresponding values for FR range from 11.8%–15.2%. Also, some inter-city variation can be seen in how well the VIs explain the variation in IS fractions (Table 3). For the polynomial regression models the explained IS variation ranges from 71% in the worst case to 87% in the best case. Corresponding values are 68%–87% and 72%–83% for the linear models and the FR method, respectively.

Table 3. Linear and second order polynomial regression models, R² and RMSE between reference IS and NDVI/SAVI, and linear regression models and RMSE between reference IS fractions and FR. The improved performance of the polynomial regression models are validated and found to be statistically significant.

			Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	Linear regression	a	−129	−136	−164	−149	−147	−174	−182	−176
		b	129	140	156	140	140	152	148	142
	R2		0.78	0.68	0.71	0.83	0.73	0.74	0.75	0.71	0.74
	RMSE		13.9	16.2	12.6	13.5	15.9	14.3	15.8	14.3	14.6
	Polynomial regression	a	−118	−227	−191	−38	−178	−162	−129	−170
		b	−15	104	34	−106	25	−19	−56	−33
		c	106	83	108	130	104	118	120	115
	R2		0.79	0.74	0.75	0.84	0.77	0.85	0.76	0.73	0.77
	RMSE		13.3	14.5	11.9	13.4	14.9	13.8	15.5	13.9	13.9
SAVI	Linear regression	a	−227	−208	−246	−244	−264	−245	−321	−297
		b	115	120	135	124	131	135	151	143
	R2		0.80	0.70	0.75	0.87	0.79	0.79	0.75	0.72	0.77
	RMSE		13	13.6	11.8	11.8	13.9	12.9	15.8	13.9	13.3
	Polynomial regression	a	−23	−240	−119	25	−378	−209	−194	−370
		b	−216	−76	−179	−258	−80	−129	−211	−112
		c	114	104	126	125	112	121	136	122
	R2		0.80	0.71	0.76	0.87	0.81	0.79	0.75	0.73	0.78
	RMSE		13	13.4	11.7	11.8	13.3	12.6	15.7	13.7	13.2
FR	Linear regression	a	−0.93	−1	−1	−0.92	−1	−1.1	−0.93	−0.91
		b	97.1	99.8	101.5	96.2	100.4	103.7	97.3	97.5
	R2		0.79	0.74	0.75	0.83	0.77	0.75	0.77	0.72	0.76
	RMSE		13.5	14.6	11.8	13.5	14.9	13.8	15.2	14	13.9

Table 3. Linear and second order polynomial regression models, R² and RMSE between reference IS and NDVI/SAVI, and linear regression models and RMSE between reference IS fractions and FR. The improved performance of the polynomial regression models are validated and found to be statistically significant.

			Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	Linear regression	a	−129	−136	−164	−149	−147	−174	−182	−176
		b	129	140	156	140	140	152	148	142
	R2		0.78	0.68	0.71	0.83	0.73	0.74	0.75	0.71	0.74
	RMSE		13.9	16.2	12.6	13.5	15.9	14.3	15.8	14.3	14.6
	Polynomial regression	a	−118	−227	−191	−38	−178	−162	−129	−170
		b	−15	104	34	−106	25	−19	−56	−33
		c	106	83	108	130	104	118	120	115
	R2		0.79	0.74	0.75	0.84	0.77	0.85	0.76	0.73	0.77
	RMSE		13.3	14.5	11.9	13.4	14.9	13.8	15.5	13.9	13.9
SAVI	Linear regression	a	−227	−208	−246	−244	−264	−245	−321	−297
		b	115	120	135	124	131	135	151	143
	R2		0.80	0.70	0.75	0.87	0.79	0.79	0.75	0.72	0.77
	RMSE		13	13.6	11.8	11.8	13.9	12.9	15.8	13.9	13.3
	Polynomial regression	a	−23	−240	−119	25	−378	−209	−194	−370
		b	−216	−76	−179	−258	−80	−129	−211	−112
		c	114	104	126	125	112	121	136	122
	R2		0.80	0.71	0.76	0.87	0.81	0.79	0.75	0.73	0.78
	RMSE		13	13.4	11.7	11.8	13.3	12.6	15.7	13.7	13.2
FR	Linear regression	a	−0.93	−1	−1	−0.92	−1	−1.1	−0.93	−0.91
		b	97.1	99.8	101.5	96.2	100.4	103.7	97.3	97.5
	R2		0.79	0.74	0.75	0.83	0.77	0.75	0.77	0.72	0.76
	RMSE		13.5	14.6	11.8	13.5	14.9	13.8	15.2	14	13.9

The dispersion of the observed regression lines provide evidence for the level of comparability of the IS/VI relationship between the different cities, and ultimately for the potential to develop regional quantification models (regression models which can be used for multiple cities within a larger geographical area) (Figure 5A–E). For both types of regression models, it is evident that Nice and Barcelona are characterized by a lower dynamic range of NDVI and SAVI as compared to the other cities. (Figure 5A–D). The regression models based on SAVI are more clustered for both the linear and polynomial models, indicating that SAVI may perform better (as compared to NDVI) when used as input to the regional regression models. The relationship between IS fractions and FR (Figure 5E) show similar characteristics for all the cities, with only minor variations as compared to those based on NDVI and SAVI. Again, the apparent clustering of several of the regression lines demonstrates that the IS/FR relationship appear to be more similar for urban areas located within the same terrestrial ecoregion (Figure 1).

Figure 5. (A) Linear and (B) second order polynomial regression lines between reference IS fractions and NDVI, and (C) linear and (D) second order polynomial regression lines between reference IS fractions and SAVI, and (E) linear regression lines between reference IS fractions and FR.

Figure 5. (A) Linear and (B) second order polynomial regression lines between reference IS fractions and NDVI, and (C) linear and (D) second order polynomial regression lines between reference IS fractions and SAVI, and (E) linear regression lines between reference IS fractions and FR.

All three indices provide fairly accurate estimates of subpixel imperviousness. MAE’s average is 10.8% and 10.6% for the NDVI-based estimates using a linear and second order polynomial regression model, respectively (Figure 6, Table 4). The corresponding values for the SAVI-based estimates are 10% and 9.9%. The second order polynomial models which are found to be statistically better than the linear models generally also perform slightly better, as they are able to capture the inherent non-linear behavior of the NDVI/SAVI when approaching high VI values. All the estimates are characterized by small negative MBE’s, ranging from −2.3 for the linear NDVI based model to −0.5 for the second order polynomial models (Table 4). The negative MBEs are primarily caused by the large share of pixels with high IS values within urban areas and is a result of using regression modelling for variables with physically fixed boundaries (IS fractions are between 0% and 100%). Because of the uneven distribution of pixels with high/low IS fractions, a larger number of pixels with values >100% are reduced to 100% as compared to pixels having values increased from <0% to 0%. The accuracy of the estimates based on FR (average MAE = 10.4%) is similar to that of NDVI and SAVI with only minor biases (MBE = −0.2–0.4) (Table 4). From the average MAEs and MBEs it appears that the FR method performs slightly better than regression models based on NDVI/SAVI. The inter-city variation in MAEs is found to be large as compared to the inter-model variation, i.e., the accuracy of the different models is very similar for the individual urban areas (Table 4).

Figure 6. Relationship between reference IS fractions and estimates based on NDVI/SAVI/FR. Results are shown for the second order polynomial regression models for NDVI and SAVI.

Figure 6. Relationship between reference IS fractions and estimates based on NDVI/SAVI/FR. Results are shown for the second order polynomial regression models for NDVI and SAVI.

Table 4. MAEs, MBEs and R² for IS estimates based on NDVI/SAVI and FR as compared to reference IS fractions, Lin reg = Linear regression model, Pol reg = Second order polynomial regression model, * absolute mean.

			Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	Lin reg	MAE	10.5	13.1	9.4	9.4	11.2	10.1	11.9	10.9	10.8
		MBE	−2.0	−4.2	−2.1	−0.9	−3.3	−2.6	−1.3	−2.1	2.3*
		R2	0.79	0.71	0.75	0.84	0.77	0.76	0.76	0.72	0.76
	Pol reg	MAE	10.6	12.1	9.1	9.5	11.1	10.1	11.6	10.7	10.6
		MBE	−1.2	−1.3	−0.8	−0.6	−1.5	−1.3	−0.2	−1.0	1.0*
		R2	0.79	0.73	0.75	0.84	0.77	0.76	0.76	0.73	0.77
SAVI	Lin reg	MAE	10.2	12.0	8.8	8.6	9.5	9.1	11.4	10.6	10.0
		MBE	−1.1	−0.9	−1.3	−0.4	−2.3	−2.0	−1.5	−2.0	1.4*
		R2	0.80	0.70	0.77	0.88	0.82	0.80	0.78	0.73	0.78
	Pol reg	MAE	10.2	11.8	8.8	8.6	9.4	9.1	11.3	10.4	9.9
		MBE	−1.1	−0.1	−0.9	−0.4	−1.2	−1.3	−1.0	−1.0	0.9*
		R2	0.80	0.70	0.77	0.88	0.82	0.80	0.78	0.73	0.79
FR		MAE	10.1	11.0	8.9	10.1	10.8	10.9	11.4	10.4	10.4
		MBE	0.2	0.1	−0.2	0.0	−0.1	0.0	0.4	−0.2	0.1*
		R2	0.79	0.74	0.75	0.82	0.77	0.75	0.77	0.73	0.76
		NDVI0	0.11	0.15	0.2	0.16	0.15	0.13	0.16	0.15
		NDVIs	0.85	0.89	0.83	0.83	0.83	0.80	0.71	0.67

Table 4. MAEs, MBEs and R² for IS estimates based on NDVI/SAVI and FR as compared to reference IS fractions, Lin reg = Linear regression model, Pol reg = Second order polynomial regression model, * absolute mean.

			Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	Lin reg	MAE	10.5	13.1	9.4	9.4	11.2	10.1	11.9	10.9	10.8
		MBE	−2.0	−4.2	−2.1	−0.9	−3.3	−2.6	−1.3	−2.1	2.3*
		R2	0.79	0.71	0.75	0.84	0.77	0.76	0.76	0.72	0.76
	Pol reg	MAE	10.6	12.1	9.1	9.5	11.1	10.1	11.6	10.7	10.6
		MBE	−1.2	−1.3	−0.8	−0.6	−1.5	−1.3	−0.2	−1.0	1.0*
		R2	0.79	0.73	0.75	0.84	0.77	0.76	0.76	0.73	0.77
SAVI	Lin reg	MAE	10.2	12.0	8.8	8.6	9.5	9.1	11.4	10.6	10.0
		MBE	−1.1	−0.9	−1.3	−0.4	−2.3	−2.0	−1.5	−2.0	1.4*
		R2	0.80	0.70	0.77	0.88	0.82	0.80	0.78	0.73	0.78
	Pol reg	MAE	10.2	11.8	8.8	8.6	9.4	9.1	11.3	10.4	9.9
		MBE	−1.1	−0.1	−0.9	−0.4	−1.2	−1.3	−1.0	−1.0	0.9*
		R2	0.80	0.70	0.77	0.88	0.82	0.80	0.78	0.73	0.79
FR		MAE	10.1	11.0	8.9	10.1	10.8	10.9	11.4	10.4	10.4
		MBE	0.2	0.1	−0.2	0.0	−0.1	0.0	0.4	−0.2	0.1*
		R2	0.79	0.74	0.75	0.82	0.77	0.75	0.77	0.73	0.76
		NDVI0	0.11	0.15	0.2	0.16	0.15	0.13	0.16	0.15
		NDVIs	0.85	0.89	0.83	0.83	0.83	0.80	0.71	0.67

Figure 7. (A) Average MAE and (B) average MBE for different levels of imperviousness for all the urban sub areas. NDVI Lin = IS fractions from linear regression model, NDVI Pol = IS fractions from second order polynomial regression model, SAVI Lin = IS fractions from Linear regression model, SAVI Pol = IS fractions from Second order polynomial regression model, FR = IS fractions based on FR.

Figure 7. (A) Average MAE and (B) average MBE for different levels of imperviousness for all the urban sub areas. NDVI Lin = IS fractions from linear regression model, NDVI Pol = IS fractions from second order polynomial regression model, SAVI Lin = IS fractions from Linear regression model, SAVI Pol = IS fractions from Second order polynomial regression model, FR = IS fractions based on FR.

The models were found to be most accurate for areas covered with very high levels of imperviousness (90%–100%) and were less accurate for pixels with lower values (Figure 7A–B). For areas with high IS values, the MAEs average is approximately 5% while increasing to 10%–16% for lower levels of imperviousness. A bias towards an overestimation of low values and an underestimation of high values is present for all the cities (Figure 6, Figure 7B). The bias is most pronounced for the FR-based method and the linear regression models, which clearly overestimate (MBE = 6%–12%) pixels with low levels of imperviousness and underestimate (MBE = −2%–6%) those with higher values (Figure 7B). A different pattern is found for the polynomial models where smaller MBEs are observed for all levels of imperviousness (Figure 7B). The bias is primarily caused by the inability of the regression models to adequately describe the non-linear relationship between NDVI/SAVI and IS fractions.

3.2. Spatial Transferability and Regional Regression Models

The relationship between NDVI/SAVI and reference IS for the regional models is comparable with those observed for the city-specific models with R²s of 0.71–0.79 and RMSEs of 13.5–16 (Table 5). Also here, the performances of the non-linear models are slightly better as compared to the linear models while NDVI and SAVI appear to be equally good indicators of IS fractions (Table 5).

Table 5. Linear and second order polynomial regression models, R² and RMSE between reference IS fractions and NDVI/SAVI, and linear regression models, 30 m resolution. European = sub-areas from all cities included in the regression models, Central = sub-areas from Odense, Hamburg, Strasbourg, Exeter and Vienna included in the regression models, Southern = Sub-areas from Barcelona and Nice included in the regression models.

			European	Central	Southern				European	Central	Southern
NDVI	Lin	a	−145	−151	−177	SAVI	Lin	a	−250	−245	−309
	reg	b	138	145	144		reg	b	130	130	147
	R2		0.71	0.74	0.73		R2		0.75	0.78	0.74
	RMSE		16.0	15.2	15.2		RMSE		14.6	13.7	14.9
	Pol reg	a	−100	−123	−122		Pol reg	a	−207	−178	−256
		b	−43	−22	−65			b	−139	−146	−171
		c	115	114	121			c	117	118	130
	R2		0.72	0.76	0.74		R2		0.76	0.79	0.75
	RMSE		15.6	14.6	14.9		RMSE		14.4	13.5	14.8

Table 5. Linear and second order polynomial regression models, R² and RMSE between reference IS fractions and NDVI/SAVI, and linear regression models, 30 m resolution. European = sub-areas from all cities included in the regression models, Central = sub-areas from Odense, Hamburg, Strasbourg, Exeter and Vienna included in the regression models, Southern = Sub-areas from Barcelona and Nice included in the regression models.

			European	Central	Southern				European	Central	Southern
NDVI	Lin	a	−145	−151	−177	SAVI	Lin	a	−250	−245	−309
	reg	b	138	145	144		reg	b	130	130	147
	R2		0.71	0.74	0.73		R2		0.75	0.78	0.74
	RMSE		16.0	15.2	15.2		RMSE		14.6	13.7	14.9
	Pol reg	a	−100	−123	−122		Pol reg	a	−207	−178	−256
		b	−43	−22	−65			b	−139	−146	−171
		c	115	114	121			c	117	118	130
	R2		0.72	0.76	0.74		R2		0.76	0.79	0.75
	RMSE		15.6	14.6	14.9		RMSE		14.4	13.5	14.8

The results of the spatial transferability analysis show that the city-specific models perform fairly well (comparable accuracies as for the cities where the models were developed) in estimating IS fractions for cities located in areas, which are characterized by similar vegetative and climatic conditions (Figure 8,

Table A1. MAE’s and absolute MBE’s for all possible combinations of methods and urban areas using the second order polynomial regression models. Table titles = the city where the models were developed. NDVI = NDVI models, SAVI = SAVI models, FR = models based on FR.

Oslo
		Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	MAE	10.6	13.3	9.8	9.8	10.6	10.2	13.5	13.6	11.4
	MBE	1.2	7.0	4.6	1.5	1.0	3.0	9.8	12.1	5.0
SAVI	MAE	10.2	15.7	16.0	10.1	11.6	16.9	16.8	14.7	14.0
	MBE	1.1	11.3	14.8	5.8	7.3	16.0	11.6	10.6	9.8
FR	MAE	10.1	12.5	10.3	10.1	11.1	11.5	14.2	13.8	11.7
	MBE	0.2	5.9	4.4	0.1	0.3	2.8	9.3	10.7	4.2
Hamburg
		Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	MAE	11.4	12.1	8.8	10.9	11.1	12.6	17.2	18.7	12.9
	MBE	3.1	1.3	3.5	2.5	4.7	10.2	16.3	18.5	7.5
SAVI	MAE	11.9	11.8	8.9	9.1	9.3	9.9	11.8	10.4	10.4
	MBE	8.1	0.1	2.6	3.0	2.8	4.9	0.5	1.8	3.0
FR	MAE	11.0	11.0	9.6	11.4	11.6	12.8	16.4	16.3	12.5
	MBE	6.1	0.1	1.0	6.1	5.0	7.9	14.1	14.8	6.9
Strasbourg
		Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	MAE	13.4	15.3	9.1	11.5	11.2	10.9	15.1	16.8	12.9
	MBE	3.0	6.8	0.8	3.4	0.1	6.2	13.0	16.3	6.2
SAVI	MAE	12.4	12.7	8.8	9.4	9.5	9.7	11.3	10.1	10.5
	MBE	8.5	1.1	0.9	3.0	3.8	4.1	1.9	3.6	3.4
FR	MAE	10.5	11.7	8.9	10.3	10.7	11.5	15.1	16.0	11.8
	MBE	2.2	3.2	0.2	1.8	2.6	6.8	13.2	14.9	5.6
Vienna
		Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	MAE	10.4	13.2	9.9	9.5	10.6	10.1	13.8	13.7	11.4
	MBE	0.3	6.3	4.6	0.6	0.6	3.0	10.2	12.4	4.7
SAVI	MAE	10.4	12.9	11.6	8.6	9.7	12.5	13.8	11.7	11.4
	MBE	4.1	4.9	8.6	0.4	1.8	9.9	5.2	4.2	4.9
FR	MAE	10.3	12.4	9.5	10.1	10.7	11.1	14.2	14.5	11.6
	MBE	0.3	5.4	2.8	0.0	0.5	4.3	10.8	12.6	4.6
Exeter
		Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	MAE	12.9	15.2	9.4	11.1	11.1	10.1	13.7	15.0	12.3
	MBE	3.7	8.2	3.2	4.1	1.5	4.1	10.9	14.2	6.2
SAVI	MAE	11.8	14.9	10.4	10.2	9.4	12.2	12.0	10.3	11.4
	MBE	3.2	4.6	5.8	2.8	1.2	9.5	3.8	1.7	4.1
FR	MAE	10.3	12.6	9.7	10.1	10.8	11.1	14.0	14.2	11.6
	MBE	0.2	5.9	3.5	0.4	0.1	3.7	10.2	12.0	4.5
Odense
		Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	MAE	17.2	19.7	12.4	14.7	13.1	10.1	12.0	12.8	14.0
	MBE	10.4	14.6	9.1	10.5	7.1	1.3	6.4	10.6	8.7
SAVI	MAE	13.6	12.2	8.4	10.0	10.2	9.1	11.6	10.8	10.7
	MBE	11.1	3.9	1.8	5.9	6.4	1.3	4.9	6.3	5.2
FR	MAE	11.6	15.0	11.3	10.7	11.6	10.9	12.7	12.6	12.0
	MBE	4.8	10.4	7.4	4.8	4.0	0.0	6.6	9.1	5.9
Nice
		Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	MAE	20.0	24.0	18.3	17.9	16.2	12.3	11.6	10.4	16.3
	MBE	15.9	21.2	17.3	15.7	13.4	8.7	0.2	4.4	12.1
SAVI	MAE	13.1	16.2	9.9	11.2	9.8	11.9	11.3	10.3	11.7
	MBE	4.0	2.3	2.6	2.4	0.8	7.8	1.0	1.5	2.8
FR	MAE	18.8	22.7	16.7	17.0	15.3	12.2	11.4	10.6	15.6
	MBE	15.2	20.0	15.5	14.9	12.3	7.5	0.4	4.4	11.3
Barcelona
		Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	MAE	23.5	29.1	25.3	22.0	20.5	16.9	13.0	10.7	20.1
	MBE	20.6	27.2	25.1	20.6	19.4	15.7	5.9	1.0	16.9
SAVI	MAE	13.3	17.0	10.7	12.0	10.0	13.2	11.8	10.4	12.3
	MBE	1.8	4.7	5.0	4.6	1.5	10.1	3.7	1.0	4.0
FR	MAE	22.2	27.2	22.4	20.8	19.0	15.4	12.4	10.4	18.7
	MBE	19.7	25.5	22.1	19.6	17.6	13.5	4.6	0.2	15.4
Average
		Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	MAE	14.9	17.8	12.9	13.4	13.0	11.7	13.7	14.0	13.9
	MBE	7.3	11.6	8.5	7.3	6.0	6.5	9.1	11.2	8.4
SAVI	MAE	12.1	14.2	10.6	10.1	9.9	11.9	12.6	11.1	11.6
	MBE	5.2	4.1	5.2	3.5	3.2	7.9	4.1	3.8	4.6
FR	MAE	13.1	15.6	12.3	12.6	12.6	12.1	13.8	13.5	13.2
	MBE	6.1	9.6	7.1	6.0	5.3	5.8	8.7	9.9	7.3

). In general, the estimates based on the SAVI models are most accurate while the NDVI based models perform the worst (Table A1). The average MAEs are 11.6%, 13.9% and 13.2% for the SAVI, NDVI and FR models, respectively. A similar pattern is found with regards to the MBEs with average biases of 4.6%, 8.4% and 7.3% for the models based on SAVI, NDVI and FR models. From this, SAVI appears to be a more comparable indicator of IS fractions across cities as compared to NDVI and FR, and ultimately a more suitable index in relation to the construction of regional quantification models. In general it is seen that models developled for Hamburg, Vienna, Odense, Exeter and Strasbourg appear to be very similar, and therefore perform well in quantifying IS fractions independent of where the model was developed. Conversely the models developed for Oslo do not perform well in estimating IS fractions for any of the other cities (Figure 8, Table A1).

Figure 8. MAEs and absolute MBEs for all possible combinations of urban areas for the second order polynomial regression models based on SAVI. Chart title = the city where the model were developed. NDVI, SAVI and FR are the averages of the different indices for the second order polynomial regression models.

Figure 8. MAEs and absolute MBEs for all possible combinations of urban areas for the second order polynomial regression models based on SAVI. Chart title = the city where the model were developed. NDVI, SAVI and FR are the averages of the different indices for the second order polynomial regression models.

The accuracy of the European and Central models are characterzised by relatively low MAE’s (<≈12%) and MBE’s (<≈6%) for most cities and for both the NDVI and SAVI-based models (Figure 9,

Table A2. MAE’s and absolute MBE’s for IS fraction based on NDVI/SAVI as compared to reference IS fractions for the second order polynomial regression models, European = sub-areas from all cities are included in the regression model, Central = sub-areas from Odense, Hamburg, Strasbourg, Exeter and Vienna are included in the regression model, Southern = Sub-areas from Barcelona and Nice are included in the regression model.

European (all)
		Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	MAE	12.9	15.9	10.5	11.1	11.3	9.8	12.8	13.3	12.2
	MBE	5.0	10.1	6.3	5.3	3.5	1.4	8.6	11.8	6.5
SAVI	MAE	11.7	13.6	9.6	9.4	9.2	10.9	11.6	10.0	10.8
	MBE	5.7	2.0	3.7	0.0	1.0	7.1	1.4	0.5	2.7
Central (Odense, Exeter, Hamburg, Vienna, Strasbourg)
		Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	MAE	11.5	13.8	9.0	10.3	10.7	10.3	14.5	15.6	12.0
	MBE	1.2	6.1	2.0	1.7	0.2	5.3	12.1	14.9	5.4
SAVI	MAE	11.9	12.5	9.0	9.1	9.2	10.1	11.6	10.1	10.4
	MBE	7.5	0.3	2.6	2.1	2.5	5.4	0.3	1.8	2.8
South (Nice, Barcelona)
		Oslo	Hamburg	Strasbourg	Vienna	Exeter	Odense	Nice	Barcelona	Average
NDVI	MAE	21.6	26.3	21.5	19.7	18.1	14.2	12.1	10.3	18.0
	MBE	18.3	24.1	21.0	18.0	16.3	12.1	3.1	1.6	14.3
SAVI	MAE	13.0	16.5	10.3	11.4	9.8	12.5	11.5	10.3	11.9
	MBE	3.0	3.6	4.0	3.4	0.4	9.0	2.5	0.1	3.2

). Expectedly the models (European and Central) perform the worst for Oslo, as it is located in a location with vegetative conditions that differs from the majority of the cities which are used to create the models. Conversely, the Southern model is most accurate for Nice and Barcelona (and to some degree Strasbourg and Exeter), while performing worse for the other urban areas (Figure 9). As for the spatial transferability analysis the regional models based on SAVI perform the best (lowest average MAE’s and MBE’s) as compared to those based on NDVI (Figure 9, Table A2).

Figure 9. MAEs and absolute MBEs using the regional second order polynomial regression models based on NDVI and SAVI. European = sub-areas from all cities are included in the regression model, Central = sub-areas from Odense, Hamburg, Strasbourg, Exeter and Vienna are included in the regression model, Southern = Sub-areas from Barcelona and Nice are included in the regression model. NDVI and SAVI are the averages of the indices for the second order polynomial regression models.

Figure 9. MAEs and absolute MBEs using the regional second order polynomial regression models based on NDVI and SAVI. European = sub-areas from all cities are included in the regression model, Central = sub-areas from Odense, Hamburg, Strasbourg, Exeter and Vienna are included in the regression model, Southern = Sub-areas from Barcelona and Nice are included in the regression model. NDVI and SAVI are the averages of the indices for the second order polynomial regression models.

4. Discussion

4.1. Accuracy Assessment

The results of the of the accuracy assessment, with MAE’s ≈ 10% (Table 4) for the eight urban sub-areas, are comparable with the findings of several previous studies [10,17,31,32,33]. By developing regression models between Landsat TM Tasselled cap greenness and observed IS fractions [17], standard errors ranging from 7.7%–15.9% for 14 sub-areas in Minnesota, US were obtained. It was here concluded that Landsat greenness mapping provides an effective means for estimating imperviousness and changes herein over large geographical areas at a low cost. Similarly, through the adaptation of a Landsat ETM based regression tree modelling approach [31], we found MAEs ranging from 8.8%–11.4% for various combinations of input variables for three study areas in US. At a larger scale, Parece and Campbell 2013 observed MAEs ranging from 2.1%–27% for 11 census tracts in the city of Roanoke using multiple classification approaches based on various combinations of Landsat TM bands. In a European context, and through the application of classification approaches and SMA using SPOT and Landsat TM/ETM data, similar accuracies (MAE’s ≈ 10%) were obtained by [32,33].

The improved ability of the polynomial models to describe the non-linear relationship between VIs and IS fractions is expected to be the primary reason for the relatively small biases of the non-linear models as compared to the linear and FR based approaches. The FR based estimates were found to clearly overestimate low IS fractions and underestimate values for areas which are dominated by high percentage IS (Figure 7B). Using SMA on Landsat ETM data [34] likewise found model overestimation of IS fractions in less developed areas while underestimating high values, while [35] found a very clear over/underestimation of low/high IS values for three cities in Germany using a support vector regression on Landsat ETM imagery. In the current study, it was expected that the presence of boundary NDVI values for complete IS cover (NDVI₀) and zero imperviousness (NDVI_s) would have reduced the bias significantly, which does not seem to be the case. Including an inverse calibration approach to the FR method, similar to that of the NDVI/SAVI regression models, could potentially have reduced this bias, but this was not tested.

As in the current study, [34] found their model to be most accurate when applied for highly developed areas. The better performance of all of the models in the current study (Figure 7A) for very high IS fractions is partly a consequence of the location and characteristics of the urban sub-areas, which are all positioned in within major urban areas and therefore characterized by a high degree of imperviousness (Table 2). As a consequence hereof, the regression models are developed based on a relatively large number of pixels very high IS fractional values (with little or no vegetation) as compared to the number of pixels with lower and moderate values.

The relatively low variability in accuracies across geographical locations (Table 4) confirms the general applicability of the method, and confirms that the relationship between VIs and IS fractions in urban areas. The minor variation in performance that exists between the different cities is, at least partly, caused by variations in the time-lag between the acquisition dates for the high resolution images and the Landsat 8 imagery. The reference data are based on images for the period 2009–2012, whereas all Landsat 8 images are from 2013, which implies a time-lag of one to four years between the acquisitions dates. Any land cover changes having occurred within this period will generally decrease the accuracy of the models. That said, as most urban development occurs on an inter-annual or even decadal scale, and since the most extensive urban development often takes place at the urban fringes [36], it can be assumed that the time difference only moderately influences our analyses. In general, the methods performed the best for Odense and worst for Hamburg, Nice and Barcelona (Table 4). There is a time-lag of four years (2009–2013) for Hamburg while the images for Odense are only one year apart which might explain part of the variation in accuracy. Additionally, the spatial resolution of the imagery used to create the reference IS fractions is higher for Odense (0.1 m) as compared to the other cities (0.3–0.4 m), which is likely to have influenced the accuracy of the reference IS data for Odense and thus ultimately the overall accuracy. A lower dynamic range in NDVI and SAVI for natural vegetation in southern Europe may be affecting the accuracy for Nice and Barcelona. Shadows from tall buildings and trees covering pervious surfaces reduce the NDVI/SAVI signal as measured by the sensor, and consequently reduces the overall accuracy of the models. Visual inspection of the high resolution images suggests that shadows are not a major feature for most of the urban areas covered in the study. Similarly, trees covering IS (e.g., trees planted on sidewalks or near streets) are a common feature in most major cities in Europe. As growing season images are used in the current study, the effect of tree crowns covering IS is expected to reduce the accuracy for all the urban areas and using all of the methods. None of the urban sub-areas were characterized by a high degree of bare soils, hence this is not considered to have had a major impact on the observed variation in accuracies between the different cities. The presence of a strong linear relation between FR and reference IS fractions (Figure 4, Table 3) confirms the assumption that bare soils are only to a lesser degree present within the examined cities, and that urban areas in Europe are in general characterized primarily by either IS or vegetation.

In combining the findings of this study with the results of previous research efforts that are based on earlier Landsat sensors (especially Landsat 5 and 7), it is evident that IS can be mapped solely from Landsat imagery at the same level of accuracy for the past 30–40 years. Consequently, this enables simple and systematic analyses of historical changes in small scale imperviousness for most cities.

4.2. Spatial Transferability and Regional Regression Models

The results of the spatial transferability analysis clearly indicate that city specific VI based models can be used to quantify IS fractions for other urban areas without losing high levels of accuracy. The VI/IS relationship is found to be similar for cities characterized by comparative vegetative and climatic conditions (Table 3, Figure 5) and cross validation of the developed models (Figure 9) show equivalent results with relatively low MAEs and MBEs for a number of different combinations of city-specific models and urban sub-areas. For this purpose, SAVI was found to perform the best, followed by the non-linear NDVI based regression models. As compared to NDVI, SAVI reduces the influence of variations in soil background color/building materials and consequently may improve the inter-city comparability of the regression models, which arguably could be the cause for the better performance of the models based on SAVI. Also, it may be argued that the high level of comparability between the city-specific models could be, at least partly, a result of using MVCs instead of single date Landsat images.

The results of the accuracy assessment for the regional regression models are promising in relation to developing models, which can be applied for multiple cities over a larger geographical area (Figure 9, Table A2). The accuracy of the regional regression models (SAVI: MAE’s = 10.8%, 10.4%, 11.9% and MBE’s = 2.7%, 2.8%, 3.2%) were found to reduce the accuracy only marginally as compared to the results of the individual city-specific models. This demonstrates that VIs, and especially SAVI, show similar characteristics for cities located in areas with similar vegetative and climatic conditions, and that single regression models based on Landsat VIs (SAVI in particular) can be applied uniformly for multiple urban areas with adequate precision. However, if regional regression models are to be developed and applied consistently it is required that they are built on a somewhat larger sample as compared to what has been included in the current study. The development of such models would enable simple and resource efficient quantification of small scale IS fractions accessible to a much wider audience as compared to what is currently the case.

5. Conclusions

As major European urban areas are almost exclusively characterized by a combination of IS and green vegetation, information on vegetation cover from remote sensors can be utilized to provide accurate and cost-efficient estimates of the quantity and spatial distribution of IS and changes herein. From a regression modelling approach, the accuracies of Landsat 8 NDVI, SAVI and fractional vegetation cover (FR) in quantifying IS fractions were examined for eight cities in Europe. The three methods were found to perform almost equally well for the individual cities (MAE’s ≈ 10%, MBEs < 2%), while SAVI was found to be a superior index for mapping of multiple cities using regional models. At the same time, only minor variations in the accuracy were observed for the different cities, which suggest that the method can be applied for urban areas at other geographical locations as long as the requirement of consisting primarily of IS and vegetation cover is fulfilled. In most cases considered in this study, we find that the VI based models are readily transferable with only limited or acceptable loss of precision. Additionally, the accuracy of the regional regression models indicates that it is possible to develop quantification models, which can be used for measuring IS for multiple cities within a larger geographical area. Using vegetation indices to map urban IS is most effective where areas of bare soil is only marginally present and thus the method should be used with caution for cities in, e.g., developing countries. The development of accurate methods for estimating changes in small scale IS fractions for urban areas within a larger geographical context will allow for resource efficient and systematic assessment of long term changes in urban land cover, which are useful for a wide range of applications, including projections of future urban land cover change, environmental impacts of urbanization, and for analysis of the importance of land cover changes for the exposure of cities towards the impacts of climate extremes, such as flooding and heat waves.