The use of combined Landsat and Radarsat data for urban ecosystem accounting in Canada

2.1Urban characterization for ecosystem accounting

Urban ecosystem accounting has considered the use of different levels of detail in the delineation of urban ecosystem assets and types – from the detailed level of the individual green assets, such as rows of street trees or green roofs, to larger patches of urban area with common characteristics to which ecosystem services can be attributed. The focus, in this study, is this latter perspective: analyzing urban morphological structure at the landscape level to define urban ecosystem types. The challenge in doing so lies in the spectral complexity of the urban environment [11]. The urban environment is composed of various natural (e.g., vegetation, bare soil, bare rock and sand) and artificial features (e.g., buildings, roads and urban open space) constructed with different kinds of materials. Many natural surfaces are spectrally similar to the artificial features. Remote sensing analysis should therefore consider not only the spectral characteristics of different materials, but also contextual information [2].

Assessing the ecosystem services generated by urban ecosystem types at landscape level requires that urban areas be subdivided into smaller patches that maintain structural diversity while being regrouped within homogeneous areas [12]. For example, urban heat island studies have developed local climate zones (LCZs), a generic, climate-based typology of urban and natural landscapes that integrates surface structure (height and density of buildings and trees) and surface cover (pervious or impervious) [13]. The classes are local in scale, climatic in nature, and zonal in representation. LCZs are defined as “regions of uniform surface cover, structure, material, and human activity that span hundreds of meters to several kilometres in horizontal scale” and have a “characteristic screen-height temperature regime” [13]. The physical properties of all zones are measurable and nonspecific as to place or time. A generic urban classification scheme, which addresses the constructional characteristics of a surface (i.e., the land cover and not the usage), facilitates the classification of urban land cover into classes representing real world features with an organized system that enables the comparison between urban classes and their components [2].

The Stewart and Oke LCZs classification consists of 17 standard classes, of which 15 are defined by surface structure and cover, and two by construction materials and anthropogenic heat emissions [13].

2.2Earth observation sensors

Satellite earth observation uses remote-sensing technologies including optical, infrared and radar sensors to collect information about physical and other characteristics of the Earth, such as land cover, vegetation, temperature and aerosols. The choice of sensors for the identification of detailed urban ecosystem types should consider and balance between the sensor’s technical characteristics such as spatial, spectral, radiometric and temporal resolution, and feasibility, so that the resulting data are of sufficient quality and are useful for ecosystem accounting. Desirable characteristics include high to very high spatial resolution; good coverage of the electromagnetic spectrum, including number and width of spectral bands; radiometric accuracy and reliability (radiometric resolution is determined by the sensitivity to the magnitude of the electromagnetic energy e.g. 8, 16, 32 bits); and suitable temporal scales, including historical coverage and revisit times. These characteristics will also have an impact on the practicality and efficiency of the processing.

The Landsat series22 satellites are the most common optical Earth observation (EO) data sources for land cover mapping, including for urban, peri-urban and rural areas [14]. Zhang et al. described Landsat data as efficient and cost effective for mapping urban areas in Canada [15]. This is particularly relevant for a large country like Canada, where population centres – areas with a population of at least 1,000 and a population density of 400 persons or more per square kilometre [16] – cover an area of approximately 18,000 km2 (or 1.5% of Canada landmass), distributed across the 1,210,519 km2 (or 12% of Canada landmass) that correspond to the 2016 population ecumene [17].

The Landsat series provides nearly continuous coverage since the early 1970s. Landsat-8 is the latest satellite in the series, providing imagery since 2013. The seven spectral bands (multispectral) at a spatial resolution of 30 m, coupled with the 15 m visible band (panchromatic), offer sufficient spatial and spectral information to map major urban land use classes, such as residential, recreational and industrial/commercial/institutional. Another benefit is that the United States Geological Survey (USGS) provides analysis ready data (ARD) for Landsat 4–8, resulting in a highly accurate surface reflectance image, which significantly reduces the pre-processing required by users [14].

The Global Human Settlement Layer (GHSL) produced by the Joint Research Centre (JRC) of the European Commission, used Landsat data for mapping built-up areas worldwide. These datasets contain a multi-temporal information layer of built-up presence derived from Landsat image collections from 1975 to 2014 [18, 19, 20]. The GHSL-built dataset was combined with a population grid to produce the degree of urbanisation. The degree of urbanisation is described by the classes: 1) cities, 2) towns and suburbs, and 3) rural areas [21, 22]. However, a preliminary review of the accuracy of the GHSL built-up area in Canada has identified some issues [5]. In some cases, the product struggles to accurately delimit the urban extent or cannot differentiate the urban landscape into separate urban classes. However, although GHLS overestimates the area of settlements in some places and does not capture the built-up area in others, it remains a useful tool as a basis for change analysis and comparison of built-up areas.

The Sentinel-2 satellites from the European Space Agency (ESA) Copernicus programme can be used to complement Landsat-8 imagery for urban area mapping and monitoring [22, 23], although the satellites (A and B) were only launched recently (2015 and 2017) reducing the capacity to do historical comparisons. The twin satellites have 13 spectral bands, with spatial resolution of 10 m, 20 m and 60 m in the visible, near infrared and shortwave infrared bands of the electromagnetic spectrum. They offer a revisit time of five days.

Very high resolution images (< 5 m) can map detailed urban features and can therefore provide useful additional information to understand the heterogeneity of urban areas that is not otherwise available using Landsat or Sentinel. However, such datasets are costly and they require a more significant use of computing resources for analysis compared to Landsat data, especially when covering a large area. Also, with very high resolution images, the heterogeneity of urban areas is amplified by the illumination of surfaces exposed to the sun and the resulting shadows create noise in the signal measured on the image [24].

Synthetic aperture radar (SAR) has been recognized as an effective tool for urban analysis, as it is less influenced by solar illumination or weather conditions (e.g., cloud cover) compared to optical or infrared sensors [25]. SAR sensors can acquire information in different polarizations (these are related to the orientation of the electromagnetic field), which helps provide a more complete description of the target. SAR data have been increasingly used for urban land use and land cover (LULC) classification [26, 27, 28]. The difficulty in detailed urban mapping using SAR data is mainly due to the complexity of the urban environment – its many different kinds of natural and built objects and materials, orientations, shapes, sizes, etc. complicate the interpretation of SAR images [25, 29, 30, 31].

The integration of optical and polarimetric SAR data is seen as a promising approach to improve urban impervious surface mapping [29, 32], because it helps resolve issues with the diversity of urban land covers and the spectral overlap between covers. The Global Urban Footprint (GUF) project from the German Aerospace Agency (DLR) demonstrated that SAR data can be used with optical data to delineate urban areas. SAR data usefully captures built environments, since these exhibit strong scattering properties as a result of double bounce effects and direct backscattering from vertical building structures [33, 34, 35, 36]. The GUF project, like most studies about urban mapping using SAR data was, however, limited to the identification of the urban extent and relative built-up density.

The approach developed by DLR for mapping urban areas using Sentinel-1 data (SAR) has also been tested with success on Radarsat-2 images [37]. In addition, previous comparative studies on the combined use of multispectral optical and polarization SAR data to identify urban areas have demonstrated that polarization data improved the accuracy of results [28, 38, 39].

The Canadian Radarsat-2 https://mdacorporation. com/geospatial/international/satellites/RADARSAT-2/ satellite has a C-band synthetic aperture radar (SAR) with multiple polarization modes, including a fully polarimetric mode in which HH, HV, VV and VH polarized data are acquired. The spatial resolutions range from 1 m to 100 m. The radar signal in different polarizations determines physical properties of the observed object and the use of polarimetric parameters can provide a more complete characterization of the target. Polarimetric target decomposition has become the standard method for the extraction of parameters from polarimetric SAR data [39]. Once decomposed, the physical properties of the target are extracted and interpreted through the analysis of the simpler responses [31, 40, 41, 42].

Touzi multiresolution decomposition is a method that considers the target’s structure to generate the decomposition [41, 43]. It uses a scattering-vector model to represent each coherency eigenvector in terms of unique target characteristics. Each coherency eigenvector is uniquely characterized by five independent parameters. Scattering type is described with a complex entity, whose magnitude (alpha_s) and phase (phi) characterize the magnitude and phase of target scattering. The helicity (tau) characterizes the symmetric-asymmetric nature of the target scattering. The angle (psi) is the target orientation. To complement the decomposition, discriminators can be used to provide the intensity of the target’s response. The algorithm developed by Touzi [40] derives the maximum and minimum degrees of polarization and total intensity of the scattered wave. Touzi decomposition and discriminators have been shown to provide better classification for urban areas with high variability and complexity, compared to other decomposition methods [31].

2.3Classification approach

Depending on the resolution of the data provided by the sensor, pixels may contain more than one type of land cover. This is particularly a problem in heterogeneous areas like cities, where a given 30 m pixel may contain built structures, impervious surfaces, and different types of vegetation [11]. Conventional pixel-based classifiers, such as maximum likelihood classification (MLC), cannot effectively handle the mixed-pixel problem in urban areas. In the MLC classification, the distribution of each class is assumed to come from a normal distribution. The probability that a given pixel belongs to a specific class is calculated using the mean vector and the covariance matrix. Each pixel is then assigned to the class that has the highest probability (the maximum likelihood) [11].

Moreover, in the pixel-based approach, the proportion of the signal coming from the surrounding area is ignored. Alternative approaches that use geographic object-based image analysis (GEOBIA) with machine learning classifiers provide better results [24].

The GEOBIA approach provides a geographic information system-like functionality for classification that integrates not only the spectral but the spatial characteristics during the image segmentation process. The result is individual areas with shape and spectral homogeneity referred to as segments or objects [44]. This method can make use of varied data sources including multi-sensors, geospatial datasets and any spatially-distributed variable (elevation, slope, population density) [11]. GEOBIA uses a two-step approach. First, the segmentation process defines homogeneous objects and second, classifies these objects into a classification scheme based on spectral, spatial and contextual information. Objects are assigned to land cover classes representing real-world features (e.g., open low-rise, compact medium-rise, sparsely built, dense trees, and paved), instead of a somewhat arbitrary pixel structure [45]. Object-based methods, which make use of contextual information to improve mapping accuracy, are increasingly employed [11, 29, 46].

A comparison of pixel-based and object-based approaches for identifying urban land-cover classes has demonstrated that the object-based classifier produced a significantly higher overall accuracy. In addition, segmentation procedures used in GEOBIA can be applied following a hierarchical structure where classes are defined at an appropriate scale [29, 47].

Nonparametric classification approaches using machine learning algorithms such as artificial neural network, decision tree (DT), random forest (RF) and support vector machine (SVM) are widely used in land cover classification [30, 46]. Niu, X. and Y. Ban used the object-based SVM for the classification of Radarsat-2 polarimetric parameters in order to produce an urban land cover classification [30]. The results showed that SVM can map detailed urban land cover classes with high accuracy (90%); it provides homogeneous mapping results with preserved shape details; and outperforms other land cover classification approaches in a complex urban environment with limited training samples [30].

Machine learning classifiers were used to map seven land cover classes in the United Kingdom with Landsat data [48]. Results from this study suggest that the RF classifier performs equally well to SVM in terms of classification accuracy and training time. This study also concluded that fewer user-defined parameters are required by the RF classifier than for SVM and that they are easier to define. However, the performance of the classifiers was highly dependent upon several factors such as the training set design (i.e., sample selection, composition, purity and size), even though the nonparametric algorithm does not require normally distributed training data [11], as well as input imagery resolution and landscape heterogeneity [49].

Finally, other studies have demonstrated that the integration of indices such as Normalized Difference Vegetation Index (NDVI) and Normalized Difference Water Index (NDWI) into the segmentation and classification processes can improve the results for urban characterization [50, 51].

2.4Change detection

Images acquired on different dates rarely capture the landscape surface in the same way due to many factors including illumination conditions, view angles and meteorological conditions [52]. Change detection methods that use the pixel as the change unit are negatively affected by these factors. Problems with pixel change detection are usually caused by the complexity of change conditions such as small spurious changes, the accuracy of image registration, and shadows resulting from different viewing angles, which can be dominant in urban areas. Pixel-level comparison for change detection therefore requires further analysis to attach a meaning to the change observed [53]; for example, to assess a decrease of forest canopy, a geospatial analysis is required to define forest patches. By contrast, the GEOBIA approach transforms spectral and spatial properties of the image into meaningful objects, thereby creating features that have real-world meaning (e.g., buildings) [52, 53]. Object-based change detection reduces the problems associated with assessing pixel-level changes.

Object-based analysis offers great potential for identifying and characterizing LULC change processes, especially when combined with multi-temporal analysis. Segmentation-generated objects from different dates often vary geometrically, even though they represent the same geographic features. Instead of separately segmenting these images, the use of multi-temporal object change detection takes advantage of all multi-temporal states of the scene [54]. Specifically, temporally sequential images are combined and segmented together, producing spatially corresponding change objects [52]. For example, Schneider et al. used the NDVI generated from Landsat time series from 1985 to 2015 to detect urban extents, identify the conversion types and the date the change occurred [55].

3.1Study area

Milton is a town located in Southern Ontario, Canada, and is part of the Toronto census metropolitan area. Between 2001 and 2011, Milton was the fastest growing municipality in Canada, with a 71.4% increase in population from 2001 to 2006, and another 56.5% increase from 2006 to 2011. In 2016, Milton’s population was 110,128 and it is forecasted to grow to 228,000 by 2031 [56]. While most of the development is suburban in nature, larger industrial lots are being developed, and commercial and institutional facilities have been built as well. The built-up patterns in Milton are fairly representative of suburban growth within the Toronto CMA, which represents 44% of the population of Ontario and 17% of the total population of the country.

3.2Satellite data and pre-processing

Radarsat-2 fine quad pol wide images33 in different angles were acquired for Milton for the months of July to October 2019 to test the multi-polarization information. All images were converted in Sigma Nought format and then orthorectified. The pixel size of all images, even if they were acquired in different angles, was resampled at 8 m for better comparison. A Lee adaptive filter with a window of 7 by 7 was applied on the images to reduce the speckle effect [57]. The Touzi multiresolution decomposition [41] was performed using Geomatica software [57] and created an image of 15 decomposition parameters. For each of the primary, secondary, and tertiary eigenvectors, the orientation angle (psi), dominant eigenvalue (lambda), alpha_s angle, phase, and helicity (tau) are computed.

These parameters characterize the properties of scattering by computing the proportion and type of the scattering mechanism for all features in the image. For example, the Alpha S polarimetric parameter measures the double bounce response resulting from the interaction of the radar signal with a corner reflector created by two or more smooth, flat surfaces oriented at a 90-degree angle to each other (most commonly an impervious surface and a building). Four Touzi discriminators [40] were also generated for measuring the intensity of the signal returned to the sensor [57]. The computed discriminators include the Touzi anisotropy, the degrees of minimum and maximum polarization response, and the difference between the maximum and minimum response. Table 1 provides the list of the parameters used for the assessment.

Landsat analysis ready data (ARD) images were selected from the USGS Earth Explorer platform portal for census years 2001, 2006, 2011, and 2016. Imagery was also downloaded for the most recent year (2019) to have the most up-to-date classification possible. Multiple images were downloaded for each year to cover seasonal variability from spring to fall. The NDVI was calculated for each image.

3.3Training and validation data selection and classification scheme

The 2016 Census of Population geographic boundary for the population centre of Milton [58] was used to delineate the urban core and buffered by an arbitrary distance of two kilometres to include the newly built-up areas in the peri-urban area. A buffer distance of two kilometres was found to capture the majority of new built-up areas that occurred outside of the population centre between 2016 and 2019. The 2019 Radarsat-2, Landsat ARD and the buffered population centre polygon for Milton were loaded into a project using the commercial software eCognition V9.4.44 A multiresolution segmentation, based on the 2019 Landsat NDVI from spring, summer and fall, was generated to create meaningful objects, i.e., different urban ecosystem types. Identification of training classes was performed through visual interpretation of Landsat data and high resolution imagery available through Google Earth and on-site visits to assign a class to each polygon.

The approach used a classification scheme based on the local climate zone (LCZ) classes developed by Stewart and Oke for the urban ecosystem types (Table 2) [13]. For each LCZ class, 75% of the polygons were randomly selected for training. The remaining 25% were kept for validation. A total of 723 polygons were selected for training and 237 for validation (Table 3). This proportion was chosen to balance between what is statistically sound and what is practicably attainable.

3.4Object-based classification

The decision tree (DT), random forest (RF), and support vector machine (SVM) supervised classification models were evaluated. Each classifier was trained and validated, and the accuracy of each classification was compared. Radarsat-2 polarimetric parameters Alpha S, Dominant Lambda, Tertiary Lambda and Maximum Intensity, and Landsat spectral bands and NDVI were used to create classifier statistics from the same set of training samples. The accuracy of the three classifications was assessed using the same set of validation polygons.

3.5Ecosystem assets calculation

Urban ecosystem accounting requires the delineation of individual ecosystem assets that make up the entire fabric of urban areas (e.g., vegetation, buildings, impervious surface, bare soil and water).

The total area and the percentage of vegetation coverage by each LCZ class were calculated using a binary map of vegetation based on the NDVI from Landsat imagery. All available Landsat NDVI imagery was acquired for 2019 and the maximum NDVI value for each pixel was selected. A binary vegetation mask was generated by manually specifying a threshold for the maximum NDVI value observed for 2019.

The road coverage was derived from the Road Network File (RNF) from Statistics Canada, a digital representation of Canada’s national road network [59]. All road line features were buffered by four metres55 to create road polygons in order to calculate the road area of each LCZ class.

The building coverage was drawn from Microsoft’s Building footprint dataset [60] dated March 2019. This dataset includes the footprint of nearly 12 million buildings, covering all provinces and territories in Canada. The total area and percentage coverage of the building footprints were calculated for each LCZ class.

3.6Change detection

Object-based change detection analysis captures change by comparing the values of objects over several dates. The annual mean maximum NDVI values for 2001, 2006, 2011, and 2016 census years, and 2019, were compared to identify the type of landscape change and timing of its occurrence for each built-up image object from the 2019 segmentation. To calculate the annual maximum NDVI value for a given year, all available Landsat NDVI images were acquired and the highest NDVI value for each pixel was selected. The NDVI values within each built-up object were averaged to obtain the mean maximum NDVI value. Assessing the change directly from the NDVI value reduces the error introduced by post classification change analysis, where misregistration and misclassification can lead to false changes.

The Landsat-8 sensor acquires data in two thermal bands, in addition to the seven spectral bands. The USGS systematically processes the thermal bands to create a surface temperature product, which is available for all ARD image dates. A single surface temperature image was selected for each year of the analysis. The individual images were chosen to coincide with typical summer temperatures and with daily temperature readings from two nearby weather stations. Local surface temperature is dependent on a number of factors such as elevation, vegetation, topography, solar radiation and prevailing winds.

Surface temperature was calculated for each object by averaging the individual pixel values within the object. The maximum surface temperature value for each date, as measured at the weather stations, was then subtracted from the mean value of each object to adjust for differences between dates. The final surface temperature values for each object represent the deviation from the maximum surface temperature observed at the weather stations.