Multiscale mapping of local climate zones in Tokyo using airborne LiDAR data, GIS vectors, and Sentinel-2 imagery

ABSTRACT

Multisource remote sensing and geographic information system (GIS) data have contributed powerfully to the large-scale automated mapping of local climate zones (LCZs). However, the accessibility of high-resolution height data, the applicability of standard thresholds to local contexts, and the dependence of mapping scales have limited LCZ classification studies. In this study, we combined airborne LiDAR data, Sentinel-2 imagery, and GIS vector (buildings and roads) data to develop a multiscale automated LCZ classification scheme in the 23 special wards of Tokyo. Based on the optimized thresholds of seven LCZ properties, GIS-based LCZ mapping was implemented using fuzzy logic classifiers at the block scale and at different grid-cell scales (100 m–1000 m). In addition to assessing accuracy using reference samples, multidate thermal infrared data (Landsat-8 and ASTER data) were used to understand the LCZ-LST (land surface temperature) relationship at multiple scales. The results showed that the overall accuracies of LCZs could be significantly increased by threshold optimization at all scales. Significant differences in LCZs and LSTs among different mapping units were observed. The highest overall accuracy was greater than 80% at the 100-m grid-cell scale. As the size of grid cells increased, the overall accuracy of LCZ classification decreased. For each LCZ, the mean daytime/nighttime LST exhibited more variation by date than by scale. This study provides a promising picture of GIS-based LCZ mapping and LCZ-LST relationships at multiple scales.

KEYWORDS:

• Local climate zones
• Urban form
• Multisource data fusion
• Fuzzy logic
• Threshold optimization
• Land surface temperature

1. Introduction

1.1. Concept of local climate zone (LCZ) scheme

More than half of the world’s population lives in urban areas, and such areas account for the majority of the world’s economic activity, energy use, and energy-related emissions (Seto et al. 2017; United Nations, Department of Economic and Social Affairs, Population Division 2019; Zhu et al. 2019). Along with the increasing amount of human activity that accompanies urban development (United Nations, Department of Economic and Social Affairs, Population Division 2021), the energy balances at and near the surface are altered in urban areas, typically resulting in higher air and surface temperatures in urban areas than in surrounding rural areas (i.e. the urban heat island phenomenon) (Voogt and Oke 2003; Oke et al. 2017; Masson et al. 2020).

Considerable previous studies have used conventional land use/cover classification maps to investigate the urban heat island effect (Zhao et al. 2017; Deilami, Kamruzzaman, and Liu 2018; Xian et al. 2021). However, not all classification systems use consistent or climate-related properties to define classes (Stewart and Oke 2012). As a result, a “local climate zone” (LCZ) classification system was developed by Stewart and Oke (2012) to study the thermal performances of different land cover types at the local scale to reflect the diversity of urban areas, replacing the traditional “urban‒rural dichotomy.” The standard LCZ scheme consists of 17 basic classes (i.e. built LCZ classes [LCZs 1‒10] and land cover LCZ classes [LCZs A – G]) based on surface structure, cover, material, and human activity differences. The LCZ framework provides the value ranges of 10 LCZ properties for each of the 17 standard LCZ classes, relating to the surface structure (e.g. the height of roughness elements, sky view factor, aspect ratio, and terrain roughness class), surface cover (e.g. the building surface fraction, impervious surface fraction, and pervious surface fraction), surface material (e.g. the surface admittance and surface albedo), and human activity (e.g. the anthropogenic heat output) (Stewart and Oke 2012; Stewart, Oke, and Krayenhoff 2014).

1.2. LCZ mapping approaches

The LCZ classification system has developed rapidly and been broadly applied in studies of urban morphology, urban climate, and urban sustainability (Bechtel et al. 2019b; Taubenböck et al. 2020; Eldesoky, Gil, and Pont 2021; Aslam and Rana 2022; Xia et al. 2022). High-quality LCZ classification maps are key basic data used to support urban sustainable development and can provide reference values for the improvement of urban form with respect to the thermal environment. The development of remote sensing (RS) and geographic information system (GIS) technologies has enabled automated LCZ mapping on a large scale. In the field of LCZ classification, the performance of classification methods varies depending on the regional context, data used, and quality of the samples. Mainstream LCZ mapping methods can generally be divided into the following three types: (1) RS-based methods, (2) GIS-based methods, and (3) combined methods. A good overview of recent research can be found in the referenced review papers (Jiang et al. 2021; Lehnert et al. 2021; Ma et al. 2021; Quan and Bansal 2021; Aslam and Rana 2022; Huang et al. 2023).

Among these LCZ mapping methods, RS-based methods have the potential for large-scale mapping, for example, on regional and global scales (Demuzere et al. 2019; Demuzere, Bechtel, and Mills 2019; Demuzere et al. 2020, 2022; Rosentreter, Hagensieker, and Waske 2020; Zhu et al. 2022). Most notably, the World Urban Database and Access Portal Tools (WUDAPT) project was dedicated to developing an LCZ map database and a map generator for major cities in the world using freely available RS data (Bechtel et al. 2015, 2019a; Ching et al. 2018; Demuzere, Kittner, and Bechtel 2021). RS-based methods, such as supervised classification methods, strongly rely on highly subjective training samples (Bechtel et al. 2017; Yoo et al. 2019) and RS data. Currently, multisource RS data used in RS-based methods are limited mostly to pixel-level stacking (e.g. La, Bagan, and Yamagata 2020; Chen et al. 2021a, 2021b). In addition, these data often still lack information about the heights of surface objects. Moreover, most machine learning (including deep learning) methods used for LCZ classification are still black-box models with little interpretability. Unlike RS-based methods, GIS-based methods usually allow for the generation of LCZ classification results with high accuracies, especially for built LCZ classes (Shi et al. 2018). GIS-based methods rely on the calculation of properties (also called parameters or indicators) that are closely related to LCZs and are interpretable; therefore, these methods have high requirements for the input datasets.

1.3. LCZ mapping units

LCZs characterize relatively homogeneous urban surfaces (10²–10⁴ m) that share similar energy budgets (Stewart and Oke 2012; Stewart, Oke, and Krayenhoff 2014). However, various studies have not been in agreement on the mapping units used. To date, the mapping units used for LCZ classification have usually been grid-cell units and parcel units. Grid cells, the most common mapping units, enable the user to perform integrated data analyses from multiple sources and to achieve a balance between the feasibility of computations and the need for details (Bagan and Yamagata 2012, 2014; Chen et al. 2021a). Parcel units, defined as spatial units bounded by line elements, can usually be presented as plots/lots, blocks/islets, triangulated irregular networks (TINs) and skeletons (Rodler and Leduc 2019). Given the complexity and heterogeneity of urban surfaces, the spatial resolution of the data and mapping units can affect the LCZ classification results. However, only a few studies have tested the scale sensitivity of LCZ mapping using grid cells (Kotharkar and Bagade 2018; Hu et al. 2019; Mouzourides et al. 2019; Zhou et al. 2022).

1.4. Challenges in GIS-based LCZ mapping

Usually, GIS-based LCZ mapping can be applied using the algorithms shown in Table 1. The rule-based classifier is a common decision-making algorithm based on “if-then” rules and the value ranges of LCZ properties (i.e. thresholds). LCZ properties are critical inputs for GIS-based LCZ classification. Not all LCZ properties, however, can be calculated quantitatively for GIS-based LCZ classification. The increasing availability and quality of high-spatial-resolution RS and GIS data have the potential to allow for the quantitative calculation of LCZ properties. Highly accurate two-dimensional (2D) and three-dimensional (3D) information is very beneficial in this context. However, the input data required in GIS-based LCZ classification tasks are not always available, especially highly accurate and fine-scale height information. The availability of light detection and ranging (LiDAR) technology provides a promising opportunity to address this problem. Recently, a few studies have integrated airborne LiDAR data to GIS-based LCZ mapping methods (Bartesaghi Koc et al. 2017, 2018; Du et al. 2020; Zhao et al. 2019, 2020). However, these studies have used classifiers with a deterministic nature (i.e. the classifier returns the same output for a given input; e.g. the rule-based classifier or LULC type-based reclassifier). These methods are unable to address uncertainties that are usually introduced in the data collection, LCZ property calculation, and threshold optimization processes.

Table 1. Types of common algorithms for GIS-based LCZ mapping.

Algorithm	References
Rule-based classifier	(Mitraka et al. 2015; Nassar, Blackburn, and Whyatt 2016; Bartesaghi Koc et al. 2017, 2018; Jin et al. 2020)
Combination of a rule-based classifier and a land use/land cover (LULC) type-based reclassifier	(Wicki and Parlow 2017; Kotharkar and Bagade 2018; Wang et al. 2018b; Zheng et al. 2018; Quan 2019; Zhao et al. 2019, 2020; Du et al. 2020; Chang et al. 2021; Cilek and Cilek 2021; Wu, Liu, and Wang 2022)
Rule-based classifier with fuzzy membership	(Lelovics et al. 2014; Unger, Lelovics, and Gál 2014; Skarbit, Gál, and Unger 2015; Skarbit et al. 2017; Geletič and Lehnert 2016; Geletič, Lehnert, and Dobrovolný 2016; Geletič et al. 2019; Quan et al. 2017; Wu et al. 2018)
Fuzzy logic classifier	(Estacio et al. 2019; Muhammad et al. 2022)
K-means clustering classifier	(Wang and Ouyang 2017)
LULC type-based reclassifier	(Alexander and Mills 2014; Oliveira, Lopes, and Niza 2020)
Naive Bayes classifier	(Hammerberg et al. 2018)
Random forest classifier	(Hu et al. 2019; Zhou et al. 2022)
Other classifiers	(Hidalgo et al. 2019)

Among these classifiers, fuzzy logic is promising in dealing with uncertainty (Zadeh 1988; Fonte et al. 2019). The main structure of the fuzzy logic classifier includes membership, aggregation, and decision functions (Quan and Bansal 2021). Based on the membership functions and value ranges of LCZ properties, the membership degree of each LCZ property can be calculated for each LCZ class. The composite membership degree for each LCZ class is obtained through the aggregation function, and the best-fitting LCZ class is assigned to the corresponding mapping unit through the decision function. By combining fuzzy logic with GIS-based LCZ classification, there is promising potential to address the “holes” (i.e. mapping units not belonging to any LCZ) and “overlaps” (i.e. mapping units being assigned to more than one LCZ) that may be caused by the combination of the LCZ classes and LCZ properties.

To date, GIS-based methods have not reached uniformity in their use of thresholds. Due to data accessibility and method specification factors, there is considerable diversity in the inclusion and calculation methods of the LCZ properties used in relevant studies (Quan and Bansal 2021). Some studies have strictly followed the thresholds defined by Stewart and Oke (2012) for classification (e.g. Lelovics et al. 2014; Unger, Lelovics, and Gál 2014; Wang et al. 2018b; Estacio et al. 2019; Zhao et al. 2019; Jin et al. 2020). However, the original thresholds are still not fully universal for major cities in the world. Previous studies have compared their newly calculated thresholds with the original thresholds and found that several of them did not match exactly (Leconte et al. 2015; Wang et al. 2018a; Zheng et al. 2018; Hidalgo et al. 2019; Zhao et al. 2019, 2020; Zhou et al. 2020a, 2020b). Moreover, some studies have adjusted some of the thresholds according to the local context (e.g. Geletič and Lehnert 2016; Nassar, Blackburn, and Whyatt 2016; Bartesaghi Koc et al. 2017, 2018; Zheng et al. 2018; Du et al. 2020; Chang et al. 2021; Cilek and Cilek 2021; Muhammad et al. 2022), thereby achieving a trade-off between generalization and localization.

Different from RS-based methods, GIS-based methods have not reached uniformity in accuracy assessments (Jiang et al. 2021; Quan and Bansal 2021). Most studies have examined the differences in air temperature from in situ measurements or land surface temperature (LST) observed from thermal infrared sensors (e.g. airborne thermal infrared sensors, ECOSTRESS, ASTER, Landsat, and MODIS) among LCZ classes. Some studies using GIS-based methods have performed accuracy assessments using reference samples. These studies have reported the overall accuracy of LCZ classification (also called the agreement between the LCZ classification results and reference samples), for example, 72% for Hong Kong (Wang et al. 2018b), 81% for five southern European cities (Oliveira, Lopes, and Niza 2020), 81%–90% for several Central European cities (Geletič and Lehnert 2016; Geletič, Lehnert, and Dobrovolný 2016; Geletič et al. 2019), >84% for the main urban areas of Nanjing (Hu et al. 2019), >90% for Beijing (Quan 2019), 92% for Berlin (Muhammad et al. 2022), and>93% for the main urban areas of Xi’an (Zhou et al. 2022). However, only a few studies have combined both validation methods (Hu et al. 2019; Quan 2019; Geletič et al. 2019; Oliveira, Lopes, and Niza 2020; Zhou et al. 2022).

1.5. Main objectives of this study

Although various GIS-based LCZ mapping methods have been developed, limited comprehensive research has been conducted on LiDAR data integration, uncertainty handling, multiscale mapping units, and multiple accuracy assessments. Against the above background, we combined airborne LiDAR data, large-scale building vector data, large-scale road vector data, and Sentinel-2 optical imagery to generate LCZ maps for the 23 special wards of Tokyo (Tokyo 23-Ku) based on a fuzzy logic classifier. The main objectives of this study are (1) to investigate the effect of the threshold optimization of LCZ properties on the LCZ classification results, (2) to explore the effects of different mapping units (blocks and multiple grid cells) on the LCZ classification results, and (3) to evaluate the LCZ classification quality by analyzing the multiscale relationships between LCZs and multidate LSTs using Landsat-8 and ASTER thermal infrared data, in addition to the confusion matrices obtained from the reference samples. In this study, the term “LCZ classification results” refers to the output; the term “LCZ classification” refers to the classification process; and the term “the LCZ classification quality” refers to the accuracy of the results.

2. Materials

2.1. Study area

Tokyo is the center of the Tokyo Metropolitan Area, the largest metropolitan area in Japan. Tokyo consists of 23 special wards (Tokyo 23-Ku), 26 cities, five towns, and eight villages. As a windy coastal city, Tokyo has a subtropical monsoon climate, with an average annual temperature of 16.5°C, annual precipitation of 1874 mm, and average annual relative humidity of 70%, based on a 2019 summary of meteorological statistics (https://www.toukei.metro.tokyo.lg.jp). Tokyo is densely populated, with numerous buildings, a variety of building types, and an irregular street network. The population of Tokyo 23-Ku reached approximately 9.73 million in 2020, corresponding to approximately 69% of Tokyo’s population (https://www.toukei.metro.tokyo.lg.jp; https://www.stat.go.jp/data/nenkan/index1.html).

The study area boundary (the green line in Figure 1) was selected by combining the administrative boundaries (gray lines in Figure 1) of Tokyo 23-Ku in 2020 and the actual area covered by the data. The study area (139°34′–139°55′E, 35°32′–35°49′N) covers an area of approximately 659.28 km² and contains a variety of building types as well as various land use and land cover types. The study area has rolling terrain to the west and relatively flat terrain to the east (Figure 1). The vast majority of the ground in the study area has been replaced by artificial features.

Figure 1. Location of the study area (outlined by the green line) and its digital terrain model (DTM) derived from airborne LiDAR data. The gray lines indicate administrative boundaries.

2.2. Data and preprocessing

The details of the data used in this study are shown in Table 2. All data were projected and/or transformed to JGD_2000_Japan_Zone_9 (WKID: 2451; Authority: EPSG), which is the Japan Plane Rectangular Coordinate System IX with Japanese Geodetic Datum (JGD) 2000. Most of the data used for LCZ classification were concentrated approximately 2016, except for the airborne LiDAR data. However, inconsistencies in the acquisition dates of the datasets have little effect on the LCZ classification results because the heights of surface objects in Tokyo 23-Ku did not change notably at the local scale during the data acquisition interval. The percentage of new buildings (by number) increased by less than 1.3% from 2001 to 2017, based on the calculation method proposed by Chen et al. (2022). In addition, the land use/cover changes of the natural classes in Tokyo 23-Ku were small and negligible at the local scale (Figure 122).

Table 2. Summary of data used in this study.

Data type	Data	Date (yyyy/mm/dd, local time)	Spatial resolution/ scale	Usage	Source
Remote sensing	Airborne LiDAR data	2001–2002	4–5 points/m^2	Classification	Geospatial Information Authority of Japan (GSI)
	Sentinel-2 MSI Level 1C imagery	2015/01/01–2016/12/31	10 m, 20 m	Classification	Google Earth Engine (GEE) platform (https://earthengine.google.com/)
	Landsat-8 collection-2 level-2 daytime surface temperature products	2016/03/01 10:15 (Spring) 2016/03/17 10:15 (Spring) 2016/05/04 10:15 (Spring) 2016/07/07 10:15 (Summer) 2016/10/27 10:16 (Autumn) 2016/12/30 10:16 (Winter)	100 m	Verification	United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/)
	ASTER level-2 AST-08 nighttime surface temperature products	2011/03/26 21:48 (Spring) 2012/12/25 21:48 (Winter) 2013/08/31 21:42 (Summer) 2013/10/09 21:48 (Autumn) 2014/12/24 21:42 (Winter) 2015/02/01 21:48 (Winter) 2017/06/23 21:42 (Summer)	90 m	Verification	Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov/products/ast_08v003/)
GIS	Zmap-TOWNII building footprint product	2017	Polygon vector (1:2500)	Classification	The residential map database from the ZENRIN Company, Japan (https://www.zenrin.co.jp/product/category/gis/basemap/zmaptown/index.html)
	Road vector data	2016	Polygon vector (1:2500)	Classification	The basic land use status survey database of the Tokyo Metropolitan Development Bureau
		2013	Line vector (1:2500)	Division of block units	The standard national digital road map (DRM) database of the Japan DRM Association (https://www.drm.jp/)
	LCZ reference samples	2016	Polygon vector	Verification	Google Earth

During the 2001–2002 period, aerial flight projects in Tokyo were conducted to collect high-spatial-resolution LiDAR data under cloudless, windless, and dry conditions. With a point density of approximately 4–5 points/m², LiDAR point cloud data have an uncertainty of less than 25 cm in the vertical direction. The LiDAR data were preprocessed as follows. (1) The LiDAR point clouds from the first returned laser pulse were interpolated into Delaunay-based TINs to generate a digital surface model (DSM) (Figure 111a) at a grid-cell size of 0.5 m. (2) LiDAR point clouds from the last returned laser pulse were filtered to obtain point clouds representing the ground. (3) The “ground” point clouds were interpolated into the Delaunay-based TINs to generate a digital terrain model (DTM) (Figure 111b) at a grid-cell size of 0.5 m. (4) A normalized digital surface model (nDSM) (Figure 111c) was generated by subtracting the DTM from the DSM.

All large-scale GIS vector data (Figures 111d, e, and f) were checked for geometry and identical records. The 164 scenes of Sentinel-2 MSI Level 1C images (top-of-atmosphere reflectance) (Drusch et al. 2012) were used to smooth out variations in shadows and natural objects. The final image (Figure 111h) was obtained by preprocessing using the Google Earth Engine platform, including by applying sensor-invariant atmospheric correction (https://eartharxiv.org/repository/view/1034/), cloud removal (using the cloud mask band “QA60”), and median extraction. In addition, it should be noted that because surface temperature data are scarcer during nighttime than during the daytime, we collected nighttime surface temperature products during the 2011–2019 period to reflect seasonal characteristics.

2.3. LCZ reference samples

LCZ reference samples were collected for 15 LCZ classes over the study area (Figures 111h and 2). First, sites were randomly selected within the study area (Congalton and Green 2019). Next, the LCZ polygons enclosing the sites were delineated manually and digitally through the visual interpretation of high-spatial-resolution imagery collected in 2016 using Google Earth (Demuzere, Kittner, and Bechtel 2021; Demuzere et al. 2022). Finally, Google Earth imagery, nDSM data, and other ancillary data (i.e. large-scale land use data, airborne aerial images, and street view images) were used to assign LCZ labels to each polygon. It should be noted that LCZs 7 (lightweight low-rise) and F (bare soil or sand) were omitted. Although single-story wood buildings (heights: 2 m–4 m) and bare soil/sand exist within the study area, these elements did not constitute LCZ classes at the local scale because they occurred only in small parts of the study area.

Figure 2. (a) Summary of LCZ reference samples and (b) the number of pixels corresponding to the LCZ reference samples at multiple grid-cell scales (100 m–1000 m).

3. Methods

The overall workflow of this study is visually summarized in Figure 3. We used two types of mapping units for LCZ classification: grid-cell units and block units. At the local scale, we generated grid cells of different sizes based on the boundary of Tokyo 23-Ku. The size of the grid cells ranged from 100 m to 1000 m at intervals of 100 m. For the block units, we delineated urban blocks within the boundary of Tokyo 23-Ku using road centerline data.

Figure 3. The overall workflow of this study. Triangulated irregular network, TIN; digital surface model, DSM; digital terrain model, DTM; normalized digital surface model, nDSM; sky view factor, SVF; bottom of atmosphere, BOA; all mapping parameters are defined in Section 3.

3.1. Calculation of LCZ properties for classification

3.1.1. Height of roughness elements (HRE)

The nDSM raster data (0.5-m pixel spacing) extracted from LiDAR data were used to calculate the height of roughness elements (HRE), which refers to the average heights of surface objects within the mapping units (Figure 3). To calculate the HRE more accurately, the ground area in the nDSM raster data must be filtered out. To do so, we extracted pixels with nDSM values less than 0.25 m and united these pixels with the road polygons. The resulting polygons were assumed to be ground (ground mask polygons), and the ground in the study area consists mainly of paved areas, water surfaces, grasses, etc. To obtain the nDSM for nonground areas (nonground nDSM), ground mask polygons were used to filter the initial nDSM raster data, and the nonground nDSM was aggregated into mapping units by using the weighted average (by area percentage) rule.

3.1.2. Sky view factor (SVF)

In terms of geometry, the sky view factor (SVF) measures the percentage of the visible portion of the sky within the hemisphere at a given point on the ground (Stewart and Oke 2012; Oke et al. 2017). This factor was computed using the DSM raster data (0.5-m pixel spacing) and Relief Visualization Toolbox (Zakšek, Oštir, and Kokalj 2011). Based on previous studies (Gál et al. 2009; Zakšek, Oštir, and Kokalj 2011; Chen et al. 2012; Unger, Lelovics, and Gál 2014), the number of horizon search directions was set to 32, and the maximum search radius was set to 200 m. The resulting SVF (2-m pixel spacing) is shown in Figure 111g. To facilitate an accurate LCZ classification, we obtained the ground SVF by filtering out nonground areas from the SVF raster data by using ground mask polygons (Figure 3). Finally, the ground SVF was aggregated into the mapping units by using the weighted average (by area percentage) rule.

3.1.3. Aspect ratio (AR)

The aspect ratio (AR), also known as height/width (H/W), is the mean height-to-width ratio of street canyons, tree spacing, and spaces between buildings and trees (Stewart and Oke 2012). The AR was computed in this study by using ground mask polygons as follows (Zhao et al. 2019): (1) mapping units were intersected with the ground mask polygons; (2) the distance of the deepest point within each intersected portion from the closest surrounding edge, namely, the thickness, was computed (Figure 133); (3) the thickness was aggregated into the mapping units by using the averaging rule; and (4) the AR for each mapping unit was obtained by dividing the HRE by its thickness.

3.1.4. Building surface fraction (BSF), impervious surface fraction (ISF), and pervious surface fraction (PSF)

The building surface fraction (BSF), also known as the building density, refers to the ratio of the building footprint area to the area of the mapping unit. We computed the BSF by using the building footprint polygons as follows: (1) the mapping units were intersected with the building footprint polygons; (2) for each mapping unit, the percentage of each intersected portion within the mapping unit was calculated by dividing the area of each intersected portion by the area of the corresponding mapping unit; and (3) the percentage of each intersected portion was aggregated into corresponding mapping units by using the summation rule.

The impervious surface fraction (ISF) is mainly associated with paved roads and exposed rocks. Because rocks in the study area are rarely exposed, we computed the ISF using road polygons. Even though the roads actually contain a small amount of vegetation, we ignored the effect of vegetation on the LCZ classification at the local scale. The ISF calculation procedure was similar to that of BSF.

Pervious surfaces are composed mainly of water, vegetation, and bare soil. Generally, the ground surface is assumed to be a combination of buildings, impervious surfaces, and pervious surfaces. Therefore, the PSF was derived from the following equation: $PSF=1-BSF-ISF$.

3.1.5. Surface albedo (SA)

For a surface, the surface albedo (SA) is the ratio of the amount of reflected solar radiation to the amount of received solar radiation. We computed SA using the Sentinel-2 data (Figure 3) as follows: (1) the surface broadband albedo (10-m pixel spacing) was retrieved from six bands of Sentinel-2 (bands 2, 3, 4, 8, 11, and 12) with the narrow-to-broadband conversion coefficients developed by Bonafoni and Sekertekin (2020); then, (2) it was aggregated into the mapping units using the weighted average (by area percentage) rule.

3.2. Fuzzy logic classification

The original thresholds were further optimized (Figure 3) based on the statistics of the LCZ reference samples (Figure 4), taking into account the specific context of Tokyo. Both the optimized and original property thresholds corresponding to each LCZ class are shown in Table 3. Because there were deviations between the original thresholds (Table 3) and the statistical values (Figure 4), we designed two LCZ classification schemes: one based on the original thresholds (Scheme 1, S1) and one based on the optimized thresholds (Scheme 2, S2).

Figure 4. Box plots of seven LCZ properties in all LCZ reference samples. The white circle and the black vertical line within each box represent the mean and median values, respectively. The LCZ class numbers/letters/colors correspond to the classes defined in Figure 2.

Table 3. Threshold values of seven LCZ properties used for fuzzy logic classification. The optimized thresholds are in bold, and the original thresholds are in parentheses. The LCZ class numbers/letters correspond to the classes defined in Figure 2.

LCZ Class	HRE (m)	SVF (–)	AR (–)	BSF (%)	ISF (%)	PSF (%)	SA (–)
LCZ 1	>25	0.2–0.6 (0.2–0.4)	>2.00	>40 (40–60)	<60 (40–60)	<20 (<10)	0.10–0.20
LCZ 2	10–25	0.3–0.9 (0.3–0.6)	0.75–2.00	>40 (40–70)	<50 (30–50)	<30 (<20)	0.10–0.20
LCZ 3	3–10	0.2–0.9 (0.2–0.6)	0.75–1.50	>40 (40–70)	<50 (20–50)	<40 (<30)	0.10–0.20
LCZ 4	>25	0.5–0.7	0.75–1.25	10–40 (20–40)	20–40 (30–40)	40–50 (30–40)	0.12–0.25
LCZ 5	10–25	0.5–0.8	0.30–0.75	10–40 (20–40)	10–50 (30–50)	30–50 (20–40)	0.12–0.25
LCZ 6	3–10	0.6–0.9	0.30–0.75	20–40	10–50 (20–50)	40–70 (30–60)	0.12–0.25
LCZ 7	2–4	0.2–0.5	1.00–2.00	60–90	<20	<30	0.15–0.35
LCZ 8	3–10	>0.7	0.10–0.30	>30 (30–50)	<50 (40–50)	<20	0.15–0.25
LCZ 9	3–10	>0.8	0.10–0.25	10–20	<20	60–80	0.12–0.25
LCZ 10	5–15	0.6–0.9	0.20–0.50	20–30	20–40	40–50	0.12–0.20
LCZ A	3–30	<0.8 (<0.4)	>1.00	<10	<10	>90	0.10–0.20
LCZ B	3–15	0.5–0.9 (0.5–0.8)	0.25–0.75	<10	<10	>90	0.10–0.25 (0.15–0.25)
LCZ C	<2	0.7–0.9	0.25–1.00	<10	<10	>90	0.10–0.30 (0.15–0.30)
LCZ D	<1	>0.9	<0.10	<10	<10	>90	0.15–0.25
LCZ E	<0.25	>0.9	<0.10	<10	>90	<10	0.15–0.30
LCZ F	<0.25	>0.9	<0.10	<10	<10	>90	0.20–0.35
LCZ G	<0.25	>0.9	<0.10	<10	<10	>90	0.02–0.10

The fuzzy logic classification was achieved as follows. Consider a set (universe of discourse) $G=\left\{{\mathit{g}}_1,{\mathit{g}}_2,\dots ,{\mathit{g}}_m\right\}$ consisting of m mapping units to be classified. Each mapping unit ${\mathit{g}}_i=\left(g_{i,1},g_{i,2},\dots ,g_{i,n}\right)$ is a vector consisting of n properties, where $g_{i,j}\left(1\le i\le m,1\le j\le n\right)$ is the value of the j-th property at the i-th mapping unit. We denote classes by $k\left(1\le k\le s\right)$, where s is the total number of classes.

The s classes of mapping units can be expressed as s fuzzy subsets of the set G. A fuzzy set $C_k$ of G is a set of ordered pairs:

(1) $C_k=\left\{\left({\mathit{g}}_i,{\mu }_{C_k}\left({\mathit{g}}_i\right)\right)|{\mathit{g}}_i\in G\right\}\left(1\le k\le s\right),$(1)

where ${\mu }_{C_k}\left({\mathit{g}}_i\right)$ is the degree of membership of ${\mathit{g}}_i$ in $C_k$ and ${\mu }_{C_k}:G\to \left[0,1\right]$ is the membership function of $C_k$ in G.

Each fuzzy subset $C_k$ is characterized by n properties. Let $C_{k,j}$ be each fuzzy subset (covering the j-th property) in a given fuzzy set $C_k$. A fuzzy set $C_{k,j}$ is also a set of ordered pairs:

(2) $C_{k,j}=\left\{\left.\left(g_{i,j},{\mu }_{C_{k,j}}\left(g_{i,j}\right)\right)\right|1\le i\le m\right\}\text{break}\hspace{0.17em}\left(1\le k\le s,\hspace{0.17em}1\le j\le n\right),$(2)

where ${\mu }_{C_{k,j}}\left(g_{i,j}\right)$ is the degree of membership of $g_{i,j}$ in $C_{k,j}$.

Based on the threshold ranges in Table 3, we used three types of trapezoidal membership functions (Table A1). There are clear ranges for some properties, including SVF (0–1), BSF (0–100), ISF (0–100), PSF (0–100), and SA (0–1). However, for properties with no explicit ranges (HRE and AR), the boundary thresholds of their membership functions were set depending on the maximum and minimum values of the actual data.

The thresholds (Table 3) and membership functions (Table A1) of the properties enabled us to calculate ${\mu }_{C_k}\left({\mathit{g}}_i\right)$ as follows:

(3) ${\mu }_{C_k}\left({\mathit{g}}_i\right)=\underset{n}{\overset{j=1}{\wedge }}{\mu }_{C_{k,j}}\left(g_{i,j}\right)=\underset{1\le j\le n}{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{{\mu }_{C_{k,j}}\left(g_{i,j}\right)\right\}.$(3)

Finally, the class $t\in \left\{1,2,\dots ,s\right\}$ of ${\mathit{g}}_i$ was identified as follows:

(4) ${\mu }_{C_t}\left({\mathit{g}}_i\right)=\underset{s}{\overset{k=1}{\vee }}{\mu }_{C_k}\left({\mathit{g}}_i\right)=\underset{1\le k\le s}{\mathrm{m}\mathrm{a}\mathrm{x}}\left\{{\mu }_{C_k}\left({\mathit{g}}_i\right)\right\}.$(4)

To include more contextual information, we adopted a 3 × 3 majority filter for the LCZ results at multiple grid-cell scales (100 m–1000 m). We did not implement post processing for the LCZ results at the block scale due to the irregular shapes, different sizes, and complex neighborhoods of the blocks.

3.3. Accuracy assessments

We evaluated the LCZ results for multiple mapping units in terms of confusion matrices based on LCZ reference samples. The overall accuracy (OA), user’s accuracy (UA), and producer’s accuracy (PA) were derived from the confusion matrices (Congalton and Green 2019).

Originally, and most ideally, each LCZ should exhibit distinct screen-height air temperatures (Stewart and Oke 2012; Stewart, Oke, and Krayenhoff 2014). Subsequently, it has been suggested that each correctly classified LCZ class should demonstrate its specific surface thermal environment (Geletič, Lehnert, and Dobrovolný 2016). Therefore, the LCZ results obtained under S2 were further evaluated by using LSTs from multiple dates. To explore the differences in the mean LSTs among LCZ classes, we carried out a one-way analysis of variance (ANOVA). First, we evaluated assumptions by using normality tests (histogram comparisons, Q – Q plots, and Kolmogorov‒Smirnov tests) and a test of homogeneity of variances (Levene’s test). All the test designs (Kolmogorov‒Smirnov tests and Levene’s test) met the two assumptions at the 99.9% confidence level (i.e. a significance level of 0.001; p value<0.001). When the one-way ANOVA F test indicated a statistically significant difference between LSTs (p value<0.001), we chose Tamhane’s T2 post hoc test to perform pairwise multiple comparisons based on the results of the above tests. Finally, detailed insight into the results of the pairwise multiple comparisons was obtained. For each given LCZ class, the percentages (0–100%) of the number of other LCZs that were significantly (p < 0.05) different from the given LCZ in terms of the mean LST (also called “hits”) were calculated.

4. Results

4.1. Comparison between S1 and S2

The LCZ results at multiple scales under S1 and S2 are presented in Figures 5, 144, and 155. Table 4 shows the corresponding OAs. The percentages of the area occupied by each LCZ class are shown in Figures 6a,b, and the differences in the PAs and UAs for the different LCZ maps are shown in Figures 6c-f.

Figure 5. LCZ maps at the 100-m grid-cell and block scales obtained using S1 (original thresholds) and S2 (optimized thresholds). The LCZ class numbers/letters correspond to the classes defined in Figure 2.

Figure 6. (a, b) Percentages, (c, d) producer’s accuracies (PAs), and (e, f) user’s accuracies (UAs) of the LCZ classes obtained at multiple scales using S1 (original thresholds) and S2 (optimized thresholds). The LCZ class numbers/letters correspond to the classes defined in Figure 2.

Table 4. Overall accuracies (OAs) (%) of the multiscale LCZ results obtained using S1 (original thresholds) and S2 (optimized thresholds).

Scheme	Grid cells (m)
	100	200	300	400	500	600	700	800	900	1000	Blocks
S1	46.28	44.43	39.07	35.85	31.16	32.06	29.21	29.14	29.76	28.06	39.01
S2	80.34	71.81	65.91	58.59	52.84	48.46	42.92	37.63	39.19	36.51	66.88

From a qualitative point of view (Figures 5, 144, and 155), the spatial distribution of LCZs obtained using S2 was notably different from that obtained using S1. Using S1, LCZ A (dense trees) in central Tokyo was misclassified as LCZ G (water). The OAs of the LCZ results at all scales were higher for S2 than for S1 (Table 4). Compared to the results obtained using S1, the OA was improved by 34.06% when using S2 at the 100-m grid-cell scale; this was the largest increase obtained among all scales.

Compared to S1, LCZ 3 (compact low-rise) had the largest increase in area at all scales when using S2, followed by LCZs 2 (compact mid-rise), 1 (compact high-rise), A, and B (scattered trees) (Figures 6a,b). Conversely, the areas of LCZs 6 (open low-rise), 9 (sparsely built), and G decreased at all scales when using S2. Figure 7 depicts the different LCZ amounts that remained unchanged and changed from S1 to S2 at the 100-m grid-cell scale. As shown in Figure 7, at the 100-m grid-cell scale, the largest transitions were those from LCZ 6 to LCZ 3 (approximately 129.8 km²) and from LCZ 9 to LCZ 6 (approximately 81.0 km²) after using S2.

Figure 7. Transitions in the LCZ classes from S1 (original thresholds) to S2 (optimized thresholds) at the 100-m grid-cell scale (unit: km²). The LCZ class numbers/letters correspond to the classes defined in Figure 2.

As shown in Figures 6c,d, after using S2, the greatest improvement in the PA at most scales was associated with LCZ 3, followed by LCZs A, 1, and 2. When using S2, the PAs of LCZs 6 and 5 (open mid-rise) increased at the block scale and within the 500-m grid-cell scale. Using S2, PA increases were observed for LCZs B and C (bush, scrub) at the 100-m grid-cell and block scales. After using S2, the PAs of LCZs 10 (heavy industry), D (low plants), E (bare rock or paved), 8 (large low-rise), G, and 4 (open high-rise) changed slightly at most scales. Overall, the PA improvements decreased with increased grid-cell sizes for LCZs 3, A, 1, 6, and 5.

As shown in Figures 6e,f, after using S2, the greatest UA improvements at most scales were associated with LCZ 3, followed by LCZs A, 6, and 5. The UAs of LCZs B and 2 increased at the 100-m grid-cell and block scales using S2. When using S2, the UAs of LCZs 8, C, G, and 9 increased at the block scale and within the 300-m grid-cell scale. After using S2, the UAs of LCZs 10, D, and E changed slightly at most scales. Overall, the UA improvement decreased as the grid-cell size increased for LCZs 3, A, 6, 5, and G.

The confusion matrices for LCZ maps at the 100-m grid-cell and block scales under both schemes are presented in Figure 8. Compared to S1, the confusions among built LCZ classes (LCZs 1–10) and among land cover LCZ classes (LCZs A – G) were substantially reduced at the 100-m grid-cell and block scales after using S2. The confusion between the built LCZs and land cover LCZs was slightly reduced at the 100-m grid-cell and block scales using S2. After using S2, at the 100-m grid-cell scale, large confusion reductions between LCZs 3 and 6; between LCZ G and LCZs C, B and A; between LCZs 8 and 9; and between LCZs 1 and 2 were obtained. At the block scale, large reductions in the confusion among LCZs G, C, and A were observed, as well as between LCZs B and G, LCZs 8 and 9, LCZs 3 and 6, and LCZs 5 and 9.

Figure 8. Confusion matrices for the LCZ maps obtained at 100-m grid-cell and block scales using S1 (original thresholds) and S2 (optimized thresholds). The LCZ class numbers/letters correspond to the classes defined in Figure 2.

4.2. LCZ classification at multiple scales using the optimized thresholds (S2)

As presented in the second column of Figure 5, at the 100-m grid-cell and block scales, the central region of Tokyo 23-Ku was identified as LCZ 1 (compact high-rise), surrounded by LCZs 2 (compact mid-rise), 3 (compact low-rise), and 6 (open low-rise). LCZs 4 (open high-rise), 5 (open mid-rise), and 9 (sparsely built) were scattered within a “radial ring” LCZ pattern. LCZs 1 and 2 were clustered not only in the central area of Tokyo 23-Ku, but also in the surrounding subcity areas. At the 100-m grid-cell and block scales, LCZ 6 occupied the largest area, followed by LCZs 3 and G (water) (Figure 6b).

As the size of the grid cells increased gradually to 1000 m, the basic composition of the “radial ring” LCZ pattern did not change (Figure 155). However, this pattern was more pronounced at the coarser grid-cell scales. As shown in Figure 6b, as the size of grid cells increased, the areas of LCZs 2, 1, 3, C (bush, scrub), and A (dense trees) increased, and the areas of LCZs 6, G, E (bare rock or paved), 9, 5, D (low plants), and B (scattered trees) decreased. As the size of the grid cells increased, the classes with larger area changes were LCZs 2, 6, G, 1, and 3. Additionally, a few classes (LCZs 8 (large low-rise), 10 (heavy industry), and 4) disappeared at coarser grid-cell scales (≥600 m).

For LCZ classification using S2, the highest OA (80.34%) was achieved at the 100-m grid-cell scale, whereas the lowest OA (36.51%) was obtained at the 1000-m grid-cell scale (Table 4). The OA of the LCZ classification at the block scale was lower than that at the 100-m grid-cell scale, which was 66.88%. Overall, the OA of the corresponding LCZ classification decreased as the size of grid cells increased.

As shown in Figures 6d,f, the PAs of LCZs G, A, E, D, 3, and 6 were all greater than 82% at the 100-m grid-cell scale, indicating that these classes had low omission errors. In contrast, at the 100-m grid-cell scale, the PAs of LCZs 10, 8, 4, B, and 5 were all less than 37%, exhibiting high omission errors. At the 100-m grid-cell scale, the UAs of LCZs D, G, E, 3, and 6 exceeded 80%, whereas those of LCZs 10, C, 9, and 8 were all less than 35%. LCZs G, D, E, 3, and 6 all had relatively high PAs and UAs at this scale, whereas those of LCZs 10 and 8 were both relatively low. At the block scale, the classes with relatively high PAs (>66%) were LCZs G, A, 3, C, and 6; the classes with relatively high UAs (>67%) were LCZs 3, G, E, 1 and 2. At the block scale, only LCZs G and 3 had PA and UA scores that were both above 73%. Overall, the PA and UA of each LCZ class (except for LCZ 10) showed a decreasing trend as the size of the grid cells increased, with some differences in the magnitude of each decrease.

Overall, using S2, all of the confusions were more pronounced at the block scale than at the 100-m grid-cell scale (Figure 8). At the 100-m grid-cell scale, confusion occurred mainly between LCZs B and C, LCZs 3 and 8, LCZs 1 and 2, LCZs 9 and 10, and LCZs 6 and 9. At the block scale, major confusion existed between LCZs D and G, LCZs 10 and G, LCZs E and G, LCZs B and C, LCZs 1 and 2, LCZs 5 and 8, LCZs 5 and 6, LCZs 3 and 6, LCZs 6 and 9, and LCZs 9 and A.

4.3. Relationships between LCZs and multidate LSTs

Overall, the mean LST (both daytime and nighttime) for each LCZ class was more variable by date than by scale (Figure 9a). Figure 9b shows the percentage of the number of “hits” (i.e. the number of times a specific LCZ exhibits significant differences in the mean LST compared to the other LCZ classes) in pairwise multiple comparisons. Statistically significant differences in the mean LSTs between most LCZ classes were more pronounced at the block scale and within the 200-m grid-cell scale for most dates (Figure 9b). Among grid cells larger than 200 m, the “hits” for most LCZ classes exhibited different variations in terms of both the scale and date.

Figure 9. (a) Mean LSTs for each LCZ class at multiple scales using S2 (optimized thresholds); (b) percentages of “hits” showing significant differences (p < 0.05) in the mean LSTs for each LCZ class as a result of multiple comparison tests using S2. The LCZ class numbers/letters correspond to the classes defined in Figure 2. Dates are in the format year/month/day. The numbers in parentheses in the LCZ classes represent the sequential codes of the mapping scales. As an example, the number “1” represents “100 m,” the number “10” represents “1000 m,” and the number “11” represents “blocks”.

Generally, at all scales, the mean daytime LSTs for each LCZ class were ranked in descending order as follows: July 7 (summer), May 4 (spring), October 27 (autumn), March 17 (spring), March 1 (spring), and December 30 (winter). At the block scale, the mean daytime LSTs of LCZ 3 (compact low-rise) were highest for all dates, followed by LCZ 6 (open low-rise). On March 17 (spring), May 4 (spring), and July 7 (summer), the mean daytime LSTs of LCZ 3 were highest at most grid-cell scales, followed by LCZ 6. On March 1 (spring), October 27 (autumn), and December 30 (winter), the mean daytime LSTs of LCZ E (bare rock or paved) were highest at most grid-cell scales.

At all scales, the mean nighttime LSTs for most LCZ classes were relatively high on August 31 (summer), June 23 (summer), and October 9 (autumn). In spring (March 26) and winter (December 25, December 24, and February 1), the mean nighttime LSTs of LCZ G (water) were highest at all grid-cell scales. In summer (August 31 and June 23), the mean nighttime LSTs of LCZ E were highest at all grid-cell scales, followed by those of LCZ 1 (compact high-rise). In autumn (October 9), the mean nighttime LSTs of LCZ 1 were highest at most grid-cell scales, followed by those of LCZ E. At the block scale, the highest mean LSTs during nighttime were observed in LCZ 2 (compact mid-rise) for summer (August 31 and June 23) and in LCZ 4 (open high-rise) for spring (March 26) and winter (December 25, December 24, and February 1).

We also examined the relationships between 100-m grid-cell LCZs and multidate LSTs (Figure 10). In general, there were differences in the statistical distributions of LSTs (both daytime and nighttime) among LCZ classes, regardless of the date. For the LCZ compact classes (LCZs 1–3), the mean daytime LSTs were ranked in descending order as LCZs 3, 2, and 1; the mean nighttime LSTs were ranked in descending order as LCZs 1, 2, and 3. For LCZ open classes (LCZs 4–6), the mean daytime LSTs were ranked in descending order as LCZs 6, 5, and 4; the mean nighttime LSTs were ranked in descending order as LCZs 4, 5, and 6. Regardless of the date, the mean LST (both daytime and nighttime) was higher for LCZ 3 (compact low-rise) than for LCZ 6 (open low-rise). The mean LST was lower for LCZ 1 (compact high-rise) than for LCZ 4 (open high-rise) only on December 30 (winter daytime), March 26 (spring nighttime), and February 1 (winter nighttime). The mean LST was lower for LCZ 2 (compact mid-rise) than for LCZ 5 (open mid-rise) only on March 1 (spring daytime), March 17 (spring daytime), October 27 (autumn daytime), and December 30 (winter daytime).

Figure 10. Violin density plots of LSTs for each LCZ class at the 100-m grid-cell scale using S2 (optimized thresholds). The white circles indicate the mean, and the horizontal white lines indicate the median. The LCZ class numbers/letters correspond to the classes defined in Figure 2. Dates are in the format year/month/day.

As shown in Figure 10, regardless of the date, the mean daytime LSTs of LCZ E (bare rock or paved) were highest among the land cover LCZ classes (LCZs A – G). For dates excluding December 30 (winter), the mean daytime LST of LCZ G (water) was lowest among the land cover LCZ classes. On December 30 (winter), the mean daytime LST of LCZ A (dense trees) was lowest, followed by that of LCZ C (bush, scrub). On March 1 (spring), March 17 (spring), and October 27 (autumn), the mean daytime LSTs of LCZ A (dense trees) were the second lowest among the land cover LCZ classes. The mean daytime LSTs of LCZ D (low plants) were the second lowest among the land cover LCZ classes on May 4 (spring) and July 7 (summer). In spring (March 26) and winter (December 25, December 24, and February 1), among the land cover LCZ classes, the mean nighttime LSTs of LCZ G were highest, whereas the mean nighttime LSTs of LCZs D and C were lowest. In summer (August 31 and June 23) and autumn (October 9), among the land cover LCZ classes, the mean nighttime LSTs of LCZ E were highest, whereas the mean nighttime LSTs of LCZs G and D were lowest.

5. Discussion

5.1. Effects of threshold values of LCZ properties

As expected, performing LCZ classification after threshold optimization yielded satisfactory results. This revealed the potential of combining fuzzy logic and threshold optimization for GIS-based LCZ classification in urban areas. When using LCZ properties to classify different cities with different regional climates, it seems reasonable to adjust the threshold values of the LCZ properties according to the specific urban form (Geletič and Lehnert 2016; Nassar, Blackburn, and Whyatt 2016; Bartesaghi Koc et al. 2017, 2018; Perera and Emmanuel 2018).

We found that the SVF values calculated using DSM data were higher than the original threshold values, especially for LCZs 1 (compact high-rise), 2 (compact mid-rise), 3 (compact low-rise), and A (dense trees). The SVFs of the LCZ compact classes (LCZs 1–3) were higher than the original thresholds, probably because of the heterogeneity of building heights and building densities. Typically, LCZ A and LCZ G (water) are quite different, but they were confused in some areas (e.g. the imperial palace and moat in central Tokyo) under S1. This was because of the deviation of the SVF threshold values for LCZ A. The biased SVFs of LCZ A were consistent with the results reported by Zhao et al. (2019) because LiDAR technology can obtain information only on the tree canopy, while the ground is mostly obscured. In addition, the SVF results calculated using DSM data were affected by the search radius. The influence of the SVFs calculated at different search radii on the LCZ classification results should be explored in the future. Due to the complexity of and variability in the geometry of urban surface objects (especially buildings), it is still difficult to calculate AR in a standard way. The ARs in most LCZ classes were considerably different from the original threshold, probably because of differences in the calculation methods. The potential of the approach used to calculate AR in this study for LCZ classification should be further investigated.

However, some challenges remained even after optimizing the threshold. In particular, the classification accuracies of LCZ 10 (heavy industry) under S2 were still very low at all scales. In Tokyo 23-Ku, LCZ 10 was scattered and distributed over a small area, making it difficult to identify. Another potential reason for this uncertainty may be that the common thresholds for LCZ 10 are still difficult to define, which may be a conceptual problem (Geletič and Lehnert 2016). This problem may also exist for LCZ 8 (large low-rise) in the study area. It is worth considering that these two LCZ classes (i.e. LCZs 8 and 10) may be best distinguished by their thermal, radiative, and metabolic properties rather than by their geometric and surface cover properties. In addition, threshold optimization may not work well on relatively coarse grid cells. On coarse grid-cell scales, threshold optimization did not lead to much improvement in PA and UA for some LCZ classes in this study, probably because of the inherent scale effect of LCZ classification (Bechtel et al. 2019a).

For the LCZ results with the highest classification accuracy, the confusions between LCZ classes persisted (Figure 8); this may be explained by the heterogeneous mixture of actual land cover and difficulty in characterizing the 2D/3D spatial arrangement and configuration of urban surface objects. Considering the heterogeneity and complexity of the buildings in Tokyo 23-Ku, there may not be clear boundaries between different LCZs in the real situation. Even though LCZs are defined as uniform surfaces, each LCZ in Tokyo 23-Ku is not actually an absolutely homogeneous region. For example, in the central area of Tokyo, high-rise and mid-rise buildings alternate densely. Using LCZ subclasses can help to address this situation (Kotharkar and Bagade 2018; Perera and Emmanuel 2018; Zhou et al. 2020a, 2020b), but there are still challenges in delineating LCZ subclasses using universal LCZ properties. A major challenge associated with GIS-based LCZ identification is the variability in LCZ properties within a single LCZ class (Figure 4). For example, the land cover LCZs (LCZs A – G) in urban areas are typically characterized by higher percentages of impervious and building surfaces than those in rural areas.

5.2. Effects of different mapping units

In this study, both grid cells (100 m–1000 m) and blocks were regarded as a kind of scene complex with specific semantic features at their corresponding scales. The spatial distribution, area percentage, and classification accuracy of the LCZ results differed among different mapping units, revealing the inherent scale dependency of the LCZs and the edge effects of the mapping units. Among the grid-cell scales, the 100-m grid cells showed the most complicated, informative, and fine-scale LCZ details; the 1000-m grid cells showed the simplest, most abstract, and coarsest LCZ information. As expected, the LCZ classification accuracy gradually decreased with increasing grid-cell size, consistent with the findings of previous studies (Hu et al. 2019; Zhou et al. 2022). The apparent differences in the classification accuracies of LCZs at different mapping units indicated that the reference samples used for verification also suffered from scale effects. The shape, size, and location of the reference samples impact their LCZ properties, and these samples are limited by human experience and expert knowledge. Ideally, performing an accuracy assessment by using reference samples that match the scale of each type of mapping unit may address this issue. However, characteristics such as the shape, size, and location of reference samples as well as their seasonal variations need to be considered. Moreover, the selection of reference samples also requires a great deal of subjective expert knowledge, time, and resource costs (Bechtel et al. 2017; Zhu et al. 2020; Xu et al. 2021). Therefore, assessing the accuracy of LCZ classification results using LSTs may become another potential method.

As expected, relatively fine grid cells lead to the fragmentation of the LCZ classification results. As the size of grid cells increases, the regular edges of LCZs become progressively clearer, and the area within the grid cell becomes more heterogeneous and mixed; this is the main limitation of coarse grid cells (Kotharkar and Bagade 2018; Mouzourides et al. 2019). In general, the surface objects contained within large grid cells are abundant. When these various surface objects are aggregated into coarse grid cells, some LCZ classes characterized by more dispersion or smaller area shares may not be dominant in relatively large grid cells, causing some details to be omitted. For this reason, in most cases, LCZ open classes (i.e. LCZs 4–6) are more difficult to classify than LCZ compact classes (i.e. LCZs 1–3) at the same building height level.

Compared to the LCZ results obtained at the 100-m grid-cell scale, the LCZs at the block scale followed irregular boundaries more intuitively but with a lower accuracy. The LCZ classification results obtained at the block scale were limited by the road centerline data that delineated urban blocks. Compared to urban areas with well-developed road networks, the delineation of blocks in natural areas with less developed road networks was less accurate (e.g. the areas close to Tokyo Bay and Haneda Airport). Notably, some LCZ E (bare rock or paved) areas in our LCZs at the block scale exhibited small linear features, mainly because the road centerline data used in this study differentiated directions on the same road, resulting in the presence of more than one centerline on some roads. At the block scale after optimizing the thresholds, a few areas that should have been classified as LCZ G (water) were classified as LCZs A (dense trees) and B (scattered trees) (e.g. the river closest to the central area in Figure 5). This may have been because of the shapes and sizes of irregular block units. These block units were usually larger than the actual coverage of the water bodies. As a result, these block units, which should be classified as LCZ G, were misclassified by including some other surface objects.

5.3. LST differentiation of LCZs

Typically, the inter-LCZ variability in multidate (including daytime and nighttime) LSTs revealed that individual LCZs exhibited specific surface thermal environments associated with their surface characteristics (Stewart and Oke 2012; Stewart, Oke, and Krayenhoff 2014). For example, in terms of the LCZs, the negative relationship between the building heights and daytime LSTs may have been affected by aerodynamic differences and shadows (Li et al. 2011; Zhao et al. 2014; Manoli et al. 2019; Chen et al. 2022). Conversely, the positive relationship between building height and nighttime LST may have been related to heat storage release (Guo, Wu, and Schlink 2021). In addition, variability in multidate (including daytime and nighttime) LSTs occurred within each LCZ class (i.e. the intra-LCZ variability in multidate LSTs), indicating the potential influence of surrounding heterogeneous environments. The relationships between LCZs and multidate LSTs presented in this study were compatible with the context (subtropical monsoon climate and coastal sea and land breezes) in which Tokyo 23-Ku is located.

Similar to the evaluation results obtained based on reference samples, the evaluation results obtained using the multidate LSTs exhibited the inherent scale variability in LCZs. Small grid cells and blocks are likely good representatives of LCZs, as they are used in most studies, e.g. in the WUDAPT method (Bechtel et al. 2019a). There is some agreement between these two types of accuracy assessments. For example, the LCZ classes with low numbers of “hits” (e.g. LCZs 10 (heavy industry), 8 (large low-rise), 4 (open high-rise), and B (scattered trees)) also had low PAs and UAs and occupied small areas (Figures 6 and 9b). It seems feasible to perform an accuracy assessment using a combination of LCZ reference samples and multidate LSTs, especially in cases of unbalanced LCZ reference samples. Combining multiple thermal infrared sensors to form LST time series is an appropriate means to gain a more complete understanding of the LCZ-LST relationship. However, satellite-acquired LSTs still suffer from limited available data (especially at night), retrieval uncertainty, and thermal anisotropy (Berger et al. 2017).

5.4. Limitations and future research directions

One of the major limitations is the LCZ concept itself, as it is defined at a horizontal scale of hundreds of meters to several kilometers (Stewart and Oke 2012; Stewart, Oke, and Krayenhoff 2014). The 100-m scale, the default scale for most RS-based LCZ mapping methods (e.g. the WUDAPT method (Bechtel et al. 2015, 2019a; Ching et al. 2018; Demuzere, Kittner, and Bechtel 2021)), lies between the local scale and microscale in climate research. For GIS-based LCZ mapping, the vast majority of studies have been conducted at scales between 100 m and 1000 m. Because a certain range of appropriate scales is allowed in the LCZ concept (Bechtel et al. 2015), other scales except for 100 m–1000 m have been investigated in several studies. Recently, RS-based LCZ mapping has been attempted at relatively fine spatial resolutions (e.g. 50 m (Yoo et al. 2020), 30 m (Verdonck et al. 2017), and 10 m (La, Bagan, and Yamagata 2020; Rosentreter, Hagensieker, and Waske 2020; Chen et al. 2021b)). Additionally, a few GIS-based LCZ maps have been constructed at the 50-m grid-cell scale (Oliveira, Lopes, and Niza 2020). Therefore, the effect of our workflow for automated LCZ classification at relatively fine grid-cell scales needs to be further investigated. Moreover, the spatial sensitivity of the geographic locations of grid cells (Zheng et al. 2018) used for LCZ mapping needs to be further explored.

The LCZ results obtained using S2 indicated that the LCZ properties and their corresponding thresholds can be further optimized, for example, by incorporating more precise data sources and property calculations (e.g. anthropogenic heat release). Although the fuzzy logic classifier can help address uncertainties, the performance of the classifier is limited directly by the LCZ properties and their corresponding thresholds (Quan and Bansal 2021). In the future, more localized property definitions and corresponding thresholds will also contribute to GIS-based LCZ mapping, but this will require more reference samples and prior knowledge. A comprehensive characterization of the 2D/3D spatial arrangement and configuration of surface objects will contribute to a greater understanding of LCZ classification. Given the heterogeneity and complexity of urban forms and landscapes, the subclasses of LCZs need to be considered. In addition, automated components need to be introduced in threshold optimization.

LCZ maps have become increasingly valuable in terms of urban form and urban climate. From a policy perspective, climate change adaptation strategies for cities are becoming increasingly important. In high-density cities, the identification of LST differences for each LCZ class provides insight into the prioritization and targeting of such mitigation and adaptation strategies. LCZ maps can help support urban planning, policy-making, and climate-responsive designs. Future work could aim to explore the application of LCZ maps in relatively wide fields (including crosscutting phenomena), such as the urban heat island effect (e.g. Bechtel et al. 2019b; Oliveira et al. 2022; O’Malley and Kikumoto 2022), thermal comfort and heat stress (e.g. Geletič et al. 2018; Verdonck et al. 2018; Wang et al. 2021; Wu, Liu, and Wang 2022), energy consumption (e.g. Yang et al. 2020), carbon emissions (e.g. Wu et al. 2018), and human health (e.g. Brousse et al. 2019).

6. Conclusions

In this study, a fuzzy logic classifier was developed for LCZ mapping based on multisource RS (airborne LiDAR and Sentinel-2) and GIS (building vectors and road vectors) data in the 23 special wards of Tokyo. We examined the effects of threshold optimization and mapping units (i.e. multiple grid cells (100 m–1000 m) and blocks) on GIS-based LCZ mapping. Multisource thermal infrared data (i.e. Landsat-8 and ASTER data) were used to enhance the understanding of the relationships between multiscale LCZs and multidate LSTs.

Utilizing complementary multisource remote sensing and GIS data to classify LCZs more effectively is important for studying the urban form and surface thermal environment. Our study highlighted the potential of multisource RS and GIS data for GIS-based LCZ mapping. In terms of LCZs, the 3D information provided by LiDAR and the 2D information provided by GIS vectors are particularly important. In GIS-based LCZ mapping, it is necessary to optimize the original thresholds to achieve satisfactory results for localization. Considering the uncertainties in GIS-based LCZ mapping, using a fuzzy logic classifier is a promising method for creating a feasible and universal framework. The LCZ maps and their corresponding multidate LSTs varied among the different mapping units, and 100-m grid cells appeared to be well-characterized for GIS-based LCZ mapping. A combination of accuracy assessments using reference samples and multidate LSTs was found to be effective. The reference samples can reveal LCZ accuracies, but they are also affected by scale effects. The multidate LSTs are able to reflect the unique surface thermal characteristics of each LCZ class.

Our LCZs indicated a “radial ring” form of Tokyo 23-Ku, changing from the inside to the outside in the order of LCZs 1 (compact high-rise), 2 (compact midrise), 3 (compact low-rise), and 6 (open low-rise). The urban form defined by LCZs is crucial for sustainable urban development. Therefore, the datasets used and LCZ classification results obtained in our present study need to be further explored for use in research on urban morphology and urban climate.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are available from the corresponding author, H. Bagan, upon reasonable request. The data are not publicly available due to legal restrictions.

Additional information

Funding

The work was supported by the  National Key R&D Program of China [2022YFE0119500] and by the Science and Technology Commission of Shanghai Municipality, China [22010503600].

Appendices

Figure 111. (a) Digital surface model (DSM), (b) digital terrain model (DTM), and (c) normalized digital surface model (nDSM) (0.5-m pixel spacing) derived from airborne LiDAR data; (d) building footprint polygons; (e) road polygons; (f) road centerlines; (g) the sky view factor (SVF) (2-m pixel spacing) computed from the DSM; (h) Sentinel-2 MSI bottom-of-atmosphere (BOA) reflectance (R/G/B = bands 4/3/2, 10-m spatial resolution) and LCZ reference samples.

Figure 122. (a) Land use map of Tokyo 23-Ku in 1997 (100-m mesh data); (b) land cover map of Tokyo 23-Ku in 2001 (30-m spatial resolution); (c) land use map of Tokyo 23-Ku in 2016 (version 1, 100-m mesh data); (d) land use map of Tokyo 23-Ku in 2016 (version 2, 50-m mesh data); (e) land use map of Tokyo 23-Ku in 2016 (version 3, vector data). [Sources: (a), (c), and (d) are from the National Land Information Division, National Spatial Planning and Regional Policy Bureau, Ministry of Land, Infrastructure, Transport and Tourism of Japan; (b) is from Bagan and Yamagata (2012); and (d) is from the Tokyo Metropolitan Development Bureau].

Table A1. Descriptions of three types of trapezoidal membership functions.

Form	Function
	${\mu }_{C_{k,j}}\left(g_{i,j}\right)=\left\{\begin{array}{c}\hspace{0.17em}0,g_{i,j}\le a,\\ \frac{g_{i,j}-a}{b-a},a<g_{i,j}<b,\\ \hspace{0.17em}1,b\le g_{i,j}\le c,\\ \frac{d-g_{i,j}}{d-c},c<g_{i,j}<d,\\ \hspace{0.17em}0,g_{i,j}\ge d.\end{array}\right.$
	${\mu }_{C_{k,j}}\left(g_{i,j}\right)=\left\{\begin{array}{c}\hspace{0.17em}1,g_{i,j}<a,\\ \frac{b-g_{i,j}}{b-a},a\le g_{i,j}<b,\\ \hspace{0.17em}0,g_{i,j}\ge b.\end{array}\right.$
	${\mu }_{C_{k,j}}\left(g_{i,j}\right)=\left\{\begin{array}{c}\hspace{0.17em}0,g_{i,j}\le a,\\ \frac{g_{i,j}-a}{b-a},a<g_{i,j}\le b,\\ \hspace{0.17em}1,g_{i,j}>b.\end{array}\right.$

Figure 133. Schematic illustration of the calculations of the thickness (red). A single grid cell (light blue) was intersected with the ground mask polygon (light gray).

Figure 144. LCZ maps at grid-cell scales of 200 m, 300 m, 400 m, 500 m, 600 m, 700 m, 800 m, 900 m, and 1000 m derived using S1 (original thresholds). The LCZ class numbers/letters correspond to the classes defined in Figure 2.

Figure 155. LCZ maps at grid-cell scales of 200 m, 300 m, 400 m, 500 m, 600 m, 700 m, 800 m, 900 m, and 1000 m derived using S2 (optimized thresholds). The LCZ class numbers/letters correspond to the classes defined in Figure 2.