MULTIMODAL, MULTITEMPORAL, AND MULTISOURCE GLOBAL DATA FUSION FOR
LOCAL CLIMATE ZONES CLASSIFICATION BASED ON ENSEMBLE LEARNING
Naoto Yokoya∗,†,‡ , Pedram Ghamisi†,‡ , and Junshi Xia∗
∗Department of Advanced Interdisciplinary Studies, the University of Tokyo, Japan
Remote Sensing Technology Institute (IMF), German Aerospace Center (DLR), Germany‡Signal Processing in Earth Observation (SiPEO), Technische Universität München (TUM), Germany†ABSTRACTThis paper presents a new methodology for classification oflocal climate zones based on ensemble learning techniques.
Landsat-8 data and open street map data are used to extract
spectral-spatial features, including spectral reflectance, spectral indexes, and morphological profiles fed to subsequentclassification methods as inputs. Canonical correlation forests
and rotation forests are used for the classification step. The final classification map is generated by majority voting on different classification maps obtained by the two classifiers using
multiple training subsets. The proposed method achieved an
overall accuracy of 74.94% and a kappa coefficient of 0.71 in
the 2017 IEEE GRSS Data Fusion Contest.
Index Terms— Local climate zones (LCZs), canonicalcorrelation forests, rotation forests, morphological profiles.
1. INTRODUCTION
In recent years, there has been a growing interest in the mapping of local climate zones (LCZs) for urban temperaturestudies [1]. LCZs comprises 17 classes based on propertiesof 3D surface structure (e.g., height and density of buildings
and trees) and surface cover (e.g., previous or impervious).
The classification of LCZs is a challenging task due to a large
intra-class variability of spectral signatures caused by theregional variations of vegetation and artificial materials.
The LCZ classification scheme has been successfullydeveloped via the World Urban Database and Access Portal Tools (WUDAPT1 ) to collect information regarding theform and function of cities worldwide [2]. The contributorsof WUDAPT create training data using high-resolution images obtained by Google Earth. The LCZ classification is
performed on Landsat data using random forests [3] in a System developed for Automated Geoscientific Analyses (SAGAGIS). Although the scheme has been well established, it still
This research has been partly supported by Japan Society for the Promotion of Science (JSPS) KAKENHI 15K20955 and 16F16053, and Alexandervon Humboldt Fellowship for postdoctoral researchers.
1 http://www.wudapt.org/

requires an enormous amount of human efforts, particularly
for the specific step of creating training data.
To overcome this limitation, we promote innovative classification algorithms, which have high generalization ability
and transferability. Therefore, the Image Analysis and DataFusion Technical Committee (IADFTC) of the IEEE Geoscience and Remote Sensing Society (GRSS) organized thecontest of the LCZ classification in 2017, aiming at accelerating progress on data fusion methodologies for multitemporal,
multisource, and multimodal remote sensing data.
We present here a new multimodal data fusion methodology for the LCZ classification based on ensemble learningtechniques. Canonical correlation forests (CCFs) [4] and rotation forests (RoFs) [5] are applied to spectral-spatial features
extracted from satellite images and semantic layers. Our results achieved an overall accuracy of 74.94%, which ranked
first in the contest among more than 800 submissions.
The remainder of the paper is organized as follows. Section II describes the data sets provided by the fusion committee for the contest. Section III introduces our proposed
methodology. Experimental results and discussion are presented in Section IV. Section V concludes the paper.
2. DATA SETS
The data sets comprise multitemporal, multisource, and multimodal data for nine cities: five cities (i.e., Berlin, Rome,Paris, Sao Paulo, and Hong Kong) for training and four cities(i.e., Amsterdam, Chicago, Madrid, and Xi’an) for testing.
For each city, the following three kinds of image data havebeen provided for the contest.
• Landsat data2 : Eight multispectral bands (i.e., visible,short and long infrared wavelengths) resampled at aground sampling distance (GSD) of 100 m with multitemporal (i.e., 2-6) acquisitions.
• Sentinel-2 data: Nine multispectral bands (i.e., visible,vegetation red edges and short infrared wavelengths)
resampled at a 100 m GSD for a single acquisition.
2 Courtesy of the U.S. Geological Survey

Ŷ1
Ŷ2
Ŷ3
Ŷ4
Ŷ5
Ŷ6
Ŷ7
Ŷ8
Ŷ9
Ŷ 10
Ŷ 11
Ŷ 12
Ŷ 13
Ŷ 14
Ŷ 15
Ŷ 16
Ŷ 17(a) Berlin(b) Rome(c) Paris(d) Sao Paulo(e) Hong KongFig. 1. Color composite images (top) and ground truth (bottom) for five training cities: (a) Berlin, (b) Rome, (c) Paris, (d) Sao
Paulo, and (e) Hong Kong.

Table 1. LCZs and numbers of pixels for training and test.
No.

ColorName

Train

Test

12
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17



















Compact high-riseCompact mid-riseCompact low-riseOpen high-riseOpen mid-riseOpen low-rise
Lightweight low-rise
Large low-rise
Sparsely built
Heavy industry
Dense trees
Scattered trees
Bush, scrub
Low plantsBare rock or paved
Bare soil or sand
Water

1642
6103
5738
2098
4759
8891
0
4889
1156
449
17716
2819
1741
14457
323
503
8561

242
4892
1522
2270
2255
8265
0
11230
1072
920
3170
4528
1284
12994
1104
391
4454

• Open Street Map (OSM) data3 : Three layers of “buildings,” “land-use,” and “water” at a 5 m GSD.
For the training cities, ground-truth samples of the LCZ
classes were provided in several areas of each city in the formof raster images at a 100 m GSD. Table 1 summarizes classnames and colors of the LCZ classes and the numbers of pixels for training and test. Fig. 1 shows the color composite
images created from the Landsat-8 data as well as the groundtruth maps for the training cities.
Among the data mentioned above, only Landsat and OSM
data have been used in this work. We used Landsat data but
not Sentinel-2 data to exploit long infrared wavelengths. By
using the available time-series of the Landsat data, one cantake temporal-spectral variability into consideration to effectively train the classifiers. It should be noted that the originalLandsat-8 images at a GSD of 30 m were downloaded via
Amazon Simple Storage Service (Amazon S3) for all images
provided in the contest. All eleven bands were used as input,
which is detailed in the next section.
3 Data c OpenStreetMap contributors, available under the Open Database

Licence http://www.openstreetmap.org/copyright3. METHODOLOGY
Fig. 2 illustrates the flowchart of our algorithm developed forthe classification of LCZs in the framework of the 2017 IEEE
GRSS Data Fusion Contest. The algorithm is mainly composed of four steps: preprocessing, feature extraction, classification, and postprocessing as detailed below. Due to the factthat the number of available data sets and the corresponding
training samples provided by the fusion committee was high,a particular emphasis in the proposed framework was dedicated to fast, automatic, yet effective approaches to achieve
accurate results in an acceptable CPU processing time.
3.1. Preprocessing
Atmospheric correction was performed on the original Landsat8 images using ATCOR-2/3 version 9.0.0 with the haze removal option. The Landsat-8 images were upsampled at aGSD of 10 m using bicubic interpolation. The subareas usedin the contest have been extracted using sub-pixel precision
image matching based on phase correlation at a GSD of 100
m. All OSM images have been normalized between 0 and
1 and spatially downgraded to a GSD of 10 m to reduce thecomputational cost in the subsequent processing steps.
3.2. Feature Extraction
The Landsat-8 images as well as the OSM images have been
investigated to extract input features suitable for the classification of LCZs. A total number of 44 features, which arecomposed of spectral reflectance, spectral indexes, and spatialfeatures, have been extracted at a GSD of 100 m as describedin detail below:• Mean and standard deviation have been calculated foreach patch of 10 × 10 pixels for all bands of the 10
m-GSD Landsat-8 data (22 features).

Original
Landsat-8

Atmospheric
Correction &Upsampling

10 m GSD
Landsat-8

CCFs
FeatureExtraction

100 m GSD
Features

Spectral reﬂectance,Indexes, MPs, etc.

OSM

RoFs

SpatialFiltering

Majority
VotingSpatialFiltering

Majority
VotingDecision
FusionOutputFig. 2. Flowchart of the classification algorithm.
• Mean and standard deviation have been calculated foreach patch of 10 × 10 pixels for three indexes obtained
from the 10 m-GSD Landsat-8 data, namely, the normalized difference vegetation index (NDVI), the normalized difference water index (NDWI), and the baresoil index (BSI) (6 features). The advantage of using
these indexes has been already shown for LCZ classification [2].
• Mean has also been calculated for each patch of 10×10pixels for the 10 m-GSD OSM images of “buildings,”“land-use,” and “water” (3 features).
• Morphological profiles (MPs) [6] composed of opening and closing by reconstruction have been used to extract spatial information from the 10 m-GSD NDVI andOSM “building” images. To do so, a circular structuring element with threshold values of 1, 2, and 3 havebeen taken into account (14 features).
3.3. Classification
Two ensemble learning methods, namely CCFs4 [4] andRoFs [5], were used for the classification step. For rotation offeature axes before constructing decision trees, CCFs adopt
canonical correlation analysis between features and labels,
whereas RoFs use random splits of features and principal
component analysis. The superior performance of these twoensemble classifiers have been already proven in the remotesensing community in terms of classification accuracy, CPUprocessing time, and generalization capability [7].
The number of trees was set to 20 for both methods withreference to the studies reported in [7]. 15 different training data sets were prepared. The first ten sets were createdby splitting the whole training data provided by the fusioncommittee into ten subsets whereas the other five sets wereobtained by randomly extracting the same number of trainingsamples (i.e., 500) for all classes. In this way, one can increase the diversity of the forests, which plays an importantrole to boost the classification performance of ensemble andmulti-classifiers. Although some classes include fewer than
500 training pixels, the total number of training samples for
all classes is over 500 since each city has multiple Landsat-8images. It should be noted that some Landsat-8 images that
include clouds were not used for the training step. Both CCFand RoF were built on each training data set, resulting in 30
4 The source code is available at https://bitbucket.org/twgr/ccf

Table 2. OA and kappa for two versions of training data.
Sampling of training data
Uniform
Uniform & imbalance-corrected

OA (%)
70.29
74.94

Kappa
0.66
0.71

different forests in total. Classification was performed on selected cloud-free images for each city using CCFs and RoFs.
3.4. Post-processing
Spatial filtering was applied to all classification maps using a
3 × 3 median filter to reduce the labeling uncertainty and salt
and pepper appearance of labeled pixels. The final classification map was obtained using majority voting on 15 × N classification maps obtained by each ensemble learning method,where N is the number of Landsat-8 images for each city. For
the final classification maps, the results of CCFs were used for
Amsterdam, Chicago, and Xi’an, whereas those of RoFs wereused for Madrid.
4. EXPERIMENTAL RESULTS AND DISCUSSION
Fig. 3 shows the LCZ classification maps for the test cities
obtained by the proposed algorithm. Our result achieved an
overall accuracy (OA) of 74.94% and a kappa coefficient of0.71. Table 2 shows OA and kappa for two versions of ourresults. The first one was obtained by using the first ten training data sets with uniform sampling (i.e., different numbersof training samples for different classes). The second one
demonstrates our final result obtained by merging the resultswith the other five training data sets, which include the same
number of training samples for all classes (i.e., imbalancecorrected sampling). By integrating the classification results
obtained using the additional five subsets of training data, theclassification accuracy of minor classes has been increased,leading to the improvement of the overall result.
Table 3 presents the confusion matrix with the producer’saccuracy (PA) and the user’s accuracy (UA) for the final result to further analyze the performance of our algorithm. The
classification accuracies of classes 6, 8, 11, 14, 16, and 17 are
relatively high, achieving over 75%. Many of those classes
include a high number of training samples as shown in Table 1. For instance, classes 6, 11, 14, and 17 have more than
8000 pixels for training. These results demonstrate that highŶ1
Ŷ2
Ŷ3
Ŷ4
Ŷ5
Ŷ6
Ŷ7
Ŷ8
Ŷ9Ŷ 10Ŷ 11
Ŷ 12
Ŷ 13
Ŷ 14
Ŷ 15
Ŷ 16Ŷ 17(a) Amsterdam(b) Chicago(c) Madrid(d) Xi'anFig. 3. LCZ maps for four test cities: (a) Amsterdam, (b) Chicago, (c) Madrid, and (d) Xi’an.
profiles [8]) will also be an interesting line of research.

Table 3. Confusion matrix.

Prediction

112
34568
910
111213141516
17234568

Ground truth
910111213141516

17167 151 151122002400000102 3033 321 97 166 41
160300000000181 624 125
6230 375 04600084670012 350 18 1243 121 528108115010090469 202 242 1540 1223 740151140111001343 168 250 252 6432 1004 86 85
65 227 12330 2501659 294 149 290 154 89 9055 0 625 2566051 429 15800002013 124 18
1101007130
110001515305 19252103107 16
2300000614 28 20 2690 201 11207
2030220650
700190 3230 131 454 300 0601120766
31002391 940
57110015718215596 929 6
177 233 190 11788 29 26
2800000000000000000016000450
9000178 294
6010301
110002001032 39 4272

UA (%)
43.72
82.44
35.90
65.04
39.67
69.19
79.56
8.61
15.49
84.54
73.78
62.67
86.45
—
65.48
97.33

6. ACKNOWLEDGEMENT
The authors would like to thank the following institutions, initiatives, and companies for acquiring, providing, and preparing the data used in this study: U.S. Geological Survey, European Space Agency, OpenStreetMap, World Urban Databaseand Access Portal Tools (WUDAPT), GeoWIKI, IEEE GRSS
Data and Algorithm Standard Evaluation website, Ticinum
Aerospace S.r.l. (Italy), and IEEE GRSS IADTC.

PA (%) 69.01 62.00 40.65 54.76 68.29 77.82 80.63 1.68 10.00 84.86 71.33 73.21 90.72 0.00 75.19 95.91

7. REFERENCES

generalization ability and transferability of classifiers can beobtained by feeding sufficient training data to the algorithm.
On the other hand, the classification accuracies of classes
3, 4, 9, 10, and 15 are poor (less than 60%). There are twodifferent reasons: (1) For classes 10 and 15, the number oftraining samples is relatively limited, which makes it impossible to learn intra-class variability in an effective way. As aresult, it is difficult to increase the generalization ability andtransferability of the classifiers. (2) Classes 3, 4, and 9 are
misclassified into classes whose spectral-spatial features aresimilar. For instance, the following class pairs are confused:classes 2 and 3, classes 4 and 5, classes 9 and 14. To distinguish these pairs, a digital surface model is a key data source.

[1] I. D. Stewart and T. R. Oke, “Local climate zones for urban temperature studies,” Bulletin of the American Meteorological Society, vol. 93, no. 12, pp. 1879–1900, 2012.

5. CONCLUSION[6] J. A. Benediktsson, J. A. Palmason, and J. R. Sveinsson, “Classification of hyperspectral data from urban areas based on
extended morphological profiles,” IEEE Trans. Geos. RemoteSens., vol. 43, no. 3, pp. 480–491, 2005.

We proposed a new methodology for LCZ classification basedon ensemble learning techniques. CCFs and RoFs are applied
to classify spectral-spatial features (e.g., spectral reflectance,
spectral indexes, and MPs) extracted from the Landsat-8 andOSM data. Our algorithm achieved 74.94% OA in the 2017
IEEE GRSS Data Fusion Contest.
Our future work will be on analysis of feature importance. Further study on investigating the impact of using
recent advances in spatial feature extraction (e.g., extinction[2] B. Bechtel, P. J. Alexander, J. Böhner, J. Ching, O. Conrad,
J. Feddema, G. Mills, L. See, and I. Stewart, “Mapping local
climate zones for a worldwide database of the form and function of cities,” ISPRS Inter. Jour. Geo-Information, vol. 4, no. 1,pp. 199–219, 2015.
[3] L. Breiman, “Random forests,” Mach. Learn., vol. 45, no. 1,pp. 5–32, 2001.
[4] T. Rainforth and F. Wood, “Canonical correlation forests,” ArXiv
Preprint ArXiv:1507.05444.
[5] J. J. Rodriguez and L. I. Kuncheva, “Rotation forest: A new classifier ensemble method.,” IEEE Trans. Pat. Anal. Mach. Intel.,vol. 28, no. 10, pp. 1619–1630, 2009.

[7] J. Xia, N. Yokoya, and A. Iwasaki, “Hyperspectral image classification with canonical correlation forests,” IEEE Trans. Geos.
Remote Sens., vol. 55, no. 1, pp. 421–431, 2017.
[8] P. Ghamisi, R. Souza, J. A. Beneiktsson, X. X. Zhu, L. Rittner,and R. Lotufo, “Extinction profiles for the classification of remote sensing data,” IEEE Trans. Geos. Remote Sens., vol. 54,
no. 10, pp. 5631–5645, 2016.