Sustainable use of riverine systems and riparian habitats are directly affected by changing land use patterns  . Modeling land use patterns is an important technique for the projection of alternative pathways into the future    . Geographic Information Systems (GIS) combined with satellite remote sensing has varied application and has been recognized as a powerful tool in detecting land use and land cover change (LULC)    . Satellite data is cost effective and the information obtained from them can be used as inputs to build land use and land cover datasets  . Spatial representation of the LULC change is essential for regional planners and management  . To elucidate the optimal use of land and to provide input data for watershed models, it is necessary to have information on existing LULC change patterns  .
Geostatistics deals with problems pertaining to spatial serial data, mapping and interpolation of the data on a statistical platform that are related to a time analysis  . It has an ability of distinguishing the continuous nature of LULC and is able to detect random variations during modeling, dependent on the spatial correlation within the ecosystem  . Prediction using sample points is carried out by the spatial behavior and spatial distribution of parameters to minimize the error while doing any type of image classification  . Inverse Distance Weighting and Splines are deterministic interpolation methods to analyze change in land use patterns but these methods tend to oversimplify the results, as the spatial autocorrelation of the data is not considered  . A geostatistical method is usually preferred where sample data points can be transformed into continuous surfaces to understand the spatial autocorrelation within the data  . The parameters used for any analysis can be aggregated from pixels to object class representation using image segmentation  .
Kriging Interpolation is a very popular geostatistical method   . Ordinary Kriging estimates the mean as a constant in the searching neighborhood  . The Kriging technique has recently become very common for analyzing spaceborne data  . The values of unsampled locations are estimated by Kriging models using weighted averaging of the known sampled locations, which provide a correlation among the neighboring values that can be modelled as a function of the geographical distances between each location across the study area within the variogram  . Global and local information in predictions can be obtained from Kriging, but the ability of the variogram in describing spatial dependence is a function of the quality and quantity of the data samples  . According to a study by  , exponential models are often best-fitted semivariogram models as they use the weighted least-squares method.   and  introduced the semivariogram to remote sensing and discovered that the parameters of the variogram can be directly related to a feature in an image. The primary assumption of a geostatistical analysis when assuming spatial continuity is that samples that are located close to each other are similar than samples that are far apart  . This variation in geographic data or the spatial relation can be analyzed from a semivariogram model. An ideal semivariogram has associated features such as the lag, nugget, range and sill. The direction and distance are commonly referred to as the lag, the nugget is variability at zero distance and represents sampling and analytical errors, the range of influence in a semivariogram designates the extent beyond which autocorrelation between sampling sites is very less or zero and the sill represents the variability of spatially independent samples  .
An effective way of mapping vegetation and analyzing LULC change is using the Tasseled-Cap transformation  -  . The Tasseled-Cap transformation is a conversion of the original bands of an image into a new set of bands with defined interpretations that are useful for vegetation mapping  . The term Tasseled-Cap comes from the shape of the plot of the data that resembles a cap. A Tasseled-Cap transform is performed by taking “linear combinations” of the original image bands―similar in concept to principal components analysis   . Tasseled-Cap reduces the volume of the data without any loss of information and its spectral features are directly related to land features    . This transformation in remote sensing is the conversion of the readings in a set of data into composite values that is the weighted sums of separate data readings  . One of these weighted sums measures roughly the brightness or greenness of each pixel in the scene  .  reported that the composite values are linear combinations of the values of the separate data readings, but some of the weights are negative and others are positive. The composite values represent the degree of greenness of the pixels or the degree of yellowness of vegetation or perhaps the wetness of the soil   . Usually there are just three composite variables listed within a remote sensing interface. The Tasseled-Cap transformation of Landsat thematic mapper (TM) consists of six multispectral features, all of which could be potentially differentiated in terms of stability and change in a multitemporal dataset    . The first three features, which are brightness index, greenness index, and wetness index, respectively usually account for the most variation in a single-date image  .  analyzed Landsat data for environmental studies and found the Tasseled Cap transformation to be a consistent indicator of assessing forest change as it captures Shortwave Infrared (SWIR). Tasseled-Cap’s Greenness Index (TCGI) would ideally identify forest cover but it is less sensitive to any topographic effect  .  led a study to distinguish old growth and mature forests in the Pacific Northwest using Landsat datasets. In their study, the Tasseled-Cap brightness index did not separate old growth and mature forests due to their sensitivity to topography but the greenness and the wetness index were able to identify forest disturbances.
The primary objectives of this study were: 1) to apply Ordinary Kriging Interpolation technique to smooth TCGI values, extracted from 30 m to 60 m spatial resolution Landsat images in order to analyze spatio-temporal transformations; 2) to apply change detection techniques to the interpolated prediction maps to yield the intensity of the LULC change.
2. Site Description
The Pipestem Creek watershed, 8-Digit Hydrologic Unit Code (HUC) (10160002) sub-basin is approximately 257,178 hectares covering parts of 4 counties (Foster, Kidder, Stutsman, and Wells) in the Missouri Region―James Sub-Region of North Dakota. Of the 257,178 hectares, Stutsman County contains 65%, Wells 22%, Foster 8%, and Kidder 5%   . North Dakota’s land base is mostly a prairie land i.e. it mostly consists of grasslands  . Most of the forests in this region are found scattered along the river bank or in urban patches   . Pipestem Creek starts from the Pipestem Dam downstream to its confluence with the James River, which is about 9.01 kilometers  . The mean annual precipitation is from 330 to 559 millimeters and mean elevation ranging from 390 to 780 meters. The type of soil found at this location is Williams-Bowbells loams which is a well-drained, non-saline clay loam with calcium carbonate of about 20%  . Figure 1 shows the study area map and its location in North Dakota.
3. Materials and Methodology
3.1. Image Processing
Image classification and processing was done on a remote sensing platform― ENVI®4.5. Six Landsat images (Table 1) covering the study site were downloaded from the Global Land Cover Facility  . The images were acquired by
Figure 1. Map of study area―Pipestem Creek watershed in North Dakota.
Table 1. Landsat time series scenes used in the study.
*Path/row of the MSS image is listed according to Worldwide Reference System-1 (WRS-1) while those of TM and ETM+ are according to WRS-2.
different sensors (MSS, TM, and ETM+) and were from six different time periods, as listed in Table 1. The Landsat images were processed by applying a dark object subtraction and then converting the image digital numbers to reflectance values. Dark object subtraction was applied to remove shadows, scattering and electrical gains within the datasets, e.g.  . This was done to obtain a sound quantitative analysis of the images. Reflectance values were used to calculate several vegetation indices for each image subset. These include the Normalized Difference Vegetation Index (NDVI) and Tasseled-Cap indices. NDVI was calculated in ENVI®4.5 using the equation NDVI = (NIR − Red)/(NIR + Red)  . Band 3 was used as red and band 4 was used as near IR to generated NDVI and the output datasets were saved as floating point data type. Tasseled-Cap was calculated in ENVI®4.5 using the Transform tool where the reflectance images of years 1976 to 2015 were used as inputs. The output image generated four bands― Brightness index, Greenness index (TCGI), Wetness index and a null or Non- index. These individual bands were displayed and linked to acquire regions of interest (ROI) representing forested areas. 30 training sites were acquired. Spectral separability analysis was performed using ROI Separability tool on NDVI and TCGI, incorporating mean and standard deviation values of extreme classes in each scene, to analyze the most suitable index for differentiating between forested areas and non-forested areas using methodologies of  and  . The ROI statistical results displayed univariate statistics such as minimum value, maximum value, mean, standard deviation among other values. Since the resolution of MSS and TM/ETM+ are different, for the MSS image, a 3 × 3 pixels window was selected and for the TM/ETM+ images, 6 × 6 pixels window was selected. The window values were averaged to generate new pixel values using raster calculator in ArcMap® 10.4.1. Thus, the resolution of the images was reduced by factors of 3 and 6, to specifically fit the MSS and TM/ETM+ images respectively.
3.2. Geostatistical Analysis
To better understand the spatial structure of the imagery on a given date and location, an empirical semivariogram model such as Ordinary Kriging was applied to the sampled data. The six transformed datasets corresponding to the six years were imported to ArcMap® 10.4.1 where they were clipped to the watershed boundary shapefile. Geostatistical Analyst was used to perform Ordinary Krigingon each dataset for each model type: Gaussian, Exponential, J-Bessel, K-Bessel, Circular, and Spherical. To consider the model that best fitted the study, certain parameters were considered―1) Cross-validation scatter plot where the measured and predicted values were compared by using the difference between them, 2) Mean estimation error where the difference between the estimated and the known point values were considered, and 3) Mean standardized squared estimation error. These parameters were generated in ArcMap® 10.4.1 using Geostatistical Analysis. Based on these parameters, among other Kriging models, the Exponential model was found to be the best fit for this study. A smoothing factor of 0.2 was applied to the search neighborhood type for all the datasets. The scatterplots derived from the Exponential model for each datasets is shown in Figure 2. The
Figure 2. Cross validation scatterplots of the Exponential model of Pipestem Creek watershed datasets for years 1976 to 2015 where (a) 1976, (b) 1991, (c) 2000, (d) 2005, (e) 2011, and (f) 2015.
resultant interpolation images are shown in Figure 3. Isotropic distribution was assumed in all cases, similar to the study by  .
Figure 3. Kriging Interpolation maps of Pipestem Creek watershed based on TCGI or GI (Greenness Index) values for years 1976 to 2015.
4. Results and Discussion
4.1. Image Processing Analysis
Figure 4 is the visual representation of TCGI for the Pipestem creek watershed. The greener areas represent forest cover and the light green to white areas
Figure 4. TCGI images of Pipestem Creek watershed for years 1976 to 2015.
represents non-forested areas to barren areas. TCGI is high for the years 1976 to 1991 and it gradually decreases from 2000 to 2015. The images generated for NDVI (Figure 5) for years 1976 to 2015 showed similar results to TCGI when compared
Figure 5. NDVI images of Pipestem Creek watershed for years 1976 to 2015.
visually, but the spectral separability analysis generated low standard deviation for TCGI which is indicative of data clustering around the mean, implying data reliability. TCGI efficiently determined each class on what it resembled most in the base image. TCGI was originally designed to examine vegetation properties, its advantage lies in its ability to compare different sensors with different spectral bands, as it subsets different spectral bands to one normalized layer of TCGI values  .
4.2. Semivariogram Analysis
All six TCGI images were used for the geostatistical analysis. First, empirical semivariograms for the six periods were established. The rationale for using a semivariogram model was the similarity in spatial structure of most of its variables, gradually increasing or decreasing as a function of the increasing distance from the river until the boundary of the watershed, and the typical shape of the variogram. In the current case, the presence or absence of sill may be an indicator of presence or absence of forested areas. The level of sill may be related to the level of spatial correlation within the watershed. Therefore, if semi-va- riance reaches its maximum point (sill), beyond that, the data may not be correlated. The range may be the defining boundary of the watershed since it incorporates all the pixel values within the image that are correlated.
Results of the cross-validation analyses of the semivariogram models are presented in Tables 2-6. For an ideal model, the Mean prediction error should be near 0 (this investigates bias), Root Mean Square (RMS) prediction error should be small, average standard error should be close to RMS error, Mean-standar- dized prediction error should be near 0 and RMS standardized prediction error should be near 1, indicating that the estimated prediction uncertainty is consistent  . In Table 2, the RMS ranged from 0.043 to 0.088, average standard error ranged from 0.036 to 0.084, which is not very close to the RMS values. Mean-standardized prediction error values ranged from 0.041 to 0.057. RMS standardized prediction error values were greater than 1 for years 1976 to 1991, ranging from 1.690 to 1.721. This model did not prove the ability to reproduce the observed values accurately. In Table 3, the RMS ranged from 0.008 to 0.031,
Table 2. Semivariogram parameters for Gaussian Geostatistical Simulation models fitting the Tasseled Cap Greenness Index (TCGI) products for the pipestem creek watershed (Nugget = 0; Lag = 1000 m).
Table 3. Semivariogram parameters for Spherical Geostatistical Simulation models fitting the Tasseled Cap Greenness Index (TCGI) products for the pipestem creek watershed (Nugget = 0; Lag = 1000 m).
Table 4. Semivariogram parameters for Circular Geostatistical Simulation models fitting the Tasseled Cap Greenness Index (TCGI) products for the pipestem creek watershed (Nugget = 0; Lag = 1000 m).
Table 5. Semivariogram parameters for J-Bessel Geostatistical Simulation models fitting the Tasseled Cap Greenness Index (TCGI) products for the pipestem creek watershed (Nugget = 0; Lag = 1000 m).
average standard error ranged from 0.002 to 0.063, which is not very close to the RMS values. Mean-standardized prediction error values ranged from 0.001 to 0.027. RMS standardized prediction error values were greater than 1 for years 1976 to 2011, ranging from 1.057 to 1.501. This model also did not prove the ability to reproduce the observed values accurately. In Table 4, the RMS ranged from 0.017 to 0.038, average standard error ranged from 0.002 to 0.054, which is not very close to the RMS values. Mean-standardized prediction error values ranged from 0.001 to 0.033. RMS standardized prediction error values were greater than 1 for years 1976 to 2015, ranging from 1.591 to 1.681 which was not considered to be a best fit. In Table 5, the RMS ranged from 0.001 to 0.013, average standard error ranged from 0.012 to 0.073, which is not very close to the RMS values. Mean-standardized prediction error values ranged from 0.001 to 0.013. RMS standardized prediction error values were greater than 1 for some datasets, ranging from 0.772 to 1.613. This model did not prove the ability to reproduce the observed values accurately. In Table 6, the RMS ranged from 0.001 to 0.128, average standard error ranged from 0.010 to 0.042, which is not very close to the RMS values. Mean-standardized prediction error values ranged from 0.005 to 0.033. RMS standardized prediction error values were greater than 1 for years 1976 to 2015, ranging from 1.001 to 1.672 which implied that this model was not found to be a good fit.
In Table 7, the RMS ranged from 0.017 to 0.038, average standard error ranged from 0.016 to 0.063, which is quite close to the RMS values. Mean-stan- dardized prediction error values ranged from 0.007 to 0.041. RMS standardized prediction error values were closer to 1 ranging from 0.591 to 0.681. The slope coefficient was very close to unity and the intercept coefficient was very close to zero, proving the ability of the chosen exponential model to reproduce the observed values   . The least-squares measure of fit was used, incorporating exponential model, as shown in Figure 3.
All semivariograms were processed with 20 lags of 1000 m each. Because of the irregular distribution of forests and non-forests, data values exactly separated by 1000 m could not be expected, thus the range of 1000 m to 20,000 m was selected.
Table 6. Semivariogram parameters for K-Bessel Geostatistical Simulation models fitting the Tasseled Cap Greenness Index (TCGI) products for the pipestem creek watershed (Nugget = 0; Lag = 1000 m).
Table 7. Semivariogram parameters for Exponential Geostatistical Simulation models fitting the Tasseled Cap Greenness Index (TCGI) products for the pipestem creek watershed (Nugget = 0; Lag = 1000 m).
Lag values were determined by trial and error process to optimize the above- mentioned criteria. Ordinary Kriging interpolation maps were produced based on the exponential models with the parameters presented in Table 7. Figure 2 shows the scatterplots generated for cross validation of the Exponential model. These graphs show the spatial correlation of the data. Most of the data values lie along the line in the scatterplots indicating how closely related the data is. During ground truthing and field observations, no evidence was found to support an anisotropic pattern, as in  study, that may explain the direction of forest reduction or the direction in which the agricultural lands are increasing. So, isotropic distribution was assumed in all cases. Figure 3 depicts the final results for the distribution of the TCGI values for the six periods. The dark-red areas in the images are related to forested areas. The surrounding light red and yellow belts represent a mixed zone where forested areas and non-forested areas overlay each other or create a stable spectral balance. The zone colored by blue tones is considered to be non-forested areas that include mostly agricultural land. Forested areas are more concentrated in the 1976 image and is gradually seen to reduce for the rest of the images through 2015.
TCGI was selected to describe the spatial surface patterns since it produced the best contrast in terms of separability among the spectral indices. It produced the best contrast in terms of separability among all examined spectral indices. TCGI was able to compare between the different sensors with different spectral bands, as it reduced their different spectral bands to one normalized layer of TCGI values. The semivariance analysis was found to be a suitable method for gaining insight to the spatial structure present in the imagery for a given date and location. The similarity between the shape of the semivariogram and the directional change of the environmental variables is a logical reason for using this method. The Kriging interpolation technique using the Exponential Geostatistical Simulation model was used as a smoothing filter in which each pixel was being replaced with the solution for the semivariogram equation (exponential model in the current case) calculated from all other pixels in the image. As a result, it reduced spatial errors and fine scale variability and helped to better identify the transition from forest to non-forested areas.
This research is supported primarily by US Forest Service, award #13-DG- 11010000-004, and CFDA Cooperative Forestry Assistance, CFDA Number 10.664, award year FY 2016, NDSU Project FAR0026072 and the ND View Scholarship 2016. The research was carried out in a facility that was funded by NSF EAR Division of Earth Sciences Award #: 0963486. This work was conducted at the NDSU Environmental Geomechanics Research Facility.
 Rozario, P.F., Oduor, P., Kotchman, L. and Kangas, M., (2016) Quantifying Spatiotemporal Change in Land Use and Land Cover and Assessing Water Quality: A Case Study of Missouri Watershed James Sub-Region, North Dakota. Journal of Geographic Information System, 8, 663-682.
 Lambin, E., Turner, B., Geist, H., Agbola, S., Angelson, A., Bruce, J., Coomes, O., Dirzo, R., Fischer, G., Folke, C. and George, P. (2001) Our Emerging Understanding of the Causes of Land Use and Cover Change. Global Environmental Change, 11, 261-269.
 Ehlers, M., Jadkowski, M.A., Howard, R.R. and Brostuen, D.E. (1990) Application of SPOT Data for Regional Growth Analysis and Local Planning. Photogrammetric Engineering and Remote Sensing, 56, 175-180.
 Méaille, R. and Wald, L. (1990) Using Geographical Information System and Satellite Imagery within a Numerical Simulation of Regional Urban Growth. International Journal of Geographical Information System, 4, 445-456.
 Treitz, P.M., Howarth, P.J. and Gong, P. (1992) Application of Satellite and GIS Technologies for Land Cover and Land Use Mapping at the Rural-Urban Fringe: A Case Study. Photogrammetric Engineering and Remote Sensing, 58, 439-448.
 Piccini, C., Marchetti, A. and Francaviglia, R. (2014) Estimation of Soil Organic Matter by Geostatistical Methods: Use of Auxiliary Information in Agricultural and Environmental Assessment. Ecological Indicators, 36, 301-314.
 Harris, P.M. and Ventura, S.J. (1995) The Integration of Geographic Data with Remotely Sensed Imagery to Improve Classification in an Urban Area. Photogrammetric enginEering and Remote Sensing, 61, 993-998.
 Eldeiry, A.A. and Garcia, L.A. (2010) Comparison of Ordinary Kriging, Regression Kriging, and Cokriging Techniques to Estimate Soil Salinity Using LANDSAT Images. Journal of Irrigation and Drainage Engineering, 136, 355-364.
 Robinson, T.P. and Metternicht, G. (2006) Testing the Performance of Spatial Interpolation Techniques for Mapping Soil Properties. Computers and Electronics in Agriculture, 50, 97-108.
 Wang, L., Sousa, W.P. and Gong, P. (2004) Integration of Object-Based and Pixel-Based Classification for Mapping Mangroves with IKONOS Imagery. International Journal of Remote Sensing, 25, 5655-5668.
 Huang, C., Wylie, B., Yang, L., Homer, C. and Zylstra, G. (2002) Derivation of a Tasselled Cap Transformation Based on Landsat 7 At-Satellite Reflectance. International Journal of Remote Sensing, 23, 1741-1748.
 Mishra, U., Lal, R., Slater, B., Calhoun, F., Liu, D. and Van, M.M. (2009) Predicting Soil Organic Carbon Stock Using Profile Depth Distribution Functions and Ordinary Kriging. Soil Science Society of America Journal, 73, 614-621.
 Miller, J., Franklin, J. and Aspinall, R. (2007) Incorporating Spatial Dependence in Predictive Vegetation Models. Ecological Modelling, 202, 225-242.
 Woodcock, C.E., Strahler, A.H. and Jupp, D.L.B. (1988) The Use of the Variogram in Remote Sensing: I. Scene Models and Simulated Images. Remote Sensing of Environments, 25, 323-348.
 Woodcock, C.E., Strahler, A.H. and Jupp, D.L.B. (1988) The Use of the Variogram in Remote Sensing: II. Real Digital Images. Remote Sensing of Environments, 25, 349-379.
 Jin, S. and Sader, S.A. (2005) Comparison of Time Series Tasseled Cap Wetness and the Normalized Difference Moisture Index in Detecting Forest Disturbances. Remote Sensing of Environment, 94, 364-372.
 Bauer, M.E., Burk, T.E., Ek, A.R., Coppin, P.R., Lime, S.D., Walsh, T.A., Walters, D.K., Befort, W. and Heinzen, D.F. (1993) Satellite Inventory of Minnesota Forest Resources. National Aeronautics and Space Administration, Washington DC, 87-106.
 Cohen, W.B. and Spies, T.A. (1992) Estimating Structural Attributes of Douglas-Fir /Western Hemlock Forest Stands from Landsat and SPOT Imagery. Remote Sensing of Environment, 41, 1-17.
 Cohen, W.B., Spies, T.A. and Fiorella, M. (1995) Estimating the Age and Structure of Forests in a Multi-Ownership Landscape of Western Oregon, U.S.A. International Journal of Remote Sensing, 16, 721-746.
 Collins, J.B. and Woodcock, C.E. (1996) An Assessment of Several Linear Change Detection Techniques for Mapping Forest Mortality Using Multitemporal Landsat TM Data. Remote Sensing of Environment, 56, 66-77.
 Dymond, C.C., Mladenoff, D.J. and Radeloff, V.C. (2002) Phenological Differences in Tasseled Cap Indices Improve Deciduous Forest Classification. Remote Sensing of Environment, 80, 460-472.
 Fiorella, M. and Ripple, W.J. (1995) Analysis of Conifer Forest Regeneration Using Landsat Thematic Mapper Data. Photogrammetric Engineering and Remote Sensing (American Society for Photogrammetry and Remote Sensing), 59, 1383-1388.
 Fiorella, M. and Ripple, W.J. (1995) Determining Successional Stage of Temperate Coniferous Forests with Landsat Satellite Data. Photogrammetric Engineering and Remote Sensing (American Society for Photogrammetry and Remote Sensing), 59, 239-246.
 Franklin, J.F., Spies, T.A., Van, P.R., Carey, A.B., Thornburgh, D.A., Berg, D.R., Lindenmayer, D.B., Harmon, M.E., Keeton, W.S., Shaw, D.C. and Bible, K. (2002) Disturbances and Structural Development of Natural Forest Ecosystems with Silvicultural Implications, Using Douglas-Fir Forests as an Example. Forest Ecology and Management, 155, 399-423.
 Skakun, R.S., Wulder, M.A. and Franklin, S.E. (2003) Sensitivity of the Thematic Mapper Enhanced Wetness Difference Index to Detect Mountain Pine Beetle Red-Attack Damage. Remote Sensing of Environments, 86, 433-443.
 Watkins, T. (2005) The Tasseled Cap Transformation in Remote Sensing.
 Crist, E.P. and Cicone, R.C. (1984) A Physically-Based Transformation of Thematic Mapper Data—The TM Tasseled Cap. IEEE Transactions on Geoscience and Remote Sensing, 22, 256-263.
 Crist, E.P., Laurin, R. and Cicone, R.C. (1986) Vegetation and Soils Information Contained in Transformed Thematic Mapper Data. Proceedings of International Geoscience and Remote Sensing 86 Symposium, Zurich, 1986, 1465-1470.
 Cohen, W.B., Fiorella, M., Gray, J., Helmer, E. and Anderson, K. (1998) An Efficient and Accurate Method for Mapping Forest Clearcuts in the Pacific Northwest Using Landsat Imagery. Photogrammetric Engineering and Remote Sensing, 64, 293-300.
 Cohen, W.B. and Fiorella, M. (1998) Comparison of Methods for Detecting Conifer Forest Change with Thematic Mapper Imagery. In: Lunetta, R.S. and Elvidge, C.D., Eds., Remote Sensing Change Detection: Environmental Monitoring Methods and Applications, Taylor & Francis Press, 89-102.
 Mondal, A., Khare, D., Kundu, S., Mondal, S., Mukherjee, S. and Mukhopadhyay, A. (2016) Spatial Soil Organic Carbon (SOC) Prediction by Regression Kriging Using Remote Sensing Data. The Egyptian Journal of Remote Sensing and Space Science, 20, 61-70.
 Yuan, D., Elvidge, C.D. and Lunetta, R.S. (1998) Survey of Multispectral Methods for Land Cover Change Analysis. In: Lunetta, R.S. and Elvidge, C.D., Eds., Remote Sensing Change Detection: Environmental Monitoring. Methods and Application, Taylor and Francis Ltd., 21-39.
 Yeh, A.G.O. and Li, X. (1997) An Integrated Remote Sensing and GIS Approach in the Monitoring and Evaluation of Rapid Urban Growth for Sustainable Development in the Pearl River Delta, China. International Planning Studies, 2, 193-210.
 Madurapperuma, B., Rozario, P., Oduor, P. and Kotchman, L. (2015) Land Use and Land Cover Change Detection in Pipestem Creek Watershed, North Dakota. International Journal of Geomatics and Geosciences, 5, 416-426.
 Rozario, P.F., Oduor, P., Kotchman, L. and Kangas, M. (2017) Transition Modeling of Land-Use Dynamics in the Pipestem Creek, North Dakota, USA. Journal of Geoscience and Environment Protection, 5, 182.
 Haugen, D.E., Harsel, R., Bergdahl, A., Claeys, T., Woodall, C.W., Wilson, B.T., Crocker, S.J., Butler, B.J., Kurtz, C.M., Hatfield, M.A. and Barnett, C.H. (2013) North Dakota’s Forests 2010. U.S. Department of Agriculture, Forest Service, Northern Research Station, Washington DC, 52.
 Global Land Cover Facility (2017) Landsat Imagery.
 Kaufman, Y.J. and Remer, L.A. (1994) Detection of Forests Using MID-IR Reflectance—An Application for Aerosol Studies. IEEE Transactions on Geoscience and Remote Sensing, 32, 672-683.
 Lasaponara, R. (2006) Estimating Spectral Separability of Satellite Derived Parameters for Burned Areas Mapping in the Calabria Region by Using SPOT Vegetation Data. Ecological Modelling, 196, 265-270.
 Friedland, C.J., Joyner, T.A., Massarra, C., Rohli, R.V., Trevino, A.M., Ghosh, S., Huyck, C. and Weatherhead, M. (2016) Isotropic and Anisotropic Kriging Approaches for Interpolating Surface-Level Wind Speeds across Large, Geographically Diverse Regions. Geomatics, Natural Hazards and Risk, 1, 1-18.
 Weng, Q. (2001) A Remote Sensing? GIS Evaluation of Urban Expansion and Its Impact on Surface Temperature in the Zhujiang Delta, China. International Journal of Remote Sensing, 22, 1999-2014.
 Westmoreland, S. and Stow, D.A. (1992) Category Identification of Changed Land—Use Polygons in an Integrated Image Processing/Geographic Information System. Photogrammetric Engineering and Remote Sensing, 58, 1593-1599.