Multitemporal Land Use and Land Cover Classification from Time-Series Landsat Datasets Using Harmonic Analysis with a Minimum Spectral Distance Algorithm

Abstract

An understanding of historical and present land use and land cover (LULC) information and its changes, such as urbanization and urban growth, is critical for city planners, land managers and resource managers in any rapidly changing landscape. To deal with this situation, the development of a new supervised classification method for multitemporal LULC mapping with long-term reliable information is necessary. The ultimate goal of this study was to develop a new classification method using harmonic analysis with a minimum spectral distance algorithm for multitemporal LULC mapping. Here, the Jiangning District of Nanjing City, Jiangsu Province, China was chosen as the study area. The research methodology consisted of two main components: (1) Landsat data selection and time-series spectral reflectance reconstruction and (2) multitemporal LULC classification using HA with a minimum spectral distance algorithm. The results revealed that the overall accuracy and Kappa hat coefficients of the four LULC maps in 2000, 2006, 2011, and 2017 were 97.03%, 90.25%, 91.19%, 86.32% and 95.35%, 84.48%, 86.74%, 80.24%, respectively. Further, the average producer accuracy and user accuracy of the urban and built-up land, agricultural land, forest land, and water bodies from the four LULC maps were 92.30%, 90.98%, 94.80%, 85.65% and 90.28%, 93.17%, 84.40%, 99.50%, respectively. Consequently, it can be concluded that the newly developed supervised classification method using harmonic analysis with a minimum spectral distance algorithm can efficiently classify multitemporal LULC maps.

1. Introduction

By definition, “land use refers to what people do on the land surface (e.g., agriculture, commerce, settlement) while land cover refers to the type of material present on the landscape (e.g., water, crops, forest, wetland, humanmade materials such as asphalt)” [1]. At present, the speed, magnitude, and scale of human activities in changing the Earth’s landscape are the most significant in human history. In particular, urban development leads to the replacement of natural landscapes with impermeable materials, alteration of the biophysical environment, and transformation of the land surface energy process [2]. An understanding of historical and current land use and land cover (LULC) status and change information, urbanization, and urban growth are critical for city planners, land managers, and resource managers in any rapidly changing environment because LULC changes will cause alterations in environmental conditions [3,4,5]. A modification in LULC leads to alterations in nature, such as green cover destruction and water resource pollution [6]. When a LULC change occurs due to urbanization along a city boundary, it increases the size of the city [7]. Accelerated urban growth and LULC changes exert pressure on the natural environment and human welfare and have become a global concern [8]. Several studies on LULC changes and their impacts have been conducted worldwide from multiple dimensions using satellite remote sensing and GIS technology [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]. All of these studies require time-series datasets that are mostly derived from Earth observation satellites to classify multitemporal LULC maps.

In principle, LULC maps that are classified using multiple single-temporal dates frequently lack deficiency information about the LULC change process and its transformations (e.g., its rate of change). On the contrary, multitemporal LULC mapping can overcome these shortcomings [24,25]. At the same time, the availability of Earth observation satellite data with an extensive range of temporal and spatial resolutions continues to increase, which guides the advanced methods to produce frequent and accurate LULC maps. Still, it is unlikely that there is only one method for LULC classification using time-series satellite data [25]. Gómez et al. [25] reported that there are two approaches for generating multitemporal LULC information from time-series datasets. Firstly, independent LULC maps are produced for each time interval (i.e., for every month or year) [26,27]. Secondly, multitemporal LULC maps are created using a base map as a reference condition with change information for multiple dates [28,29].

In recent decades, available time-series remotely sensed datasets, such as the NOAA-Advanced Very High Resolution Radiometer (AVHRR) [30,31,32,33,34], Terra/Aqua-Moderate Resolution Imaging Spectroradiometer (MODIS) [35,36,37,38,39], SPOT Vegetation [40,41,42], Landsat [43,44,45,46], and Sentinel-2 [47,48] datasets, have been increasingly used by remote sensing scientists and researchers in multitemporal LULC classification and mapping, LULC change detection, and vegetation monitoring through adopting harmonic wave models with novel classification methods. In particular, Landsat datasets have become a promising resource for multitemporal classification and change detection because they have a long-time period (over 40 years), better spatial resolution (30 meters), better spectral resolution (3 visible, 1 near infrared (NIR) and 2 shortwave infrared (SWIR) bands), and shorter revisit cycle (an 8 day revisit cycle when different Landsat sensors overlap) [49]. These studies can provide better characteristics of LULC type with dynamic changes for any given time. For example, Zhu and Woodcock [45] developed a Continuous Change Detection and Classification (CCDC) algorithm for change detection and classification.

Meanwhile, numerous time-series reconstruction models have been developed in different application fields and regions [50,51,52,53,54]. A series of algorithms based on discrete Fourier analysis [55,56,57], also known as Harmonic Analysis (HA), which is one of the most extensively used models to fit time-series data because it can identify and remove noise (e.g., cloud, cloud shadow, snow, and Landsat 7 strips) [58,59] and fit a time-series curve with fewer input data by temporal interpolation [60]. For example, Weng, Lu, and Schubring [58] examined the capability of HA for cloud removal and reconstructed a 10-day composite AVHRR data. Julien, Sobrino, and Verhoef [31] assessed vegetation changes by combining LST and NDVI data with the removal of clouds. Additionally, frequency domain information (e.g., the magnitude and phase of harmonic components) obtained from HA can be applied to quantify vegetation phenology and classify cropland [32,33,34,61].

To produce highly accurate LULC maps for any given time similar to the CCDC approach [45], we propose a new supervised multitemporal LULC classification workflow using harmonic analysis and a minimum spectral distance algorithm with time-series Landsat 5, 7, and 8 datasets between 2000 and 2017. The specific research objectives were (1) to retrieve available Landsat data between 2000 and 2017, to reconstruct time-series spectral reflectance using the time-series model, and (2) to develop a supervised classification method using harmonic analysis with a minimum spectral distance algorithm for multitemporal LULC mapping. The developed method was applied to classify multitemporal LULC maps from 2000, 2006, 2011, and 2017 for the Jiangning District of Nanjing City, China.

2. Materials and Methods

2.1. Study Area

Jiangning, as the largest district among the eleven districts of Nanjing City, Jiangsu Province, China was chosen as the study area. Jiangning covers an area of 1587 km² or about 24% of Nanjing City (6,598 km²) and is situated between 118°28’-119°7’E and 31°37’-32°7’N. The climate of Jiangning is humid and subtropical, and its landform is exemplified by hillocks, low mountains, hills, plains, rivers, and lakes (Figure 1).

The selected study area represents the rapid urbanization and economic growth of Nanjing City. According to reports of the Ministry of Housing and Urban-Rural Development (MOHURD) of the People’s Republic of China [62] and the Nanjing Municipal Bureau of Statistics [63], the population in Nanjing City increased from 3.1 million persons in 2000 to 5.9 million persons in 2016, and the urban and built-up area expanded from 201.40 km² to 773.79 km² in the same period. Meanwhile, the gross regional product of the city has grown from 107.35 Billion Yuan in 2000 to 1,015.30 Billion Yuan in 2016. These statistics indicate the rapid urbanization and economic growth of Nanjing City. Consequently, many agricultural lands in the city were transformed into urban and built-up land due to population growth and rapid urbanization.

2.2. Datasets

All available Landsat 5/7/8 images of Level-1 products for Worldwide Reference System-2 with Path 120 and Row 38 between 2000 and 2017, with a total of 673 scenes, were downloaded through the USGS online portal (www.earthexplore.usgs), as shown in Figure 2. A total of 229, 340, and 104 scenes were taken from Landsat 5, 7, and 8, respectively. Further, very high spatial resolution images from Google Earth for 2000, 2006, 2011, and 2017, used as the ground reference data, were downloaded for the training area selection via a maximum likelihood classifier (MLC) and an accuracy assessment of multitemporal LULC maps using harmonic analysis with a minimum spectral distance algorithm (Figure 3).

2.3. Research Methodology

The research methodology workflow (input, process, and output) consisted of two main components: (1) Landsat data selection and time-series spectral reflectance reconstruction and (2) multitemporal LULC classification using HA with a minimum spectral distance algorithm (Figure 4).

2.3.1. Landsat Data Selection and Time-Series Spectral Reflectance Reconstruction

Landsat data selection and time-series spectral reflectance reconstruction were implemented in four separate steps, as follows:

Cloud Cover Assessment

In this step, we considered the percent of total images and the percent of the total clearly observed pixels by percent of cloud cover to reduce the number of scenes and costs and to retain the highest percentage of total clearly observed pixels for the data analysis.

In practice, the metadata for cloud cover percentage from the USGS web portal was observed, recorded, and calculated for the available Landsat 5/7/8 images between 2000 and 2017. The percentage of total images and total clearly observed pixels, according to 10 interval classes (≤10, ≤20, ≤30, ≤40, ≤50, ≤60, ≤70, ≤80, ≤90, and ≤100) of cloud cover (%), was calculated using Equations (1) and (2):

$$P_{image,i}=\frac{N_{image,i}}{N_{image}}*100\%$$

(1)

$$P_{pixel,j}=\frac{{{\displaystyle \sum }}_0^j{N}_{image,j}*\left(1-j\right)*M}{N_{image}*M}*100\%$$

(2)

where i is the cloud cover (%) that is less than or equal to the interval class (10, 20, …, 100), $P_{image,i}$ is the percentage of total images, $N_{image,i}$ is the number of images, and $N_{image}$ is the number of total images. j is the cloud cover (%), $P_{pixel,j}$ is the percentage of total clearly observed pixels, and $M$ is the number of pixels in the image.

Clearly Observed and Contaminated Pixel Recognition

In this step, clearly observed and contaminated pixels of Landsat images were first recognized by using the QA band, and contaminated pixels were then identified using the HA model.

(1) Clearly observed and contaminated pixel recognition based on the QA band.

Typically, each pixel in the QA band contains information related to the terrain, radiometric saturation, cloud, and cloud shadow. Users can apply the information from the QA band of Landsat 5, 7, and 8 data products to recognize the clearly observed and contaminated pixels using Landsat QA tools [64].

In this study, per-pixel filtering with the Landsat QA tools of USGS [64] was applied to examine the contaminated pixel coverage according to information on the QA band. If the coverage of the contaminated pixels of an image in the study area was greater than 90%, the image was removed.

(2) Contaminated pixel recognition base on the HA model.

The developed algorithm by Zhu and Woodcock [45] was adopted to recognize contaminated pixels in this study. In practice, the surface reflectance values of the green band and SWIR1 band of the Landsat data were first transformed into a time-series model (HA), and the actual Landsat observations were then compared with the corresponding derived harmonic waves to recognize the contaminated pixels using the two conditions in Equations (3) and (4), as suggested by Zhu and Woodcock [45].

$$\rho \left(blue, x\right)-\widehat{\rho }{\left(blue, x\right)}_{RIRLS}>0.04$$

(3)

$$\rho \left(SWIR1, x\right)-\widehat{\rho }{\left(SWIR1, x\right)}_{RIRLS}<-0.04$$

(4)

where x is the Julian date, ρ(blue, x), and ρ(SWIR1, x) are the observed values of Landsat at Julian date x, and $\hat{\rho }$(blue, x)_RIRLS and $\hat{\rho }$(SWIR1,x)_RIRLS are the predicted values of Landsat at Julian date x.

Conversion of Digital Numbers to Spectral Reflectance

Radiometric correction is a crucial step to ensure the homogeneity of the data for change detection [65]. In operation, the USGS uses the LEDAPS and the LaSRC systems to generate standard surface reflectance products [66,67].

In this study, the Level 1 product of Landsat was converted to its spectral reflectance using Equation (5), as suggested by [66,67].

$$\rho {\lambda }^{\prime }=M_{\rho }*Q_{cal}+A_{\rho }$$

(5)

where $\rho {\lambda }^{\prime }$ is the spectral reflectance, M_ρ is the reflectance multiplicative scaling factor, A_ρ is the reflectance additive scaling factor, and Q_cal  is the digital number value.

Time-Series Spectral Reflectance Reconstruction

The clearly observed and contaminated pixels of the Landsat images using the QA band and HA model (with value 0 or 1) from Step (2) were used to create a spatiotemporal cube using the reshape function in the MATLAB software. Likewise, the derived six spectral reflectance bands (with values 0 to 1) from Step (3) were also used to create a spectral reflectance spatiotemporal cube in the same manner. Then, information on the clearly observed and contaminated pixels and the spectral reflectance of the six bands were converted from the spatial dimension to the time dimension using the reshape function (the data were read pixel by pixel from the spatiotemporal cubes) in the MATLAB software. After that, two-time dimension datasets were combined to remove the contaminated pixels for the time-series spectral reflectance reconstruction (see Figure 5). The derived time-series spectral reflectance data will be further applied to classify and map multitemporal LULC.

2.3.2. Multitemporal LULC Classification Using HA with the Minimum Spectral Distance Algorithm

Under this component, five main steps were separately implemented, as follows.

Stable Pixels of the LULC Type Extraction

The different LULC types (urban and built-up land (U), agricultural land (A), forest land (F), and water bodies (W) for 2000, 2006, 2011, and 2017) were first classified using the MLC in the ERDAS Imagine software. In practice, the common training areas of four LULC types from different years were first selected using visual interpretation with ancillary data from Google Earth. Then, the selected training areas were separately applied to classify the LULC type from different Landsat images using the MLC. Then, four classified LULC maps (2000, 2006, 2011, and 2017) were simultaneously superimposed to identify the common areas of each LULC type from four different years via an overlay analysis using the ESRI ArcMap software. The derived output presents the stable pixels of each LULC category between 2000 and 2017, as well as the LULC changes during this period.

Next, sample points from each stable LULC type were randomly selected by using the create random points function in the ESRI ArcMap software. The sample points of each LULC type were further transformed into harmonic function curves and used to construct standard harmonic curves for the six spectral bands in the next step.

Harmonic Function Curve Transformation and Standard Harmonic Curve Construction

The time-series curve of each data point is expressed as the sum of a series of cosine or sine waves; each wave is determined by a different amplitude and phase [68]. These continuous amplitudes and phases are summed to produce a complex curve [32]. Since the cosine or sine wave is a typical periodic change curve, HA can be used to simulate the periodic change of spectral reflectivity.

To obtain the harmonic curve characteristic of different LULC types, the selected stable pixels for each LULC type between 2000 and 2017 were first transformed into a spectral harmonic curve using Equation (6), modified from Zhu and Woodcock [45].

$${\mathrm{y}}_i=a_i+b_{i}t+A_i cos\left(\frac{2\pi }{T}t-{\phi }_i\right)$$

(6)

where t is the Julian date, i is the i^th Landsat band, T is the number of days per year (T= 365), a_i is the coefficient of the intercept value, b_i is the coefficient of the slope value, A_i is the coefficient of the amplitude value, ${\phi }_i$ is the coefficient of the phase value, and y is the reconstructed reflectance value at Julian date t. In this study, typical spectral bands among Landsat 5, 7, and 8 (namely Blue, Green, Red, NIR, SWIR1, and SWIR2) were selected for harmonic curve function transformation.

Next, the median values of the four fitted coefficients from all sample points of each LULC type for the six spectral bands were extracted, and the derived median values were then used to construct the standard harmonic function curve of each LULC type for the six spectral bands using Equation (6).

Spectral Distance Measurement and Probability Calculation

In principle, spectral distance can be measured by simple subtraction of the spectral reflectance value between the standard harmonic function curve of each LULC type (U, A, F, and W) and an unclassified pixel at a specific time point. However, the spectral distance values of four LULC types for an unclassified pixel have different range values over time and cannot be directly applied to compare different LULC categories for a specific time point due to their different scales. Consequently, the spectral distance value of each LULC type was normalized into the range [0,1].

In practice, Equations (7) and (8) were first applied to calculate the maximum and minimum spectral distance between the standard harmonic function curves of each LULC type (U, A, F, and W) and an unclassified pixel at each specific time point. Then, the maximum and minimum values from any LULC type at the same specific time point were further used to calculate the normalized spectral distance between the standard harmonic function curve of each LULC type (U, A, F, and W) and an unclassified pixel at the same specific time point using Equations (9)–(12), respectively. In this way, the lowest normalized spectral distance between any LULC type and an unclassified pixel (i.e., when the spectral distance equals 0) will provide the highest likelihood of being a specific LULC type (i.e., the probability equals 1).

$$MA{X}_{ij}=max\left(\left|R_{Uij}-R_{X_{ij}}\right|,\left|R_{Aij}-R_{X_{ij}}\right|,\left|R_{Fij}-R_{X_{ij}}\right|,\left|R_{Wij}-R_{X_{ij}}\right|\right)$$

(7)

$$MI{N}_{ij}=min\left(\left|R_{Uij}-R_{X_{ij}}\right|,\left|R_{Aij}-R_{X_{ij}}\right|,\left|R_{Fij}-R_{X_{ij}}\right|,\left|R_{Wij}-R_{X_{ij}}\right|\right)$$

(8)

$$Normalized spectral distanc{e}_{Uij}=\frac{\left|R_{Uij}-R_{X_{ij}}\right|-MI{N}_{ij}}{MA{X}_{ij}-MI{N}_{ij}}$$

(9)

$$Normalized spectral distanc{e}_{Aij}=\frac{\left|R_{Aij}-R_{X_{ij}}\right|-MI{N}_{ij}}{MA{X}_{ij}-MI{N}_{ij}}$$

(10)

$$Normalized spectral distanc{e}_{Fij}=\frac{\left|R_{Fij}-R_{X_{ij}}\right|-MI{N}_{ij}}{MA{X}_{ij}-MI{N}_{ij}}$$

(11)

$$Normalized spectral distanc{e}_{Wij}=\frac{\left|R_{Wij}-R_{X_{ij}}\right|-MI{N}_{ij}}{MA{X}_{ij}-MI{N}_{ij}}$$

(12)

where $R_{Uij}$, $R_{Aij}$, $R_{Fij}$, and $R_{Wij}$ are the spectral reflectance from the standard harmonic function curves of U, A, F, and W, respectively, at time point i using band j. $R_{X_{ij}}$ is the spectral reflectance of an unclassified pixel at time point i using band j.

Next, the normalized spectral distance of an unclassified pixel to any LULC type was applied to calculate the probability of an unclassified pixel being a specific LULC type (U, A, F, and W) using Equations (13)–(16), respectively.

$$Probabilit{y}_{Uij}=1-Normalized spectral distanc{e}_{Uij}$$

(13)

$$Probabilit{y}_{Aij}=1-Normalized spectral distanc{e}_{Aij}$$

(14)

$$Probabilit{y}_{Fij}=1-Normalized spectral distanc{e}_{Fij}$$

(15)

$$Probabilit{y}_{Wij}=1-Normalized spectral distanc{e}_{Wij}$$

(16)

where $Probabilit{y}_{Uij}$, $Probabilit{y}_{Aij}$, $Probabilit{y}_{Fij}$ and $Probabilit{y}_{Wij}$ are the probability of an unclassified pixel being U, A, F, and W, respectively, at time point i using band j.

Multitemporal LULC Classification

In principle, the probability of an unclassified pixel being a specific LULC type from Equations (13)–(16) at a specific time point from one spectral band can be separately applied for multitemporal LULC classification.

In this study, all six standard spectral bands of Landsat 5, 7, and 8, including the blue, green, red, NIR, SWIR1, and SWIR2 bands, were applied for the multitemporal LULC classification. In this way, the average probabilities of an unclassified pixel being a specific LULC type (U, A, F, and W) from among the six spectral bands were first separately calculated using Equations (17)–(20). Then, the average probabilities of an unclassified pixel being a specific LULC type (U, A, F, and W) from among the four LULC types were compared to identify the highest value; the corresponding LULC type that provides the highest probability was then assigned to an unclassified pixel at a specific time point.

$$Average probabilit{y}_{Ui}={\displaystyle \sum }_{j=1}^{6}Probabilit{y}_{Uij}/6$$

(17)

$$Average probabilit{y}_{Ai}={\displaystyle \sum }_{j=1}^{6}Probabilit{y}_{Aij}/6$$

(18)

$$Average probabilit{y}_{Fi}={\displaystyle \sum }_{j=1}^{6}Probabilit{y}_{Fij}/6$$

(19)

$$Average probabilit{y}_{Wi}={\displaystyle \sum }_{j=1}^{6}Probabilit{y}_{Wij}/6$$

(20)

where, $Average probabilit{y}_{Ui}$, $Average probabilit{y}_{Ai}$, $Average probabilit{y}_{Fi}$, and $Average probabilit{y}_{Wi}$ are the probability of an unclassified pixel being U, A, F, and W, respectively, at time point i using six bands.

After multitemporal LULC classification based on average probability, post-classification processing was applied to remove unexpected errors because differences in the acquisition and atmospheric conditions of time-series datasets can create spectral shifts [69]. In this study, the mode function with a moving window (1 x 9 in size) of spatiotemporal filtering under the MATLAB software was used to eliminate the unexpected errors.

Accuracy Assessment

Accuracy assessment of LULC maps has always been a difficult task because of the lack of large-scale, high temporal, and spatial resolution maps for reference. For the study of change detection using time-series Landsat data, the best verification data are the Landsat data themselves, as suggested in [70]. In this study, very high spatial resolution images from Google Earth (see Figure 3) and Landsat images were used as ground reference information for accuracy assessment. Meanwhile, the number of sample points was estimated based on a multinomial distribution, and a stratified random sampling technique was applied to allocate the sample points for the ground reference test information, as suggested in [71,72]. In this way, the overall accuracy (OA), producer’s accuracy (PA), user’s accuracy (UA), and Kappa hat coefficient (KHAT) were reported.

Furthermore, the four stable LULC types (U, A, F, and W) and the LULC change areas from the classified LULC maps in 2000, 2006, and 2011, and 2017 were reclassified into two groups: changed and unchanged areas for change detection accuracy assessment using OA, PA, and UA, as suggested by Congalton and Green [71].

3. Results

In this research, four Landsat images for 2000 (12 June 2000), 2006 (31 July 2006), 2011 (29 July 2011), and 2017 (18 May 2017), from time-series Landsat 5, 7, and 8 datasets between 2000 and 2017, were selected to classify multitemporal LULC data using HA with a minimum spectral distance algorithm under the MATLAB software environment.

3.1. Selection of Landsat Data According to Cloud Cover Assessment

Figure 6 shows the percentage of images and clearly observed pixels as a cumulative histogram over the Nanjing City scene (Path 120 and Row 38) from the available Landsat 5/7/8 datasets between 2000 and 2017. Based on this information, if images with cloud cover ≤10% are selected, about 26% of the total images (673 scenes) can be used. These images include about 50% clearly observed pixels; therefore, about 50% of clearly observed pixels are ignored.

In this study, we chose all images with a cloud cover ≤90%, or about 75% of the total images (673 scenes), since these images include about 99% clearly observed pixels, while only about 1% of clearly observed pixels are ignored. This selection can decrease the number of scenes from 673 scenes (or 100%) to 505 scenes (or 75%) of the total scenes. However, seven scenes from Landsat 5 and three scenes from Landsat 8 were removed due to poor image quality and incorrect information on cloud cover. Subsequently, 168 images from Landsat 5, 251 images from Landsat 7, and 76 images from Landsat 8, with a total of 495 images, were selected to identify clearly observed and contaminated pixels in the study area based on the QA band and HA model.

3.2. Reselection of Landsat Data Using the QA Band and HA Model

After cloud coverage assessment, images with less than or equal to 90% contaminated pixels in the study area were reselected using the QA band and HA model. As a result, 121 images from Landsat 5, 201 images from Landsat 7, and 66 images from Landsat 8, with a total of 388 images, were finally selected and used in this study. Additionally, the digital values of the clearly observed and contaminated pixels of all chosen images are 0 and 1, respectively. Examples of the distribution of clearly observed and contaminated pixels from Landsat 5 and 7 are displayed in Figure 7. In Figure 7a, most of the contaminated pixels are cloud and cloud shadow, whereas the contaminated pixels in Figure 7b include cloud, cloud shadow, and gaps.

3.3. Spectral Reflectance Data

The final selection of Landsat data, with a total of 388 scenes, was converted to spectral reflectance for data analysis. Examples of the spectral reflectance data from Landsat 5 and 7 are displayed in Figure 8. As a result, the existing gaps in the Landsat 7 data, as shown in Figure 8b, are considered here as contaminated pixels and are ignored in the data analysis.

3.4. Time-Series Spectral Reflectance Reconstruction

By combining the spatiotemporal cube of the recognized clearly observed and contaminated pixels (with value 0 or 1) with the spatiotemporal cube of the derived spectral reflectance data of six bands (with values from 0 to 1) using the reshape function of the MATLAB software, the contaminated pixels are removed from original time-series spectral reflectance data of the six bands.

Figure 9a shows the original time-series spectral reflectance data of the NIR band (with a value from 0 to 1) from one pixel among the 388 dates. This observation reveals that some pixel values deviate from their normal reflectance range due to contamination, such as cloud/cloud shadows. Meanwhile, Figure 9b shows the time-series spectral reflectance data of the NIR band from the same pixel after the removal of contaminated values on some dates. This reveals that the reflectance values oscillate up and down with temporal variations. In this study, the reconstructed time-series spectral reflectance data of 388 images serve as the primary input data for multitemporal LULC mapping using HA with the minimum spectral distance algorithm.

3.5. Stable Area of LULC Type Extraction

The stable areas of four LULC types and change areas between 2000 and 2017 from four LULC maps in 2000, 2006, 2011, and 2017 are summarized in

Table 1. Area and percentage of the stable land use and land cover (LULC) type and its changes between 2000 and 2017 using maximum likelihood classifier (MLC) algorithm.

LULC Type	Area in km2	Percent (%)
1. Stable urban and built-up land	62.05	3.91
2. Stable agriculture land	335.4	21.14
3. Stable forest land	129.59	8.17
4. Stable water bodies	49.82	3.14
5. LULC change between 2000 and 2017	1,009.87	63.64
Total	1586.74	100

and displayed in Figure 10. As a result, the stable areas of U, A, F, and W during 2000 to 2017 are 62.05 km² (3.91%), 335.40 km² (21.14%), 129.59 km² (8.17%), and 49.82 km² (3.14%), respectively. Most of the stable areas of urban and built-up land are located in an old core area of the Jiangning district, while stable areas of agricultural land are randomly distributed in lowland areas far from the core area of the district. Meanwhile, most of the stable areas of forest land are situated in hilly and mountainous regions. In contrast, most of the stable areas of water bodies are located on the Yangtze River and the Qinhuai River.

On the other hand, the unstable area of LULC types was 1,009.87 km² (63.64%). This finding indicates that LULC changes during the three studied periods (2000–2006, 2006–2011, and 2011–2017). The unstable areas of LULC types spread outward from the old core area of Jiangning through the agricultural land.

The distribution of random sampling points of U (9,395 pixels), A (9,893 pixels), F (9,711 pixels), and W (9,181 pixels) from stable LULC areas between 2000 and 2017 is displayed in Figure 11.

3.6. Spectral Harmonic Function Curve of LULC Types

The spectral harmonic function curves of four LULC types that were transformed from the sample points of the stable pixels between 2000 and 2017 are displayed in Figure 12. Since the harmonic function curves of each LULC type come from nearly 10,000 sample points, the harmonic function curves of each LULC type from each spectral band appear in a ribbon form, and the harmonic function curves of different LULC types overlap among the LULC types. For example, the harmonic function curve of forest land overlaps with that of water bodies in the NIR band, while the harmonic curve of the water bodies overlaps with that of the other LULC types in the green band. Thus, it remains difficult to determine the similarities among LULC types and to compare the derived harmonic function curve with the spectral polyline of an unclassified pixel. Consequently, it is necessary to find a standard harmonic function curve that can represent different LULC types in various bands.

3.7. Standard Harmonic Function Curve of Each LULC Type

Table 2. Median values for the coefficients of the samples of the four LULC types for the six spectral bands.

Spectral Features	LULC Type	Coefficients
		Intercept	Slope (E-04)	Amplitude	Phase
BLUE	U	0.1143	0.0114	™0.0321	0.4298
	A	0.0996	0.0092	™0.0225	0.5379
	F	0.0916	0.0063	™0.0210	0.4925
	W	0.1103	0.0136	™0.0287	0.3571
GREEN	U	0.1039	0.0020	™0.0325	0.4572
	A	0.0881	0.0045	™0.0228	0.6110
	F	0.0753	0.0027	™0.0206	0.5527
	W	0.1028	0.0105	™0.0293	0.3148
RED	U	0.1020	™0.0014	™0.0314	0.5205
	A	0.0777	0.0024	™0.0163	0.9323
	F	0.0597	™0.0014	™0.0126	0.8894
	W	0.1001	™0.0019	™0.0305	0.1774
NIR	U	0.1368	0.0161	™0.0626	0.3530
	A	0.1613	0.0515	™0.0836	0.2996
	F	0.1548	0.0619	™0.1045	0.3088
	W	0.0685	™0.0214	™0.025	0.2702
SWIR1	U	0.1250	0.0102	™0.0512	0.4175
	A	0.1103	0.0387	™0.0360	0.4323
	F	0.1032	0.0230	™0.0426	0.4420
	W	0.0198	0.0003	™0.0089	0.7526
SWIR2	U	0.0947	0.0083	™0.0364	0.5518
	A	0.0641	0.0200	™0.0166	0.9975
	F	0.0500	0.0088	™0.0151	1.0717
	W	0.0122	0.0002	™0.0059	0.8611

shows the median values of the four fitted coefficients of each LULC type for the six spectral bands between 2000 and 2017, while the constructed standard harmonic function curves of four LULC types from six spectral bands are displayed in Figure 13. These data reveal that there are no overlapping strips among LULC types, making it is easy to differentiate the harmonic function curves of LULC types.

3.8. Spectral Distance Measurement

Figure 14 shows the spectral distance between the standard harmonic function curves of the four LULC types and the polyline of one unclassified pixel for similarity visualization. As a result, it is easy to observe that the spectral reflectance of the unclassified pixel as a polyline (black color) is similar to the standard harmonic function curves of agricultural land and forest land before 2008. However, this curve is also similar to the standard harmonic function curves of urban and built-up land after 2008. As a result, the LULC type of the unclassified pixel between 2000 and 2008 should be agricultural land or forest land that was transformed into urban and built-up land after 2008.

3.9. Probability of an Unclassified Pixel Being a Specific LULC Type

The probability of an unclassified pixel being a specific LULC type (U, A, F, and W) between 2000 and 2017, based on the six spectral bands, is displayed in Figure 15. This figure shows the variation in the probability of an unclassified pixel being a specific LULC type between 2000 and 2017 and indicates the probability of an unknowing pixel being a specific LULC type at a particular time point. For example, if we select the blue band to classify multitemporal LULC data, most of the LULC types before 2000 were forest land or agricultural land, but most of the LULC types after 2008 were urban and built-up land.

Meanwhile, the average probability of an unclassified pixel being a specific LULC type (U, A, F, and W) between 2000 and 2017, based on the six spectral bands, is presented in Figure 16. This figure shows the variation of the average probability of an unclassified pixel for each LULC type from the six spectral bands between 2000 and 2017 at a specific time point. As a result, the unclassified pixel belongs to different LULC types at different time points. Before 2008, this pixel alternatively shows a high average probability to appear in forest land and agricultural land, but after 2008, this pixel displays a high average probability to appear in urban and built-up land. This finding indicates that the LULC type of this pixel can be classified at a specific time point according to the highest average probability value of the corresponding LULC type.

3.10. Multitemporal LULC Classification and Mapping

Figure 17a shows the LULC type of the unclassified pixel between 2000 and 2017 at specific dates of the 388 Landsat data before post-classification processing. As a result, the LULC type of the pixel alternately belongs to forest land, agricultural land, and urban and built-up land between 2000 and 2008. It also frequently belongs to agricultural land during this period. In contrast, this pixel belongs to urban and built-up land after 2008. However, this pixel belongs to water bodies at three dates (as shown by the red oval). Based on our prior knowledge, it is difficult for the LULC type to undergo two sharp changes in a short period. The LULC type of this pixel at the three-time points indicates unexpected errors like salt and pepper noise [1].

On the contrary, Figure 17b shows the final results after post-classification processing for the unclassified pixel; this figure reveals that the unexpected errors from water bodies in Figure 17a have been removed. Accordingly, the classified LULC type of this pixel was agricultural land before 2008 and was urban and built-up land after 2008.

The results of the selected multitemporal LULC maps in 2000, 2006, 2011, and 2017 from time-series Landsat datasets between 2000 and 2017 using HA with the minimum spectral distance algorithm are presented in Figure 18 and

Table 3. Area and percentage of LULC classification using HA with the minimum spectral distance algorithm for four different years.

LULC Type	2000		2006		2011		2017
	Area (km2)	%	Area (km2)	%	Area (km2)	%	Area (km2)	%
U	147.43	9.29	259.17	16.33	354.22	22.32	444.85	28.04
A	1,145.08	72.17	911.97	57.47	867.33	54.66	691.66	43.59
F	226.15	14.25	350.63	22.10	309.41	19.50	338.78	21.35
W	68.08	4.29	64.97	4.09	55.77	3.52	111.45	7.02
Total	1586.74	100.00	1586.74	100.00	1586.74	100.00	1586.74	100.00

.

As a result, the percentage of urban and built-up land extends from 9.29% in 2000 to 28.04% in 2017, the proportion of forest land increases from 14.25% in 2000 to 21.35% in 2017, and the percentage of water bodies also increases from 4.29% in 2000 to 7.02% in 2017. On the contrary, the percentage of agricultural land decreases from 72.17% in 2000 to 43.59% in 2017. These results indicate a dramatic increase in urban and built-up land due to urbanization in the study area from 2000 to 2017. The expansion of urban and built-up land in the three different periods (2000–2006, 2006–2011, and 2011–2017) is displayed in Figure 19.

This result reveals that areas of urban and built-up land between 2000 and 2006 expanded from the old core Jiangning district, and newly developed areas occurred near the Yangtze River and Lukou international airport, covering an area of 111.74 km². Meanwhile, areas of urban and built-up land between 2006 and 2011 expanded from the old core Jiangning district in the northern part of the district, and newly developed areas appeared near the Yangtze River and Lukou town, covering an area of 95.05 km². In the meantime, areas of urban and built-up land between 2011 and 2017 were expanded by the conversion of agricultural land in urban and built-up land, covering an area of 90.63 km². This finding is consistent with the statistical report of MOHURD [62], who reported that urban and built areas in the entirety of Nanjing City increased from 201.40 km² in 2000 to 773.79 km² in 2016 due to urbanization.

3.11. Accuracy Assessment of LULC Map

Results of the accuracy assessment of LULC maps in 2000, 2006, 2011, and 2017 are displayed in

Table 4. Accuracy assessment of the classification LULC map in 2000 for 12 June 2000.

LULC Type	U	A	F	W	Row Total
U	69	6	0	1	76
A	2	403	2	0	407
F	0	2	98	3	103
W	0	0	0	50	50
Column Total	71	411	100	54	636
PA	97.18%	98.05%	98.00%	92.59%
UA	90.79%	99.02%	95.15%	100.00%
OA	97.03%
KHAT	95.35%

,

Table 5. Accuracy assessment of classification LULC map in 2006 for 31 July 2006.

LULC Type	U	A	F	W	Row Total
U	109	3	2	0	114
A	1	319	4	4	328
F	0	45	96	3	144
W	0	0	0	50	50
Column Total	110	367	102	57	636
PA	99.09%	86.92%	94.12%	87.72%
UA	95.61%	97.26%	66.67%	100.00%
OA	90.25%
KHAT	84.48%

,

Table 6. Accuracy assessment of the classification LULC map in 2011 for 29 July 2011.

LULC Type	U	A	F	W	Row Total
U	137	6	1	1	145
A	13	287	8	2	310
F	0	17	106	8	131
W	0	0	0	50	50
Column Total	150	310	115	61	636
PA	91.33%	92.58%	92.17%	81.97%
UA	94.48%	92.58%	80.92%	100.00%
OA	91.19%
KHAT	86.74%

and

Table 7. Accuracy assessment of classification LULC map in 2017 for 18 May 2017.

LULC type	U	A	F	W	Row Total
U	142	31	0	4	177
A	31	228	7	6	272
F	1	4	130	2	137
W	0	1	0	49	50
Column Total	174	264	137	61	636
PA	81.61%	86.36%	94.89%	80.33%
UA	80.23%	83.82%	94.89%	98.00%
OA	86.32%
KHAT	80.24%

, respectively. These tables reveal that the OA of the four LULC maps varied from 86.32% in 2017 to 97.03% in 2000, and the KHAT of the four LULC maps varied from 80.24% in 2017 to 95.35% in 2000. Here, the average OA and KHAT of the four classification maps are 91.20% and 86.70%, respectively.

Further, the PA varies from 81.61% for urban and built-up land in 2017 to 99.09% for urban and built-up land in 2006. Meanwhile, the UA varies from 66.67% for forest land in 2006 to 100% for water bodies in 2000, 2006, and 2011. In the meantime, the average value of the PA for the four LULC types in 2000, 2006, 2011, and 2017 is 96.46%, 91.96%, 89.81%, and 85.80%, respectively, and the average value of the UA of the four LULC types in the same years is 96.24%, 89.89%, 92.00%, and 89.24%, respectively.

In addition, the OA and KHAT of the LULC maps in 2000 are higher than those of the other LULC maps in 2006, 2011, and 2017. This finding indicates the urban development stage of the study area between 2000 and 2017. In 2000, urban and built-up land was concentrated in the core area of the city, and most of the new construction was located near the core area of the town. Most LULC types were distributed with a large homogeneous patch, and each LULC type was easy to classify. On the other hand, many agricultural lands and forest lands were converted into urban and built-up land, and the developed areas in 2006, 2011, and 2017 were more dispersed from the core area of the city, resulting in a large number of small heterogeneous patches. These phenomena affected the thematic accuracy of LULC maps in 2006, 20011, and 2017 due to the limitation in the spatial resolution of the Landsat data.

3.12. Accuracy Assessment of Change Detection

Figure 20 shows the stable areas of four LULC types and the changed areas between 2000 and 2017 from the four LULC maps for 2000, 2006, 2011, and 2017 by using the newly developed algorithm. As a result, the spatial distribution of the stable LULC type and its change from the newly developed algorithm is similar to the MLC algorithm (see Figure 10). Additionally, the area and percentage of stable LULC and its changes are summarized in

Table 8. Area and percentage of the stable LULC type and its change between 2000 and 2017 using the newly developed algorithm.

LULC Type	Area in km2	Percent (%)
1. Stable urban and built-up land	66.97	4.22
2. Stable agriculture land	389.88	24.57
3. Stable forest land	142.22	8.96
4. Stable water bodies	38.10	2.40
5. LULC change between 2000 and 2017	949.57	59.84
Total	1586.74	100

. Accordingly, the derived areas and percentages from both algorithms are slightly different (see Table 1). Nevertheless, the newly developed algorithm can be easily applied to classify multitemporal LULC maps and detect LULC changes at any specific time.

The result of the change detection accuracy assessment is summarized in

Table 9. The accuracy assessment of change detection between 2000 and 2017.

	Changed areas	Stable areas	Row Total
Changed areas	247	53	300
Stable areas	6	196	202
Column Total	253	249	502
PA	97.63%	78.71%
UA	82.33%	97.03%
OA	88.25%

. As a result, the OA for change detection is 88.25%. Also, the PA for the changed areas is 97.63%, and the PA for the stable areas is 78.71%, while the UA for the changed areas is 82.33%, and the UA for stable areas is 97.03%.

4. Discussion

4.1. Procedure for Landsat Image Selection

As a result of the Landsat selection using cloud cover assessment alongside the clearly observed and contaminated pixel recognition using the QA band and HA model, we can reduce the number of Landsat 5/7/8 scenes between 2000 and 2017 for the time-series data analysis from 673 to 388 scenes (a reduction of about 44.74%). This procedure can not only reduce the number of images and the time for data analysis but can also recognize the clearly observed and contaminated pixels. However, users are required to observe, record, and calculate the percentage of total images and the percentage of total clearly observed pixels according to the cloud cover information of each scene provided by the USGS web portal.

4.2. Semi-Automatic Time-Series Spectral Reflectance Reconstruction

Time-series spectral reflectance data, which is constructed semi-automatically by the multiplication between the spatiotemporal cube of the recognized clearly observed and contaminated pixels (with values of 0 or 1) and the spatiotemporal cubes of the derived six spectral reflectance bands (with values from 0 to 1) using the reshape function of the MATLAB software can remove contaminated pixels from original time-series spectral reflectance data. This process can be quickly implemented, but users should understand basic programming with the MATLAB software.

4.3. Technique for the Extraction of a Stable Area for Each LULC Type

In this study, stable areas for each LULC type between 2000 and 2017 were extracted by overlay analysis based on the classified LULC maps in 2000, 2006, 2011, and 2017 using MLC. Even if some errors in the classified LULC types remain for a particular year, the error of the stable area for each LULC type is minimal since the stable areas of each LULC category are extracted from the multitemporal LULC maps. However, LULC classification with MLC requires a great deal of time to select suitable training areas with a normal distribution. In the future, the extraction of a stable field for each LULC type should be an automatic or semi-automatic process. Currently, many researchers have developed training data automation as an alternative to manual training area selection. Huang et al. [73] applied a dark object concept to automatically generate a training area for mapping forest cover changes with the advanced support vector machines algorithm.

4.4. Supervised Classification Approach for Multitemporal LULC Mapping Using Harmonic Analysis with a Minimum Spectral Distance Algorithm

The developed supervised classification approach for multitemporal LULC mapping using HA with a minimum spectral distance algorithm is a semi-automatic process under the MATLAB software environment and requires little human interference. This approach can efficiently classify multitemporal LULC maps, offering with highly accurate results based on the standard harmonic function curve for each LULC type obtained from the stable areas of each LULC type from multiple years.

Moreover, this approach overcomes the limitations of image selection by considering individual pixels from the available Landsat data. For example, the scan line gaps of Landsat 7 are regarded as contaminated pixels, but available clearly observed pixels could also be used (see Figure 5). Also, this approach does not require the relative normalization of each image for multitemporal LULC classification and change detection because phenological and solar angle differences do not affect classification under the time-series model. Furthermore, this approach applies spatiotemporal filtering with a mode function under post-classification processing to the remove unexpected errors caused by contaminated pixels (e.g., clouds and cloud shadows) or meteorological conditions (e.g., flooding and drought) (see Figure 17). More specifically, this approach can quickly generate a LULC change map for any period and can provide “from–to” change information as a post-classification comparison change detection algorithm, which is frequently utilized for change detection [3,74,75,76,77].

However, this approach has some limitations. Firstly, this approach demands images with high temporal frequencies to ensure accuracy. Secondly, this approach requires a large amount of data storage space (about 500 GB). Finally, the standard harmonic function curve of each LULC type is constructed based on the permanent LULC type over a long period (2000–2017). Thus, LULC types with small areas (such as grassland, wetland, and bareland) might be misclassified into four LULC types because stable areas of LULC types with small patches are difficult to identify for an extended period.

4.5. Accuracy Assessment of LULC Maps

An OA value of more than 85% in LULC maps, as reported in Section 3.11, can provide acceptable results, as suggested by Anderson et al. [78]. Likewise, a KHAT value greater than 80% represents an almost perfect agreement between the classified LULC map and the ground reference data [79]. Further, all the average PA and UA values of different LULC types in 2000, 2006, 2011, and 2017 are higher than 85%. Thus, the developed classification method using HA with a minimum spectral distance algorithm using time-series Landsat datasets is adequate for multitemporal LULC classification.

Furthermore, the overall derived accuracy in this study is comparable with that of other studies that applied time-series Landsat datasets to classify multitemporal LULC maps with a specific algorithm. For example, Gebhardt et al. [80] applied 135 Landsat scenes to classify a series of seven maps of Mexico between 1993 and 2008 by using a decision tree algorithm, which provided an overall accuracy of about 76%. Zhu and Woodcock [45] applied their developed CCDC algorithm with an RF classifier to classify multitemporal land cover maps from time-series Landsat datasets (1982–2011) in coastal New England, United States, achieving an OA of about 90%. Gounaridis et al. [18] applied an RF classifier with time-series Landsat datasets (1991–2016) to classify LULC maps in Attica, Greece and attained an overall accuracy varying from 90.5% to 93.5%. Lu et al. [15] applied Landsat images from 1990 to 2015 at five-year intervals to classify impervious surface and non-impervious surface areas using a linear spectral mixture analysis of the six selected metropoles in the coastal and inland metropoles in 2015, achieving an overall accuracy varying from 94% to 95%. Mi et al. [14] applied an RF classifier with time-series Landsat datasets (1987–2017) to detect the LULC changes in a mining area, achieving an average OA of about 84%. Buitre et al. [13] applied a support vector machine with time-series Landsat datasets (1987–2016) to classify mangroves, non-mangroves, seawater, and clouds in the Philippines and attained an average OA of about 84%.

5. Conclusions

In this study, a new supervised classification method for multitemporal land use and land cover classification was successfully developed using harmonic analysis with a minimum spectral distance algorithm under the MATLAB environment through systematic Landsat selection and time-series spectral reflectance reconstruction by converting the space domain to the time domain. The derived overall accuracies and Kappa hat coefficients of the land use and land cover maps of the four selected years (2000, 2006, 2011, and 2017) were greater than 80%, and the average producer accuracy and user accuracy for different land use and land cover types from the four selected years were also more than 85%.

Therefore, it can be concluded that the newly developed supervised classification method using harmonic analysis with a minimum spectral distance algorithm can be efficiently used to classify and map multitemporal land use and land cover and to detect its changes from time-series Landsat datasets with high reliability for the information. However, this developed classification method should be examined in other areas to determine its spatial and temporal transferability with a more detailed set of classes. Nevertheless, the presented workflow of the research methodology can be used as a guideline for software developers for (semi) automatic land use and land cover classification and mapping.