An Automatic and Operational Method for Land Cover Change Detection Using Spatiotemporal Analysis of MODIS Data: A Northern Ontario (Canada) Case Study

Abstract

Mapping and understanding the differences in land cover and land use over time is an essential component of decision-making in sectors such as resource management, urban planning, and forest fire management, as well as in tracking of the impacts of climate change. Existing methods sometimes pose a barrier to the effective monitoring of changes in land cover and land use, since a threshold parameter is often needed and determined based on trial and error. This study aimed to develop an automatic and operational method for change detection on a large scale from Moderate Resolution Imaging Spectroradiometer (MODIS) data. Super pixels were the basic unit of analysis instead of traditional individual pixels. T² tests based on the feature vectors of temporal Normalized Difference Vegetation Index (NDVI) and land surface temperature were used for change detection. The developed method was applied to data over a predominantly vegetated area in northern Ontario, Canada spanning 120,000 sq. km from 2001–2016. The accuracies ranged between 78% and 88% for the NDVI-based test, from 74% to 86% for the LST-based test, and from 70% to 86% for the joint method compared with manual interpretation. Our proposed method for detecting land cover change provides a functional and viable alternative to existing methods of land cover change detection as it is reliable, repeatable, and free from uncertainty in establishing a threshold for change.

1. Introduction

The complexity of the Earth’s dynamics requires continuous monitoring of its data to be able to predict future events and respond accordingly. It is vital to identify changes in land use and land cover (LULC) and atmospheric conditions as they impact the global climate. Monitoring forest and agricultural land is sometimes carried out to assist in land use planning or for tracking the exploitation of natural resources. Our objective for LULC change analysis is to examine the impacts of climate change from land management schemes. The area under study is in Northern Ontario, Canada. The region is expected to see an increase in land use change due to the opportunities the warmer weather brings. With the different activities that will occur, greenhouse gas emissions (GHG) and stored organic carbon will be displaced. This study is part of a larger work which will assess the amounts of soil organic carbon and GHG’s from the land cover changes which occur in the region.

In addition, the information on the land cover and land use change can be factored into critical decisions made by governments. There are two common approaches to change detection. In the first category, the objective is to detect whether a change has occurred or not. The outcomes are always dualistic: “Change” or “no change.” In the second category, the nature of the change is also determined, and thus the solution provides “from-to” information. For studies focusing on large-scale monitoring, “change/no-change” detection is often carried out. In addition, results from “change/no-change” detection could be used as the first step to identify areas for further detailed “from-to” change detection. The “change/no-change” detection was the focus of this study. A detailed change detection based on high spatial resolution imagery will be carried out for the areas and for the periods where/when changes occurred.

Advances in remote sensing have permitted the development of automated or semiautomated change detection methods [1,2]. The thresholding method of change analysis establishes a threshold beyond which a change is assumed to have occurred. For instance, image differencing uses a threshold for grey levels of the image to note pixels which have undergone land cover changes [3,4]. The image differencing can be performed through algebraic image differencing, vegetative index differencing, image ratioing, regression, or change vector analysis [5,6,7,8].

Specifically, Kleynhans et al. [4] applied a threshold method to detect where change had occurred, and obtained a change detection accuracy of 89%. Lu et al. [9] applied the Breaks for Additive and Seasonal Trend (BFAST) model to the spatiotemporal region to observe change for a period of 12 years. The feature used was the temporal enhanced vegetation index. The highest producer’s accuracy obtained on the changed pixels was just over 50% at a 5% significance level. Another study used a threshold method and statistical analysis [10]. For change detection, the authors considered the pixel neighborhood, where the size of the neighborhood again was chosen by a trial-and-error method.

This current work aims to improve on a change detection method which examines pixel-by-pixel differences. In those methodologies, the analyst must choose the threshold to identify where change has occurred or not. For that method, the analyst assigns the minimum or maximum threshold value to decide if change in land cover has occurred or not. Even though these pixel-based methods have successful accuracies, a drawback to the method is having to choose an appropriate threshold, finding a threshold is typically a trial-and-error process. The shortcoming of this process is not knowing which threshold is appropriate for the analysis. Furthermore, some pixel-based image differencing methods produce false detection due to registration errors [10]. As a result, further research is required for addressing the abovementioned issues.

A change detection technique is chosen according to the application. For the application of large-scale monitoring of binary change, our proposed method is appropriate and operational. In literature, there are various features which have been exploited for land cover/land use change detection, including spectral or spatial features of remotely sensed images, the texture, or vegetation indices. This work used land temperature and a vegetation index as features to detect change.

The Normalized Difference Vegetation Index (NDVI) is a well-used vegetation index because of its success in detecting vegetation, the ease with which it is calculated and interpreted, and its accessibility from various satellite data. The NDVI is known to be the primary vegetation index for global operational applications [11]. The NDVI was used in our work to monitor change in the study areas. The monitoring was done on an intra-annual basis and an interannual basis. The other feature considered in this work is the Land Surface Temperature (LST). Studies have shown the impact that land cover and land use changes have on LST [12,13,14,15]. Depending on the type of land cover present, the surface temperature in an area could vary [16]. Similarly, when the land undergoes a change in land cover or land use type, it could experience a simultaneous change in LST. Majumder et al. [17] noted that LST could decrease in an area which increases in agriculture. Zhang and Liang [18] used China’s land use/cover datasets (CLUDs) to study the impact of land cover change on surface temperature. The database contains images of 1km raster size having 25 land cover classes. These classes are further grouped into cropland, grassland, woodland, unused land, water, and urban areas. The study conducted by Zhang and Liang [18] showed that land cover change from cropland to built-up areas increased the LST, while a change from forest to cropland reduced the surface temperature.

Although a multitude of studies have examined the relationship between NDVI, LST, and land cover/land use change, few have aimed to use LST change to predict or express a difference in land cover. This paper, therefore, proposes a model to predict an area of land cover change from both LST and NDVI values.

Other statistical change detection methods have been proposed in the literature. In the study by Waylen et al. [19], interannual changes were detected over a 25-year period from an AVHRR NDVI dataset. The development of the timeseries persistence analysis was based on a null hypothesis which stated that the values of NDVI are normally distributed and serially independent [19]. The metrics used were measures of directional persistence, relative directional persistence, and massive persistence of NDVI change. The first two metrics only consider the direction of change, while the latter also accounts for the magnitude of change. The use of the massive persistence test requires a separate classification to obtain the statistical control pixels which did not experience change in land cover or land use. Each persistence method compares a monthly NDVI to a user-selected benchmark NDVI value (directional persistence) or to the monthly NDVI observation of the previous year (relative directional persistence and massive persistence method) to determine if there has been a positive or negative change. Monte Carlo simulations were used to compute the variables’ distributions. Teng et al. [3] used two hypothesis-test-based land cover change detection methods: The bivariate joint distribution method and the conditional distribution method. The former method takes a group of no-change pixels (identified from a post-classification method) to estimate the mean and covariance matrix of each land use class. This distribution could vary with land use classes. The conditional distribution similarly uses parameters which are class-dependent and need to be estimated from the no-change pairs of pixels of individual land use classes. Testing the methods using SPOT-4 images, the methods had accuracies ranging from 92% to 94%. For bivariate joint distribution, the highest overall accuracy was achieved at the 5% level of significance, while the highest overall accuracy was achieved at the 10% level of significance for the conditional distribution method. Although these statistical methods did not require a threshold for change detection, they are reliant on a prior land use classification.

In this study, we propose an automatic change detection process. Our method uses a hypothesis test as the basis for determining where land cover change has occurred. To the best of our knowledge, many methods of change analysis use threshold methods, but a hypothesis-based test is scarcely used for such analyses. The developed method was validated based on data collected over a large area in Northern Ontario, Canada.

2. Materials and Methods

2.1. Study Area and Data

The study was conducted over Ontario’s Great Claybelt region located in northeast Ontario, Canada. The region spans 120,000 sq. km, centered at ($49.476 °\mathrm{N},82.283 °\mathrm{W}$). Two areas within the Claybelt region were selected for validation: The town of Hearst at the westernmost part of the region, and the town of Cochrane in the southeast portion of the Claybelt, as illustrated in Figure 1. Hearst covers an area of 15,000 sq. km, while the town of Cochrane covers 539 sq. km [20]. Although large spatially, these areas are sparsely populated. For instance, Hearst has a population density of 51.1 persons per sq. km, and Cochrane’s population density is 9.9 persons per sq. km [20].

This study examined change within 3 periods: From 2001 to 2005, from 2005 to 2009, and from 2010 to 2016. Initially, to study the land use in the study area, the maps from Agriculture and Agri-Foods Canada (AAFC) were reviewed. Figure 2, obtained from AAFC datasets, demonstrates one of the challenges in characterizing land usage in Northern Ontario: Datasets are decadal, and they often do not extend over the entire region. The lack of reliable, accessible, and continuous geospatial information in Northern Ontario has resulted in expensive and protracted solutions, or broad portrayals of the region. This example—with no data coverage at the topmost portion of the Claybelt—highlights the urgent need to develop suitable methodologies for detecting and mapping land cover dynamics in Northern Ontario for forest and agricultural lands. These, in turn, can be used to study the changes in soil carbon stock and greenhouse gas emissions.

Based on existing land cover maps such as the Ontario Far North Land Cover Compilations [21,22], the Cochrane map indicates that, in 2016, the town was comprised of mostly mixed forest in the central and northern parts, while evergreen needleleaf coniferous forests were in the northwestern part of the town. The land cover types categorized as ‘savannah trees’ or ‘woody savannahs’ were found mostly in the west-central to south part of the town. (The savannahs correspond to fen, bog, and trees land covers in the Ontario Far North Land Cover map. They will henceforth be referred to as fen/bog in this writeup.) That year, Cochrane comprised mostly of fen/bog, grassland, and deciduous broadleaf trees in the central part of the town (see Figure 3). The mixed forest covered most of the northern part of the town, and fen/bog and mixed forest had the highest population in the south. The land cover map of Cochrane also distinctly shows an area of cropland in the central part of the town. Given the moderate resolution of the map, this is significant, indicating that there is a vast amount of farmland in that portion of Cochrane.

That swath of cropland in central Cochrane was a change from the town’s landscape a decade and a half prior. In 2001, Cochrane was mostly fen/bog from the central part of the town to the south, with some grassland and deciduous trees. The contrast is evident from the land cover maps, which show that a large portion of the fen/bog areas from 2001 had been converted to grassland and cropland in the central part of the town by 2016. Furthermore, the deciduous trees’ population has increased whereas the fen/bogs have decreased.

There have also been changes in land cover in the Hearst forest area over the 16 years examined. For instance, in 2001, the southernmost part of the town was covered by fen/bogs, grassland, and mixed forest. In 2005, that portion of the town was almost entirely fen/bog and mixed forest only. By 2010, the mixed forest had taken up most of the area where the fen/bog had been, and some of the fen/bogs were replaced by deciduous trees. In 2016, very little of the fen/bog remained, and the area was predominantly mixed forest with a little of deciduous trees. A similar significant change is also noticed in the northeast part of Hearst. In 2001 the northeast part of Hearst began as a mix of coniferous trees, mixed forest, and some fen/bogs. By 2005, some of the coniferous trees had been replaced by fen/bogs and grassland. In 2010, there was even more dominance of the fen/bogs in the areas where there was previously mixed forest and coniferous trees, and by 2016, most of the mixed forest was replaced by larger portions of grassland and fen/bogs. The number of coniferous trees did not appear to change much in that part of the town between 2010 and 2016. These land cover changes were observed from the land cover maps produced for these years.

The remotely sensed data employed in this study were obtained by MODIS. The NDVI and LST images were taken from Google Earth Engine. This web platform hosts publicly available satellite imagery from around the globe which can be manipulated for a variety of purposes. The MODIS Terra data products used for the change analysis were the 16-day 250 m resolution product for NDVI (MCD13Q1) [23], the 500 m resolution surface reflectance (MOD09A1.006), and the 8-day 1 km resolution product for LST (MOD11A2) [24]. Linear interpolation was performed for missing data, interpolating both temporally and spatially. An Image Collection with a total of 23 images was used for the months’ NDVI maps, and for the LST, an Image Collection of 4 images was used for each month. The data in Google Earth Engine have typically undergone preprocessing. The analysis was performed with the Google Earth Engine code editor [25], ArcMap 10.7 [26], and MATLAB [27]. It is worth mentioning that though NDVI was used because of its historical success as a vegetation monitoring measure, the developed method should work for other products, such as the Enhanced Vegetation Index (EVI), Leaf Area Index (LAI), and Fraction of Photosynthetically Active Radiation (FPAR), which will be explored in future work.

Image preprocessing was carried out to remove the unwanted artefacts such as those created by atmospheric interferences. The MODIS Surface Reflectance product at 500 m resolution was estimated to be the surface reflectance at ground level. It was corrected for atmospheric conditions such as gasses, aerosols, and Rayleigh scattering. Similarly, the MODIS NDVI data from Google Earth Engine consists of atmospherically corrected bi-directional surface reflectance which have been masked for water, clouds, heavy aerosols, and cloud shadows. The LST Level 3 products used in the simulations had cloud-contaminated LSTs removed already. Another element of image preprocessing was the smoothing of input data to remove noise. In this work, we used a Savitzky-Golay filter as the smoothing technique. Other smoothing methods considered were the Best Index Slope Extraction, logistic and double logistic smoothing, asymmetric Gaussian, and mean- and median-moving average. The premiere method did not change or smooth any of the NDVI data that it was applied to and thus was not used further. Of the latter 2 methods, the moving median gave a better representation of our data compared to the moving mean. The rest of the filtering methods tested also produced a good degree of smoothing to the input data. Of all methods, the Savitzky-Golay achieved the best result, using the greatest number of datapoints while still maintaining the general trend of the plots. A window size of 5 was chosen as our proposed approach obtains its optimal performance with that window size. Figure 3 is an example of the effect of smoothing on the NDVI data in 2010. The rest of this section describes the approach taken to create our land cover change detection method and the procedure used for assessing the accuracy of the method.

2.2. Methods

As mentioned earlier, the change detection in this study was performed over 16 years, beginning in 2001, using NDVI and LST data to assess the change in land cover roughly every 5 years. Super pixels were used as the basis of analysis for detecting change instead of the traditional pixels-based method. They were relatively small in size to ensure the areas covered by them were homogeneous. Alternatively, an arbitrary neighborhood around each pixel could be used as samples, for which each observation (pixel value) may not follow the same random distribution. Moreover, the arbitrary neighborhood may also be computationally expensive. To determine the size of the super pixels, the spatial variations of the areas of interest exhibited in the MODIS data were analyzed using a semivariogram. The input to the semivariogram was the first principal component of the MODIS reflectance. The super pixels were obtained based on the first 3 principal components of MODIS imagery covering the optical bands. The statistical test to determine if significant change had occurred was performed on the individual super pixels of the NDVI and LST maps. Below is the step-by-step process followed in our methodology.

Principal Component Analysis (PCA) [28] was performed on a 2001 image of the surface reflectance with the first 7 MODIS bands. The PCA showed that more than 98% of the data could be explained by the first 3 components.

To characterize the spatial variation of the study area, a semivariogram was first calculated based on the first principal component. This allowed us to determine the dominant object at each site so that a homogenous area could be found. The semivariogram provided information on the appropriate size for the super pixels. Figure 4 illustrates a semivariogram and how the sill and range are determined. For example, in Figure 4, the plateau of semivariance values is reached at a range of about 104 pixels. Whereas the sill explains the variance structure of the data, the range shows the farthest lag distance at which observations are correlated. Hence, there is no spatial correlation for the observations separated by more than 104 pixel units.

An image can be segmented such that each segment represents a collection of about 104 pixels having a similar characteristic such as pixel value or pixel intensity. The pixels in this segment are termed super pixels. The super pixels can be compared to determine if a change has occurred in the environment.

The super pixel package used in this study to perform the image segmentation was the Simple Linear Iterative Clustering (SLIC) package. A characteristic of the SLIC is that it makes smooth, regular-sized super pixels in smooth regions and highly irregular super pixels in textured regions [29]. Image segmentation was carried out using the first 3 principal components.

The result from the image segmentation performed on the Hearst site is shown in Figure 5 as an example. The unit for the super pixel was obtained from a 2005 LST image.

Following the principle of ergodicity [30], pixel values within each super pixel were considered as samples in the hypothesis test. Ergodicity is important because we do not have a large number of samples from which to compute sample averages. For a random process which is ergodic, only 1 sample function is required. In our remote sensing application, we applied the spatial property of the observation to exhibit randomness. Hence, using super pixels which are homogenous over a spatial area to represent our data, we drew on the ergodic feature, which allowed us to use the population for the statistical test.

Hotelling’s T² test [31] was selected in this method as the T² test can compare multivariate data. Thus, where univariate data can be compared using a t-test, multivariate data can be compared using the T² test. The T² test, like the t-test, compares samples using means. In Hotelling’s case, the multivariate mean is used because it compares the means of 2 or more samples. Thus, the T² statistic can be considered as comparing a point in space defined by the means of all the individual variables.

When X₁, …., X_N are k-component vectors, and the normal distribution has a vector mean µ and a covariance matrix Σ with the sample mean defined as

$$\overline{x}=\left(\frac{1}{N}\right){\displaystyle \sum }_{i=1}^N{X}_i$$

(1)

and the sample covariance matrix as

$$S=\frac{1}{n-1}{\displaystyle \sum }_{i=1}^n\left(X_i-\overline{x}\right){\left(X_i-\overline{x}\right)}^{\prime }$$

(2)

then the generalized T² statistic is given as

$$T^2={\left(\overline{x}-{\mu }_0\right)}^{\prime }{\left(\frac{S}{n}\right)}^{-1}\left(\overline{x}-{\mu }_0\right)$$

(3)

where ${\mu }_0$ is a known vector, and there are n − 1 degrees of freedom.

The paired Hotelling’s T² test statistic is given by the expression:

$$T^2=n{\overline{x}}^{\prime }{S}^{-1}\overline{x}$$

(4)

The T² test is more advantageous than a t-test because the Type I error rate is better controlled in the former than in the latter. Another benefit of the T² test is that it considers the relationship between the different variables, thereby the differences between groups can be obtained. The assumptions for this statistic are that the samples have normal distributions, they are independent, and if there are only 2 samples, then they have equal variance-covariance matrices.

The test was performed to confirm if change had occurred in an area. Thus, the null hypothesis statement is: No land cover change has occurred in these areas in a 5-year period. The hypothesis was tested using the difference in mean values of the parameters under study for the 5 years. The null hypothesis would be rejected for the alternative hypothesis if the sample data used for the study showed the null hypothesis to be false. That is, if change had occurred, the null hypothesis would be rejected. The null hypothesis implies that the mean LST and the mean NDVI remained the same. That is, their population mean vectors are equal.

For both the NDVI and LST maps, the monthly mean of the parameter was obtained for each super pixel. Then, a vector of these values was created in each year under study. For instance, when comparing 2005 records to 2009, a 12 × 1 vector of the differences of means was created. Also, for each segment, the covariance of both years’ data was found, creating a 12 × 12 matrix.

The range from the semivariogram was used to determine the maximum number of pixels which should form a super pixel. Using Equation (4), the T² statistic was found for each super pixel. Thereafter, the Hotelling test was performed at a significance level of 5%. The algorithm returns the value of 0 if the test indicates that an insignificant change has occurred (that is, the null hypothesis should be accepted). Similarly, it returns a value of 1 if the mean vector difference indicates that a significant difference has occurred—that is, the null hypothesis is rejected.

An extension of the proposed change detection method was explored whereby the NDVI and LST data of each year were combined into 1 matrix, and their land cover change was predicted using the algorithm described here. This extension entailed normalizing the NDVI and LST data before making the comparison between the years (hereafter referred to as the joint method). More details on this approach will be presented in the Section 4 and Section 5.

The proposed method is well suited to forest areas. Croplands may produce results of more occurrences of land cover change depending on the crop senescence. Thus, the next part of this study will be important in agricultural areas, identifying the “from-to” changes. Our method does not emphasize the magnitude nor type of change. Instead, it identifies areas of change so we may concentrate on those areas to later identify the type of land cover change to calculate carbon amounts.

To quantitatively estimate the performance of our change detection method, an assessment of the accuracy was performed using higher resolution imagery. In total, 50 random samples (super pixels) were taken over the study area for the years 2001, 2005, 2009, 2010, and 2016. For these years, the sample locations were overlaid on maps of the area at Landsat-level (30 m) resolution. Thereafter, the 50 locations of the Landsat-level maps were compared via visual interpretation. The “change/no-change” verification was performed by visual interpretation of 2 independent researchers. In a few instances, there was disagreement in the manual interpretation of whether change had occurred or not in the sample sites. The discrepancy in those cases was resolved by accepting the option which represented the predominant status, e.g., if the sample location was mostly unchanged, then it was classified as unchanged. Finally, the verification results were compared to the output from the statistical test.

3. Results

Following the statistical tests, the program output an image of the study area, with 0 representing a no change area and 1 representing a changed super pixel. The areas where there was a zero value were made transparent, and changed areas were given a grey color. The results in the form of binary maps are shown in Figure 6, Figure 7 and Figure 8 below.

Figure 6 illustrates the differences in areas where the LST test detected change and where the NDVI test detected change between 2001 and 2005. Overall, the LST detected more areas of significant change than the NDVI test. The most noticeable differences were observed in the north central part of Hearst and in the mid-eastern part as the NDVI test returned less significant change in land cover than the LST in those areas. The parts which were more similar were in the western and central part of the town.

As shown in Figure 7, the land cover change detected between the LST and NDVI tests was more similar in the west and central portions of the town. The change detection results were dissimilar mainly in the north central part of the town, where the NDVI test predicted mostly no significant change in that area, whereas the LST predicted that significant change had occurred in that area between 2005 and 2009.

There was good agreement in the change detection results between the NDVI and LST for the period 2010 to 2016 (Figure 8). Apart from the north central part of the town, and an area in the central west part of the town, the change detection maps produced from NDVI and the LST indicate similar land cover change occurring. In the areas of dissimilarity noted, the NDVI indicated no significant change while the LST had predicted a significant change in land cover had occurred.

An accuracy assessment was conducted on each of the NDVI and LST change maps produced, comparing areas of “change/no-change” with the high-resolution maps—maps of the area at 30 m resolution, which is the Landsat-level resolution—for the years under study.

Table 1. Confusion matrix representing change detection accuracy of NDVI-based test: 2005–2009.

	Reference Data
Classified Data	Changed	Unchanged	Totals	User’s Accuracy (%)
Changed	21	4	25	84
Unchanged	8	17	25	68
Total	29	21
Producer’s Accuracy (%)	72	81	Overall Accuracy: 76% Kappa Coefficient: 0.52

,

Table 2. Confusion matrix representing change detection accuracy of LST-based test: 2005–2009.

	Reference Data
Classified Data	Changed	Unchanged	Totals	User’s Accuracy (%)
Changed	24	3	27	89
Unchanged	5	18	23	78
Total	29	21
Producer’s Accuracy (%)	83	86	Overall Accuracy: 84% Kappa Coefficient: 0.68

and

Table 3. Confusion matrix representing change detection accuracy of NDVI-based test: 2010–2016.

	Reference Data
Classified Data	Changed	Unchanged	Totals	User’s Accuracy (%)
Changed	21	5	26	85
Unchanged	1	23	24	96
Total	22	28
Producer’s Accuracy (%)	96	85	Overall Accuracy: 88% Kappa Coefficient: 0.76

provide examples of the results achieved. The accuracies obtained varied from 74% to 88%. For instance, the overall accuracy of the change predicted by NDVI between 2001 and 2005 was 78%, while the LST had a 74% accuracy. Similarly, the accuracy of the LST change detected was 84% for the period 2005–2009 (see Table 2). The average overall accuracy of the change detection process was 81%.

The overall accuracy of the NDVI change map was 76% with a kappa coefficient of 0.52 in the period 2005 to 2009. The user’s accuracy for changed areas was 84%, while the user’s accuracy for the unchanged areas was 68%. Similarly, the producer’s accuracy was 72% for the changed areas because eight areas which had undergone land cover change were incorrectly classified as being unchanged. The producer’s accuracy was better for the unchanged areas as only four of the unchanged areas were predicted by the statistical test as being changed (see Table 1).

Table 2 below provides the accuracy results from the same timeframe when LST was the feature on which the change assessment was made.

The LST change map for the change between 2005 and 2009 yielded an overall accuracy of 84%. Similar to the NDVI test for that period, more errors of omission were found in the changed areas than errors of commission. However, for the unchanged areas, most of those chosen for accuracy assessment were correctly classified as unchanged.

Table 3 displays another example of the results from the change detection test. When NDVI was used to predict land cover change, the overall accuracy was 88% with a kappa coefficient of 0.76.

In the detection period 2010–2016, there was less confusion in the NDVI change maps pertaining to the changed areas. Only one area was incorrectly classified as unchanged (see Table 3). This resulted in a high producer’s accuracy for the changed areas.

When the joint method of NDVI and LST together were used to predict land cover change, the accuracies were in the same range as the standard (our proposed) method using individual parameters. The average overall accuracy using this approach was 80%.

Table 4. Accuracy results from NDVI- and LST-based (joint) test: 2010–2016.

	Reference Data
Classified Data	Changed	Unchanged	Totals	User’s Accuracy (%)
Changed	20	2	22	91
Unchanged	5	23	28	82
Total	25	25
Producer’s Accuracy (%)	80	92	Overall Accuracy: 86% Kappa Coefficient: 0.72

below illustrates an example of the land cover change prediction results from the joint method.

The resulting map from the joint method had a lower producer’s accuracy than the user’s accuracy of the changed areas, signifying a higher omission error than commission error (Table 4). The lower producer’s accuracy was evident in the changed category because five changed super pixels were identified as unchanged, while only two super pixels which were unchanged were incorrectly classified as changed. The overall accuracy from this joint method was 86% for the period 2010–1016.

To examine the possible change in the study areas exclusive of land cover type, the variations in NDVI and LST were observed from their monthly maps. Figure 3, Figure 4, Figure 5 and Figure 6 illustrate an example of how the maps of the two parameters were used to obtain a cursory perception of the change in land cover in the study areas.

The NDVI maps created in Figure 9 show the difference in NDVI for the Hearst area between 2001 and 2010. In January 2001, there was more healthy vegetation in the northern parts than there was in January 2010. However, in March 2010, the vegetation cover was much higher (on average 0.4) than in the same month in 2001 where the average NDVI value was 0.17. Similar NDVI values were observed for May between 2001 and 2010. The NDVI in November was similar in the 2 years, and the southern part of Hearst had a slightly higher amount of healthy vegetative structure in 2001 than in 2010. Their average NDVI values were 0.4 in 2001 and 0.45 in 2010.

Similarly, Figure 10 is a representation of the NDVI changes which occurred in Cochrane between 2001 and 2010. Figure 11 and Figure 12 display the differences in temperature maps for the towns.

Figure 11 shows that the land surface temperature was quite different in the 2 years under comparison: January and March 2010 had a higher average temperature than in 2001, particularly in the east, south, and northwest. For example, the average temperature in Hearst was 6 °C in March 2010, whereas it was −3 °C the same month in 2001. The maps show that the temperatures were similar in May, with 2010 showing about 2–3 °C higher temperature, particularly in the east and north of the town. November 2001 had a higher average temperature than the same month in 2010. The average temperature in 2001 was 3.5 °C, while that in 2010 was −2.4 °C.

Figure 12 presents the differences in average surface temperature for Cochrane. As shown in Figure 12, 2001 and 2010 had similar average surface temperatures in January and November. Everywhere in Cochrane, March 2010 had higher average temperature than in 2001. The differences in temperature were minor, with an average of 2.8 °C and −2.9 °C, respectively. May 2010 also showed a 2 °C higher average temperature than May 2001, particularly in the central portion of the town.

Figure 10, Figure 11 and Figure 12 represent a cross-section of the spatial distribution of NDVI and LST in Hearst and Cochrane over the examined years.

4. Discussion

There is a multiplicity of ways carbon is introduced into an ecosystem. Our project, of which this paper is one aspect, seeks to quantify the carbon and GHGs linked to land use conversion in farmlands and forests. Studies have shown that land cover change and land use change can promote global climate change [32]. When there are disturbances to these environments, there is a corresponding effect on the carbon levels. For example, there could be a decrease in ecosystem carbon when trees are logged and removed from the forest floor. Similarly, when the farmland undergoes change to become a woodlot or is abandoned, there could be a net carbon increase. Between 2007 and 2016, the total emissions from forestry and other land use change (e.g., afforestation, deforestation, wood harvest, peatland burning) amounted to 5.2 (±2.6) Gt CO₂ per y, whereas the total worldwide emissions from other sources such as emissions from industrial sources, aviation, and waste amounted to 33.9 (±1.8) Gt CO₂ per y in the same period [33]. This paper was the first phase of the project: Identifying areas which had undergone some land cover change. Thereafter, these changed areas will be examined more closely to determine the type of change to quantify the net carbon emission from each land use change.

The NDVI has long been used as an indication of the presence or absence of healthy vegetation. Through remote sensing observations, the vegetation canopies and characteristics of the landscape are estimated. The LST is also an important indicator of the changes in an environment or landscape. For example, if an area was barren land and then converted to being built up, or if it was bare and then changed to being highly vegetated with copious tree coverage, the average LST of the area would be different between the two states. The LST has been shown to be affected by land cover or land use changes, whether induced by natural causes or human activity [34], with bare ground being most affected by changes in LST [34]. These two parameters—the NDVI and LST—can be used in monitoring tools of global climate change effects. Textural features are commonly used and are more effective when higher spatial resolution imagery is used. In this study, the main objective was to develop an operational method. With a low spatial resolution imagery involving broad cover types, the seasonal changes in foliage and temperature are more effective measures which could be reflected by temporal NDVI and surface temperature. In future research, other indices will be investigated.

As mentioned in the Introduction, this study presents an operational method to determine where land cover change has occurred. The proposed method for determining the land cover change that took place in a given period uses an algorithm we developed. It applies the Hotelling T² statistical test to determine when a significant difference between two sets of data has occurred. An analyst only needs to define the required significance level of the T² test. This is easier for the user and provides more consistent results than existing methods of change detection which necessitate choosing a threshold for determining change, where that choice is made by trial-and-error. Image processing techniques of Principal Component Analysis and image segmentation were performed on the input data. In our study, we subdivided the input images into super pixels. The same super pixels obtained based on the first three principal components of MODIS imagery covering the optical bands were used for both NDVI and LST. It is worth mentioning that the same super pixels were used for both NDVI and LST. There is the possibility that using the same super pixels might not be “ideal” for LST analysis. It is assumed that the super pixels obtained from the high spatial resolution optical imagery would capture the variations exhibits in the low spatial resolution thermal imagery. The image processing also entailed Savitzky-Golay data smoothing. Although other methods of smoothing were explored, similar results were obtained.

It should be noted that our proposed method is not a classification effort. Therefore, performance comparisons were made to land cover change detection methods which only considered a change or no change, as our method does. An assessment of the accuracy of the method was performed to estimate the agreement with the reference data, the Landsat-7 imagery of the study areas. The tests resulted in overall accuracies ranging from 74% to 88% for the three time periods under study. This result is consistent with similar studies of land cover change detection using MODIS data [1,3,10,35]. Although the LST has a spatial resolution of 1km while the NDVI data has a resolution of 250 m from the MODIS satellite, the accuracy results within the study period are comparable for both parameters.

The joint method compared well to the standard change detection method proposed. Apart from the first change period, the overall accuracy of the joint method was roughly equal to or better than the average accuracy obtained for the standard method which uses individual features. For the 2001–2005 change period, the accuracy of the joint method was 70% compared to the standard’s average of 78%. The joint method had an 84% accuracy in change detection between 2005 and 2009 compared to an average of 80% for the standard methodology, and the joint method had an 86% accuracy for the change detected between 2010 and 2016. This is about equal to the average accuracy of 87% using the standard method. In general, the predicted change results are consistent in each of the detection periods. Upon examining the results further, it was observed that most of the discrepancy in change predicted by the two features occurred in areas where the cover type at the end of the detection period was mixed forest. For instance, several errors of omission occurred where the cover type had changed from deciduous treed to mixed forest. In the study, the joint method provided similar results to the proposed standard method. Before running the statistical test, we checked if there was a significant correlation between the LST and the NDVI—if the LST can be a predictor to characterize the vegetative state of the land under study. We observed that the NDVI and LST parameters were well correlated. Testing the correlation at different areas of the Great Clay Belt during the 16-year study period, we obtained correlation coefficients ranging from 0.85 to 0.94. In situations of high correlation between parameters, the joint method could be a suitable substitute, though our standard change detection method is being proposed here.

As part of the post-calculation results assessment, we compared our results to another method in the literature. Thus, change vector analysis (CVA) was performed. A change vector can be calculated between two vectors in n-dimensional space over the same geographical location within a time period [36]. Figure 13, Figure 14 and Figure 15 provide examples of the change amount and direction within a detection period.

Figure 13 illustrates the change calculated in Cochrane between April 2005 and April 2006. The CVA was done in ArcMap. The figure illustrates a moderate change throughout most of the area, except in the southernmost part of the town where there was less change. Similarly, comparing the change between November 2005 and 2006, there was change mostly between 0.04 and 4.6. The direction of change is also illustrated in Figure 13, which shows there was mostly a positive shift in the feature between the 2 years.

Similarly, the CVA on NDVI maps from Hearst from 2001 to 2005 and from 2009 to 2016 is illustrated in Figure 14 and Figure 15. These figures show that the landscape changed or shifted in various parts of the Hearst forest area over the periods, with the change ranging between 10 to 60 in the first period and averaging about 10 in the second period. The direction of change was mostly negative in the eastern part of the Hearst forest area between 2009 and 2016, while it was mostly positive in the Hearst region between 2001 and 2005. Overall, the CVA maps indicate most of the Cochrane and Hearst region experienced change during the detection period, though of varying magnitude.

We observed that our proposed statistical method found several super pixels which had experienced change in the years under study. This might have to do with the geography of the study areas: Both towns examined are heavily vegetated. Furthermore, there are forest areas in the Great Clay Belt which are continuously logged. Figure 16 illustrates some of the disturbances which regularly occur in the Claybelt area, with the figure displaying some events which occurred between 2001 and 2009. With the various activities in the region which could impact the landscape, the changes we observed in our analysis is to be expected. They are also borne by the visual differences observed in Figure 9, Figure 10, Figure 11 and Figure 12. As the figures show, there are differences between the surface temperature and vegetation index between 2001 and 2010. Thus, Figure 16 gives a possible explanation as to why several areas of change were detected. As a quick check to confirm the changes in vegetated areas in Hearst, the Hansen et al. [37] method of detecting forest loss and forest gain was used to estimate the yearly forest area changes over the study period. It was seen that there was approximately 18,000,000 m² loss between 2006 and 2007, and between 2008 and 2009, there was 8,000,000 m² loss. In addition, 26,000,000 m² forest loss occurred between 2010 and 2011, and approximately 12,000,000 m² forest gain was estimated in the town between 2014 and 2015. These figures again confirm the expectation of observing significant change from the statistical test method.

Assessing “change/no change” in agricultural lands may require a different approach. These lands may experience changes more than once in the year following the pattern of the growing season. Thus, one may want to demarcate those areas to be investigated separately. For instance, the vegetation index would measure the variation in greenness to not only determine if and where agricultural land changed during the year. It can also determine if, for instance, an area of land was converted from a forested habitat to an agricultural area. These types of land conversions will be of great interest when studying land use change for the purpose of determining the carbon or GHGs emitted from the conversion.

5. Conclusions

In this study, we proposed an automatic method to detect binary land cover change from a statistical test. Using the NDVI and LST parameters of MODIS data, we determined where significant land cover change had occurred in the study areas. The study areas were analyzed within super pixels as necessitated by the statistical test. A high average accuracy was obtained from this study even though the data used for validation were of a much higher spatial resolution than MODIS’s moderate resolution. Detailed change detection based on high spatial resolution imagery [38,39] will be carried out for the areas and for the periods where/when changes occurred.

The contribution of this study was to find an operational yet automatic method to determine where change in land cover has occurred. The only parameter an analyst might want to adjust is the statistical level of significance, based on their application. Additionally, using NDVI and LST, which are readily available data from satellite imagery, we present a method that can automatically detect land cover change without much other influence or input from the user. The platform used for processing the data, Google Earth Engine, is also easily accessible. Thus, our proposed method is highly operational and can find wide applicability in the monitoring of land cover change. Since our proposed method of change detection is suitable for large-scale monitoring, it can be used on the entire Great Claybelt area, again using the super pixel classification.

This study is also significant in that NDVI and LST were used together to predict land cover changes. To the best of our knowledge, such experiments are not widely seen in literature, though the features have been used individually for that purpose. The result of this study indicates that it is a valuable area to continue developing, since its outcomes are similar to the existing and established change detection methods.

An extension was made to the proposed standard change detection method. In that scenario, the NDVI and LST were used together as input parameters. The comparison of results indicated that the joint method could yield similar accuracies to the standard method. Further investigations should be made to ascertain if the joint method could be a reliable alternate to our standard model of change detection. In addition, future studies will examine the land cover and land use change in greater detail, specifically the progressive transformation of the land cover (as a “from-to” study). The “change/no change” study is important for large-scale monitoring as was performed in this work. However, opportunities exist to extend the study to explore how agricultural land has changed over the recent decades, and to examine how the land cover has changed during the agricultural seasons. In addition, as a continuation of this study, we will perform multiscale studies of the region by employing other remotely sensed data at higher spatial resolutions to observe land cover change.