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

: 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.


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].
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.

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).

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).

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): where i is the cloud cover (%) that is less than or equal to the interval class (10, 20, …, 100), , is the percentage of total images, , is the number of images, and is the number of total images. j is the cloud cover (%), , is the percentage of total clearly observed pixels, and 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].

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]. * where is the spectral reflectance, M is the reflectance multiplicative scaling factor, A is the reflectance additive scaling factor, and Q 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.

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].
where t is the Julian date, i is the i th Landsat band, T is the number of days per year (T= 365), ai is the coefficient of the intercept value, bi is the coefficient of the slope value, Ai is the coefficient of the amplitude value, 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).
, , , where , , , and are the spectral reflectance from the standard harmonic function curves of U, A, F, and W, respectively, at time point i using band j.
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. 1 where , , and 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.

6
(17) where, , , , and 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 largescale, 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].

Results
In this research, four Landsat images for 2000 ( were selected to classify multitemporal LULC data using HA with a minimum spectral distance algorithm under the MATLAB software environment. 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.

Selection of Landsat Data According to Cloud Cover Assessment
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.

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  7(a), most of the contaminated pixels are cloud and cloud shadow, whereas the contaminated pixels in Figure 7(b) include cloud, cloud shadow, and gaps.

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 8(b), are considered here as contaminated pixels and are ignored in the data analysis.

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 9(a) 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 9(b) 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.  Table 1 and displayed in Figure 10 On the other hand, the unstable area of LULC types was 1,009.87 km 2 (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.

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.    Table 2 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.   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.

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.

Multitemporal LULC Classification and Mapping
Figure 17(a) 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 17(b) shows the final results after post-classification processing for the unclassified pixel; this figure reveals that the unexpected errors from water bodies in Figure 17(a) 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.  Figure 18 and Table 3.        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. 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.

Accuracy Assessment of Change Detection
The result of the change detection accuracy assessment is summarized in Table 9. 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%.

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.

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.

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.

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].

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  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  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%.

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.