Abstract
This paper describes an approach for combining Landsat and Radarsat satellite images to generate national statistics for urban ecosystem accounting. These accounts will inform policy related to the development of mitigation measures for climatic and hydrologic events in Canada. Milton, Ontario was used as a test case for the development of an approach identifying urban ecosystem types and assessing change from 2001 to 2019. Methods included decomposition of Radarsat images into polarimetric parameters to test their usefulness in characterizing urban areas. Geographic object-based image analysis (GEOBIA) was used to identify urban ecosystem types following an existing classification of local climate zones. Three supervised classifiers: decision tree, random forest and support vector machine, were compared for their accuracy in mapping urban ecosystems. Ancillary geospatial datasets on roads, buildings, and Landsat-based vegetation were used to better characterize individual ecosystem assets. Change detection focused on the occurrence of changes that can impact ecosystem service supply – i.e., conversions from less to more built-up urban types. Results demonstrate that combining Radarsat polarimetric parameters with the Landsat images improved urban characterization using the GEOBIA random forest classifier. This approach for mapping urban ecosystem types provides a practical method for measuring and monitoring changes in urban areas.
Introduction
Urban areas are among the most complex manifestations of man’s impact on the environment [1]. They consist of a wide array of heterogeneous materials with diverse properties and multi-faceted interactions. Combinations of buildings (e.g., low- and high-rises), impervious surface covers (e.g., roads and parking lots), vegetation (e.g., parks and sports fields), bare soil (empty lots and unattended garden plots) and water (e.g., wetlands and streams) are fundamental components of the urban ecosystem [1, 2]. More research is needed on the impact of land conversion on ecosystem assets and services in urban areas, and on the resilience of urban areas to climate change given its pressures on ecosystems, both in Canada and elsewhere.
Statistics Canada is creating a suite of ecosystem accounts that link environmental information to socio-economic data available in the System of National Accounts and elsewhere in the national statistical system. These ecosystem accounts will contain information on both the quantity and quality of ecosystem assets and their services, in both physical and eventually, monetary terms. The System of Environmental – Economic Accounting – Experimental Ecosystem Accounting (SEEA-EEA) provides guidelines for compiling ecosystem accounts [3].
Considering the large and growing proportion of the world population living in cities, with 68% of the world population projected to live in urban areas by 2050 [4], an emerging interest concerns ecosystem accounting for urban areas. Urban areas can be considered a specific category of ecosystem, or can alternatively be understood as a combination of multiple ecosystem types.1 In both cases, urban areas produce and benefit from ecosystem services.
Statistics Canada has already developed statistical data and analysis for urban areas by integrating maps generated from satellite Earth observation (EO) data and socio-economic data to provide new insights on urban expansion and natural and semi-natural ecosystems impacted by changes in the landscape surrounding metropolitan areas [5]. This analysis made use of geospatial datasets produced by other federal departments, including Agriculture and Agri-Food Canada [6, 7] and Natural Resources Canada [8]. These datasets, although of high quality, were produced to meet the mandates and needs of those departments for tracking cropland and forests. However, they do not provide the precision needed to delineate urban areas as a specific ecosystem type that combines impervious surface, buildings, vegetation and water or to delineate intra-urban level ecosystem assets. As such, they do not meet the ‘fit for purpose’ criteria required of spatial datasets for more focused urban ecosystem accounting.
In 2018, Statistics Canada and the Canadian Space Agency (CSA) signed an agreement, under the CSA’s Government Initiatives Related Program (GRIP), to collaborate on the improvement of urban delineation and characterization in the context of Climate Change Impacts and Ecosystem Resilience (CCIER). The Earth Observation for Ecosystem Accounting in Canada (EO4EA-Can) project aims to integrate earth observation data with socio-economic information to better measure the climate change resilience of built-up areas in the context of ecosystem accounting. The main objective of this project is to map urban areas, including both the built-up structures and ‘green and blue infrastructure,’ and assess the changes in the extent and density of these types of assets, using optical and radar imagery, complemented with census and geospatial data. Green and blue ecosystem assets are of specific policy interest in urban areas since they provide ecosystem services that can improve human health and well-being and buffer against natural disturbances and extreme weather events. Maintaining these services may have a crucial role in increasing adaptive capacity to cope with climate change [9].
This case study on the integration of Landsat and Radarsat for urban ecosystem accounting is a first step towards this larger objective. Milton, Ontario has been used as a test case and subsequent test cases are being developed in other regions of the country. The development of better spatial data to support urban ecosystem accounting will help inform policy and investment decisions related, for example, to the identification and development of mitigation measures for climatic and hydrologic impacts (i.e., those relating to urban heat islands and energy consumption, air pollution, increased runoff, modified streamflow dynamics, or water quality) at a regional scale across the country; and in this way will contribute to the health, security and well-being of Canadians.
Background
The delineation and characterization of areas into a complete set of mutually exclusive and contiguous spatial units is at the foundation of ecosystem accounting. The Technical Recommendations of the SEEA-EEA [3] defines three main spatial units relevant to the production of accounts: 1) ecosystem assets, 2) ecosystem types and 3) ecosystem accounting areas. A basic spatial unit is additionally described in the recommendations; it represents fine-scaled gridded cells or polygons that underlie the delineation of ecosystem assets.
In ecosystem accounting, stocks of ecosystem assets generate flows of ecosystem services. Individual ecosystem assets can be grouped according to ecosystem types (e.g., forest, wetland, cropland, urban) that are expected to generate a similar bundle of ecosystem services [3]. Typically, ecosystem accounts are compiled and presented by ecosystem type rather than by individual ecosystem asset and they may be reported for different ecosystem accounting areas (e.g., at a national or sub-national level or for a specific ecosystem type (e.g., forest or urban). Thus, a classification describing the ecosystem types and a map showing their occurrences are essential components of ecosystem accounting that allow tracking changes over time.
If spatial data on ecological characteristics are not available to delineate the ecosystem assets, a land cover map may be used as a starting point [3]. Land cover can be defined as the “observed (bio) physical cover of the Earth’s surface,” [10] and is a synthesis of the many processes taking place on the land. It reflects land occupation by various natural, modified or artificial systems. Land cover is one of the most easily detectable indicators of human intervention on land [10].
Urban 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].
Earth 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 series2 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 km
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 (
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
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].
Classification 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].
Change 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].
Methods
In support of the creation of detailed urban maps for ecosystem accounting, the objective of this case study was to develop a classification approach to improve the accuracy of urban land cover maps generated from Radarsat-2 and Landsat Operational Land Imager (OLI) observations at 30 m. This paper describes both the mapping approach using machine learning classifiers applied on Radarsat-2 polarimetric parameters and the characteristics of the produced map.
Study 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.
Satellite data and pre-processing
Radarsat-2 fine quad pol wide images3 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.
List of Radarsat-2 polarimetric parameters
List of Radarsat-2 polarimetric parameters
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.
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.4 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.
Description of local climate zone classes found in Milton, Ontario [13]
Description of local climate zone classes found in Milton, Ontario [13]
Number, area and average area of polygons used for training and validation per LCZ class
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.
Ecosystem 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 metres5 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.
Change 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.
Results and discussion
Results show that the polarimetric parameters contributed to the classification of the local climate zone (LCZ) urban ecosystem types. The best result is obtained with the random forest (RF) algorithm, which had an overall accuracy of 74% and a kappa of 0.68 (note that a kappa value of 1 represents perfect agreement, while a value of 0 represents no agreement) (Table 4). The RF algorithm applied on the Landsat-8 combined with the Radarsat-2 image provided better results for most of the classes when compared to the Landsat-8 only classification (Table 5), for both the producer accuracy and user accuracy.
Comparison of the overall accuracy for the classifications
Comparison of the overall accuracy for the classifications
Comparison of the producer’s and user’s accuracy for Landsat-8 (L8) and Landsat-8 combined to Radarsat-2 (RS2) at class level (%)
Producer accuracy is the probability that a value in a given class was classified correctly from the point of view of the map maker. User accuracy is the probability that a value predicted to be in a certain class really is in that class. The probability is based on the fraction of correctly predicted values to the total number of values predicted to be in a class. For example, in Table 5, the producer accuracy for the ‘Compact low-rise’ class was 83% while the user accuracy was 56%. This means that even though 83% of the reference sites classified as Compact low-rise areas have been correctly identified as Compact low-rise, only 56% percent of the areas identified as Compact low-rise in the classification were actually Compact low-rise. The producer’s accuracy is a complement of the omission error, while the user’s accuracy is a complement of commission error. Errors of omission refer to reference sites that were left out (or omitted) from the correct class in the classified map. In the case of Compact low-rise, 17% of this class was omitted, while 44% of polygons identified as Compact low-rise were, in reality, other classes.
Polarimetric decomposition parameters can explain the target scattering mechanisms in different ways. However, their basic scattering mechanisms include three main classes: surface scattering, volume diffusion scattering and double bounce scattering. Table 6 explains the way in which polarimetric parameters contributed to the improvement of the classification of the LCZ urban type classes.
Scattering mechanisms associated with the polarimetric parameter for each LCZ class
The overall accuracy of the results of this approach may seem low compared to the literature, where a threshold of 85% is typically recognized as an acceptable level of accuracy [61]; however, the results should be analyzed at the class level. Moreover, even if the accuracy of one class is low, if it represents a green and blue ecosystem type that is important for defining ecosystem services, then this class can be conserved in the result. Nonetheless, the assessment of this class should be analyzed with ancillary data to avoid unacceptable errors.
Given that urban type class definitions from LCZs are based on a combination of structural and biophysical criteria observed on the ground, the use of remote sensing to map LCZs classes may introduce overlap between classes. However, this overlap can be explained when looking at the spectral response. For example, the accuracy of the ‘Paved’ class is low compared to the other classes (45% for producer accuracy and 50% for user accuracy) because of its confusion with other built-up classes that have a paved component (Table 7). This confusion is judged acceptable; confusion with a class such as water or dense trees would have been unacceptable. In this case, however, given the lower relevance of paved areas for ecosystem service supply, the class could be merged with another related class.
Error matrix of the random forest classification – Landsat-8 and Radarsat-2 for each LCZ class
Note that since the accuracy assessment was performed at the object level, the comparison is based on the count of the polygons, even though the objects do not have the same size. A comparison using the area of the objects, where each pixel is considered in the calculation, increases the user’s accuracy for the Paved class from 50% to 80%. The level of accuracy is higher for all classes when calculated with the area, but the relative comparison between classes is similar (Table 8).
Accuracy assessment calculated with the area (number of pixels) %
The RF object-based classification identified 10 homogeneous urban sub-classes based on the LCZs for Milton in 2019 at 30 m resolution (Fig. 1). All of the urban classes inside the socio-economic boundary of the buffered population centre were summed to give a total urbanized area for Milton. The total area calculated based on the RF object-based classification is 45.3 km
Areas by local climate zone class
Milton urban local climate zones – 2019.
The application of the GEOBIA segmentation and classification approach delineated Landsat images segmented into spatially continuous and homogeneous regions according to LCZ classes based on spectral, spatial and contextual information. The error matrix calculated from the validation dataset demonstrated that the combination of Landsat8 spectral bands and NDVI with the Radarsat2 polarimetric decomposition parameters provided a more accurate urban characterization than the Landsat images alone (Table 5).
The use of ancillary geospatial datasets for measuring ecosystem assets can also be used to validate whether the pixel grouping produces meaningful results that meet urban ecosystem accounting requirements. The areas of buildings, buffered roads and vegetation were calculated for each LCZ class from the geospatial datasets described in the previous section. Limited information was available for impervious surfaces (e.g., parking lots), bare soil and water. These land cover classes were included as ‘Others’ by subtraction from the total area. Figure 2 provides the ratios of the building, road and vegetation areas for each LCZ class. The ratios are compliant with the description of the LCZ classes as described in Table 2 and as estimated by Stewart and Oke [13]. Note that the percentage of vegetation in the water class is caused by the inclusion of shorelines, islands and forest canopy in the object.
Ratio of ecosystem assets by LCZ class.
In defining the segments corresponding to the LCZ classes (or urban ecosystem types), it is helpful to consider how urban ecosystem assets differ in their ability to supply ecosystem services to people. For example, different levels of vegetation will impact climate regulation services in mitigating urban heat island effects. Figure 3 shows that the vegetated urban classes experienced lower temperatures compared to the urban classes where building and road area were more predominant, based on Landsat8 surface temperature data for August 1, 2019. The highest temperatures were seen in the Compact low-rise class, which also accounted for the largest proportion of the town’s population, indicating that a greater number of people were exposed to higher temperatures (Fig. 4).
Surface temperature by LCZ for August 1, 2019.
Total population estimated by LCZ.
Changes in land cover composition are important landscape characteristics for assessing changes in ecosystem condition and functions. The object-based change detection analysis assessed the vegetation change at the landscape level (segment) using NDVI.
Figure 5 provides an example of the changes between 2001 and 2019 for one segment. The 2001 and 2006 sub-images show agriculture fields. Between 2006 and 2011, there was a conversion of agricultural fields to a combination of built-up and bare soil. Between 2011 and 2016, the entire object was converted to built-up area. Overall, 0.05 km
Changes observed in polygon 688 from Landsat series, 2001 to 2019.
Mean maximum NDVI values for polygon 688 from 2001 to 2019.
Mean temperature for polygon 688 from 2001 to 2019.
Changes in ecosystem characteristics have an impact on ecosystem conditions. For example, the surface temperature measured with Landsat8 (adjusted with ground sensors) demonstrated that the conversion of agricultural land (vegetated during the growing season) to built-up (buildings and impervious surface) correlates with an increased ground-level temperature (Fig. 7).
This study investigated Radarsat-2 polarimetric parameters combined with Landsat data for urban delineation and characterization using object-based machine learning classifiers. The approach was successfully tested in Milton, Ontario, and is being applied at two other study sites in Drummondville, Quebec, and Richmond, British Columbia to evaluate its replicability before a broader application of the method to other population centres across the country. The goal would be to eventually apply this approach to the 1,005 population centres in Canada, which represent 28.6 million people or 81% of Canada’s total population, and approximately 18,000 km
Based on the results, the following conclusions are drawn:
Geographic object-based image analysis (GEOBIA) using the local climate zone (LCZ) classes provides a meaningful description of biophysical processes that can be used to define ecosystem services (e.g., temperature increase mitigation). Although developed for urban heat island studies, the LCZ classification provides representative groupings of individual urban components within urban ecosystems. The LCZ classification is flexible enough to apply to a wide variety of urban environments. An adoption of the approach used in this study could improve ecosystem accounts in other contexts and could be explored by other countries. The approach may work best in urban centres where urban planning and zoning results in developments with homogeneous urban types. The polarimetric parameters Alpha S, Lambda 1 and 3 (Touzi decomposition) and Maximum Intensity (Touzi discriminators) contributed to the improvement of the urban characterization following the LCZ scheme. The random forest classifier provided better results than the decision tree and support vector machine classifiers. The object-based change detection tracks the change at the patch (parcel) level, which is better adjusted to the heterogeneity of urban landscapes. The object-based approach combines the spectral properties of urban features reducing the effects associated with the single pixel that are not constant through time, thus the change detected at the object level represents a real conversion of land cover. Since the classifier is dependent on the quality and quantity of training and validation datasets, the selection of areas should be considered as a critical step toward the classification process. Sampling methodology should be scientifically and statistically sound. Moreover, a good coverage of the spatial variability of the training sites can contribute to the generation of unique set of spectral signature and/or polarimetric parameters, and then facilitate transferability between study sites.
Further research is needed to explore:
Use of Sentinel-2 images with their resolution of 10 m, integrated with the Landsat images in the eCognition project, as this could improve the segmentation process, producing a more accurate and detailed delineation of the objects and facilitating the classification of the objects into LCZ classes. Use of cloud-based geospatial platforms, such as Google Earth Engine, to access imagery and enable analysis by facilitating scaling of the methodology to sites across the country. Assessment of the extent and condition of ecosystem assets through landscape metrics (e.g., fragmentation, connectivity, and mean patch size), which provide an accurate representation of the spatial arrangement by quantifying the extent and spatial characteristics at patch, classes or landscape level at different times.
The use of earth observation technologies to support official statistics remains a challenge, although methods based on earth observation for land use and land cover analysis are robust and scientifically recognized. Official statistics require that the data and documentation processes be compliant with the Quality Assurance Framework (QAF) [62] to ensure the credibility of the statistical agency. An assessment of products derived from remote sensing and earth observation data should consider the six dimensions of quality specified by the QAF including how they provide or support:
Relevant information, through use of appropriate approaches (i.e., methodology, date and resolution), adapted to the purpose of the study. Accurate results, through documentation of errors and limitations. Improved timeliness of data, with the increased availability of satellite EO information (e.g., RADARSAT Constellation Mission (RCM) [63]). Accessible EO data that is openly available through portals and platforms (e.g., UN platform and Google Earth Engine), which also provide access to open source software and scripts, in addition to big data computing capacity for implementation at national and global scale. Improved coherence and comparability, through access to complete coverage of analysis ready data (ARD) ensuring the use of calibrated data at national scale. Interpretability provided by metadata standards.
The use of satellite earth observation approaches in accordance with these dimensions, will lead to lower costs and better services, contributing to the UN SEEA-EEA, assessments of climate change resilience, and other Statistics Canada programs related to urban areas.
Footnotes
Acknowledgments
This study was supported by Statistics Canada and the Canadian Space Agency. The authors would like to thank Dr. Ridha Touzi for advice on Radarsat polarimetric parameters; Katherine Strong for her support in mapping; Mark Henry for advice on ecosystem accounting and Dave Craig for his GIS analysis. The authors are grateful to the reviewers for their thoughtful comments which contributed to the improvement of this paper.
