Urban Landscape Change Analysis Using Local Climate Zones and Object-Based Classiﬁcation in the Salt Lake Metro Region, Utah, USA

: Urban areas globally are vulnerable to warming climate trends exacerbated by their growing populations and heat island e ﬀ ects. The Local Climate Zone (LCZ) typology has become a popular framework for characterizing urban microclimates in di ﬀ erent regions using various classiﬁcation methods, including a widely adopted pixel-based protocol by the World Urban Database and Access Portal Tools (WUDAPT) Project. However, few studies to date have explored the potential of object-based image analysis (OBIA) to facilitate classiﬁcation of LCZs given their inherent complexity, and few studies have further used the LCZ framework to analyze land cover changes in urban areas over time. This study classiﬁed LCZs in the Salt Lake Metro Region, Utah, USA for 1993 and 2017 using a supervised object-based analysis of Landsat satellite imagery and assessed their change during this time frame. The overall accuracy, measured for the most recent classiﬁcation period (2017), was equal to 64% across 12 LCZs, with most of the error resulting from similarities among highly developed LCZs and non-developed classes with sparse or low-stature vegetation. The observed 1993–2017 changes in LCZs indicated a regional tendency towards primarily suburban, open low-rise development, and large low-rise and paved classes. However, despite the potential for local cooling with landscape transitions likely to increase vegetation cover and irrigation compared to pre-development conditions, summer averages of Landsat-derived top-of-atmosphere brightness temperatures showed a pronounced warming between 1992–1994 and 2016–2018 across the study region, with a 0.1–2.9 ◦ C increase among individual LCZs. Our results indicate that future applications of LCZs towards urban change analyses should develop a stronger understanding of LCZ microclimate sensitivity to changes in size and conﬁguration of urban neighborhoods and regions. Furthermore, while OBIA is promising for capturing the heterogeneous and multi-scale nature of LCZs, its applications could be strengthened by adopting more generalizable approaches for LCZ-relevant segmentation and validation, and by incorporating active remote sensing data to account for the 3D complexity of urban areas.


Introduction
Worldwide expansion and warming of urban regions present an important concern for the well-being of their residents, creating an urgent need for cost-effective strategies to track changing microclimates and inform current and future planning [1][2][3][4]. Previous research has extensively focused on urban heat island (UHI) phenomena manifested in significantly higher temperatures of cities compared to their surroundings due to unique thermal properties of urban land cover/use (LULC), building materials, human activities and other factors [5][6][7][8][9]. Increases in both ambient urban Pronounced urban development makes the Salt Lake Valley a particularly relevant place for LCZ analysis, as it is one of the fast-growing metro regions, by population, in the United States. Within the study period of 1993 to 2017, the population increased by approximately 323,000 (40%) [70], while the 2017 average annual percent change in population growth for Salt Lake County in 2017 was 1.27, compared to 0.72 for the entire USA [70]. Furthermore, the proportion of population increase due to migration is growing, which can be attributed to Utah's booming tech economy [71]. In January 2017 the Salt Lake Metro region had the lowest unemployment rate among all metro areas in the United States with greater than 1 million people [72]. As a result, the urban built up area has also expanded rapidly over the study period, primarily into the southwestern part of the region constrained by the Oquirrh Mountains bordering the West of the Valley. This has caused major changes in land uses and urban morphology within the existing urban footprint, thus providing an informative case for investigating changes in the extent and spatial pattern of LCZs. For this analysis, the study area boundary (Figure 1) was selected to capture the built-up area within Salt Lake County Utah based on 2010 US Census block group polygon geometries. To focus on the main urban extent, polygons that primarily overlaid the desired built up area were selected and combined to form the final boundary, whereas some large or irregularly shaped polygons that extended far outside the desired study area were omitted. Subsequently, object-based image classification was used to classify LCZs from the input Landsat satellite imagery separately for 2017 and 1993, followed by a post-classification change analysis, as described in the following sections ( Figure 2). Pronounced urban development makes the Salt Lake Valley a particularly relevant place for LCZ analysis, as it is one of the fast-growing metro regions, by population, in the United States. Within the study period of 1993 to 2017, the population increased by approximately 323,000 (40%) [70], while the 2017 average annual percent change in population growth for Salt Lake County in 2017 was 1.27, compared to 0.72 for the entire USA [70]. Furthermore, the proportion of population increase due to migration is growing, which can be attributed to Utah's booming tech economy [71]. In January 2017 the Salt Lake Metro region had the lowest unemployment rate among all metro areas in the United States with greater than 1 million people [72]. As a result, the urban built up area has also expanded rapidly over the study period, primarily into the southwestern part of the region constrained by the Oquirrh Mountains bordering the West of the Valley. This has caused major changes in land uses and urban morphology within the existing urban footprint, thus providing an informative case for investigating changes in the extent and spatial pattern of LCZs. For this analysis, the study area boundary (Figure 1) was selected to capture the built-up area within Salt Lake County Utah based on 2010 US Census block group polygon geometries. To focus on the main urban extent, polygons that primarily overlaid the desired built up area were selected and combined to form the final boundary, Remote Sens. 2019, 11, 1615 6 of 27 whereas some large or irregularly shaped polygons that extended far outside the desired study area were omitted. Subsequently, object-based image classification was used to classify LCZs from the input Landsat satellite imagery separately for 2017 and 1993, followed by a post-classification change analysis, as described in the following sections ( Figure 2).

Remote sensing data and image segmentation
The first step in each year's LCZ classification was generation of primitive mapping units (i.e., objects) via segmentation of a single-date Landsat image using the Segment Mean Shift tool in ArcGIS version 10.5 (ESRI Inc.) software. Such an algorithm can also be implemented in open-source GIS software, such as QGIS. Segmentation was performed using one satellite image for each time period of interest, representative of the summer-early fall time frame when atmospheric temperatures tend to be higher and annual herbaceous vegetation is typically senescent, highlighting their contrast with woody species (Table 1). For 1993, the input image was a Landsat 5 Thematic Mapper, Level 1 Precision Terrain Corrected, Tier 1 product, collected on 22 September 1993 (this date was chosen as the closest to the timing of the available digital orthophoto quad high spatial resolution reference imagery of 10 September 1993). Bands 2 (0.52-0.60 µm), 3 (0.63-0.69 µm), and 4 (0.76-0.90 µm) from this image were used in segmentation. For 2017, the input image was a Landsat 8 OLI/TIRS, Level 1 Precision Terrain Corrected, Tier 1 product, collected on 22 July 2017. Landsat bands 3 (0.533-0.590 µm), 4 (0.636-0.673 µm), and 5 (0.851-0.879 µm) from this image were used in segmentation. Electromagnetic regions of these bands correspond to the green, red, and near-infrared spectral wavelengths, expected to help differentiate between green and senescent vegetation, water and built-up/impervious surfaces. Here and later during classification, Landsat images were used at their original spatial resolution of 30 m for visible and infrared bands, and 120 m spatial resolution of thermal bands resampled to 30 m. All satellite images in this study were obtained from the U.S. Geological Survey's Earth Explorer database (https://earthexplorer.usgs.gov).
The Segment Mean Shift tool in ArcGIS software generates objects via a bottom-up, region-growing algorithm which groups the neighboring pixels together depending on the criteria for spectral similarity, spatial properties of resulting regions and minimum object size determined by a set of user-specified parameters (Table A1, Appendix; [73,74]). The parameter that controls the importance of spectral properties was set to the maximum value of 20 (on the scale from 0 to 20) because it was intended to achieve the greatest discrimination between features with similar spectral characteristics (such as vegetated classes). The parameter for relative importance of spatial detail (Table A1, Appendix) was informed by testing multiple values and visually examining the output against base map imagery, which suggested that the relatively low value of 2 (on the scale from 0 to 20) was optimal for preserving landscape element mixtures relevant to LCZ definitions.

Remote Sensing Data and Image Segmentation
The first step in each year's LCZ classification was generation of primitive mapping units (i.e., objects) via segmentation of a single-date Landsat image using the Segment Mean Shift tool in ArcGIS version 10.5 (ESRI Inc.) software. Such an algorithm can also be implemented in open-source GIS software, such as QGIS. Segmentation was performed using one satellite image for each time period of interest, representative of the summer-early fall time frame when atmospheric temperatures tend to be higher and annual herbaceous vegetation is typically senescent, highlighting their contrast with woody species (Table 1). For 1993, the input image was a Landsat 5 Thematic Mapper, Level 1 Precision Terrain Corrected, Tier 1 product, collected on 22 September 1993 (this date was chosen as the closest to the timing of the available digital orthophoto quad high spatial resolution reference imagery of 10 September 1993). Bands 2 (0.52-0.60 µm), 3 (0.63-0.69 µm), and 4 (0.76-0.90 µm) from this image were used in segmentation. For 2017, the input image was a Landsat 8 OLI/TIRS, Level 1 Precision Terrain Corrected, Tier 1 product, collected on 22 July 2017. Landsat bands 3 (0.533-0.590 µm), 4 (0.636-0.673 µm), and 5 (0.851-0.879 µm) from this image were used in segmentation. Electromagnetic regions of these bands correspond to the green, red, and near-infrared spectral wavelengths, expected to help differentiate between green and senescent vegetation, water and built-up/impervious surfaces. Here and later during classification, Landsat images were used at their original spatial resolution of 30 m for visible and infrared bands, and 120 m spatial resolution of thermal bands resampled to 30 m. All satellite images in this study were obtained from the U.S. Geological Survey's Earth Explorer database (https://earthexplorer.usgs.gov).
The Segment Mean Shift tool in ArcGIS software generates objects via a bottom-up, region-growing algorithm which groups the neighboring pixels together depending on the criteria for spectral similarity, spatial properties of resulting regions and minimum object size determined by a set of user-specified parameters (Table A1, Appendix A; [73,74]). The parameter that controls the importance of spectral properties was set to the maximum value of 20 (on the scale from 0 to 20) because it was intended to achieve the greatest discrimination between features with similar spectral characteristics (such as Remote Sens. 2019, 11, 1615 7 of 27 vegetated classes). The parameter for relative importance of spatial detail (Table A1, Appendix A) was informed by testing multiple values and visually examining the output against base map imagery, which suggested that the relatively low value of 2 (on the scale from 0 to 20) was optimal for preserving landscape element mixtures relevant to LCZ definitions. Finally, the parameter controlling segmentation affects the amount of spectral variation allowed within segments and thus minimum segment size (the same amount of local variation in more homogenous landscapes would be captured by larger segments than in heterogeneous areas). Due to potential variability in the size of LCZs, the goal of segmentation was not to delineate their complete extents, but instead to obtain object primitives that could be classified to reproduce such full extents. Thus, segmentation outcomes would ideally sufficiently capture minimum representative groupings of real-world objects and covers for different local climate zone types. To determine the appropriate minimum segment size given these objectives, 30 iterations of the segment mean shift tool were run on the input images, while varying the parameter on each iteration, from 2 to 60 by increments of 2. For each resulting segmented map, the spectral variation within each segment was summarized as the standard deviation of the near infrared band of the input image, and the mean of these within-segment standard deviations was used to summarize within-segment variation for the entire map, and to compare the output maps to one another. Scale parameter value of 20 was chosen for primitive object generation as the inflection point after which the rate of change in within-segment variation starts slowing towards saturation ( Figure 3).  Finally, the parameter controlling segmentation affects the amount of spectral variation allowed within segments and thus minimum segment size (the same amount of local variation in more homogenous landscapes would be captured by larger segments than in heterogeneous areas). Due to potential variability in the size of LCZs, the goal of segmentation was not to delineate their complete extents, but instead to obtain object primitives that could be classified to reproduce such full extents. Thus, segmentation outcomes would ideally sufficiently capture minimum representative groupings of real-world objects and covers for different local climate zone types. To determine the appropriate minimum segment size given these objectives, 30 iterations of the segment mean shift tool were run on the input images, while varying the parameter on each iteration, from 2 to 60 by increments of 2. For each resulting segmented map, the spectral variation within each segment was summarized as the standard deviation of the near infrared band of the input image, and the mean of these within-segment standard deviations was used to summarize within-segment variation for the entire map, and to compare the output maps to one another. Scale parameter value of 20 was chosen for primitive object generation as the inflection point after which the rate of change in within-segment variation starts slowing towards saturation ( Figure 3).

Image Classification
Next, the resulting segment polygons were classified into LCZs using the following process. For each time period, four cloud-free Landsat scenes (one for each season) were used (Table 1)

Image Classification
Next, the resulting segment polygons were classified into LCZs using the following process. For each time period, four cloud-free Landsat scenes (one for each season) were used (Table 1) in order to capture phenological variation in LCZs that may aid in their classification . For each primitive object in the segmentation results, basic summary statistics (minimum, maximum, mean, median and standard deviation) were calculated for each spectral band for each date. For 2017, training samples of LCZ classes (Table 2) were then selected following their definitions in [6] from ArcGIS and Google Earth high resolution imagery for 2017 (see examples of class patterns in Figure A1, Appendix A). A total of 60 segments were selected to represent each class using a random sampling technique where 100 segments within the study area were selected randomly and each was manually classified by visually interpreting their identity from reference imagery. This process was then repeated until each local climate zone class had accumulated a minimum of 60 samples. However, some classes were not prevalent enough within the study area to achieve the desired 60 samples using this methodology, so specific instances of these classes were sought out, classified manually and isolated prior to automated classification. Specifically, for LCZs 1 (compact high-rise) and 10 (heavy industry), all segments within the study area had been sufficiently identified, and therefore did not need to be included in the supervised image classification. It was discovered that classes 4 (open high-rise) and 7 (lightweight low-rise) could not be identified in the study area, while class C (brush, scrub) could not be reliably distinguished from low plants via visual interpretation alone. Therefore, the resulting candidate LCZs for classification included 2, 3, 5, 6, 8, 9, A, B, D, E, F, and G. For each of them, 50% of the previously identified samples were randomly selected for use in classifier training, and the remaining 50% were reserved for independent accuracy assessment. For 1993, training sample selection followed a slightly different process since there was no high-resolution reference imagery in Google Earth ( Table 2). The best historical (1993) data available for reference were black and white U.S. Geological Survey digital orthophoto quadrangles at 1-m resolution from September 10, 1993. Coverage of the study area by this imagery was also limited; therefore, the number of training samples for each class varied ( Table 2). The polygon data representing training sample segments was converted to csv file format for the subsequent LCZ classification.
Next, the open-source Waikato Environment for Knowledge Analysis (Weka) software [75] was used for the supervised classification of each Landsat image set into LCZs based on the identified training samples. The advantage of machine-learning methods for LCZ mapping is that, in contrast to maximum likelihood approaches, they often do not rely on rigid assumptions about class statistical distributions and can handle large numbers of discriminating features that may be useful for complex typologies. In our workflow, the Random Forest classifier recommended by WUDAPT was applied in both years. It was implemented using Weka's graphical user interface (GUI) with the specification as shown in Table A2, Appendix A. For details on how to run Weka the reader is referred to the supporting documentation at https://www.cs.waikato.ac.nz/ml/weka/documentation.html and to [75]. To maintain Remote Sens. 2019, 11, 1615 9 of 27 comparability with other studies using the WUDAPT method, a similar feature set was used, but with the addition of the Normalized Difference Vegetation Index (NDVI; normalized difference between near-infrared and red bands, [76]) and the Normalized Difference Water Index (NDWI; normalized difference between green and near-infrared bands, [77]), following previous studies [50,78]. For the 2017 period, the Landsat 8 feature set included the segment-level means of bands 2-7, 10, 11, NDVI and NDWI, while the 1993 Landsat 5 feature set included the segment-level means of bands 1-7, NDVI and NDWI.
Independent accuracy assessment was subsequently conducted for 2017 using the previously designated validation samples. Accuracy was assessed with the standard metrics (i.e., overall, user's and producer's accuracies) based on the contingency matrix comparing reference and mapped identities of the classified LCZs for validation samples [79]. In addition, a bootstrapping procedure was applied using the Python package, Scikit-learn [80] to evaluate the distributions of overall classification accuracy, kappa, overall accuracy of urban-only LCZs and overall accuracy of two aggregated classes, urban (LCZs 1-10) and natural (LCZs A-D, F and G excluding E common to both urban and natural areas), following [30]. No independent accuracy assessment was conducted for 1993 due to the lack of reference data, but the same method and similar input features (Table 1) were used here as for 2017.

LCZ Change Analysis
A post-classification change analysis was conducted by computing the geometric intersection of the 1993 and 2017 LCZ maps. Areas where different LCZs intersected were categorized as areas of change, and the specific type of change was labeled according to the previous (From) and current (To) LCZ classes. The sum area of each LCZ was calculated for each analysis year, and the change in LCZ proportional composition of the study area was analyzed. The sum areas of each specific 'from-to' change type were also analyzed for the original individual LCZs and for specific LCZ groups sharing similarity in both definitions and directions of change.
Finally, we estimated changes in two landscape metrics characterizing spatial pattern of individual LCZs: mean patch size, representing the average extent of LCZ patches in each respective period, and mean proximity index characterizing the degree of isolation among patches given their size [81]. The latter index can distinguish among scattered, smaller patches (low mean proximity index values) and clusters of larger, more aggregated patches (high mean proximity value).

Changes in the Proxy of Surface Temperature
To investigate the potential microclimatic implications of the observed LCZ changes, we used the land surface temperature (LST) datasets from the Landsat Level-2 Provisional Surface Temperature (ST) Science Products [82] for the same dates as the Landsat images used in classification (Table 1). Between two periods, these selected Landsat dates represented similar seasons and vegetation phenology stages (Table 1). Thus, to compare them, we computed per-pixel annual LST averages and then quantified changes in their mean values for each LCZ between two periods. To assess potential changes in LST distribution within the study area, we also normalized the annual average LSTs within its boundary by subtracting their mean values and dividing by standard deviation. We then evaluated how the skewness and kurtosis of the normalized value distributions changed between two periods.

Supervised Classification of LCZs
The overall accuracy of the 2017 LCZ classification was 64% with kappa coefficient of 0.61 (Table 3), while the cross-validation accuracy initially performed by Weka using the training set was 71% (with kappa equal to 0.68). Overall accuracy distribution determined by bootstrapping largely varied between these two estimates ( Figure A2 Appendix A), showing a mean value of 0.67. Furthermore, similar to the results in [30], bootstrapping produced relatively lower accuracy values for urban-only LCZs, with the minimum values reaching~0.55 ( Figure A2 Appendix A) and a mean of 0.62. In turn, simplifying the scheme to natural versus urban LCZs showed a highly successful distinction of these groups with overall accuracies exceeding 91% ( Figure A2 Appendix A), also similar to the results in [30].  The individual class accuracies in the independent assessment (Table 3) indicated several notable successes and challenges. Open Low-Rise (LCZ 6), the most prevalent class in 2017 (Figure 4), had the highest producer's accuracy at 90%. It was followed by Bare Soil (Class F) with 78%, Sparsely Built LCZ 9 (73%) and Dense Vegetation (71%). The highest user's accuracy was observed for distinct, predominantly non-built-up classes (  Table 3. Contingency matrix for LCZ classification accuracy assessment.

Mapped LCZs
Reference LCZs User's Accuracy 2 The individual class accuracies in the independent assessment (Table 3) indicated several notable successes and challenges. Open Low-Rise (LCZ 6), the most prevalent class in 2017 ( Figure  4), had the highest producer's accuracy at 90%. It was followed by Bare Soil (Class F) with 78%, Sparsely Built LCZ 9 (73%) and Dense Vegetation (71%). The highest user's accuracy was observed for distinct, predominantly non-built-up classes (Table 3): Dense Vegetation (79%), Bare Soil or Sand (78%) and Water (78%).  In contrast, a substantial amount of the classification error was associated with LCZs representing a greater intensity of development, as suggested by the off-diagonal counts in Table 3. For instance, Compact Mid-Rise (LCZ 2) and Open Mid-Rise (LCZ 5) both had the lowest producer and the lowest user accuracies of all classes (Table 3). They also showed a substantial mutual confusion with each other (Table 3). Other notable class pairs with mutual confusion were Paved (LCZ E) and Large Low-Rise (LCZ 8), likely due to their common association with flat impervious surfaces, and between Scattered Trees (LCZ B) and Low Plants (LCZ D), likely due to high prevalence of low-height vegetation in each. User's accuracy for the Open Low-Rise (LCZ 6) was relatively lower (73%), which together with its higher producer's accuracy suggested that other classes were likely to be misclassified as this LCZ (Table 3).

Change in LCZ Distributions
Spatial distribution and prevalence of LCZs in 1993 and 2017 (Figures 4 and 5) indicated a substantial degree of urbanization during this period ( In contrast, a substantial amount of the classification error was associated with LCZs representing a greater intensity of development, as suggested by the off-diagonal counts in Table 3. For instance, Compact Mid-Rise (LCZ 2) and Open Mid-Rise (LCZ 5) both had the lowest producer and the lowest user accuracies of all classes (Table 3). They also showed a substantial mutual confusion with each other (Table 3). Other notable class pairs with mutual confusion were Paved (LCZ E) and Large Low-Rise (LCZ 8), likely due to their common association with flat impervious surfaces, and between Scattered Trees (LCZ B) and Low Plants (LCZ D), likely due to high prevalence of low-height vegetation in each. User's accuracy for the Open Low-Rise (LCZ 6) was relatively lower (73%), which together with its higher producer's accuracy suggested that other classes were likely to be misclassified as this LCZ (Table 3).

Change in LCZ Distributions
Spatial distribution and prevalence of LCZs in 1993 and 2017 (Figures 4 and 5) indicated a substantial degree of urbanization during this period (Table A3 Appendix). Open Low-Rise (LCZ 6) increased from being the second largest (17.4% of the study area) to the most prevalent (35%) in 2017. This class showed the largest gains overall, adding 126.3 km 2 (88.5% of its 1993 extent), equivalent to 16.8% of the study area (Table A4 Appendix). This was consistent with the population growth observed during the study period, as LCZ 6 typically corresponded with residential suburban development. Other substantial gains in area were observed in the Paved (LCZ E) and Large Low-Rise (LCZ 8) classes, adding, respectively, 160.6% and 82.35% of their 1993 extents ( Figure 5, Table A4 Appendix). By their definition and distribution within the study area, these classes both typically corresponded with commercial and manufacturing uses, the growth of which can also be visually identified on the Local Climate Zone maps. In contrast, the two other most extensive 1993 LCZs, Low Plants (LCZ D) and Sparsely Built (LCZ 9), decreased by 37.9% and 46.1%, respectively, even though both remained among the 3 largest classes in 2017 (Table A4 Appendix). Area losses for these classes together with a 49% decrease in Bare Soil and Sand (LCZ F) indicates that all three likely represented preferred sites for new development. Among most classes that already had been developed in 1993, losses were relatively small, with the exception of Open Mid-Rise (LCZ 5), which Other substantial gains in area were observed in the Paved (LCZ E) and Large Low-Rise (LCZ 8) classes, adding, respectively, 160.6% and 82.35% of their 1993 extents ( Figure 5, Table A4 Appendix A). By their definition and distribution within the study area, these classes both typically corresponded with commercial and manufacturing uses, the growth of which can also be visually identified on the Local Climate Zone maps. In contrast, the two other most extensive 1993 LCZs, Low Plants (LCZ D) and Sparsely Built (LCZ 9), decreased by 37.9% and 46.1%, respectively, even though both remained among the 3 largest classes in 2017 (Table A4 Appendix A). Area losses for these classes together with a 49% decrease in Bare Soil and Sand (LCZ F) indicates that all three likely represented preferred sites for new development. Among most classes that already had been developed in 1993, losses were relatively small, with the exception of Open Mid-Rise (LCZ 5), which lost 48.3% of its area, and Compact Low-Rise (LCZ 3), a large part of which was converted to Open Low-Rise (LCZ 6). Both of the latter results were not intuitive, and in addition to potential LCZ transitions, they could be caused by the classification error and potential expansion of vegetation cover. For instance, growth of woody species or introduction of specific landscaping practices could affect the balance between vegetated and built-up spectral contributions, making the respective locations more similar to LCZ 6.

Changes in Spatial Pattern of LCZ Distribution
The spatial pattern of LCZ distribution also showed several notable changes during the study period. Mean proximity index increased substantially for LCZs 6, 8, and E (Figure 6a), resonating with the large net increases in their respective class areas (Figure 6b, Table A4 Appendix A). Among these, LCZ 6 in particular showed an extreme increase, indicating that its expansion was accompanied by greater local aggregation, likely due to extension of existing patches (as opposed to new leapfrog development). While some relatively isolated patches of LCZ 6 remained in the southwestern part of the study area (Figure 4), its mean patch size nearly doubled over the study period (Figure 6b).
Low-Rise (LCZ 6). Both of the latter results were not intuitive, and in addition to potential LCZ transitions, they could be caused by the classification error and potential expansion of vegetation cover. For instance, growth of woody species or introduction of specific landscaping practices could affect the balance between vegetated and built-up spectral contributions, making the respective locations more similar to LCZ 6.

Changes in Spatial Pattern of LCZ Distribution
The spatial pattern of LCZ distribution also showed several notable changes during the study period. Mean proximity index increased substantially for LCZs 6, 8, and E (Figure 6a), resonating with the large net increases in their respective class areas (Figure 6b, Table A4 Appendix). Among these, LCZ 6 in particular showed an extreme increase, indicating that its expansion was accompanied by greater local aggregation, likely due to extension of existing patches (as opposed to new leapfrog development). While some relatively isolated patches of LCZ 6 remained in the southwestern part of the study area (Figure 4), its mean patch size nearly doubled over the study period (Figure 6b).
In contrast, LCZs 9, D, and F, the classes that showed substantial net decreases in class area, all decreased in mean proximity index as well. On the 1993 map, LCZ 9 formed an aggregated and relatively large grouping of patches in the South and Southwest portions of the study area. This region is where much landscape change took place. Other LCZs, primarily LCZ 6, expanded into this region, forming a more heterogeneous patchwork of LCZs. By looking at the LCZ maps (Figure 4), LCZ 9 clearly shows change to a more fragmented pattern of smaller patches. The mean patch size of LCZ 9 also decreased by approximately 45%.

a.
b.

Changes in thermal characteristics of the urban landscape
Finally, the comparison of Landsat brightness temperature variation across the study area between 1993 and 2017 revealed a pronounced overall warming ( Figure A3 Appendix), and increases from 0.1 to 2.9 °C on average temperature of individual LCZs (Figure 7a). Pixels with net decreases in temperature represented ~15% of the study area, corresponding to small wetland features, vegetated areas, golf courses and a few urban shadow areas in high-rise neighborhoods. In 1993, the highest mean temperatures were observed for LCZ 8 (25.8 °C) followed by LCZs 2, 10 and In contrast, LCZs 9, D, and F, the classes that showed substantial net decreases in class area, all decreased in mean proximity index as well. On the 1993 map, LCZ 9 formed an aggregated and relatively large grouping of patches in the South and Southwest portions of the study area. This region is where much landscape change took place. Other LCZs, primarily LCZ 6, expanded into this region, forming a more heterogeneous patchwork of LCZs. By looking at the LCZ maps (Figure 4), LCZ 9 clearly shows change to a more fragmented pattern of smaller patches. The mean patch size of LCZ 9 also decreased by approximately 45%.

Changes in Thermal Characteristics of the Urban Landscape
Finally, the comparison of Landsat land surface temperature variation across the study area between 1993 and 2017 revealed a pronounced overall warming ( Figure A3 Appendix A), and increases from 0.1 to 2.9 • C on average temperature of individual LCZs (Figure 7a). Pixels with net decreases in temperature represented~15% of the study area, corresponding to small wetland features, vegetated areas, golf courses and a few urban shadow areas in high-rise neighborhoods. In 1993, the highest mean temperatures were observed for LCZ 8 (25.8 • C) followed by LCZs 2, 10 and E (>24 • C in each). In 2017, LCZs 8 and 10 remained the warmest on average (27.7 and 26.7 • C, respectively), followed by 2, 3, E and 10 (26.1-26.6 • C). The lowest temperatures, in turn, were observed for types A and B in both 1993 (20.4 and 22.7 • C, respectively) and in 2017 (21.5 and 22.8 • C, respectively). The highest increases in mean LST (>2 • C) were observed for open LCZs 5 and 6, compact low-rise (LCZ 3) and industrial (LCZ 10). In turn, the lowest LST changes (within 0.1-1.2 • C) corresponded to natural LCZs A-F, indicating that the background thermal conditions did not differ substantially between the two periods.  Mean values of normalized annual average LST indicated that most LCZs remained consistently above or below the study area mean values (Figure 7b). Among natural classes, vegetated (A, B) and water (G) were substantially below the average, while most of the built-up classes were above the average, especially in 2017 (Figure 7b). Two exceptions were Sparsely Built (LCZ 9) with substantial prevalence of non-built cover, and Compact High-Rise (LCZ 1) where lower-than-average LSTs could emerge due to tall building shadows at the time of the satellite overpass. Across the study area, the skewness of normalized LST distributions changed from −0.88 in 1993 to −0.95 in 2017, indicating the increase in the left tail and greater frequency of higher values in the later time period. In turn, the kurtosis of normalized LST distributions has also increased from 7.03 in 1993 to 10.05 in 2017, indicating a greater concentration of values around the mean-i.e., increased homogeneity of the landscape, which may be expected with urbanization.

The Informative Value of LCZs for Monitoring Urban Transformation
Major environmental transformations experienced by globally expanding urban regions call for a revision of traditional LULC classes to facilitate the understanding of potential processes underlying these transitions [21,31,53]. The local climate zone framework has been proposed as one of such novel, alternative typologies bridging physical characteristics of cities (and particularly their built environment) with spatial and morphological structural properties applicable to urban architecture, planning and design [9,17,21]. As the appeal of LCZs is recognized by a wider community of researchers and planners and more cities are included in the WUDAPT efforts [20,42] Mean values of normalized annual average LST indicated that most LCZs remained consistently above or below the study area mean values (Figure 7b). Among natural classes, vegetated (A, B) and water (G) were substantially below the average, while most of the built-up classes were above the average, especially in 2017 (Figure 7b). Two exceptions were Sparsely Built (LCZ 9) with substantial prevalence of non-built cover, and Compact High-Rise (LCZ 1) where lower-than-average LSTs could emerge due to tall building shadows at the time of the satellite overpass. Across the study area, the skewness of normalized LST distributions changed from −0.88 in 1993 to −0.95 in 2017, indicating the increase in the left tail and greater frequency of higher values in the later time period. In turn, the kurtosis of normalized LST distributions has also increased from 7.03 in 1993 to 10.05 in 2017, indicating a greater concentration of values around the mean-i.e., increased homogeneity of the landscape, which may be expected with urbanization.

The Informative Value of LCZs for Monitoring Urban Transformation
Major environmental transformations experienced by globally expanding urban regions call for a revision of traditional LULC classes to facilitate the understanding of potential processes underlying these transitions [21,31,53]. The local climate zone framework has been proposed as one of such novel, alternative typologies bridging physical characteristics of cities (and particularly their built environment) with spatial and morphological structural properties applicable to urban architecture, planning and design [9,17,21]. As the appeal of LCZs is recognized by a wider community of researchers and planners and more cities are included in the WUDAPT efforts [20,42], the need to better understand the opportunities offered by different classification approaches (and possibilities to improve them with new methodological advances) is becoming even more pertinent.
A closer look at the transitions among individual LCZs in the study area ( Figure 5) highlights the regional tendency towards suburbanization at the expense of lateral low-density development. For example, the largest net increase of the low-density built-up areas (LCZ 6) and relatively low (17%) conversion of their 1993 area to any other class (Table A4 Appendix A) suggest that it is not an intermediate stage of development, but rather an ultimate stage. Our results further suggest that much of the new development associated with commercial/industrial land uses (LCZs 8 and E) emerged primarily from direct open space (Figures 4 and 5, Table A3 Appendix A) and in some cases sparsely built (LCZ 9) rather than existing low-density development. These patterns are consistent with the broader tendencies in the rapidly urbanizing Wasatch Front of north-central Utah encompassing our study area [83,84]. Developed areas here are projected to expand by 48-80% by 2030, while expected increases in population density range only between 15-25 people/ha, suggesting a prevalence of low-density peripheral urban growth [83]. Similar to urbanizing regions in other parts of USA and Canada [85], such lateral expansion occurs at the expense of agriculturally productive lands and poses multiple threats to natural habitats and ecologically sensitive areas in the region [83,84].
Our results further show that LCZs provide a useful framework to characterize urban transformation [36] and elucidate microclimate-relevant aspects of urban heterogeneity [27,31,50]. Differences in spatial extents and patterns of LCZs between 1993 and 2017 ( Figure 4) highlighted the increases in relatively warmer LCZs 8, E and 6 and losses of cooler types A and 9 ( Figure 5, Table A4 Appendix A). Relative differences in mean brightness temperature among these LCZs (Figure 7) were consistent with earlier studies [29,37,45], although in our case type D (low plants) appeared among the warmer types likely due to dry and warm state of herbaceous natural vegetation in the dry-summer climate of the study area. The overall character of LCZ change showed large increases in LCZs 6, 8, and E, and decreases in LCZs 9, D, and F, which collectively indicate increases in landscape structural and morphological features contributing to potential heat island phenomena. These results were similar to an LCZ change analysis conducted in Bogor City, Jakarta [36], in which an overall trend of increasing area among the built LCZs and decreasing area among the vegetated LCZs was observed along with increases in land surface temperature.
Our results further suggest that potential microclimatic transitions indicated by LCZ changes could be amplified by shifts in their spatial configuration ( Figure 6), particularly, consolidation and expansion of LCZs 6, 8 and E, and reduction and fragmentation of LCZs 9, D and F. Hypothetically, not all development transitions should elevate urban heat, particularly when encroaching onto relatively hot open space. For example, Bare Soil (LCZ F) tends to be a relatively warm LCZ in arid and semi-arid cities [45], and thus loss of area and fragmentation of patches of LCZ F may be expected to produce some amount of local cooling if new development adds shading structures, tree and shrub vegetation or irrigation, which may be expected for suburban residential land management. In our case, 35% of the area lost by LCZ F over the study period was replaced by LCZ 6, a class expected to be relatively cooler due to open nature of development and greater vegetation cover (Figure 7), while~5.6 km 2 also converted from LCZ F to LCZ 9, another relatively cooler LCZ.
Notably, however, the overall shifts in brightness temperatures in our focal area do not indicate any substantial local cooling ( Figure A2 Appendix A) with the exception of a few irrigated golf course areas. Furthermore, the increase of average temperatures for all LCZs regardless of the amount of their change (Figure 7, Figure A3 Appendix A) suggests that the whole region has experienced an increase in surface temperature between the compared periods. The magnitudes of this increase should be interpreted with caution because they could be influenced by ambient climate warming and the differences between two Landsat instruments and short-term weather fluctuations despite the use of top-of-atmosphere brightness temperature products and their averages over 3-year summer periods. Nevertheless, the observed region-wide warming and shifts in normalized LST distributions were consistent with the expected potential warming due to the expansion of urban area, leading to local temperature increases even in neighborhoods not experiencing changes in urban morphology per se [13,86]. For instance, modeling of summer climate indices based on LCZs in Brno, Czech Republic [86] reported an increase in summer warming indicators for all LCZs under different climate change scenarios, especially the extreme ones.
Importantly, these results imply that microclimates may change without substantial changes in urban morphology defining LCZs, due to synergistic effects of broader-scale urban configuration and increasing ambient and surface temperatures. Changes in urban energy consumption and anthropogenic heat sources could also contribute to thermal shifts without an apparent change in urban morphology [19,87]. These considerations should be taken into account by spatial land use planning and policy to strategically plan for local cooling via landscape design measures to mitigate the broader warming introduced by urban development and promote more sustainable management of water and energy resources [16,19,[87][88][89][90]. As many LCZ analyses to date have been performed as single-year, the feedbacks between LCZ-specific microclimates and the broader urban environment represent an important direction for future work.

Successes and Challenges in Object-Based LCZ Classification
Our results indicate that, overall, object-based image analysis is appealing for classification of Local Climate Zones particularly as high spatial resolution imagery becomes incorporated in such efforts; however, the workflow does not necessarily outperform the successes of previous pixel-based analyses and needs more refinement to become a generalizable framework. Although the overall LCZ classification accuracy of 64% for the full range of considered classes was relatively low, it was comparable with several other previous studies performing supervised classifications of LCZs, such as 56% in [46], 64% in [35] and 67.7% in [42]. Furthermore, most of the classification error was a cumulative outcome of small misclassifications among individual class pairs (Table 3), while the sizable pool of 12 candidate classes and their inherent similarities contributed to greater likelihood of such misclassifications. The observed mutual confusions among LCZs were also common in other studies, particularly the prominent mutual confusion between Compact Mid-Rise (LCZ 2), Large Low-Rise (LCZ 8), and Paved (LCZ E). For instance, [21] noted the misclassification between Large Low-Rise (LCZ 8) and Paved (LCZ E), while [44] reported confusion between the compact built zones and Large Low-Rise (LCZ 8).
Our results further suggest that these confusions may result from both vertical (i.e., building height and rise) and horizontal (i.e., openness and spatial arrangement) similarities among LCZs (Table 3). Mutual confusions among built-up classes with different levels of building height indicate the limitations of passive remote sensors to capture the differences in rise, while misclassifications among compact and open LCZs suggest potential similarities in relative contributions of their land covers despite the expected differences in configuration. Similarly, substantial confusion among 'sparse' LCZs 9, B and D (Table 3) places varying sparseness of woody plant cover as an important challenge for both LCZ mapping and inference of microclimatic properties in less developed and sparsely built areas. This evidence implies that complexity of 3D urban structure is an important contributor to LCZ mapping uncertainties regardless of the classification method. Future work should thus more rigorously explore the utility of datasets representing vertical variability in urban surface, such as novel metrics of building density, tree cover and broader local heterogeneity developed from LiDAR [28] and radar [49] active remote sensing products.
Some of the observed classification successes were likely attributed to the advantages of the OBIA framework, particularly, the high producer's accuracy for a heterogeneous-by-definition Open Low-Rise (LCZ 6) ( Table 3). Using image objects as multi-pixel mapping units may allow capturing complex spatial patterns composed of both constructed elements and various types of open and vegetated spaces between them under the same semantic identity. As open low-rise development appears to be one of the key transformations in the study area's urban sprawl, these results suggests that the object-based framework may help address the definitional complexity of LCZs as an inherent challenge in their mapping. Our results also suggest that mapping LCZs with potentially larger amounts of vegetation (e.g., A, B, D, 9) can benefit from the OBIA approach as well by helping to separate them from more developed LCZs that may include vegetated pixels as part of, e.g., backyard and street spaces. However, primitive objects generated from moderate-resolution data might not be sufficiently sensitive to nuanced variation in low-density woody cover or presence of sparsely built elements.
Differentiating among such LCZs with confidence may ultimately require higher-resolution imagery [46,48,67], where, again, objects may offer greater flexibility for summarizing local heterogeneity and complex LCZ boundaries compared to square block mapping units [67]. For current and future analyses, such opportunities are emerging with the accumulation of open-access Sentinel-2 imagery and Harmonized Landsat-Sentinel-2 product [91]. Moderate-resolution imagery remains appealing for broad-scale and historical LCZ mapping due to the availability of thermal products and unparallelled archive length. However, emerging urban studies are increasingly uncovering the microclimatic impacts of fine-scale LULC patterns and urban structure [15,16,63,64], which necessarily rely on higher detail of landscape characterization. As such, mechanistic understanding of LCZs and their microclimatic properties would necessarily invoke the applications of high spatial resolution data [51], even if thermal information is not available at the same spatial scales. As the benefits of greater spatial detail come at the cost of greater local spectral noise and confounding heterogeneity, the advantages of OBIA [59,69] and its earlier applications in urban mapping [53,56] may be increasingly useful to LCZ analyses in the future.
The limitations encountered in our OBIA implementation for LCZ mapping call for more research in three important directions. First, OBIA performance, as well as its added cost and logistical demand, should be compared against the established pixel-based workflows, particularly WUDAPT [20,21,30]. The relatively straightforward implementation and universality of the latter are important advantages when urban structure can be well represented by 100-150 m pixel units, as suggested by high classification accuracy in some of the previous studies [21,43,44]. In contrast, OBIA's possibility to use more flexible and variable mapping units in complex landscapes comes at a cost of extra steps required to parameterize the segmentation and limited generalizability, which may not be equally necessary among different cities. Thus, to provide a more informed guidance about method selection, future research should perform comparative assessments of the alternative workflows in landscapes with different levels of structural and morphological complexity, using high-quality reference data.
Second, the role of segmentation and strategies to optimize the choice of its parameters need to be identified more explicitly, because the degree of spectral heterogeneity in LCZs is expected to vary due to their definitions [58,92]. Although the Segment Mean Shift tool is able to generate variable sizes and properties of primitive objects, the contrasts between more homogeneous LCZs and classes with inherent configurational complexity present the risks for both under-segmentation (i.e., covering more than one visually discernable LCZ) and over-generalization (delineating single homogeneous cover signatures, such as building rooftops, rather than capturing local scale structures). Furthermore, the flexibility of variable-size objects comes at a cost of difficulties to anticipate specific properties of mapping units, which may need to be chosen according to city-specific landscape properties. Thus, a more careful discrimination should be taken in future analyses to choose training objects that are relatively large and "homogeneous" with respect to LCZ definition. Incorporating spatial units that are likely to influence the pattern of development associated with hypothesized LCZ structure-such as land parcels, zoning units or census units sharing similar population and housing densities-may also facilitate this task.
Third, there remains an important uncertainty about which spectral, textural and other features may be most effective at differentiating LCZs given the complexity of their typology (Table 2). Phenological characteristics of vegetation may further contribute to such uncertainties; for example, in dry-summer climates such as in our study area, herbaceous vegetation may be spectrally similar to trees in green stage and to bare soil in senescent stage. More detailed future investigations would thus be useful both at the segmentation stage (since the Segment Mean Shift tool in ArcGIS takes as input 3 user-specified image bands at a time) and at the stage of classification to identify the most useful discriminating features and optimal date(s) [47].

Other Key Lessons and Future Research Directions
Opportunities and challenges revealed by our analysis suggest several other important future research directions to make LCZ applications more robust and transferable across global regions. There is an obvious need to refine and standardize the validation of LCZ mapping outcomes, because different accuracy assessment strategies in previous studies make it difficult to draw generalized lessons about method performance. Spatial dimensions of testing units may also affect accuracy assessment; for instance, a study with classification accuracy of 64% similar to ours [35] commented that even independent validation pixels may be problematic as LCZs are supposed to be at least 1 km 2 in size. To avoid a potential bias due to proximity of training and test samples, in [44], independent validation polygons were purposefully chosen from training areas instead of randomly selected pixels; the resulting overall accuracies ranged from 79.6 to 90.2% in different cities. Lessons from a global LCZ transferability assessment [50] highlight the importance of both training and test sample selections and suggest the need to consider sample sources broader than an individual focal city. A unique challenge in validating object-based LCZ classification thus becomes a decision as to what types of reference objects might perform as the most representative, and what criteria might be especially important in guiding the initial segmentation process to capture their LCZ relevance.
The latter challenge, in turn, highlights another interesting and important research need: developing a stronger understanding of the 'characteristic' spatial scales of different LCZs to better inform the minimum mapping unit selection (or segmentation scale in OBIA workflows). Both the original definitions of LCZs and the insights from our results make it obvious that a single segmentation scale or moving window size might not fully accommodate their differences, which has been long recognized in complex landscapes [55,92]. However, in the context of LCZ mapping, it further translates into the need to understand scales that effectively represent thermal outcomes of specific patterns. This understanding is currently limited both by challenges to discriminate among LCZs despite the inclusion of the thermal remote sensing datasets [93] and by the dynamic behavior of urban temperatures that may shift even when morphological aspects of LCZs remain stable. Future studies should therefore more explicitly investigate the fundamental relationships between local spatial patterns and their thermal outcomes across ranges of urban scales beyond individual pixels or multi-pixel square neighborhoods.
Together, these issues also raise important questions about the interpretation and uncertainty of the LCZ boundaries. Urban microclimates are often continuous, scale-dependent [17,20,22,39] and difficult to discretize due to various processes contributing to atmospheric circulation and energy transfer [6]. In contrast, the morphology of LCZs treats urban landscapes as mosaics of discrete semantic entities, e.g., buildings, trees, roads or their groupings, with detectable limits [6]. However, with moderate to coarse spatial resolution units, such as 30 m Landsat pixels or larger sizes recommended by WUDAPT [21,30,42], the boundaries between urban entities and LCZs themselves are likely to be diluted. As such, it can be difficult to interpret the mapped LCZ boundaries based on the local gradients of temperature or other microclimatic parameters. This uncertainty is further augmented by both the shorter-term variation of urban temperatures (e.g., diurnal, seasonal or interannual) and longer-term shifts (e.g., due to continuing urban warming). As a result, spatially explicit validation of LCZs and their boundaries from the urban microclimate perspective remains an important and interesting challenge for the future research.
Finally, an important research need is testing the potential of complementary remote sensing datasets to elucidate structural, configurational and microclimatic differences among LCZs. The 'ideal' datasets for characterizing urban 3D complexity, such as LiDAR point clouds or waveform datasets, are not yet available globally as systematically acquired and publicly available products. However, advances in moderate-resolution satellite radar systems, including publicly available Sentinel-1 instruments, may provide useful information about volume and complexity of urban structures to complement the passive optical remote sensing datasets even if thermal bands are not available [43]. Very few studies to date, however, have attempted such multi-source analyses in the context of LCZs, and the utility of radar inputs was not uniformly obvious among their results [43,93]. The earlier mentioned harmonization of Landsat and Sentinel-2 products [91] and emerging higher-resolution open-access land surface temperature datasets [94] would further support such analyses at global scales and repeated time frames.

Conclusions
This study classified Local Climate Zones [6] in the Salt Lake Metro Region (USA) and assessed their change between 1993 and 2017 using object-based image analysis with 30-m resolution Landsat data. Our results show that urbanization of the study area during this period was accompanied by several notable LCZ transformations, particularly the expansion of open low-rise development common in residential neighborhoods and large low-rise and paved LCZs common in commercial/industrial areas. Both of these types grew primarily at the expense of non-developed or sparsely built areas and suggested that suburbanization and lateral sprawl, rather than intensification of existing development, were the key players in the observed transformations. While hypothetically the observed transitions could be expected to produce both local warming and local cooling (e.g., due to conversion of bare soil to low-density developed areas with greater green vegetation cover), the mean values for Landsat-based surface temperature increased for all LCZs by 0.1-2.9 • C, showing the warming tendency across the vast majority of the region. While our results cannot reliably disentangle urbanization from ambient climate warming, new anthropogenic energy sources and potential confounding effects of satellite data availability, they highlight an important need to investigate thermal dynamics of LCZs which might occur in growing cities without pronounced changes in the morphology of their neighborhoods.
From a practical perspective, our findings indicate that OBIA applied to moderate-resolution satellite imagery shows potential to facilitate LCZ classification and change analysis by capturing local heterogeneity of urban landscapes within objects as mapping units. However, as most of the classification error resulted from the confusions among developed LCZs with varying horizontal density of landscape elements and vertical height of vegetation and built-structures, several important directions should be pursued in future research. These directions include: (1) identifying characteristic spatial scales of different LCZ patterns to more effectively guide segmentation of images into objects as mapping units; (2) testing additional feature sets as discriminating attributes in LCZ classification, particularly, measures of local variance and texture [44,47,66]; and (3) testing the utility of complementary datasets, particularly active remote sensing data, to support LCZ classification by providing additional information on 3D structure and heterogeneity of urban landscapes [43,93].
In sum, our findings concur with previous studies that Local Climate Zones show potential as a valuable framework for landscape change analysis as they can provide a more detailed, climate-relevant perspective of urban change when compared with traditional land use and land cover classification schemes. Our results corroborate the potential of the open-access Landsat imagery to support the analysis of urban microclimatic transformation consistently with the goals of the WUDAPT project and similar initiatives and encourage LCZ classifications for past time periods and more detailed longer-term change assessments. Figure A1. Examples of spatial patterns of different LCZs from high-resolution reference data.