Detecting LandCover Change using Stochastic Simulation Models and Multivariate Analysis of Multi-Temporal Landsat Data for Cass County , North Dakota

Understanding forest transiting at wildland-urban interfaces offers a glimpse into the effect of anthropogenic activities that may threaten biota.We examined forest conversion from 2006 to 2011 at urban-wildland fringes in Cass County, North Dakota. Grid data from the National Agricultural Statistic Service, published by USDA, was used as preliminary inputs to ascertain land-use and land-cover dynamics. Markovian transition probabilities were derived for each pair of years from 2006 to 2011. These transition probabilities were further subjected to multivariate analysis to detect forest change in one-year time steps. From this study, pairwise combinations of years yielded two distinct statistical groups. The first group comprised of seven pairs of year combinations displaying high transition probability of unchanged forest (0.54 ≤ Pff ≤ 0.68), while the second group comprised of eight pairs of year combinations and showed a low transition probability of unchanged forest (0.26 ≤ Pff ≤ 0.37). A third group displayed comparatively high transition probabilities of forest transiting to non-forest (0.26 ≤ Pfnf ≤ 0.36), such as forest to row crops, with an increasing trend over time. We also generated the forest cover in relation to soil characteristics. We can surmise that forest cover at poorly drained soils showed a higher distribution, which could be due to the unsuitability of this soil for crop cultivation. The results of this study on how land-cover has changed in Cass County for the last six years could be used by policy makers and forest managers in applying BMPs (Best Management Practices).


Introduction
Land-use and land-cover change (LULCC) is a dynamic, widespread and accelerating process.LULCC is mainly driven by natural phenomena and anthropogenic activities, which in turn have driven changes that impact natural ecosystems (Ruiz-Luna & Berlanga-Robles, 2003;Turner & Ruscher, 2004).Anthropogenic-induced conversion of semi-natural managed forest into non-forested agricultural land is currently driven by global demand for cash crops.Intensification of agricultural land usage has been identified as one of the most common forms of land-cover modification (Matson et al., 1997).An example is within the eastern United States.Agricultural expansion has caused a decrease in the forest cover by more than 4% (Drummond & Loveland, 2010).
To design appropriate polices for sustainable development, it is fundamental to understand how technical innovations and agricultural intensification in forest frontier areas influences land-use at forest margins (Maertens et al., 2006).Accurate and up-to-date land-cover change information is necessary to understand these influences.The land-cover information is essential for developing long-term strategies for managing priority landscapes, delineating and prioritizing forest resources, assessing risks, conditions and trends in forests.Secondarily, this information can be used to assess the environmental consequences such as climate change, biodiversity conservation, ecosystem assessment, and environmental modeling (Giri et al., 2005).
Adhoc forest resource mapping can be applied as a powerful tool to identify forest resource threat patterns and also to manage natural forest resources sustainably.Remote sensing has the capability of capturing LULC data, extracting the LULCC information from satellite data requires effective and automated change detection techniques (Roy et al., 2002).LULCC mapping and detection of change are of paramount importance to planners, geographers, environmentalists and policy makers for the sustainable development of natural resources (Abbas et al., 2010).The basic premise for using remote sensing data to detect change is that the process can identify change between two or more dates that is uncharacteristic of normal variation (Kasereka, 2010).Modeling land-use and land-cover change has become popular since it addresses when, where and why land-use and land-cover change occurs (Baker, 1989;Lambin, 1997;Theobald & Hobbs, 1998).
Markov Chain Monte Carlo simulation models of LULCC can offer a glimpse in understanding the changes of landscape (Oduor et al., 2012) over a time scale and can then extend those patterns into the future allowing for prediction (Brown et al., 2000;Clarke & Gaydos, 1998).Markovian models are one of the most widely used approaches for predicting forest cover change with discrete parameter space (Acevedo et al., 1996;Hall et al., 1991;Logofet & Lesnaya, 2000;Wu et al., 2006).A Markov chain consists of a system with a series of changes from one state to another over time, which is measured in discrete intervals (Kemeny & Snell, 1976).The Markov chain models are employed by ecologists to study the vegetation dynamics and the succession process, to evaluate conservation intervention and to simulate forest dynamic transformation patterns (Mondal & Southworth, 2010;Yuliang et al., 2008;Leps, 1988).
This study is primarily based on geo-spatial modeling and application of stochastic methods to address year-to-year stratified changes in land-cover with a special emphasis on forestry in Cass County, North Dakota, which is a rapidly growing metropolitan area.Our main goal was to estimate transition probabilities for forest to non-forest conversion using random-walk stochastic processes and Markovian models.This pixel-wise transition model approach can be used to predict forest-cover changes and ongoing factors associated with changing forest management at landscape scales in Cass County.
The research questions addressed in this paper were: We hypothesized that the demand for agricultural crops, such as soybean, and population growth changed the forest land into agricultural crops from 2006 to 2011 in Cass County.In addition, we hypothesized that forest cover changed based on the suitability of soil types for growing crops.For instance, land-cover associated with well-drained soil type would be more likely converted to agricultural cropland than poorly-drained soils.

Study Area
Cass County is located in southeastern North Dakota in the Red River valley bounded by 46° N and 47° N latitude and 96° W and 97° W longitude (Figure 1).The Red River of the North establishes Cass County's eastern border, separating it from Minnesota.The county experienced a great population growth at the beginning of the 20 th century, averaging 19.4% growth for the first 30 years of the century and another noticeable spike was a 21.6% (149,778 in population count) increase in 2010 (US Census Bureau, 2010).Cass County contains over 457, 699 ha of total land area, covering a nearly square area roughly 71 km west to east by 68 km north to south.The county has over 52,000 parcels of land extending 454, 867 ha (Cass County Government, 2005).The extent of forest cover in Cass County isestimated to be about 6, 691ha (http://ndfsdss.ndsu.nodak.edu/maps/map_data/County/Cass.pdf).The average annual temperature for Fargo is 41.5° F with an average yearly precipitation of 53.82 cm (V.Godon & N. Godon, 2002).Four general land-use classes (agricultural, rural non-farm, small city and metropolitan area) can be identified within Cass County.The majority of land is used for crop production, specifically soybeans, wheat, and barley.

Land-Cover Maps
We generated a GIS database based on the analysis of a sequence of satellite images dating from 2006 to 2011.National Agricultural Statistics Service (NASS) datasets originally derived from classified LandSat ETM + satellite imagery for Cass County were imported into ArcMap-ArcInfo ® 9.3.Theneach dataset was exported into ArcMap-ArcInfo ® 9.3 to obtain grid datasetswith 30 m spatial resolution.Raster reclassification was performed where minor NASS classes were collapsed into the major categories, for example, class 22 (durum wheat) and class 21 (barley) were reclassified as grains hay seeds (Oduor et al., 2012).The total number of classes was 8, namely, (i) row crops (ii) grains, hay, seeds (iii) other crops (iv) idle cropland/fallow (v) grass, pasture, non-agriculture (vi) forest (vii) urban/developed and (viii) water.Some of these original Landsat ETM+ images were classified using both NASS classes and NLCD classes (National Land Cover Database).We collapsed those classes also into eight major classes.For instance, class 141 (deciduous forest), class 142 (evergreen forest) and class 143 (mixed forest) were reclassified as forest.

Markov Model
A finite time step first-order Markov process was generated as a probability matrix(t) representing the mutual pixel wise transition from one category defined or redefinedas aforementioned.The basic equation used was, with time homogeneity transition probability expressed as,  , solved as a set of transition matrices where p ij are the elements in the matrix of transition probabilities (Wu et al., 2006;McCauley, 2007).Transition probabilities were derived using SemGrid (a freely available program).This was done by first converting each derived raster grid to ASCII files, which are the generic input file format for SemGrid.Transition probabilities and pertinent stochastic parameters for each selected pair of years, for example, 2006-2007, 2006-2008 and 2006-2009, were then derived by running stochastic analyses on each pairing.The transition probabilities are calculated using a Markovian random process algorithm, which factors in corresponding pixel centroids and pairwise relationships.Brown et al. (2000) subsidiary two state variation of transition probabilities was also adopted for this study to determine probabilities of lumped areas converting from forest to non-forest cover.

Multivariate Analysis
In order to trend results generated from Markov chain analyses, we introduced a new method to illustrate the in-situ variation between each pair of transition states using ordination and classification techniques.The pixels that changed for each land-use and land-cover category for each pair of years were subjected to multivariate analysis using PC-ORD (1999).Cluster Analysis and Canonical Correspondence Analysis (CCA) were performed to identify the trends and relationships between pixel change in each land-use and land-cover class and time steps (Madurapperuma 2013;Madurapperuma et al., 2011).The dendrogram was constructed using the distance between objects index measured by Bray-Curtis, and groups were linked by the nearest neighbor method (Bray & Curtis, 1957).The Bray-Curtis method was adopted due to its retention of sensitivity in more heterogeneous data sets and because it weights fewer outliers compared to Euclidean distance (Roberts, 1986).
An ordination technique was used to reveal subtle intrinsic patterns that may not be readily evident in the dataespecially on a multidimensional scale.The primary data matrix for the CCA was assembled using pixel change in LULC of forest to non-forest (P fnf ) and non-forest to forest (P nff ) for each time step.The secondary data matrixwas generated using the forests' extent in relation to underlying soil characteristics, population and housing density data.
Although the time steps (rows) of the two matrices are equal, the variables (columns) of the two matrices differ.The eight variables for the primary data matrix were: row crops; grains, hay, seeds; other crops; idle cropland; grass, pasture, non-agriculture; forest; urban/developed; and water.The six variables for the secondary data matrix were forest acreage, moderately well drained to somewhat poorly drained soil, moderately well drained soil, poorly drained soil and somewhat poorly drained soil, population density and housing density.
The two data matrices were overlaid (the joint plot scale was set to 0.10) through canonical corresponding analysis and the results given by the ordination diagram.The results of the joint plots in the ordination diagram show the relationship between soil characteristics and LULC types through a diagram of radiating lines.The angle and length of each line indicates the direction and strength of the relationship (McCune & Grace, 2002).Bray-Curtis distance matrices created by the cluster analysis were used to separate groups in the ordination diagrams (Bray & Curtis, 1957).

Soil-Forest Associations
We investigated the spatial distribution of forests with respect to soil types.Forest acreage in each time period was derived from the classified Landsat ETM+ image using a raster-to-vector conversion algorithm.Soil data was obtained from the Soil Survey Geographic (SSURGO) Database of the Natural Resources Conservation Service.Soil types grouped by soil characteristics for example well drained soils, were overlain with the forests layer to tease out any relationships between soil types and forest vibrancy.The forest acreage in each of the soil types (for example moderately well drained soil (MWD), moderately well drained to somewhat poorly drained soil (MWPD), somewhat poorly drained soil (SPD) and poorly drained (PD) was derived (Table 1).The data gathered in the study was further analyzed using a Canonical Correspondence Analysis (Figure 4).We examined the ordination diagram from this analysis in relation to the groups identified in the cluster analysis (Figure 3).As the diagram showing axis 1 and 2 (Pearson correlation 0.98 & 0.96 respectively) gave the best separation of land-use and land-cover data (Figure 4), we therefore adopted it to gauge interrelationships among variables.Table 2 shows the forest to non-forest transition (P fnf ) probabilities for pairwise comparisons of images according to the groups separated from the cluster analysis (Figure 3) during the period 2006 to 2011 for Cass County.Table 2 shows the forest to non-forest transition (P fnf ) probabilities for pairwise comparisons of images according to the groups separated from the cluster analysis (Figure 3) during the period 2006 to 2011 for Cass County.

Land-Use/Cover Change
Results showed that unchanged forest (P ff ) transition probabilities (0.5400-0.6806) were high for the first cluster group.In contrast, the transition probabilities for other crops and idle cropland were low.Of the three main axes in the ordination diagram (Figure 6), Axes 1 and 2 were selected to describe the inherent variation among variables due to the high Pearson correlation (0.99 & 0.95 respectively).The first group, which is similar to the cluster analysis, is separated at the right-most side of Axis 1.The first group is closely associated with forest, MWD soil and SPD soil and human population.The second group separates at the left-most side of Axis 1 and is associated with water and grass, pasture, non-agriculture land-use.
The NFF transition probabilities, where pixel changing from forested to non-forested, of each land-use class for 2006-2011 are listed in Table 3.The first and the second group showed very low transition probability values for all land-use classes except forest.The transition probabilities for unchanged forest were much greater in the first group (0.5400≤ P ff ≤ 0.6806) than the second group (0.2634 ≤ P ff ≤ 0.3659).Key to transition states: 1 row crops, 2 grains, hay, seeds, 3 other crops, 4 idle cropland, 5 grass, pasture, non-agriculture, 6 forest, 7 urban/developed, 8 water.

Forest Soil Characteristics
The relationship between soil pattern and forest cover in each year is shown in Figure 7. Forest cover occurrence in MWD soil and PD soil was greater than forest cover in MWPD soil and SPD soil.PD had the highest forest cover (58%), followed by WMD (25%), SPD (15%), and MWPD (2%).
Figure 7. Forest acreage in relation to soil characteristics for the study site

Land-Use/Cover Change
The multi-temporal satellite image analysis of land-use and land-cover change in Cass County showed that more land was used for cultivation of row crops from 2006 to 2011.Conversely, less land was used for cultivation of grains, hay, seeds and other crops from 2006 to 2011.Interestingly, forest cover increased slightly in 2007, but overall remained consistent from 2006 to 2011.Neau et al. (2011) reported that the major changes occurring across the U.S. landscape were the conversion of lands from traditional small grain production to corn and soybean production.North Dakota is an agricultural state and soybeans are a major row crop cultivated in this region (VanWechel et al., 2003).According to an economic census in North Dakota, the Red River Valley area, which includes Cass County and Richland County, was the major region that cultivated soybean (Bangsund & Leistritz, 1999).The demand for soybean oil for production of biodiesel could explain the increase in acreage of soybeans in Cass County (VanWechel et al., 2003).In contrast, spring wheat production in Cass County declinedfrom 2006 to 2011, which may be due to a delay in snow melt and major river flooding.These conditions are unfavourable for the cultivation of spring wheat.
North Dakota and Minnesota are the main regions, which cultivate sugarbeets in the US, but from 2003 to 2010, the industry decreased planted acreage.The decrease was greater in North Dakota (18%) than in Minnesota (11%).The decline in sugarbeet production associated with high expenditure for processing and marketing activities (Bangsund & Leistritz, 2004).

Multivariate Analysis
We derived first-order Markov chain models, which can serve as an indicator of the direction and magnitude of land-cover change in the future as well as a quantitative description of change in the past (Weng, 2002).An important aspect of change detection is to determine which land-use class is changing and how it is changing (Shiferaw, 2011).This information reveals both the desirable and undesirable changes and classes that are relatively stable overtime, which is useful for management decisions.
Two distinct transition states were observed for the entire time step.One transition state comprised of 2006 -2009 images compared to images 2010 and 2011, which formed a main group in the cluster and ordination diagrams.The forest acreage in 2010 and 2011 is more or less similar (Figure 2) and therefore forest to non-forest transition for those combinations of years make a distinct cluster.For these combinations of years, forest cover transitioned to rowcrops (Table 2).
Cass County is a metropolitan area with a 21.6% population increase recorded from 2000 to 2010 (US Census Bureau, 2010).Even though the population increased over the decade, forest cover did not change considerably.The reasons could be attributed to an increased interest by landowners including farmers in community forestry programs and also change in the main economic base of the city from farming to other white collar jobs.For instance, the agricultural community together with the forest service in North Dakota established windbreaks and shelterbelts (Kotchman, 2010), which have worked to increase forest cover.In addition, Cass County is a soil conservation district and therefore the county is required to implement natural forest, urban and community forestry conservation practices.
Also, according to the Agassiz Lake resource surveys (Lake Agassiz resource conservation and development area plan council, 2008), in Cass County, changing land-use patterns to the production of small grains and row crops caused severe erosion problems.Wind erosion continues to be the most serious conservation problem occurring in cropland areas.Windbreaks are therefore established in agricultural settings in order to protect crops.This may further promote forest acreage remaining unchanged with diversified cropping.

Forest Soil Characteristics
Understanding the relationship between LULC practices and soil characteristics is important for controlling factors associated with deforestation or reforestation.Therefore, many plant geographers are interested in the relationship between vegetation and soil characteristics.For example, Iverson (1988) used GIS techniques in assessing relationships between vegetation types and the soil characteristics to compare the soil and landscape attributes of Illinois with its historic vegetation, current land-use, and patterns of land-use change.
In our study, historic forest change in terms of soil characteristics was analyzed.
Forest occurrence in PD soil showed a higher distribution, and this soil type was unsuitable for cropland because of its texture, that is, silty clay, silty clay loam, and sandy loam, which may hold excess water.As Hopkins et al. (2012) suggested, saturated soil contributed to higher soil salinity, which prevented planting in the Dakotas.
Moreover, he reported that the Cass County soil contains high exchangeable sodium percentage at shallower depths, which apparently is undesirable for crop cultivation.
The forest occurrence at the MWD soil is also noticeable.The main soils associated with this class were silt loam, loam, sandy loam, and silty clay loam.Ideally, these soils are suitable for agricultural crops and therefore the forestlands may be confined to shelterbelts and/or windbreaks.
Finally, lands suitable for agriculture seem to be characterized by SPD soil and MWPD soil.The main soils associated with these groups were loam, clay loam, and silt loam.Merrill et al. (2013) reported crop productivity results for dry pea, spring wheat and maize in North Dakota, which showed that loam/clay loam, and sandy loam soils effected for crop productivity.For example, spring wheat yields in the loam/clay loam soil were significantly higher than sandy loam soil, while the opposite result was recorded for maize.Therefore, high nutrient content and water retention capacity of those soils contributed to growth and yield of such crops.

Conclusion
This research adopts concepts of random-walk stochastic processes and multivariate analysis to investigate land-cover change from a spatial evolution perspective in Cass County, North Dakota.In general, row crops showed an increasing trend, while grains, hay, seeds and other crops showed a declining trend from 2006 to 2011.
Of the cultivated crops, soybean showed a tremendous increase, while spring wheat and sugarbeet showed a declining trend.Interestingly, forest cover showed a stable trend.Although Cass County is one of the major metropolitan areas in North Dakota with an increasing population, our results showed that this increasing population is not resulting in a depletion of forest lands.
Multivariate analysis resulted in two distinct groups: one group with unchanged forest transition probabilities that were comparatively high (0.54 ≤ P ff ≤ 0.68), and the second group with low transition probabilities of unchanged forest (0.26 ≤ P ff ≤ 0.37).
Our results showed that the forests in Cass County are characterized by poorly drained soil, which is apparently unsuitable for growing crops due to water logging conditions.Overall, these findings may be useful to managers and foresters in North Dakota to assess their stewardship and resource management program effectiveness within Cass County.
(i) How has land-cover changed in Cass County from 2006 to 2011? (ii) What are the primary drivers of land-use and land-cover change in Cass County?(iii) Is there a relationship between soil characteristics and forest cover change?

Figure 1 .
Figure 1.A detailed map of the study area

Figure 2
Figure2shows the changes of land-use and land-cover types in Cass County from 2006 to 2011.According to our landsat-based analysis, row crops showed an increasing trend (R 2 = 0.93) of cultivation (Figure2).Soybean was the main row crop cultivated in Cass County with 515,704 acres planted in 2006, increasing to 598,929 acres in 2011.In contrast, grains, hay, seeds showed a decreasing trend from 2006 to 2011 (R 2 = 0.80).Spring wheat was the leading crop of grains, hay, seeds cultivated at 183,029 acres in 2006 and decreased to 107,642 acres in 2011.

Figure 2 .
Figure 2. Land-use land-cover change in Cass County from 2006 to 2011 in North Dakota

Figure 3 .
Figure 3. Dendrogram for the probabilities of forest to non-forest transition using Bray Curtis similarity values for Cass County for the 2006 to 2011 periodTransition probabilities for 15 pairs of year combinations were used to produce the cluster analysis.

Figure 4 .
Figure 4. Canonical corresponding analysis ordination for the probabilities of forest to non-forest transition using pairs of year combinations, land-use classes and secondary data such as population density, housing number and forest acreages in each soil types.Joint plot was performed at the 0.01 cut off levels.Two to four letter codes represent the variables such as population (Popu), housing number (Hous), somewhat poorly drained soil (SPD), moderately well drained (MWD), poorly drained (PD) and moderately well drained to somewhat poorly drained soil (MWPD).
The second group consisted of eight pairs of years.Of those pairs, image 2007 compared to images 2010-2011, image 2008 compared to images 2010-2011, and image 2009 compared to images 2010-2011 showed comparatively high transition probabilities (0.3572-0.3659) for unchanged forest compared to the transition probabilities for image 2006 compared to images 2010-2011 (0.2634-0.2832).In contrast, the transition probabilities for row crops for image 2006 compared to images from 2010-2011 were higher than other land-use and land-cover categories (0.3579 and 0.3022).

Figure 5
Figure 5 shows the results of the cluster analysis for the P nff transition from 2006 to 2011.Two distinguishable groups were identified, which was similar to the previous cluster analysis.The first group clearly separated at 92% similarity and included image 2006 compared to images analyzed for years 2007-2009, image 2007 compared to images analyzed for years 2008-2009, image 2008 compared to image 2009, and image 2010 compared to image 2011.The subgroup of the first cluster, which separated at 98% similarity, included image 2006 compared to images 2008-2009.The second group separated at 88% similarity and included the 2006 image compared to images analyzed for years 2010-2011, image 2007 compared to images 2010-2011, image 2008 compared to images 2010-2011, and image 2009 compared to images 2010-2011.The sub group of the second main cluster separated at 96% similarity and included image 2006 compared to images 2010-2011.

Figure 5 .
Figure 5. Dendrogram for the probabilities of non-forest to forest transition using Bray Curtis similarity values for Cass County for the 2006-2011 period.Transition probabilities for 15 pairs of year combinations were used to make the cluster analysis

Figure 6 .
Figure 6.Canonical corresponding analysis ordination for the probabilities of non-forest to forest transition using pairs of year combinations, land-use classes and secondary data such as population density, housing number and forest acreages in each soil types Joint plot was performed at the 0.01 cut off levels.Two to four letter codes represent the variables such as population (Popu), housing number (Hous), somewhat poorly drained (SPD)soil, moderately well drained (MWD) soil, poorly drained (PD)soil and moderately well drained to somewhat poorly drained (MWPD) soil.

Table 1 .
Change in forest cover in relation to different soil types in Cass County

Table 3 .
Land use transitional probabilities (P nff ) of each sub-period for Cass County, North Dakota