MULTI-LAYER PERCEPTRON – MARKOV CHAIN BASED GEOSPATIAL
ANALYSIS OF LAND USE AND LAND COVER CHANGE: A CASE STUDY OF
STONEY CREEK WATERSHED, BC, CANADA
by
Lei Shen
B.E., College of Science & Technology, North China Electric Power University, China, 2016

THESIS SUBMITTED IN PARTIAL FULFILLMENT OF
THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
IN
NATURAL RESOURCES AND ENVIRONMENTAL STUDIES
(ENVIRONMENTAL SCIENCE)

UNIVERSITY OF NORTHERN BRITHISH COLUMBIA
April 2019
© Lei Shen, 2019

ABSTRACT
This thesis study analyzed the land use and land cover (LULC) changes in Stoney Creek
Watershed, BC, Canada using the combination of remote sensing, GIS and modeling
approaches. The Object-Based Image Analysis (OBIA) tool in PCI Geomatica 2017 software
was applied to generate unsupervised classification LULC maps using Landsat TM and OLI
images of the years 1986, 1999 and 2016. Various band ratio were computed to improve
different classification results. Esri ArcMap 10.5 was used to produce all the LULC maps for
subsequent modeling. A modeling method using Multi-layer perceptron (MLP) neural network
and Markov Chain (MC) was performed to predict LULC changes in 2026, using hard and soft
prediction results. The outcomes of this study could provide valuable information of LULC
patterns and dynamics for supporting both environmental and economic development in this
area.

i

Table of Contents
ABSTRACT ............................................................................................................................... i
Table of Contents .......................................................................................................................ii
List of Figures ........................................................................................................................... iv
List of Tables ............................................................................................................................ vi
Acknowledgement ...................................................................................................................vii
CHAPTER 1

INTRODUCTION .................................................................................... 1

1.1

Background .......................................................................................................... 1

1.2

Objectives and significance of study ................................................................... 3

1.3

Structure of the thesis .......................................................................................... 4

CHAPTER 2
2.1

LITERATURE REVIEW .......................................................................... 6
Land use and land cover (LULC) changes .......................................................... 6

2.1.1

Rationale of land use and land cover change ............................................. 6

2.1.2

Remote sensing (RS) and Geographic Information Systems (GIS) ........... 7

2.1.3

Satellite image analysis .............................................................................. 9

2.1.4

Landsat image data .................................................................................. 11

2.2

Land use and land cover change prediction models .......................................... 14

2.3

Summary ............................................................................................................ 16

CHAPTER 3
3.1

METHODOLOGY ................................................................................. 18
Overview of the Stoney Creek Watershed (SCW)............................................. 18

3.1.1

Location and extent .................................................................................. 18

3.1.2

Stoney Creek ............................................................................................ 20

3.1.3

Sub-basins ................................................................................................ 22

3.2

Reconnaissance survey ...................................................................................... 23

3.3

Data selection and collection ............................................................................. 23

3.4

Identification of land use and land cover categories ......................................... 25

3.5

Landsat image processing and analysis ............................................................. 31
3.5.1

Image pre-processing ............................................................................... 31
ii

3.5.2

Image classification ................................................................................. 38

3.5.3

Image post-processing ............................................................................. 45

3.6

Land use and land cover change modeling ........................................................ 47
3.6.1

Change analysis ....................................................................................... 47

3.6.2

Transition potential modeling .................................................................. 48

3.6.3

Change prediction .................................................................................... 53

CHAPTER 4

RESULTS AND DISCUSSIONS ........................................................... 56

4.1

Landsat image results and discussions .............................................................. 56

4.2

Accuracy assessment results and discussions .................................................... 63

4.3

LULC change analysis ....................................................................................... 65

4.4

LULC analysis for sub-basins ........................................................................... 68
4.4.1

Stoney Creek (SC) sub-basin ................................................................... 68

4.4.2

Tachick Creek (TC) sub-basin ................................................................. 73

4.4.3

West Nulki Stream (WNS) sub-basin ...................................................... 76

4.4.4

South Nulki Stream (SNS) sub-basin ...................................................... 78

4.4.5

Corkscrew Creek (CC) sub-basin ............................................................ 80

4.4.6

Mid Cork Stream (MCS) sub-basin ......................................................... 82

4.4.7

West Cork Stream (WCS) sub-basin ........................................................ 84

4.5

Land change modeling results and discussions ................................................. 86
4.5.1

Results of transition sub-models .............................................................. 86

4.5.2

Results of Markov Chain model .............................................................. 88

4.5.3

Results of soft prediction and hard prediction ......................................... 91

CHAPTER 5

CONCLUSION ....................................................................................... 97

5.1

Summary of the study ........................................................................................ 97

5.2

Limitations and recommendations ..................................................................... 99

Bibliography .......................................................................................................................... 100

iii

List of Figures
Figure 1.1 Organization of the thesis ......................................................................................... 5
Figure 2.1. The basic phenomena and processes of land change science .................................. 7
Figure 2.2. Wavelength of sensors in different Landsat satellites ........................................... 14
Figure 3.1. Location of study area: Stoney Creek watershed .................................................. 19
Figure 3.2. Channel network within SCW .............................................................................. 21
Figure 3.3. Stoney Creek Watershed sub-basins...................................................................... 22
Figure 3.4. Agricultural Area class .......................................................................................... 26
Figure 3.5. Barren Land class .................................................................................................. 27
Figure 3.6. Coniferous Forest class ......................................................................................... 27
Figure 3.7. Cut Block class ...................................................................................................... 28
Figure 3.8. Deciduous Forest class .......................................................................................... 28
Figure 3.9. Mixed Forest class ................................................................................................. 29
Figure 3.10. Planted/Regrowth Forest class ............................................................................ 29
Figure 3.11. Urban Area class .................................................................................................. 30
Figure 3.12. Water class ........................................................................................................... 30
Figure 3.13. Dead Pine Trees class .......................................................................................... 31
Figure 3.14. Mosaiced 2016 Landsat image ............................................................................ 33
Figure 3.15. July 21, 1986 Landsat 5 TM image with 5-4-3 band combination ..................... 35
Figure 3.16. July 9, 1999 Landsat 5 TM image with 5-4-3 band combination ....................... 36
Figure 3.17. August 17, 2016 Landsat 8 OLI image with 6-5-4 band combination ................ 37
Figure 3.18. Segmentation results of 1986 pre-processed image ............................................ 39
Figure 3.19. Landsat image processing and analysis framework ............................................ 46
Figure 3.20. MLP neural network structure ............................................................................. 49
Figure 3.21. LULC change model framework......................................................................... 55
Figure 4.1. LULC classification map of SCW in 1986 ........................................................... 57
Figure 4.2. Area statistics for each LULC class in 1986 ......................................................... 58
Figure 4.3. LULC classification map of SCW in 1999 ........................................................... 59
Figure 4.4. Area statistics for each class in 1999 ..................................................................... 60
iv

Figure 4.5. LULC classification map of SCW in 2016 ........................................................... 61
Figure 4.6. Area statistics for each class in 2016 ..................................................................... 62
Figure 4.7. Land gain and loss analysis for period 1986-1999 (T1) ........................................ 65
Figure 4.8. Land gain and loss analysis for period 1999-2016 (T2) ........................................ 66
Figure 4.9. Net change analysis of LULC for period 1986-1999 (T1) .................................... 67
Figure 4.10. Net change analysis of LULC for period 1999-2016 (T2) .................................. 67
Figure 4.11. SC sub-basin LULC map for 1986 ...................................................................... 70
Figure 4.12. SC sub-basin LULC map for 1999 ...................................................................... 71
Figure 4.13. SC sub-basin LULC map for 2016 ...................................................................... 72
Figure 4.14. LULC change of SC sub-basin............................................................................ 73
Figure 4.15. LULC change of TC sub-basin ........................................................................... 74
Figure 4.16. TC sub-basin LULC map, (a) 1986, (b) 1999, (c) 2016 ...................................... 75
Figure 4.17. LULC change of WNS sub-basin........................................................................ 76
Figure 4.18. WNS sub-basin LULC map, (a) 1986, (b) 1999, (c) 2016 .................................. 77
Figure 4.19. LULC change of SNS sub-basin ......................................................................... 78
Figure 4.20. SNS sub-basin LULC map, (a) 1986, (b) 1999, (c) 2016 ................................... 79
Figure 4.21. LULC change of CC sub-basin ........................................................................... 80
Figure 4.22. CC sub-basin LULC map, (a) 1986, (b) 1999, (c) 2016 ..................................... 81
Figure 4.23. LULC change of MCS sub-basin ........................................................................ 82
Figure 4.24. MCS sub-basin LULC map, (a) 1986, (b) 1999, (c) 2016 .................................. 83
Figure 4.25. LULC change of WCS sub-basin ........................................................................ 84
Figure 4.26. WCS sub-basin LULC map, (a) 1986, (b) 1999, (c) 2016 .................................. 85
Figure 4.27. Transition potential map of CF to DT transition ................................................. 87
Figure 4.28. LULC hard prediction map for 2026................................................................... 91
Figure 4.29. Area comparison between 2016 and 2026........................................................... 93
Figure 4.30. Logical OR soft prediction map for 2026 ........................................................... 94
Figure 4.31. Maximum soft prediction map for 2026 ............................................................. 95

v

List of Tables
Table 2.1. Spectral features of Landsat 4-5 bands ................................................................... 12
Table 2.2. Usage of Bands for mapping .................................................................................. 13
Table 2.3. Spectral features of Landsat 8 bands ...................................................................... 13
Table 3.1. Description of the selected satellite images ............................................................ 25
Table 3.2. Introduction of driver variables for transition potential modeling ......................... 48
Table 3.3. General model information ..................................................................................... 52
Table 3.4. Forcing a single independent variable to be constant ............................................. 52
Table 3.5. Forcing all independent variables expect one to be constant .................................. 52
Table 3.6. Backwards stepwise constant forcing ..................................................................... 53
Table 4.1. Area (km2) statistics of LULC types for each study year ....................................... 62
Table 4.2. Error matrix results for 2016 image classification .................................................. 63
Table 4.3. Accuracy statistics results for 2016 image classification ........................................ 64
Table 4.4. LULC change analysis results for for period 1986-1999 (T1) and 1999-2016 (T2) 68
Table 4.5. Transition probability matrix .................................................................................. 90
Table 4.6. Transition areas matrix............................................................................................ 90
Table 4.7. Area statistics of 2016 and 2026 ............................................................................. 92

vi

Acknowledgement
First and foremost, my sincere gratitude goes to my supervisor, Dr. Jianbing Li, who has
unreservedly supported me at every stage of my study and research and guided me to achieve
the academic goal. I am also very grateful to my committee members Dr. Roger Wheate who
generously guided and inspired me throughout the research, and Dr. Jun Yin whose support
and company during the field work helped me a lot.
I would like to offer my gratitude to Scott Emmons, Senior GIS Lab Instructor, who is
also the instructor of my first class at UNBC. He helped and inspired me through many ways
by his kindness and profound knowledge in his field. My appreciation also goes to Ping Bai,
Lab Instructor and Anthony Jjumba, my course instructor who have provided great help to my
research.
At last, I would like to give my sincere thankfulness to my parents, who raised me and
gave me the opportunity to study abroad, and my wife, whose support is always on my back
and accompany always on my side. My thankfulness also goes to all the friends who shared
my laughs and cries throughout the time I spent outside my homeland.

vii

CHAPTER 1 INTRODUCTION
1.1 Background
Increasingly world landscapes have been influenced by anthropogenic activities, and
human development derives heavily from natural resources. Because of this, land use and land
cover (LULC) change resulting from anthropogenic activities has become an inevitable issue
and leads to environmental degradation risk around the globe (Foley, et al., 2005;
Paiboonvorachat, 2008; Stabile, 2012; Islam et al, 2018b). It was reported that half of the icefree terrestrial land in the world has been transferred to other land use and land cover types,
especially to agricultural lands (Haberl et al, 2007; Turner et al, 2007). The global LULC
changes have shown a tendency of transition from forests to agricultural lands and built-up
areas (Munoz-Villers and Lopez-Blandco, 2008; Paiboonvorachat, 2008). Deforestation and
the expansion of agricultural and built-up areas have caused various environmental issues, such
as loss of biodiversity, increased water withdrawals, decreased precipitation, restricted
hydrological cycle, degradation of soil nutrients, erosion of soils, increased greenhouse gas
emission and climate change (Turner, et al., 2007; Paiboonvorachat and Oyana, 2011; Paul,
2013).
Remote Sensing (RS) techniques have been widely utilized for LULC change detection
and dynamic observation. Geographic Information Systems (GIS) involve various functions of
spatial data management, creation, and analysis. More recently, these two sub-disciplines are
often applied in combination, providing a more effective and widely used technology for global
observing, monitoring and analyzing (Paul et al., 2012; Srivastava et al, 2013; Nguyen et al,
1

2016). By using remotely sensed data and historic maps, LULC dynamics can be mapped for
LULC change analysis, and then the information resulting from the analysis can be used for
LULC change models that are important tools for LULC management, environmental impact
assessment, and decision making etc. (Campbell, 2002; Turner et al., 2007; Paiboonvorachat,
2011; Paul, 2013). LULC models have been widely applied for studying the transition
probability among different LULC categories, a significant criterion for LULC change
prediction in modeling processes (Kamusoko, et al., 2009; Behera et al, 2012; Hua, 2017).
With abundant forestry resources and agricultural lands, Vanderhoof is located near the
geographical centre of British Columbia, Canada. It is a significant timber supply area in BC,
with several timber mills across the town (RDBN, 2014). Moreover, this area has been known
as one of the earliest agricultural settlements in BC. Due to its prolific forage production, the
Nechako River Valley and plateau have been named as the “hay basket” of the central interiornorthwest of BC. However, its capacity of agricultural production has not been fully realized
for planned and reserved agricultural lands (Agricultural Land Use Inventory, 2017). In the
area, Nechako River is the main drainage channel and is comprised of over 30 tributaries. This
study focuses on the Stoney Creek Watershed (SCW). Stoney Creek flows out of Nulki Lake
then passes through the Saik'uz First Nation Reservation, and north to Tachick Lake, before
finally, north-east towards Vanderhoof where it enters the Nechako River (Olin, 2013). It is one
of the important water resources in Vanderhoof District. However, limited studies or
information of LULC in SCW were available, posing difficulty and challenge for making
informed decision regarding its development. Thus, it is essential to examine the LULC change
and its dynamics within SCW.

2

1.2 Objectives and significance of study
This thesis study focuses on analyzing the LULC change data of SCW by using a
combined RS and GIS approach to have a better understanding of the watershed. The LULC
information obtained will provide accessible data for studies on hydrology, climate change and
food safety in SCW. Models combined with RS and GIS are used for LULC prediction, in order
to identify the scenarios of SCW in 2026. This process provides recommendations for future
environmental conservation and sustainable development planning. Finally, the methods
presented in this study can be applied to the whole Vanderhoof District in the future.

3

1.3 Structure of the thesis
This thesis is consisted of 5 chapters: Chapter1 (Introduction) introduces the background,
objectives and significance of the study. The review of existing literatures is represented in
Chapter 2 (Literature Review) together with the context explanation and method evaluation.
The detailed procedures of selected method are explained in Chapter 3 (Methodology)
including the overview of the study watershed, how satellite images were selected and used,
image processing and models for LULC change prediction. Chapter 4 (Results and Discussions)
presents the study results including processed LULC images, accuracy assessment results,
current LULC change data and LULC future modeling results and their analysis. Finally, a
summary of the study is described in Chapter 5, as well as the study limitations and suggestions
for future work. Figure 1.1 shows the organization of the thesis.

4

Chapter 1: Introduction
Introduces the
study

Study
background

Objectives and
significance

Chapter 2: Literature Review
Reviews existing literatures, explain context and choose
suitable method for the study

Chapter 3: Methodology
Overview on
SCW

Explains procedures
of selected method

Data selection
and application

LULC change
modeling

Image
processing

Chapter 4: Results and discussions
Processed
Landsat images

Presents and
analyzes results

Accuracy
assessment results

LULC change
modeling results

Current LULC
change results

Chapter 5: Conclusion
Summarizes the study, states research limitations and provide future
study recommendations

Figure 1.1 Organization of the thesis
5

CHAPTER 2

LITERATURE REVIEW

2.1 Land use and land cover (LULC) changes
2.1.1

Rationale of land use and land cover change

Land use and land cover are two different terms but are often used interchangeably by
many researchers. Land use refers to the ways how human beings utilize the land cover for
themselves and their habitat, which usually focus on the economic function of lands (e.g.
residential areas and agricultural zones) (Paul, 2013; Rawat and Kumar, 2015; Islam, 2018a).
On the other hand, land cover represents the biophysical characteristics such as vegetation, soil
and water distributed on earth’s surface. It also includes features created by anthropogenic
activities, such as settlements.
Land use and land cover have an interactive relationship. Figure 2.1 presents the basic
phenomena and processes of land change science (Turner et al., 2007). It exhibits the
relationship between the human and biophysical subsystems of land systems. It is uncertain
whether the land cover changes resulting from land use will lead to land degradation, but the
conversion of land use in the human subsystem driven by social activities can change land
cover. Moreover, these land cover changes impact natural environment and biosphere by
affecting climate, biodiversity, soil, water, radiation budgets, trace gas emissions etc. (Turner
et al., 2007; Paiboonvorachat, 2008; Rawat and Kumar, 2015; Liping et al, 2018).

6

Figure 2.1. The basic phenomena and processes of land change science (Turner et al., 2007)

All these impacts have potential to affect the long-term relationship between humans and
natural environment. Recently, the process of land use and land cover (LULC) change driven
by human activities and natural phenomena has accelerated and become a widespread
challenge. It is thus crucial to study LULC change detection for a better understanding of
landscape dynamics and future scenarios to make better land management and planning.
(Behera et al, 2012; Rawat and Kumar, 2015; Liping et al, 2018).
2.1.2

Remote sensing (RS) and Geographic Information Systems (GIS)

Remote sensing is the technique for collecting information of a surface without direct
contact. It makes use of the electromagnetic signature, reflected or diffracted by the sensed

7

objects, for the purpose of differentiating land use and land cover categories. By using this
technology, multi-spectral and multi-temporal data can be collected and converted into useful
information, for studying land patterns and dynamics (United Nations, 1986; Weng, 2001; Paul,
2013). Geographic information systems (GIS) describes geographical data collecting, storing,
processing, and presentation of geographically referenced information and spatial visions of
locations (USGS, 2007).
The combination of remote sensing and GIS has been applied and recognized widely as a
successful and effective technique in LULC change detection. Satellite remote sensing
technology is especially popular as it is supported by satellite sensors which can provide image
data with high spatial resolution and geometric precision, and by various accessible software
(Yeh and Li, 1997; Araya, 2009; Paul, et al., 2012; Stabile, 2012; Pervaiz et al, 2016).
Comparison of post-classification images and multi-date composite image change detection
are two most effective and popular ways to study change detection. Moreover, the information
of how much, where and what type of LULC change has occurred are considered as the main
results from remote sensing and GIS based change detection studies (Jensen, 1996; Weng,
2001).
For examples, Weng (2002) analyzed land use change in the Zhujiang delta of China by

combining satellite remote sensing, GIS and Markov chain modeling. Yang and Liu (2005)
applied remote sensing and GIS for coastal LULC change mapping for the Pensacola estuarine
drainage area in Mexico. By using topographic maps and multi-temporal remotely sensed data,
the LULC change in Dhaka Metropolitan area of Bangladesh from 1960 to 2005 has been
detected and monitored (Dewan and Yamaguchi, 2009). Multi-spectral satellite data was
8

applied for land use change mapping analysis in Simly watershed, Pakistan by Butt et al. (2015),
using a maximum likelihood supervised classification. Rawat and Kumar (2015) used a similar
method with Landsat image data to detect LULC change in the Hawalbagh block in India,
while emphasizing the importance of digital change detection technologies.
2.1.3

Satellite image analysis

Different land use and land cover types on the earth surface can emit and reflect a variety
of electromagnetic radiation, which can be captured by satellite sensors and then assigned
digital number (DN) values for each spectral band (Oumer, 2009; Paul, 2013). Based on the
DN values of different land use and land cover classes, all pixels in multi-spectral satellite
images can be categorized into respective LULC categories to perform image classification
(Lillesand, et al., 2008).
Various methods have been applied for satellite image classification. Yang and Liu (2005)
applied hierarchical classification, interactive classification and iterative image interpretation
for LULC change mapping. Rawat & Kumar (2015) employed supervised classification with
maximum likelihood, including Normalized Difference Vegetation Index (NDVI), Normalized
Difference Water Index (NDWI) and Normalized Difference Built Index (NDBI) channels to
support the classification. Chen et al (2006) applied a combination of NDVI, NDWI, NDBI
and NDBal (Normalized Difference Bareness Index) on the classification of urban LULC
change analysis. Most of the methods used pixel-based classification (PBC), which classifies
the images based on DN values of single pixels (Dean and Smith, 2003). It is a conventional
method and has been broadly applied as supervised and unsupervised classification in many

9

cases. However, because a pixel-by-pixel classification algorithm is applied to all available
image signals considered, the pixels with similar spectral reflectance would be grouped
together, regardless of their associated uses (Yang and Liu, 2005). Moreover, some spatial and
contextual information of image pixels are neglected, and the pixels do not represent true
geographical objects when using PBC method, thus affecting its accuracy in some cases
(Blaschke, 2010, MacLean, et al., 2013; Paul, 2013).
Compared with pixel-based classification (PBC), the object-oriented image classification
(OOIC) method classifies the satellite images by segments through creating a higher-level
image objects from the segmentation of low-level pixel-based features. While PBC classifies
by the characteristics of single pixels, OOIC generates segments and does the classification by
characteristics such as segment size, shape, texture, etc. (Benz, et al., 2003; Hay & Castilla,
2008; MacLean, et al., 2013). Compared with PBC, OOIC produces higher classification
accuracy, incorporating more information, and creating segments that are more recognizable
than pixels. It also decreases local spectral variations (Shackelford and Davis, 2003; Frohn,
2011; Chen et al, 2013).
Many previous studies used OOIC method because of its great potential for LULC image
analysis. For example, Zhou and Troy (2008) applied the OOIC approach to classify and
analyze the urban area in the Gwynns Falls watershed, Maryland, US at the parcel level. Frohn
(2011) used OOIC to map and analyze wetlands in Florida, emphasizing the advantages of this
classification method. A case study of LULC change in the Kiskatinaw River watershed in
Canada by Paul et al. (2018) also highlighted the high accuracy of the OOIC and its application
potential with land change modeling.
10

2.1.4

Landsat image data

Landsat satellites are Earth-observing satellites that have provided an archive of spacebased land remotely sensed data. The Landsat program has continuously offered the Earth’s
surface record to benefit people who work in agriculture, forestry, mapping, LULC change,
and global change research etc. since 1972 (USGS, 2012; USGS, 2016).
Landsat data has been widely applied for land use and land cover analysis (Weng, 2002;
Yang and Liu, 2005; Frohn, 2011; Zhao et al., 2012; Butt et al., 2015; Paul et al., 2018), and
the Landsat sensors have proved sensitive enough to categorise different spectral patterns
related to the LULC classes in many complex landscape conditions (Muñoz‐Villers and
López‐Blanco, 2008). In this study, Landsat 4 /5 and 8 were selected as image data sources
based on extended archive accessibility. image quality, availability and cost.
2.1.4.1 Introduction of Landsat 4-5 TM
Landsat 4 and 5 instrumented the Landsat Thematic Mapper (TM) sensor onboard from
July 1982 to May 2012. Landsat 4-5 TM collected image data with a 16-day repeat cycle
(mostly Landsat 5), and data files include seven spectral bands (Table 2.1). Bands 1-5 and 7
were collected at 30 meters resolution, while band 6 (Thermal infrared) was 120 meters then
resampled to 30 meters (USGS Landsat Archives, 2018).

11

Table 2.1. Spectral features of Landsat 4-5 bands (USGS Land Mission, 2015)

Wavelength Resolution
(micrometers) (meters)

Bands

Spectral Respond

1

Blue

0.45-0.52

30

2

Green

0.52-0.60

30

3

Red

0.63-0.69

30

4

Near Infrared (NIR)

0.76-0.90

30

5

Shortwave Infrared (SWIR) 1

1.55-1.75

30

6

Thermal

10.40-12.50

120* (30)

7

Shortwave Infrared (SWIR) 2

2.08-2.35

30

2.1.4.2 Introduction of Landsat 8 OLI and TIRS
Landsat 8 satellite carries the Operational Land Imager (OLI) and Thermal Infrared
Sensor (TIRS) launched in 2013. As with Landsat 4 / 5, Landsat 8 collects image data from the
earth with a 16-day repeat cycle, but data files consist of 11 spectral bands (Table 2.3). Figure
2.2 shows wavelength for different Landsat sensors. OLI Bands 2-7 have similar usage as
Bands 1-7 of Landsat 4-5 (Table 2.2). Moreover, the OLI sensor of Landsat 8 provides
enhancements from prior instruments, with three new spectral bands: Band 1, an ultra blue
visible channel, specifically useful for coastal and water resources studies; Band 6, a
Panchromatic channel with 15-meter resolution, useful for image sharpening; Band 9 designed
for detecting cirrus clouds contamination as a new infrared channel. The TIRS sensor was
designed to provide two thermal channels for thermal mapping and estimating soil moisture.
The data were captured with at least 100-meter resolution, then registered to OLI data with the
resolution of 30 meters (Barsi, et al., 2014; USGS Landsat Archives, 2018).
12

Table 2.2. Usage of Bands for mapping (Barsi, et al., 2014)

Bands

Usage for mapping

1

Discriminate soil from vegetation, and deciduous forest from
coniferous forest; Bathymetric mapping
Useful for emphasizing peak vegetation
Distinguish vegetation slopes
Emphasize biomass content and shorelines
Distinguish moisture content of vegetation and soil; Useful for
penetrating thin clouds
Thermal mapping; Assessing moisture content of soil
Hydrothermally altered rocks associated with mineral deposits

2
3
4
5
6
7

Table 2.3. Spectral features of Landsat 8 bands (Barsi, et al., 2014)

Bands

Spectral Respond

1
2
3
4
5
6
7
8
9
10
11

Ultra Blue (coastal/aerosol)
Blue
Green
Red
Near Infrared (NIR)
Shortwave Infrared (SWIR) 1
Shortwave Infrared (SWIR) 2
Panchromatic
Cirrus
Thermal Infrared (TIRS) 1
Thermal Infrared (TIRS) 2

Wavelength
(micrometers)
0.435 - 0.451
0.452 - 0.512
0.533 - 0.590
0.636 - 0.673
0.851 - 0.879
1.566 - 1.651
2.107 - 2.294
0.503 - 0.676
1.363 - 1.384
10.60 - 11.19
11.50 - 12.51

13

Resolution
(meters)
30
30
30
30
30
30
30
15
30
100 * (30)
100 * (30)

Figure 2.2. Wavelength of sensors in different Landsat satellites (USGS Land Mission, 2015)

2.2 Land use and land cover change prediction models
Land use and land cover models are supportive tools for analyzing the causes and
consequences of LULC changes (Verburg et al., 2004). They can be used in various fields such
as earth-atmosphere interactions, biodiversity loss, vegetation cover, forest management, urban
expansion and global environmental change (Weng, 2001; Echeverria, et al., 2008; Stabile,
2012; Kamusoko, et al., 2013; Paul, 2013). Modeling LULC change has six important features:
level of analysis, driving forces, cross-scale dynamics, temporal dynamics, spatial interaction
and neighbourhood effects, and level of integration (Verburg, et al., 2004). However, only a
few models are available to solve the causes and consequences of LULC changes, which can
be grouped into whole landscape models, distributional landscape and spatial landscape models
(Bark, 1989; Weng, 2002; Halmy et al, 2015). Unlike distributional models without spatial
details, spatial models can be manipulated quickly by using remote sensing data with
geographical referenced information (Bark, 1989; Urban and Wallin, 2017).

14

Many spatially explicit models such as Markov Chain (MC) model, Artificial Neural
Network (ANN) model, Cellular Automata (CA) model, etc. Mishra and Rai (2016) have been
incorporated with RS and GIS technologies for more effective simulation and prediction for
LULC changes, Specifically, the CA model is one of the most commonly used models in LULC
research, especially in urban development studies, in order to model complex land use change
(Verburg, et al., 2004). MC model is a stochastic model which has been adopted in many LULC
change analysis studies and researches around the world. It can examine the stochastic nature
of the LULC change data and forecast the stability of future land development (Weng, 2002).
ANN is a non-parametric approach which imitates the human biological nervous system. It can
quantify and forecast complex non-linear patterns (Pijanowski et al., 2002; Mishra and Rai,
2016). Multi-layer perceptron (MLP) is a feedforward ANN. Typically MLP consists of one
input layer, one output layer, and one or more hidden layers. It can automatically generate a
large amount of parameter values but requires less data for training (Paul, 2013). MLP has the
advantage of modeling several or all transitions at the same time (Eastman 2012).
In a specific study area, each individual model has its limitation to detect all essential
LULC change processes. Therefore, some studies combined two or even multiple models to
reduce the uncertainty of the prediction (Verburg et al., 2008; Mishra and Rai, 2016). For single
model studies, the CA model was applied by Hasbani (2008) to predict land use development
in the Calgary region, Canada. The ANN method was applied with GIS techniques by
Pijanowski et al. (2002) to analyze how influential factors affect urbanization patterns and
predict land use changes in Michigan’s Grand Traverse Bay Watershed. The MC model was
applied for land use change analysis in Zhujiang Delta, China by Weng (2002). For studies

15

with combined models, Qiang and Lam (2015) applied the combination of ANN and CA
approaches to model LULC changes in the Lower Mississippi River Basin region to support
land planning and management of the region. Freier et al. (2011) applied Markov process aided
by CA model to model rangelands considering the effects of climate change. A flexible hybrid
approach consisting of CA and Markov Chain models was used to model LULC transition
potentials, and then predict future LULC scenarios in the lower Hunter Valley, NSW (Stabile,
2012). A multi-scale and multi-model approach was developed by Verburg et al., (2008) to
analyze future LULC changes in Europe. This is an example of solving for a multitude of
changes which result from interacting processes with different operating scales and domains at
the same time. Moreover, Arekhi and Jafarzadeh (2014) predicted areas vulnerable to forest
conversion in the Zagros forest region using MLP neural network with MC model. This MLPMC modeling method was also used for land use and land cover prediction in Patna district,
India (Mishra and Rai, 2016). Similarly, modeling method with MLP neural network and MC
model was utilized by Paul et al. (2018) to predicted LULC scenarios within the study
watershed.

2.3 Summary
Because more and more anthropogenic and natural factors accelerate the LULC change
process, it is necessary to study LULC changes to have a better understanding of land patterns
and dynamics of an area. Data from Landsat satellites has been continuously provided to LULC
change researchers over the world, because of its image quality and extended archive
accessibility. The combined RS and GIS approach has been widely used in LULC change
detection to provide high quality geospatial data and flexible data processing by using various
16

software. In recent years, there is a trend to use the OOIC method for LULC image
classification which can produce higher accuracy with more incorporated information than
conventional PBC. The integrated MLP neural network and MC analysis method is a recently
developed method for quantifying and modeling spatiotemporal changes of LULC. This hybrid
approach has the advantages of both MLP and MC models as MLP has excellent ability to
model each LULC transition in the study area, while MC model is widely used for its capability
to predict future LULC changes based on transition probability matrixes over a certain period
of time. Therefore, the method using the combined RS and GIS approach to process Landsat
image data through OOIC method, and with the MLP neural network and MC hybrid modeling
is sufficient to perform geospatial analysis of LULC change and LULC change prediction.

17

CHAPTER 3 METHODOLOGY
3.1 Overview of the Stoney Creek Watershed (SCW)
3.1.1

Location and extent

The Stoney Creek Watershed is located to the south-west of Vanderhoof town, near the
geographical centre of British Columbia, Canada (RDBN, 2014) (Figure 3.1). SCW lies
between longitude 123°58’51” W to 124°25’12” W, and latitude 53°42’23” N to 54°01’32” N.
It covers a total drainage area of 497.02 km2. Tachick Lake, the biggest lake within the
watershed, lies in the north part of the SCW. Nulki Lake lies south of Tachick Lake. The Saik'uz
First Nation Reservation lies between these two lakes. Highway 16 passes through the northeast part of the study watershed. The Nechako River flows north above the SCW then east to
the city of Prince George (Figure 3.1).

18

Figure 3.1. Location of study area: Stoney Creek watershed
19

3.1.2

Stoney Creek

Stoney Creek has a total length of 20.85 km. It has a series of natural falls in the Nechako
River valley, and below these falls it flows through many different land use and land cover
types such as agricultural areas, forest areas and urban areas (Olin, 2013).
Stoney Creek lies on the south bank of the Nechako River, upstream of Riverside Park in
the Vanderhoof community, and is an important traditional area for Saik'uz First Nation. The
lower main channel of the Stoney Creek flows 5 km within the municipal boundary of the
District of Vanderhoof. The riparian zone varies widely along most of its length within the
municipal boundary (Olin, 2013). Figure 3.2 shows the main channel of Stoney Creek and its
tributaries.

20

Figure 3.2. Channel network within SCW

21

3.1.3

Sub-basins

In this thesis, the SCW area was divided into seven sub-basins including Stoney Creek,
Tachick Creek, West Nulki Stream, South Nulki Stream, Corkscrew Creek, Mid Cork Stream
and West Cork Stream (Figure 3.3) using the original watershed boundaries. The Stoney Creek
sub-basin covers the largest area of 261.41 km2 among all sub-basins, which includes the main
stream of Stoney Creek and two lakes (Tachick Lake and Nulki Lake). The other sub-basins
are Tachick Creek (36.7 km2), West Nulki Stream (39.89 km2), South Nulki Stream (19.94
km2), Corkscrew Creek (71 km2), Mid Cork Stream (30.97 km2) and West Cork Stream (37.11
km2).

Figure 3.3. Stoney Creek Watershed sub-basins
22

3.2 Reconnaissance survey
In this study, a reconnaissance field survey aimed to collect the overall information of the
land use and land cover in Vanderhoof, BC, Canada, such as historical, geographical,
humanistic, biological and economic information. Several periods of field survey have been
conducted for ground truth data collection and photos from 2017 to 2018. The obtained
information was used for image classification accuracy assessment and identifying the types
of LULC. Several sites have been chosen based on the locations where satellite images are not
easy to distinguish. A deeper insight about the SCW study area was obtained from the field
survey.

3.3 Data selection and collection
Data selection and collection should be considered based on the watershed boundary,
purpose of the study and availability of satellite images. Based on the literature review and my
prior project experiences, several rules and reasons for Landsat satellite images selection and
collection were applied:
a. The objective of the study is to capture land use and land cover changes within Stoney
Creek Watershed from the early 1980s to 2016 based on the accessibility of Landsat data.
There may be land transitions before and after 1973 as the lumber mill came to production
during the time. Because Landsat data were available after 1972, this study collects
Landsat data from 1973 to 2016;
b. Landsat 4/5 and 8 were selected for the appropriate period and time span, and with
23

minimal or no cloud cover;
c. It is important to collect images captured at about the same time of the year. Plants
experiencing seasonal changes have different levels of ability to absorb visible light and
reflect near-infrared energy. Therefore, vegetation categories (including agricultural area,
coniferous forest, deciduous forest and mixed forest) show variable appearances
seasonally. This can impact the accuracy of classification and statistic analysis. The
period when vegetation is in good health and maturity was selected to benefit the
classification.
Based on these considerations, the Landsat satellite images in the following three years
were chosen for analysis:
1) 1986 Landsat 5 TM imagery: this represents the earliest year with available spectral image
data of 30-meter resolution. The lumber mill came into production in 1973, and 1986 is
the closest date to represent the original data before LULC affected by forestry industries.
2) 1999 Landsat 5 TM imagery: the mountain pine beetle infestation occurred in the central
interior region in the late 1990s. Since 1999, significant actions have been made by the
BC government such as province-wide aerial overview surveys (Province of BC, 2016;
Ebata, 2003). Mountain pine beetle infestation has had a significant impact on the
coniferous forest and mixed forest categories.
3) 2016 Landsat 8 OLI and TIRS imagery: this can represent current status within SCW.
Landsat satellite imageries were collected from Earthexplorer and Remotepixel websites.
24

By using dates from late July to early August can effectively keep the analysis free from
seasonal variability impact. The Stoney Creek Watershed boundary I used for my research was
downloaded from iMapBC website. The whole boundary of the study area was merged with
the main Stoney Creek basin and the sub-basins around it. Those sub-basins were chosen based
on the stream sources.

Table 3.1. Description of the selected satellite images

Year

Date

Landsat

Bands Resolution(meters)

1986

July 21th

5 TM Level 1

1–5

30

1999

July 9th

5 TM Level 1

1–5

30

2016 August 17th 8 OLI & TIRS

2–6

30

For Landsat 5 image data, bands 6 and 7 were excluded. For Landsat 8 data, bands 2-6
were selected for this project.

3.4 Identification of land use and land cover categories
Information and characteristics of different land use and land cover categories were
identified based on land survey data and information gained from locals. In total 10 LULC
classes were created and named for the subsequent satellite image analysis: Agricultural Area,
Barren Land, Coniferous Forest, Cut Block, Deciduous

25

Forest, Mixed Forest,

Planted/Regrowth Forest, Urban Area, Water and Dead Pine Trees. A modified version of the
Anderson scheme has been applied as a reference for describing and identifying land use and
land cover classes. Both Level I and Level II categories have been used at the same time
(Anderson, 1976). Because of the image resolution and areas of the study area, these classes
were chosen based on their areas and the importance for the watershed.

Figure 3.4. Agricultural Area class

Agricultural Area: This class is mainly for food and fiber production. It mostly consists
of flat areas with more regular shape such as rectangles. Figure 3.4 shows cropland, which is
part of Agricultural Area class, within the study area.

26

Figure 3.5. Barren Land class

Barren Land: This class represents areas with less than one-third area vegetated or any
other type of cover, and its skill to support life is limited (Anderson, 1976). It includes some
structures such as roads, open pits etc. Figure 3.5 shows an abandoned quarry pit within the
study area.

Figure 3.6. Coniferous Forest class

Coniferous Forest (Figure 3.6): This class represents the forest areas covered mainly by
evergreen trees.
27

Figure 3.7. Cut Block class

Cut Block: This class represents the areas with forest removed after wild fire or pine beetle
infestation, or tree cutting for wood supply and industrial purposes. Figure 3.7 shows a cut
block area with grass and dead trees resulting from wild fire.

Figure 3.8. Deciduous Forest class

Deciduous Forest (Figure 3.8): This class represents the forest areas covered
predominately by trees which lose or gain leaves seasonally, and especially they lose leaves at
the end of the forest-free season or the beginning of a dry season (Anderson, 1976).
28

Figure 3.9. Mixed Forest class

Mixed Forest (Figure 3.9): This class consists of forest lands with both coniferous and
deciduous forest and neither predominates.

Figure 3.10. Planted/Regrowth Forest class

Planted/Regrowth Forest: This class represents the areas covered by young coniferous
plants, deciduous plants or mixed young plants which have naturally regrown or been planted
by human on cut block or barren land. Figure 3.10 is a typical Planted Forest area grown with
circa 4-year-old pine trees.
29

Figure 3.11. Urban Area class

Urban Area: This class represents the areas have much of the land covered by structures.
For example, cities, towns, villages, strip development along highways, transportations, etc.
(Anderson, 1976). Figure 3.11 shows a lumber mill located within Vanderhoof city area next
to highway 16.

Figure 3.12. Water class

Water: This class consists the areas mainly covered by a water body such as lakes, rivers
and creeks. Figure 3.12 is part of the main channel of Stoney Creek.
30

Figure 3.13. Dead Pine Trees class

Dead pine trees: This class represents the areas with pine trees killed due to the infestation
of mountain pine beetles. Those pine trees only show great changes since the second stage of
infestation, so only the satellite image of 2016 includes this class. Figure 3.13 shows the area
with most pine trees killed by mountain pine beetles.

3.5 Landsat image processing and analysis
Before land use and land cover modeling, processing and analysis for three satellite
images were applied by using two main software packages: PCI Geomatica 2017 and Esri
ArcGIS 10.5.1. PCI Geomatica 2017 was mostly used to conduct image pre-processing,
classification, accuracy assessment etc. ArcGIS 10.5.1 was used mainly for image postprocessing, such as format converting, mapping, statistic analysis etc.
3.5.1

Image pre-processing

On the USGS Earthexplorer website, before searching proper Landsat images, some
parameters need to be set up. In this thesis, the date of all images was selected for July - August.

31

Comparison for data quality was applied to result in ideal Landsat images for later image
classification. Landsat image data files downloaded from USGS datasets contain all bands in
TIFF format files, and an MLT text file. These bands were opened in Focus, the main platform
on PCI, through the MLT file and using Multi Spectral (MS) input option.
While choosing Landsat images, one single image in 2016 does not cover the entire study
area. The Mosaic tool in Geomatica software therefore was utilized to combine all bands of the
two adjacent images into one respectively. Figure 16 is a mosaiced Landsat image for 2016
showing false color with 6-5-4 combination bands and Adaptive enhancement. For the purpose
of generating a more accurate classification result, a larger rectangle area was used, but clipped
by the SCW boundary for later analysis. The yellow boundary is the SCW study area
surrounded by the rectangle study area (red boundary) (Figure 3.14).

32

Path 49 & Row 22

Path 49 & Row 23

Figure 3.14. Mosaiced 2016 Landsat image

Using the Clipping tool on PCI, the original Landsat images were clipped by the rectangle
boundary (shapefile). The SWIR-NIR-Red ‘false color’ composite can provide the user with a
great amount of information and color contrast, which is the 5-4-3 band combination for
Landsat 4-5 images (1986 and 1999) and 6-5-4 band combination for Landsat 8 image (2016).
[TM Band 5 has the most agricultural information.] This combination is not only widely used
33

for vegetation studies, but also useful for timber management, pest infestation etc. (James,
2001). Heathy vegetation shows in a green color and soil is mauve in this combination.
Waterbodies appear in blue color with varying depth and turbidity of water. The Adaptive
enhancement, which shows the best enhancement result, was applied for all three Landsat
images. Figure 3.15 and Figure 3.16 are the Landsat images for 1986 and 1999 with 5-4-3 band
combination and clipped to the SCW extent. Figure 3.17 shows the 2016 image with 6-5-4
band combination clipped by the SCW boundary.

34

Figure 3.15. July 21, 1986 Landsat 5 TM image with 5-4-3 band combination

35

Figure 3.16. July 9, 1999 Landsat 5 TM image with 5-4-3 band combination

36

Figure 3.17. August 17, 2016 Landsat 8 OLI image with 6-5-4 band combination

37

3.5.2

Image classification

In remote sensing, digital image classification, also known as spectral pattern recognition,
detects the spectral information for each pixel in a satellite image to cluster pixels into the same
spectral categories. Different classes contain a mosaic of pixels to composite the thematic maps
which is the classified image (Geomatica, 2017; Paul, 2013).
In this study, Object Analyst on PCI 2017 has been utilized to generate the object-oriented
image classification. It is called the object-based image analysis (OBIA) method in PCI. There
are six steps for performing OBIA by using the pre-processed images from the last step:
Segmentation, Feature extraction & calculating features, Unsupervised classification, Rulebased classification and manual editing.
3.5.2.1 Segmentation
A successful classification starts with good segmentation, which means in general most
of the segments should correspond to a single land use and land cover class. Also, the shapes
of objects should align with the boundaries of different classes in the image (Geomatica, 2017).
Moreover, in Object Analyst, combinations of different parameters including scale (SC), shape
(SP) and compactness (CP) can greatly influence the results of segmentation. Figure 3.18
shows segmentation results on the same area of the 1986 image, with various combinations of
these three parameters.

38

Figure 3.18. Segmentation results of 1986 pre-processed image

Multiple tests show that the combination (SC: 25, SP: 0.25, CP: 0.5) has the most
satisfying result for the image, as it segments most of the objects needed for later classification
and analysis. Other combinations that provide few segments or massive segments may result

39

in lower classification accuracy. The compactness was constant because it was barely changing
the segmentation results, and 0.5 has been found to be the best value for this case.
After running the Segmentation process, three vector layers were generated for each year
in .pix file format. All later processes were applied on these vector layers.
3.5.2.2 Feature extraction & calculating features
Performing feature extraction provides attributes of all segments, which are the reference
data for image object classification. In Object Analyst, two types of features were computed.
Statistical features were computed by the pixels within an object. These features include
minimum, maximum, mean and standard deviation for each of the selected image bands.
Geometrical features were computed using the geometric characteristics of objects based on
the segment boundary created during segmentation. This is the advantage of object-based
image analysis compared with the conventional pixel-based analysis. Compactness, elongation,
circularity and rectangularity were computed for geometrical features (Geomatica, 2017). All
these features were added to the attribute table of each vector layers.
The features above were used as basic bands attributes. Moreover, some extra feature
calculations were made for a better classification result by using the statistical features. In the
attribute manager of each segmentation layers, several new fields were added and new feature
layers below were computed:
1) NDVI (Normalized difference vegetation index), determines the density of green on an
area, which helps to classify live vegetation or not. It is calculated by the formula:
NDVI = (NIR — Red) / (NIR + Red) (NASA, 2000).
40

2) TM 2/3 Ratio, sharply differentiates croplands and barren lands, but not for separating
croplands, forests and water bodies;
3) TM 3/2 Ratio, for separating forests and croplands, and it is useful to distinguish broad
classes of vegetation;
4) TM 3/4 Ratio, classifies barren lands and urban area sharply, but not for defining
waterbodies, forests and croplands;
The last three ratios were applied based on the evaluation of Quinn (2001).
These features were computed for minimum, maximum, mean and standard deviation
values respectively, e.g. NDVI_Min, NDVI_Max, NDVI_Mean and NDVI_StdV.
3.5.2.3 Unsupervised classification
Before running unsupervised classification, fields of minimum, maximum, mean and
standard deviation features for band 3, 4 and 5 (4, 5 and 6 for image of year 2016) were selected.
The geometrical features and new calculated features were also included for the classification.
There were nine LULC classes for images of years 1986 and 1999, and ten LULC classes for
year 2016. Clusters were generated with twice as many as the selected LULC classes, so 20
clusters were set for the classifier. These 20 clusters could provide all the LULC classes with
better unsupervised classification results than any other numbers of cluster. Within the 20
clusters, several clusters could be classified as the selected LULC classes, because most of the
segments they covered belong to specific classes, while the rest covered segments for different
LULC classes averagely.
41

3.5.2.4 Rule-based classification
By using this rule-based classification, custom rules were created to assign class
membership to segments. Membership of a sample in a class could be determined by the criteria
selected based on the understanding of domain, data, or both (Geomatica, 2017). To edit and
improve classification results, custom classification rules were performed on classified or
unclassified segments by assigning or removing them. These two steps could simply split a
class into subclasses or combine two classes into a new one. An extra step was used to support
classification by selecting specific range of a feature. The discrimination between specific
LULC classes were performed by choosing effective features as range field. Based on some
experiment and ground truth experiences, a good range was picked by adjusting the minimum
and maximum values. For example, the cluster with both segments of barren lands and
segments of unvegetated croplands could be split by using ‘Mean_2/3 Ratio’ field as the range
field, as this feature is unique for separating barren lands and croplands. ‘Mean_3/4 Ratio’ can
improve the classification between barren lands and urban area for the cluster with both
segments of these two classes. Or by using the ‘Rectangle’ feature because of the artificial
shape of most cropland. Then adjusted both the minimum and maximum value till most of the
required segments were included, and the rule-based classification result could be satisfied for
later analysis.
3.5.2.5 Accuracy assessment
Accuracy assessment is an important step before generating a final classified image.
Usually it could be performed by using a reference map which has been field-checked or

42

contains ground truth data collected during the field survey (Weng, 2002). In this study,
because of the lack of reference map for all three years and the ground truth data was only
collected from 2017 to 2018, a training field was created for image 2016 through the tool for
editing training sites on the Object Analysis. This training field was then used in the accuracy
assessment section on the Object Analyst. By comparing unsupervised classified image and
training field created based on field survey, three accuracy reports were generated. As they
were almost the same procedures, features and rules were applied to classify the images for all
three years, the accuracy assessment of 2016 could represent the accuracy of 1986 and 1999.
To avoid the subjectivity of picking manual points, a random points layer was created
from ArcGIS software by using the points creating tool which could generate random point
layers for accuracy assessment using stratified random sampling method. In this study, the
‘Area’ field was selected as the target field, and the number of random points was set as 200
because of the extent of the study area. Therefore, those random points were created to
distribute proportionally in relative classes, which means the classes with bigger areas such as
‘Coniferous Forest’ and ‘Deciduous Forest’ covered more points than classes with smaller
areas. This method could significantly improve the objectivity of assigning segments while
creating the training field.
The random points shapefile layer was transferred from ArcGIS to Geomatica and
overlaid on the pre-processed satellite image of 2016. Segments which covered random points
were distinguished based on ground truth data and assigned to relative classes. Finally, three
accuracy assessment reports were generated:

43

a. Accuracy Statistics Report: this report shows overall accuracy, 95% confidence
intervals, kappa coefficient and producer’s and user’s accuracy for all classes.
b. Error (Confusion) Matrix Report: this report shows the matrix reference data and
classified data of all classes.
c. Sampling List Report: this report displays all assessed samples, described with
georeferenced position, classified value/name, and reference value/name.
In the accuracy assessment report, Kappa Khat method could test the data in an error
matrix better than the random, as the K coefficient estimates the agreement between a modeled
LULC contribution and ground truth (Lillesand, et al., 2008). For an error matrix with n number
of LULC categories, the equation to compute the Kappa coefficient is:
𝑝 −𝑝

0
𝑒
𝐾𝑎𝑝𝑝𝑎 = 1−𝑝
𝑒

(1)

where,
p0 = Observed agreement,
pe = Expected agreement.
Overall accuracy is the ratio between the total number of correct classifications and the
total number of classifications. User’s accuracy defines the ratio between the number of correct
classifications in a category and the total number of classifications in that category, and
Producer’s accuracy is the ratio between the number correctly identified in reference plots of
a given category and the actual number in that reference category.

44

3.5.3

Image post-processing

In this study, all the images and layers were projected to UTM zone 10 with North
America Datum (NAD) 1983. Classified images, modeling results and some maps for study
area introduction were generated as final maps using ArcGIS. Those maps were exported as
jpg image files with suitable resolution, which were helpful for result display and analysis.
Many images and tables were created for later results display and discussion, such as areas
of each class for both whole study area extent and sub-basins extent, etc.. The clipping tool in
ArcGIS software was applied to generate classified images with the extent of each sub-basin.
By working on attribute tables, new fields like ‘Areas’ could be computed. All these images
and tables have made the results discussion more visual and persuasive.
In summary, Figure 3.19 shows the flow chart of all the steps for Landsat image
processing and analysis:

45

Image pre-processing

Download Landsat images
Mosaic and clip images to
study area extent
Bands composite (5-4-3 for
1986 and 1999; 6-5-4 for 2016)

Image segmentation

Image classification
(OBIA)

Feature extraction
Unsupervised classification
Rule-based classification

Training sites editing
Accuracy assessment

Image post-processing

Manual editing and checking
Generate final maps for 1986,
1999 and 2016

Convert image format
Statistic analysis

Figure 3.19. Landsat image processing and analysis framework
46

3.6 Land use and land cover change modeling
To model LULC change, an integrated application on IDRISI Selva software (Land
Change Modeler (LCM)) is used. The LCM is used for historical LULC change analysis,
transitions potential modeling, LULC change scenarios prediction, assessing the implications
for biodiversity from specific change, and generating the evaluation for possible interventions
for maintaining ecological sustainability. It includes a set of tools for the changing assessment,
generating evaluations of gains and losses between two images, net change, persistence and
specific transitions with one-click in both maps and forms (Eastman, 2012).
3.6.1

Change analysis

The analysis for historical LULC data was performed, and graphs and maps of LULC
change were generated. The change analysis would be applied for three sets of time periods:
T1, T2 and T3. T1 represents the period from 1986 to 1999, and T2 from 1999 to 2016. These
two ranges were used for assessing the changes between different years, and creating graphs
and maps of LULC gains and losses, net change, etc.. T3 represents the time range between the
earliest (1986) and the latest time (2016). The transition potential modeling and change
prediction would be performed by using T3. During T3, the transitions with less than 5 km2
have been ignored. Some of these small changed areas were generated from irrational
transitions, e.g. transition from water to agriculture area, while some may result from
classification errors. Both types of transitions would cause negative impacts on modeling, so
those transitions were filtered and manually picked for later use. In addition, LCM is equipped
with the function of special modification to the output, so that it could mask out the transitions
47

that do not match any specific transition cases (Eastman, 2012).
3.6.2

Transition potential modeling

In the first section of the transition potential modeling tab, selected transitions from
change analysis were listed. All the transitions would be modeled one by one as sub-models.
To group these sub-models, a common name, ‘CHANGE8616’, was applied. Then in the third
part of this tab, all the driver variables were added into the model through the evaluating and
testing function. After testing it was found that some driver variable layers have different
parameters such as resolution, number of rows and columns. Therefore, these driver layers
were projected and output using image data of 1986 or 2016. In total, six driver variables were
added and listed in the table (table 3.2).
Table 3.2. Introduction of driver variables for transition potential modeling

Driver variable layers

Role

Hazard rating of mountain pine beetle This driver comprises extent and hazard rating of
infestation
pine beetle infestation.
This driver controls the general hydrology in this
Distance to major channel network
watershed.
This driver responsible for road development,
Distance to recommended permanent
built-up area development, forest clear cutting
roads
etc.
This driver determines the hydrological flow
Digital elevation model (DEM)
path i.e. overall hydrological process within a
watershed (Paul, 2013).
TWI defined as Ln(A/tanB) (Sörensen, et al.,
2006) where,
Topographic wetness index (TWI)
A = local upslope area draining through a certain
point per unit contour length, tanB = the local
slope.
Settlement reserved area and Agriculture These two can influence future built-up area and
development area
croplands development, forest harvesting etc.
48

Finally, the multi-layer perceptron (MLP) neural network technique was utilized to
perform the transition potential modeling. Non-linear neural networks can be conceived as a
complex mathematical function that is capable of converting input data to a required output.
Unlike other statistical techniques, the MLP has no pre-assumption to constrain data
contribution. It uses a back-propagation (BP) learning algorithm to model highly non-linear
functions. A typical MLP neural network consists of one input layer, one or more hidden layers
and one output layer, but generally one hidden layer is adequate for most studies. Figure 3.21
shows the MLP neural network structure (Eastman, 2012; Gardner, 1998; Lin et al., 2011)

Figure 3.20. MLP neural network structure (Eastman, 2012)

Nodes (or neurons) contained by each layer are connected to each other with unequal
connecting weights. The hidden layer nodes are significant for MLP to learn and utilize

49

interaction effects (Eastman, 2012; Paul, 2013).
For running transition sub-models, the training process used samples taken from pixels
that went through the transition being modeled, or pixels from persistence classes. 50% of the
samples were used for training and the remaining 50% were used for validation of the transition
potential model. MLP in LCM has been launched with an automatic training mode, which is
capable of monitoring and modifying the start and end learning rates by using the sample
training data (Eastman, 2012). All the parameters were assigned default values, except the
number of the hidden layer nodes. At the first time of running a sub-model, the default hidden
layer nodes value was used, then more running tests were performed. If the total accuracy and
skill score increase then the value will be doubled, otherwise, the last value is selected. Running
MLP provided the report with the aggregate accuracy and the skill measure score. The
following formula is the measure of model skill (Eastman, 2012):
S = (A – E(A)) / (1 – E(A))

(2)

where,
E(A) = excepted accuracy,
A = measured accuracy.
Moreover, the expected accuracy (E(A)) was expressed as:
E(A) = 1 / (T + P)
where,
T = the number of transitions in the sub-model,
P = the number of persistence classes in the sub-model.
50

(3)

Each sub-model was modeled with different driver variable groups. For the purpose of
getting higher modeling accuracy and skill score, it is necessary to remove the variables
without power. The backward stepwise analysis has been applied for model development
(Eastman, 2012). When running a sub-model with all the variables listed in table 3.3, each
variable was held constant one after another to figure out the one with the least effect on
modeling. If the model skill of holding a variable constant was similar with the skill of
modeling with all variables, the variable doesn’t affect the model much and it can be removed.
Using the transition from ‘Planted/Regrowth Forest’ to ‘Mixed Forest’ as an example, Tables
3.4, 3.5 and 3.6 are the sensitivity results of this transition potential sub-model by forcing
independent variables to be constant.

51

Table 3.3. General model information

Var. 1
DEM
Var. 2
Distance to major channel
Var. 3
Hazard rating of mountain pine beetle infestation
Var. 4
Distance to recommended permanent roads
Var. 5
TWI
Input layer neurons
5
Hidden layer neurons
3
Output layer neurons
2

Table 3.4. Forcing a single independent variable to be constant

Model

Skill measure

Influence order

With all variables

Accuracy
(%)
85.93

0.7185

N/A

Var. 1 constant

63.05

0.261

1 (most influential)

Var. 2 constant

81.13

0.6226

3

Var. 3 constant

82.52

0.6504

4

Var. 4 constant

78.27

0.5653

2

Var. 5 constant

85.89

0.7177

5 (least influential)

Table 3.5. Forcing all independent variables except one to be constant

Model

Accuracy (%)

Skill measure

With all variables

85.93

0.7185

All constant but var. 1

71.58

0.4316

All constant but var. 2

61.15

0.223

All constant but var. 3

49.91

-0.0018

All constant but var. 4

52.51

0.0503

All constant but var. 5

49.91

-0.0018

52

Table 3.6. Backwards stepwise constant forcing

Model

Variables included

Accuracy (%)

Skill measure

With all variables

All variables

85.93

0.7185

Step 1: var. [5] constant

[1,2,3,4]

85.89

0.7177

Step 2: var. [5,3] constant

[1,2,4]

82.36

0.6472

Step 3: var. [5,3,2] constant

[1,4]

78

0.5599

Step 4: var. [5,3,2,4] constant

[1]

71.58

0.4316

In this example, table 3.4 shows that variable 5 is the least influential. Then in table 3.6,
after running the first step of the backwards stepwise analysis, it was found that the skill
measure doesn’t change much. This means variable 5 which is the TWI layer is powerless for
the model, and it could be removed. This method was repeated till there are no such variables.
3.6.3

Change prediction

After generating all transition potential sub-models, they were applied in the Markov
chain model by using the tool for LULC change prediction on LCM. The future LULC
scenarios have been generated with the prediction date as 2026. In this procedure, the Markov
Chain determined the amount of LULC change using land cover images from the earlier date
to later date. It measured how much land is possible to transition from the later date to the
prediction date (Eastman, 2012). In the Markov Chain method, different categories were
applied as the states of a chain. The value at the previous time t+1, Xt+1, only depends on the
value at time t, Xt, and not on the process passing through Xt. It can be expressed as (Weng,
2002):

53

𝑋𝑡+1 = 𝑓(𝑋𝑡 )

(4)

If the transition probabilities are tabulated in a transition matrix, P, Xt-1 can be expressed
(Benito, et al., 2010):
𝑋𝑡+1 = 𝑋𝑡 ∙ 𝑃

(5)

The Markov Chain analysis can create a transition probabilities file, which describes the
probability of land cover change from one period to another by developing a transition
probability matrix between the earlier date and the later date. (Araya & Cabral, 2010; Houet &
Hubert-Moy, 2006).
Finally, two types of soft predictions were generated based on the present state of
transition potentials for each transition. The first one is maximum soft output, which shows a
maximum transition probability aggregation for all transition potentials. Then the ‘Logical OR’
type is calculated as (a+b – ab) for all transitions included, which can characterize a location
wanted by different transitions. The ‘a’ and ‘b’ represent the transition potential of each submodel. Along with the comprehensive soft output maps of transition probability, a hard
prediction LULC map was also created. Running the Markov Chain process provides the
quantity data. Then a multi-objective land allocation algorithm was applied to look through all
transitions and generate a list of classes both gain and lose land. While running this allocation
procedure, all the changed land of a class was allocated and then overlay to generate the output
(Eastman, 2012).
Figure 3.22 below, shows the flowchart of the LULC change model structure:

54

LCM Project
T1 (1986-1999)

LCM Project
T3 (1986-2016)

LCM Project
T2 (1999-2016)

Change

Driver
Variables

Analysis
Transition Sub-models

Net change
analysis for
T1 and T2

MLP
Transition Potentials

Figure 3.21. LULC change model framework

55

Hard prediction map
of T3

Maximum soft
prediction map of T3

Markov Chain

Logical OR soft
prediction map of T3

LULC
gain & loss
analysis for
T1 and T2

CHAPTER 4 RESULTS AND DISCUSSIONS
4.1 Landsat image results and discussions
The Stoney Creek watershed study area is 496.963 km2. In 1986, more than 75% of the
area was covered by different types of forest land, including 38.61% coniferous forest, 22.40%
mixed forest, 10.92% deciduous forest and 4.08% planted/regrowth forest. The other LULC
types include 4.85% agriculture area, 9.08% barren land, 2.00% cut block, 0.41% urban area
and 7.65% water (Figures 4.1 and 4.2).

56

Figure 4.1. LULC classification map of SCW in 1986

57

Water
Urban Area
Planted/Regrowth Forest
Mixed Forest
Deciduous Forest
Dead Pine Trees
Cut Block
Coniferous Forest
Barren Land
Agriculture Area
0

20

40

60

80

100 120 140 160 180 200 220
AREA (KM2)

Figure 4.2. Area statistics for each LULC class in 1986

Figures 4.3 and 4.4 present the image and graph results from the Landsat image analysis
of 1999. As with 1986, forest land comprised the most area of the study watershed, with
coniferous forest covering the maximum (30.83% of total area), followed by mixed forest
(22.73%), deciduous forest (14.18%). In 1999, the forest land including planted/regrowth
forest (7.10%) occupied almost 75% of the study area. Other LULC classes including
agriculture area, barren land, cut block, urban area and water covered 6.11%, 8.24%, 2.38%,
0.64% and 7.78% of the total watershed, respectively.

58

Figure 4.3. LULC classification map of SCW in 1999

59

Water
Urban Area
Planted/Regrowth Forest
Mixed Forest
Deciduous Forest
Dead Pine Trees
Cut Block
Coniferous Forest
Barren Land
Agriculture Area
0

20

40

60

80

100

120

140

160

180

AREA (KM2)

Figure 4.4. Area statistics for each class in 1999

Figures 4.5 and 4.6 exhibit image analysis output of 2016. In this period, forest land also
covered the majority of the watershed. However, a new LULC class of dead pine trees was
added in 2016 that covered 14.96% of the total area. This LULC class was mostly transferred
from coniferous forest and slightly from mixed forest. In 2016, the mixed forest occupied 21.43%
of the total area, followed by coniferous forest (18.97%), deciduous forest (13.54%), and
planted/regrowth forest (5.95%). Although it was affected by dead pine trees, forest land
comprised almost 60% of the whole study watershed. The remaining LULC features, including
agriculture area, barren land, cut block, urban area and water composed 9.62%, 6.45%, 0.83%,
0.57% and 7.69% of the SCW study area respectively. Table 4.1 shows the area statistic results
of LULC categories for all the three study years.

60

Figure 4.5. LULC classification map of SCW in 2016

61

Water
Urban Area
Planted/Regrowth Forest
Mixed Forest
Deciduous Forest
Dead Pine Trees
Cut Block
Coniferous Forest
Barren Land
Agriculture Area
0

20

40

60

80

100

120

AREA (KM2)

Figure 4.6. Area statistics for each class in 2016

Table 4.1. Area (km2) statistics of LULC types for each study year
1986
LULC Class

1999

2016

Area

Percentage

Area

Percentage

Area

Percentage

Agriculture Area

24.113

4.85%

30.352

6.11%

47.786

9.62%

Barren Land

45.108

9.08%

40.962

8.24%

32.072

6.45%

Coniferous Forest

191.898

38.61%

153.207

30.83%

94.282

18.97%

Cut Block

9.920

2.00%

11.852

2.38%

4.109

0.83%

Dead Pine Trees

0.000

0.00%

0.000

0.00%

74.339

14.96%

Deciduous Forest

54.290

10.92%

70.465

14.18%

67.306

13.54%

Mixed Forest

111.333

22.40%

112.984

22.73%

106.444

21.42%

Planted/Regrowth Forest

20.257

4.08%

35.288

7.10%

29.585

5.95%

Urban Area

2.037

0.41%

3.188

0.64%

2.835

0.57%

Water

38.006

7.65%

38.664

7.78%

38.205

7.69%

62

4.2 Accuracy assessment results and discussions
Errors of Landsat image classification can result from the lack of detailed information,
the application of medium-resolution images and other general influential factors (Paul, 2013).
For a more comprehensive discussion, results of accuracy statistics and error matrix were used
to explain and evaluate image classification output, as shown in Tables 4.2 and 4.3. The related
LULC types listed in the tables include Agriculture Area (AA), Barren Land (BL), Coniferous
Forest (CF), Cut Block (CB), Dead Pine Trees (DT), Deciduous Forest (DF), Mixed Forest
(MF), Planted/Regrowth Forest (PT), Urban Area (UA) and Water (WT).
Table 4.2 was generated by comparing the reference data (i.e. the training layer classified
by random sampling points) with the classified image data. Cells highlighted with grey shade
represent the segments that have been correctly classified, while the others, except blank cells,
represent the number of the mis-classified segments
Table 4.2. Error matrix results for 2016 image classification
Reference Data
Classified Data
WT
CF
MF
DF
PF
AA
BL
UA
CB
DT
Column Totals

WT CF MF DF PF AA BA UA CB DT Row Totals
10

47
4

31
2

1

10

51

1
35

1
17
1

19

1
11

12

63

2

1

18

18

2
15

18

5

5

8
8

1
21
24

10
49
34
23
11
21
16
5
9
22
200

By using the error matrix results, producer’s accuracy, user’s accuracy, overall accuracy,
95% confidence interval and overall Kappa value were calculated, as presented in Table 4.3.
Table 4.3. Accuracy statistics results for 2016 image classification
Overall Accuracy: 91.50%

95% Confidence Interval (87.39%

Overall Kappa Statistic: 0.90
Producer's
Class
95% Confidence
Accuracy
Name
Interval
（%）

95.62%)

User's
Accuracy
（%）

95% Confidence
Interval

Kappa
Statistic

AA

100

(97.22% 100.00%)

85.71

(68.37% 100.00%)

0.84

BL

83.33

(63.34% 100.00%)

93.75

(78.76% 100.00%)

0.93

CB

100

(93.75% 100.00%)

88.89

(62.80% 100.00%)

0.88

CF

92.16

(83.80% 100.00%)

95.92

(89.36% 100.00%)

0.95

DF

89.47

(73.04% 100.00%)

73.91

(53.79% 94.03%)

0.71

DT

87.5

(72.19% 100.00%)

95.46

(84.48% 100.00%)

0.95

MF

88.57

(76.60% 100.00%)

91.18

(80.17% 100.00%)

0.89

PF

91.67

(71.86% 100.00%)

100

(95.46% 100.00%)

1

UA

100

(90.00% 100.00%)

100

(90.00% 100.00%)

1

WT

100

(95.00% 100.00%)

100

(95.00% 100.00%)

1

Quantity Disagreement: 4.00%

Allocation Disagreement: 4.50%

In this study, the overall accuracy for image classification of 2016 is 91.50%, while the
overall Kappa statistic is 0.90. Most of the LULC classes have a larger producer’s accuracy
than overall accuracy except the classes of barren land, deciduous forest, dead pine trees and
mixed forest. The lower accuracy for these four classes can be attributed to the low Landsat
image resolution and the not very comprehensive ground truth data. However, the classification
with ≥ 90% overall accuracy is adequate for this study.

64

4.3 LULC change analysis
Based on the method stated in the previous chapter, the data of land gain and loss as well
as the net change for each LULC type have been generated through a change analysis procedure
in LCM on IDRISI.
Figures 4.7 and 4.8 represent the land gain and land loss of different LULC types for the
period of 1986-1999 (T1) and 1999-2016 (T2), respectively. Land gain results from the
transition from other classes to this class, and the opposite transition is land loss. It was
observed from Figures 4.7 and 4.8 that the mature forest land (including coniferous forest,
deciduous forest and mixed forest) occupied most of the study watershed.

Gains and losses between 1986-1999
60
40
20
0

AA

BL

CF

CB

DF

MF

PF

UA

WT

-20
-40
-60
Area Losses (km²)

Area Gains (km²)

Figure 4.7. Land gain and loss analysis for period 1986-1999 (T1)

65

Gains and losses between 1999-2016
100
80
60
40
20
0
-20
-40
-60
-80
-100
-120

AA

BL

CF

CB

DF

MF

Area Losses (km²)

PF

UA

WT

DT

Area Gains (km²)

Figure 4.8. Land gain and loss analysis for period 1999-2016 (T2)

From 1986 to 1999, coniferous forest experienced a great negative net change (i.e. greater
loss than gain) of 38.68 km2, while there is also a negative change for barren land. In contrast,
other LULC classes had more gain than loss, especially for deciduous forest and
planted/regrowth forest (Figure 4.9). The obvious changes of forest land may indicate the
forestry activities applied during this period. From 1999 to 2016, the class of dead pine trees
experienced a sharp increase in this area due to the infestation of mountain pine beetles causing
coniferous forest to lose 58.92 km2 of area. The other classes except agricultural land, also had
negative changes. An increasing trend for agricultural area since 1986 was observed, which
may be due to the changing policies and activities for agricultural development (Figure 4.10).
The slight net change of urban area indicates that there was no significant urban expansion and
industrial development. The town of Vanderhoof and the Saik'uz First Nation Reservation area
remain stable from 1986 to 2016 based on the results.

66

Net change between 1986-1999 (km²)
20
10
0

AA

BL

CF

CB

DF

MF

PF

UA

WT

-10
-20
-30
-40
-50

Figure 4.9. Net change analysis of LULC for period 1986-1999 (T1)

Net change between 1999-2016 (km²)
100
80
60
40
20
0
-20

AA

BL

CF

CB

DF

MF

PF

UA

WT

DT

-40
-60
-80

Figure 4.10. Net change analysis of LULC for period 1999-2016 (T2)
67

Table 4.4. LULC change analysis results for period 1986-1999 (T1) and 1999-2016 (T2)

LULC
Class
AA

Area Loss
(km2)
-14.62

T1
Area Gain
(km2)
20.86

BL

-33.31

29.15

-4.16

-30.80

21.92

-8.88

CF

-52.80

14.13

-38.68

-92.30

33.38

-58.92

CB

-9.85

11.78

1.93

-11.81

4.07

-7.74

DF

-28.21

44.40

16.19

-54.16

50.97

-3.19

MF

-48.69

50.32

1.63

-75.19

68.67

-6.52

PF

-3.01

18.05

15.03

-18.86

13.16

-5.70

UA

-0.19

1.34

1.15

-1.02

0.67

-0.35

WT

-0.82

1.48

0.66

-1.15

0.69

-0.46

DT

0.00

0.00

0.00

0.00

74.33

74.33

Net Change
(km2)
6.24

Area Loss
(km2)
-10.43

T2
Area Gain
(km2)
27.88

Net Change
(km2)
17.44

4.4 LULC analysis for sub-basins
Discussions on sub-basins level based on the LULC types within individual sub-basins
can provide more detailed information to have a better understanding of the LULC patterns
and dynamics of the area. It also may provide specific data for the research of certain sub-basin
in the future.
4.4.1

Stoney Creek (SC) sub-basin

SC sub-basin is located in the north part of SCW with the largest land area. It covers the
major water bodies of Stoney Creek, Tachick Lake and Nulki Lake. It also includes the Saik'uz
First Nation Reservation area and part of the town of Vanderhoof.
68

Figures 4.11, 4.12 and 4.13 show the LULC maps of the SC sub-basin area for 1986, 1999
and 2016. Based on the maps and graph (Figure 4.14), it is clear that the major LULC type in
this area was forest land where mixed forest occupied the largest land area. Coniferous forest
decreased sharply from 1986 to 2016, partly due to the infestation of mountain pine beetles,
generating the dead pine trees class covering 25.59 km2. Agriculture area increased from 20.24
in 1986 to 38.41 km2 in 2016, indicating the successful agriculture development over 30 years.
The change of deciduous forest in this area was fluctuated, and further research is needed to
explain this. The remaining features were mostly unchanged, such as the urban area.

69

Figure 4.11. SC sub-basin LULC map for 1986

70

Figure 4.12. SC sub-basin LULC map for 1999

71

Figure 4.13. SC sub-basin LULC map for 2016

72

BL

61.81
35.53

AA

25.59

30

37.45
29.11

40

27.76
26.77
23.24

50

38.41

60

42.59

55.18

70

2016

37.20
37.81
37.40

1999

80

20.24
25.96

AREA (KM2)

1986

73.55
65.75
66.06

STONRY CREEK LULC CHANGE

0

CF

DT

DF

MF

PF

2.04
3.19
2.83

CB

2.83
2.45
3.21

0.00
0.00

10

0.00
0.19
0.00

20

UA
WT
LULC CLASS

Figure 4.14. LULC change of SC sub-basin

4.4.2

Tachick Creek (TC) sub-basin

The TC sub-basin is located in the west corner of the study area. Streams in this sub-basin
flow into Tachick Lake. Figure 4.15 shows LULC change in the RC sub-basin. Figure 4.16
shows the LULC map of 1986, 1999 and 2016.

73

TACHICK CREEK LULC CHANGE
15.01

2016

AA

BL

CF

DT

0.00
0.00
0.00

2
0

8.18
4.29
5.26
1.56

4

2.14
2.87

6

5.91
5.92

8

0.00
0.00
1.74

10

8.75
7.63

12

10.56

10.94

14

11.08

1999

16

8.21

AREA (KM2)

1986

DF

MF

WT
LULC CLASS

Figure 4.15. LULC change of TC sub-basin

The forest land experienced some changes from 1986 to 2016 (Figure 4.15). As with the
SC sub-basin, coniferous forest decreased mainly because of the mountain pine beetle
infestation. The transitions from coniferous forest to the other LULC types also contributed to
this decrease, especially to the agriculture area. Agriculture area class had a 6.07 km2 growth
to 2016, due to agricultural development in this area. There are only a minimum number of
small water bodies within TC sub-basin and the area of water class showed as ‘0.00’.

74

Figure 4.16. TC sub-basin LULC map, (a) 1986, (b) 1999, (c) 2016

75

4.4.3

West Nulki Stream (WNS) sub-basin

This sub-basin covered almost half the land of the dead pine trees class with coniferous
forest representing the major land type in both 1986 and 1999. However, after suffering from
the mountain pine beetle infestation, the land area of coniferous forest class was down to only
3.96 km2 in 2016 (Figure 4.17). The land with dead pine trees is distributed on more than 60%
area of this sub-basin, which may cause negative effects on ecological and hydrological aspects.

1999

2.80
0.57
0.36

0.00
0.00

1.65
0.95
2.31

CB

DT

DF

30
25

2016

24.79

28.88

35

1986

23.49

AREA (KM2)

WEST NULKI STREAM LULC CHANGE

20

MF

0.25
0.25
0.27

CF

6.71
5.87

BL

2.53

AA

1.15
4.66
1.01

0

2.41
2.55
0.95

5

0.21
0.70
0.37

10

3.96

15

PF

WT
LULC CLASS

Figure 4.17. LULC change of WNS sub-basin

Figures 4.18 exhibits LULC maps of all three years with the WNS sub-basin extent.

76

Figure 4.18. WNS sub-basin LULC map, (a) 1986, (b) 1999, (c) 2016

77

4.4.4

South Nulki Stream (SNS) sub-basin

This sub-basin contains the stream which flows north and enters Nulki Lake at the south
part of the lake. Approximately 70% of this area was covered by mixed forest in both 1986 and
1999 (Figure 4.19). However, in 2016 most of the mixed forest land had transferred to
deciduous forest, which may indicate the forestry activities such as coniferous tree cutting; or
coniferous tress killed by mountain pine beetles in this area. Deciduous forest reached 8.57
km2 in 2016, while there was a slight increasing trend for coniferous forest.

1999

0.29
0.00
0.47

0.00
0.00
1.10

1.05
1.05

2016

CB

DT

DF

12.70
13.76

1986
16
14

8.57

12
10

AA

BL

CF

MF

0.00
0.00
0.00

0

1.65
0.06
0.62

2

0.10
0.26
0.29

4

2.52
2.54
2.09

6

1.19

5.61

8

1.63
2.26

AREA (KM2)

SOUTH NULKI STREAM LULC CHANGE

PF
WT
LULC CLASS

Figure 4.19. LULC change of SNS sub-basin

Figures 4.20 shows the LULC maps for the SNS sub-basin; the transition from mixed forest to
deciduous forest can be seen in the maps for 1999 and 2016.

78

Figure 4.20. SNS sub-basin LULC map, (a) 1986, (b) 1999, (c) 2016

79

4.4.5

Corkscrew Creek (CC) sub-basin

The CC sub-basin is in the south part of the study watershed. The stream within it flows
from the farthest south of SCW to the north and enters Nulki Lake. Coniferous forest is the
major LULC type, but it had a decreasing trend since 1986. Some of the coniferous land had
been cut or transferred to built-up areas and roads (Figures 4.21 and 4.22). This may indicate
the development of the forest industry during that time. Moreover, the coniferous forest had
decreased in 2016 because of the extent of dead pine trees. Most of the cut block area
transferred to planted or regrowth forest in 1999 and part of the forest had grown into mature
forest later in 2016. Therefore, there was a slight decrease for planted/regrowth forest, while
the mixed forest and deciduous forest increased.

3.70
1.39
8.82

2016

CB

DT

DF

9.67

20
10
5
0

1.21
2.95
0.88

15

AA

BL

CF

MF

Figure 4.21. LULC change of CC sub-basin
80

0.07
0.09
0.07

25

11.63
6.18

30

4.33

35

9.73
9.93
15.48

29.40

40

1999

0.00
0.00

45

1986

3.66
4.46
0.00

46.87
40.21

50

1.42
0.34
0.50

AREA (KM2)

CORKSCREW CREEK LULC CHANGE

PF
WT
LULC CLASS

Figure 4.22. CC sub-basin LULC map, (a) 1986, (b) 1999, (c) 2016

81

4.4.6

Mid Cork Stream (MCS) sub-basin

The MCS sub-basin is surrounded by the CC sub-basin and WCS sub-basin. The subbasin area had most forestry activities for all three years. Cut block increased from 1986 to
1999. In 2016, most of the cut block area transferred to planted/regrowth forest, while some
new cut block appeared. Similar to other sub-basins, coniferous forest had decreased sharply
because of the transitions to other LULC classes such as dead pine trees, cut block and barren
land. The cut block experienced a decrease from 1999 to 2016 as part of it had been planted or
regrown with young trees (Figures 4.23 and 4.24).

MID CORK STREAM LULC CHANGE
1999

2016

16.46

20.52

25
20

DT

DF

MF

CB

Figure 4.23. LULC change of MCS sub-basin

82

0.39
0.38
0.39

8.62
8.33
9.92
1.18
0.64
2.18

CF

0.10
0.00
2.71

BL

3.40
1.63

AA

0.00
0.00

0

0.01

5

0.16
1.53
0.11

10

6.61

7.42

15

0.00
0.22
0.01

AREA (KM2)

1986

PF
WT
LULC CLASS

Figure 4.24. MCS sub-basin LULC map, (a) 1986, (b) 1999, (c) 2016

83

4.4.7

West Cork Stream (WCS) sub-basin

The WCS sub-basin is located in the south-western part of the SCW study area. Its stream
flows south and reaches the stream of the CC sub-basin. Coniferous forest occupied the
maximum land area in this sub-basin. As it is close to the MCS sub-basin, similar effects
occurred from 1986 to 2016 within this area, e.g. the coniferous forest had a decrease trend
while the area of cut block, barren land and dead pine trees increased. The transition from cut
block to planted/regrowth forest occurred after 1999, and then part of this young forest grew
and transferred into other mature forest land (Figures 4.25 and 4.26).

WEST CORK STREAM LULC CHANGE
1999

2016

30.06
25.69

35
30

17.22

25
20

CF

CB

MF

PF

Figure 4.25. LULC change of WCS sub-basin

84

0.08
0.13
0.07

DF

5.68
7.18

DT

0.81

1.09
1.18
4.59

BL

0.90
0.00
1.18

0

0.00
0.00

5

3.16
3.22
1.64

10

4.85

15

0.98
1.19
0.36

AREA (KM2)

1986

WT
LULC CLASS

Figure 4.26. WCS sub-basin LULC map, (a) 1986, (b) 1999, (c) 2016

85

In summary, changes such as the decrease of coniferous forest and the increase of dead
pine trees from 1986 to 2016 are similar for the sub-basins because they were all affected by
the mountain pine beetle infestation before 2016. The Stoney Creek sub-basin and Tachick
Creek sub-basin which located in the north region of the SCW have rough terrain so agricultural
activities remain unchanged; however, agricultural policies such as reservation areas were
made during the period so agricultural areas increased from 1986 to 2016. Urban areas did not
have very clear boundaries as vegetation grows in and around the area so errors may happen
when separating urban areas from forest lands. The slight decrease of urban area from 1999 to
2016 may have resulted from this.

4.5 Land change modeling results and discussions
By using Land Change Modeler (LCM) in IDRISI software, LULC change analysis
between 1986 and 2016 was performed. The results of transition sub-models, Markov Chain
matrix, hard prediction and soft predictions were generated.
4.5.1

Results of transition sub-models

Based on the method mentioned in the previous chapter, effective driver variable groups
for each individual transition sub-model were selected based on the reports which generated
from running sub-models. This was accomplished by utilizing the MLP neural network on
LCM, as mentioned above. Every driver variable that had been selected for modeling process
should be positive to influence the overall accuracy and skill measure. Otherwise, the drivers
with no power were excluded from the group.
86

Running sub-models with ideal driver variable groups generated transition potential maps
for individual transitions. Using the transition from coniferous forest to dead pine trees as an
example, Figure 4.27 shows the change probability from 1986 to 2016 within the SCW study
extent.

Figure 4.27. Transition potential map of CF to DT transition

87

The map shows that most areas with more than 0.55 transition probabilities are distributed
in the middle part of the study area. In the 2016 image, the west part of the study watershed
was mostly covered by dead pine trees. Therefore, the area of dead pine trees may expand from
the west to east. The areas near dead pine trees are more likely to be affected by mountain pine
beetle infestation. Moreover, the more the position moves south, the less potential for the
transition.
4.5.2

Results of Markov Chain model

Markov Chain model determined the amount of LULC change using land cover images
from the earlier date to later date, and two matrixes were generated. Table 4.5 shows the
probability of LULC changing and persistence and Table 4.6 shows how much land is possible
to transition from the later date to the prediction date. Blank cells in both tables represent 0
value, and cells with gray shade means probability or pixel numbers of persistent area for each
class.
In Table 4.5, classes of urban area and water have very high probability with 0.917 and
0.991 respectively, meaning most of these areas can be persistent in 2026. Planted/regrowth
forest (0.765) and agriculture area (0.740) have high probability to maintain more than 70% of
their previous extent till the final period. The table tells that planted/regrowth forest is mostly
possible to transfer to mixed forest with 0.103. Besides the high probability of agriculture area
to maintain its area, barren land has 0.203 and deciduous forest has 0.167 probabilities to
transfer to agriculture area. These all represent the suitable agricultural management in this
area.
88

Around 50% of cut block will change to other LULC types in 2026. This can happen
because cut block area has been planted or regrowth to become other forest lands after decade.
Based on the probability matrix, barren land, deciduous forest and mixed forest have about 0.3
- 0.4 probabilities to maintain their extent. The probability is only 0.408 for mixed forest (Table
4.6), as it has the maximum pixels (areas) transferred from other classes through the modeling
period.
Coniferous forest has no probability to maintain its extent, as it will change mainly to
mixed forest and planted/regrowth forest with 0.433 and 0.446 probabilities respectively.
However, further research and more information about mountain pine beetle infestation areas
are needed to explain the transitions from coniferous forest to others. The 0.121 probability to
change to dead pine trees is also important because it is the most influential aspect for
coniferous forest decrease in 2016.
As mentioned before, there were no dead pine trees LULC type in the earlier date (1986).
Therefore, the LULC change determined by Markov Chain model by using the earlier date
image and later date image is not available for change prediction. Because of the lack of dead
pine trees data in the earlier image, the Markov Chain model generated the transition
probability value and the transition areas to other classes with equally divided values, which
are useless for prediction analysis. Therefore, the rows for dead pine trees in both tables are
filled as not available (N/A). Dead pine trees will not only remain its previous extent in the
final stage (2026), but also gain land which results from other transitions.

89

Table 4.5. Transition probability matrix from 2016 to 2026

Probability of changing to
Given
Class
AA
BL
CB
CF
DF
MF
PF
UA
WT
DT

AA

BL

0.740

0.178

0.203

0.399
0.013

CB

CF

0.543

DF

MF

0.056

0.027

0.002

0.088

0.014

0.044

0.000
0.167

0.055

0.011

0.023

0.019

0.247

0.011

0.053

0.030

NA

0.002

N/A

N/A

N/A

UA

WT

DT

0.232

0.015

0.003

0.001

0.058

0.057

0.024

0.001

0.306

0.433

0.446

0.303

0.459

0.293

0.408

0.006

0.012

0.103

0.765

0.083

0.000

0.002

0.004

0.000

N/A

N/A

N/A

0.000
0.002

PF

0.121
0.004
0.002

0.002
0.027

0.917
0.991
N/A

N/A

N/A

Table 4.6. Transition areas matrix from 2016 to 2026

Expected to transition to
Given
AA
BL
CB
Class
AA 353658 84949
65008 127963
BL
12067 511509
CB
CF
112539 37289
DF
24277
20018 262893
MF
3225
15740
PF
UA
611
573
WT
N/A
N/A
N/A
DT

CF

DF

MF

26531

12859

706

28108

13104

41007

0
7539
8787

UA

WT

DT

74346

4941

931

353

18514

53545

23002

566

287995

17784

18322

203944

308810

312085

434428

6282

3610

30326

226309

9
N/A

PF

2894

2365

2236

2555
7870

26091

649

1414

38

N/A

N/A

N/A

90

4984

378724
N/A

N/A

N/A

4.5.3

Results of soft prediction and hard prediction

By running a multi-objective land allocation algorithm procedure, all the changed land of
a class was allocated and overlain, then a hard prediction map was generated for 2026 (Figure
4.28).

Figure 4.28. LULC hard prediction map for 2026
91

Table 4.7 and Figure 4.29 list the area statistics for both stages (2016 and 2026). Dead
pine trees will gain areas from other classes, but it is not sure if it will face an increasing trend
because of the lack of data for hard prediction. Based on the comparison between 2016 and
2026, agricultural area shows an increasing trend in the next decade which refers to a stable
and planned agricultural development in this area. Forest land will experience a slight decrease
according to the hard prediction data, while the remaining classes have no significant change.

Table 4.7. Area statistics of 2016 and 2026

2016
LULC Class
Agriculture Area
Barren Land
Coniferous Forest
Cut Block
Deciduous Forest
Mixed Forest
Planted/Regrowth Forest
Urban Area
Water
Dead Pine Trees

2026

Area (km2)

Percentage

47.786
32.072
94.282
4.109
67.306
106.444
29.585
2.835
38.205
74.339

9.62%
6.45%
18.97%
0.83%
13.54%
21.42%
5.95%
0.57%
7.69%
14.96%

92

Area (km2)
59.492
30.773
77.487
1.809
59.565
93.031
30.687
2.846
38.205
103.129

Percentage
11.97%
6.19%
15.59%
0.36%
11.98%
18.72%
6.17%
0.57%
7.69%
20.75%

Area (km2)

120
100
80
60
40
20
0

AA

BL

CF

CB

DF

2016

MF

PT

UA

WT
DT
LULC Class

2026

Figure 4.29. Area comparison between 2016 and 2026

Similar to the hard prediction, soft prediction is also short of earlier data of dead pine trees
class. There was no sub-model for the transition from dead pine trees to other LULC types.
Therefore, in the soft prediction map, the area of dead pine trees in 2016 will be persistent and
impossible to change. Following are two types of soft prediction maps. Figure 4.30 shows the
logical OR soft prediction map, which is calculated as (a+b – ab) for all transitions. This kind
of soft prediction can describe areas with multiple transitions. The other is maximum soft
prediction that represents a maximum transition probability aggregation for all transition
potentials (Figure 4.31).

93

Figure 4.30. Logical OR soft prediction map for 2026

94

Figure 4.31. Maximum soft prediction map for 2026

95

Referring to both soft prediction maps, the middle and north region of the study area have
the highest potential to change. These transitions may be attributed to the agricultural
development and forestry activities in this area. In the central north area surrounding two lakes,
the transition probability is variable. Further south, the potential shows a decreasing trend.
However, the planted/regrowth forest extent has low transition potential, which means a high
probability to be persistent. In both soft prediction maps, water and urban area show no
probability to change.

96

CHAPTER 5 CONCLUSION
5.1 Summary of the study
This study generated the LULC change analysis and prediction in SCW by using remote
sensing, GIS and modeling methods. LULC maps with composite bands and classified maps
were generated and analyzed to provide a better understanding of LULC patterns and dynamics.
Transition potentials were then modeled as sub-models, and several LULC scenarios were
finally predicted.
Several periods of field survey were made for collecting overall information and ground
truth data in the beginning of the study. Landsat image data were collected for July – August,
1986, 1999 and 2016 and clipped to the SCW extent. The OBIA technique was utilized to
classify pre-processed LULC maps and resulted in an overall classification accuracy at 91.50%.
These three maps were then applied in LULC change analysis for two study periods, from
1986 to 1999 and from 1999 to 2016. Forest lands had been the major land cover within the
watershed from 1986 to 2016, but they had a decreasing trend especially for coniferous forest.
The LULC class of coniferous forest experienced a sharp land loss because of the infestation
of mountain pine beetles, which led to the transition from coniferous forest to dead pine trees
in map 2016. Agriculture area had a stable increasing trend (23.68 km2 land gain) for both study
periods, which may refer to the successful agricultural development activities in this watershed.
Moreover, more detailed change analysis was performed at sub-basin level, as LULC patterns
and dynamics vary from one sub-basin to another. It has been found that the dead pine trees
97

were mainly distributed in the south-western area and expanded from west to east over time.
Based on the LULC change analysis data from 1986 to 2016, transition potentials were
modeled by MLP neural network. Each transition sub-model was performed with a group of
suitable driver variables, which can improve the accuracy and skill measure. Markov Chain
model was utilized to generate future transition matrixes of 2026 by using the data resulted for
LULC change analysis from the earlier date to later date. Matrixes data were later applied for
hard prediction result generation, which represents a predicted LULC map. Transition
potentials were computed together and transferred into soft prediction results, which show the
transition probabilities calculated by Logic OR and Maximum methods. Based on the results
of soft prediction, forestry lands with higher elevations had lower transition probabilities in the
future. Areas with the highest potential to change were mainly distributed in the middle and
north part of watershed, which may refer to agricultural and forestry activities in these areas.
To summarize, this thesis study utilized one of the most widely applied geospatial
analytical method (combined RS and GIS approach) with a more recently developed image
classification method (OOIC) to generate final maps. It also performed a powerful modeling
approach (Multi-layer perceptron and Markov Chain model) for fine LULC scenarios. Results
of this study can provide accessible information for planners and researchers from LULC or
other fields, which can be supportive for both environmental and economic development in
this area.

98

5.2 Limitations and recommendations
There are several limitations for this study:
a. Free-accessed Landsat images with 30-meter spatial resolutions is impossible to classify
some features such as the main channel of Stoney Creek.
b. The lack of comprehensive historical and current ground truth data, due to the limited
accessibility of part of the watershed.
c. Restricted official information about land management for SCW limited the analysis on
LULC potential transitions.
d. The lack of ‘Dead Pine Trees’ class in earlier images meant the Markov Chain modeling
method cannot model the transitions from dead pine trees to other classes. This problem
may cause negative effects on LULC change prediction process.
For the recommendations for further study in the future, firstly, satellite images and aerial
photographs with higher resolution or larger study extent may result in more ideal LULC maps
and more detailed study. Secondly, if more comprehensive ground truth data could be applied,
more accurate and strict classification results will be created. The third recommendation is to
collect more satellite images during the period when mountain pine beetle infestation is active,
so that a better LULC change modeling can be performed. Moreover, it is useful for a further
study on when and how mountain pine beetle infestation expended in this area. Finally, more
detailed climate change studies such as climate modeling and the relationship between LULC
changes and climate changes are meaningful for SCW or larger regions.

99

Bibliography
Abuelaish, B., & Olmedo, M. T. C. (2016). Scenario of land use and land cover change in the
Gaza Strip using remote sensing and GIS models. Arabian Journal of Geosciences, 9(4),
274.
Agricultural Land Use Inventory. (2017). Vanderhoof & Electoral Area F, Regional District of
Bulkley-Nechako. British Columbia Ministry of Agriculture.
Anderson, J. R. (1976). A land use and land cover classification system for use with remote
sensor data (Vol. 964). US Government Printing Office.
Araya, Y. H. (2009). Urban land use change analysis and modeling: a case study of Setubal
and Sesimbra, Portugal. Münster, Germany: Institute for Geoinformatics, University of
Münster.
Araya, Y. H., & Cabral, P. (2010). Analysis and modeling of urban land cover change in
Setúbal and Sesimbra, Portugal. Remote Sensing, 2(6), 1549-1563.
Arekhi, S., & Jafarzadeh, A. A. (2014). Forecasting areas vulnerable to forest conversion using
artificial neural network and GIS (case study: northern Ilam forests, Ilam province,
Iran). Arabian Journal of Geosciences, 7(3), 1073-1085.
Baker, W. L. (1989). A review of models of landscape change. Landscape Ecology, 2 (2), 111133.
Barsi, J. A., Lee, K., Kvaran, G., Markham, B. L., & Pedelty, J. A. (2014). The spectral response
of the Landsat-8 operational land imager. Remote Sensing, 6 (10), 10232-10251.
Benito, P., Cuevas, J., Parra, R., Prieto, F., Barrio, J., & Zavala, M. (2010). Land use change in
a Mediterranean metropolitan region and its periphery: assessment of conservative
policies through CORINE Land Cover data and Markov models. Forest Systems, 19 (3),
315-328.
Benz, U. C., Hofmann, P., Willhauck, G., Lingenfelder, I., & Heynen, M. (2003). Multiresolution, object-oriented fuzzy analysis of remote sensing data for GIS-ready
information. ISPRS Journal of Photogrammetry and Remote Sensing, 239-258.
Blaschke, T. (2010). Object based image analysis for remote sensing. ISPRS Journal of
Photogrammetry and Remote Sensing (65), 2-16.
Brink, A. B., & Eva, H. D. (2009). Monitoring 25 years of land cover change dynamics in
Africa: A sample based remote sensing approach. Applied Geography, 29 (4), 501-512.

100

Butt, A., Shabbir, R., Ahmad, S. S., & Aziz, N. (2015). Land use change mapping and analysis
using Remote Sensing and GIS: A case study of Simly watershed, Islamabad,
Pakistan. The Egyptian Journal of Remote Sensing and Space Science, 18(2), 251-259.
Campbell, J. B., & Wynne, R. H. (2011). Introduction to Remote Sensing. Guilford Press.
Chen, J., Mao, Z., Philpot, B., Li, J., & Pan, D. (2013). Detecting changes in high-resolution
satellite coastal imagery using an image object detection approach. International Journal
of Remote Sensing, 34(7), 2454-2469.
Chen, X. L., Zhao, H. M., Li, P. X., Yin, Z. Y. (2006). Remote sensing image-based analysis
of the relationship between urban heat island and land use/cover changes. Remote Sensing
of Environment, 104(2), 133-146.
Clark Labs, (2009). The Land Change Modeler for ecological Sustainability. IDRISI Focus
Paper,
Worcester,
MA:
Clark
University.
http://www.clarklabs.org/applications/upload/Land-Change-Modeler-IDRISI-FocusPaper-pdf.
Clarke, K. C., & Gaydos, L. J. (1998). Loose-coupling a cellular automaton model and GIS:
long-term
urban
growth
prediction
for
San
Francisco
and
Washington/Baltimore. International
Journal
of
Geographical
Information
Science, 12(7), 699-714.
D BEHERA, M. D., Borate, S. N., Panda, S. N., Behera, P. R., & Roy, P. S. (2012). Modelling
and analyzing the watershed dynamics using Cellular Automata (CA)–Markov model–A
geo-information based approach. Journal of Earth System Science, 121(4), 1011-1024.
Dean, A. M., & Smith, G. M. (2003). An evaluation of per-parcel land cover mapping using
maximum likelihood class probabilities. International Journal of Remote Sensing, 24,
2905-2920.
Dewan, A. M., & Yamaguchi, Y. (2009). Using remote sensing and GIS to detect and monitor
land use and land cover change in Dhaka Metropolitan of Bangladesh during 1960–
2005. Environmental Monitoring and Assessment, 150(1-4), 237.
Dimyati, M., Mizuno, K., Kobayashi, S., & Kitamura, T. (1996). An analysis of land use/cover
change in Indonesia. International Journal of Remote Sensing, 17(5), 931-944.
Duram, M., Venancio, R., & Ribeiro, S. (2004). Influência das emoções na cognição. www. ic.
unicamp. br/~ wainer/cursos/906/trabalhos/Trabalho_E1. pdf. Acessado em, 20(05),
2013.
Dymond, J. R., Shepherd, J. D., Newsome, P. F., Gapare, N., Burgess, D. W., & Watt, P. (2012).
Remote sensing of land-use change for Kyoto protocol reporting: the New Zealand case.
Environmental Science & Policy, 16, 1-8.
101

Eastman, J. R. (2012). IDRISI Selva manual and tutorial manual version 17. Worcester (USA):
Clark University.
Eastman, J. R. (2012). Idrisi Selva tutorial. Idrisi Production, Clark Labs-Clark University.
Ebata, T. (2003, October). Current status of mountain pine beetle in British Columbia.
In Mountain pine beetle symposium: Challenges and solutions (pp. 52-56).
Echeverria, C., Coomes, D. A., Hall, M., & Newton, A. C. (2008). spatially explicit models to
analyze forest loss and fragmentation between 1976 and 2020 in southern Chile.
Ecological Modelling, 212, 439-449.
Foley, J., DeFries, R., Asner, G., Barford, C., Bonan, G., Carpenter, S., et al. (2005). Global
consequences of land use. Science (309), 570-574.
Freier, K. P., Schneider, U. A., & Finckh, M. (2011). Dynamic interactions between vegetation
and land use in semi-arid Morocco: Using a Markov process for modeling rangelands
under climate change. Agriculture, Ecosystems & Environment, 140(3-4), 462-472.
Frohn, R. C., Autrey, B. C., Lane, C. R., & Reif, M. (2011). Segmentation and object-oriented
classification of wetlands in a karst Florida landscape using multi-season Landsat-7
ETM+ imagery. International Journal of Remote Sensing, 32(5), 1471-1489.
Gardner, M. W., & Dorling, S. R. (1998). Artificial neural networks (the multilayer
perceptron)—a review of applications in the atmospheric sciences. Atmospheric
environment, 32(14-15), 2627-2636.
Geomatica. (2017). Geomatica Training Guide: Geomatica II course guide. PCI Geomatics
Enterprises, Inc..
Haberl, H., Erb, K. H., Krausmann, F., Gaube, V., Bondeau, A., Plutzar, C., ... & FischerKowalski, M. (2007). Quantifying and mapping the human appropriation of net primary
production in earth's terrestrial ecosystems. Proceedings of the National Academy of
Sciences, 104 (31), 12942-12947.
Halmy, M. W. A., Gessler, P. E., Hicke, J. A., & Salem, B. B. (2015). Land use/land cover
change detection and prediction in the north-western coastal desert of Egypt using
Markov-CA. Applied Geography, 63, 101-112.
Hasbani, J. G. (2008, January). Semi-automated calibration of a cellular automata model to
simulate land-use changes in the Calgary region. In Masters Abstracts International (Vol.
46, No. 06).
Hay, G., & Castilla, G. (2008). Geographic object-based image analysis (GEOBIA): A new
name for a new discipline. In Object based image analysis (pp. 93-112). Heidelberg,
Berlin: Springer.
102

Hobbs, R. J. (1998). Impacts of land use on biodiversity in southwestern Australia.
In Landscape disturbance and biodiversity in Mediterranean-type ecosystems (pp. 81106). Springer Berlin Heidelberg.
Houet, T., & Hubert-Moy, L. (2006). Modeling and projecting land-use and land-cover changes
with Cellular Automaton in considering landscape trajectories. EARSeL
eProceedings, 5(1), 63-76.
Hua, A. K. (2017). Application of Ca-Markov model and land use/land cover changes in
Malacca River Watershed, Malaysia. Applied Ecology and Environmental Research, 15(4),
605-622.
Islam, K., Jashimuddin, M., Nath, B., & Nath, T. (2016). Quantitative Assessment of land cover
change using landsat time series data: case of Chunati Wildlife Sanctuary (CWS),
Bangladesh. International Journal of Environment and Geoinformatics, 3(2), 45-55.
Islam, K., Jashimuddin, M., Nath, B., & Nath, T. K. (2018a). Land use classification and
change detection by using multi-temporal remotely sensed imagery: The case of Chunati
wildlife sanctuary, Bangladesh. The Egyptian Journal of Remote Sensing and Space
Science, 21(1), 37-47.
Islam, K., Rahman, M. F., & Jashimuddin, M. (2018b). Modeling land use change using
cellular automata and artificial neural network: the case of Chunati Wildlife Sanctuary,
Bangladesh. Ecological indicators, 88, 439-453.
Jensen, J. R. (1996). Introductory digital image processing: a remote sensing perspective (No.
Ed. 2). Prentice-Hall Inc.
Kamusoko, C., Wada, Y., Furuya, T., Tomimura, S., Nasu, M., & Homsysavath, K. (2013).
Simulating future forest cover change in Pakxeng District, Lao People's Democratic
Republic (PDR): Implications for sustainable forest management. Land, 2, 1-19.
Kline, J., Moses, A., Lettman, G. J., & Azuma, D. L. (2007). Modeling forest and range land
development in rural locations, with example from eastern Oregon. Urban Planning, 80,
320-332.
Lambin, E. F., Geist, H. J., & Iepers, E. (2003). Dynamics of land-use and land-cover change
in tropical regions. Annual Review of Environmental Resources, 28, 205-241.depending
on the amount
Lillesand, T. M., Kiefer, R. W., & Chipman, J. W. (2008). Remote sensing and image
interpretation (6th ed.). Hoboken, NJ: John Wiley & Sons, Inc.
Lin, Y. P., Chu, H. J., Wu, C. F., & Verburg, P. H. (2011). Predictive ability of logistic regression,
auto-logistic regression and neural network models in empirical land-use change
modeling–a case study. International Journal of Geographical Information
103

Science, 25(1), 65-87.
Liping, C., Yujun, S., & Saeed, S. (2018). Monitoring and predicting land use and land cover
changes using remote sensing and GIS techniques—A case study of a hilly area, Jiangle,
China. PloS one, 13(7), e0200493.
MacLean, M., Campbell, M., Maynard, D., Ducey, M., & Congalton, R. (2013). Requirements
for labelling forest polygons in an object-based image analysis classification.
International Journal of Remote Sensing, 34 (7), 2531-2547.
Mishra, V. N., & Rai, P. K. (2016). A remote sensing aided multi-layer perceptron-Markov
chain analysis for land use and land cover change prediction in Patna district (Bihar),
India. Arabian Journal of Geosciences, 9(4), 249.
Mishra, V. N., Rai, P. K., & Mohan, K. (2014). Prediction of land use changes based on land
change modeler (LCM) using remote sensing: a case study of Muzaffarpur (Bihar),
India. Journal of the Geographical Institute" Jovan Cvijic", SASA, 64(1), 111-127.
Mitsova, D., Shuster, W., & Wang, X. (2011). A cellular automata model of land cover change
to integrate urban growth with open space conservation. Landscape and Urban
Planning, 99(2), 141-153.
Mozumder, C., & Tripathi, N. K. (2014). Geospatial scenario based modelling of urban and
agricultural intrusions in Ramsar wetland Deepor Beel in Northeast India using a multilayer perceptron neural network. International Journal of Applied Earth Observation and
Geoinformation, 32, 92-104.
Muñoz-Villers, L. E., & López-Blanco, J. (2008). Land use/cover changes using Landsat
TM/ETM images in a tropical and biodiverse mountainous area of central‐eastern
Mexico. International Journal of Remote Sensing, 29(1), 71-93.
NASA Earth Observatory. (Aug 30, 2000). Measuring Vegetation (NDVI & EVI), Retrieved
from
the
NASA
Earth
Observatory
website:
https://earthobservatory.nasa.gov/features/MeasuringVegetation/measuring_vegetation_
2.php.
NASA. (2011). Landsat science data users handbook. Washington, D.C.: National Aeronautics
and Space Administration.
Nguyen, D. T., Iskandar, I., & Ho, S. (2016). Land cover change and the CO2 stock in the
Palembang City, Indonesia: A study using remote sensing, GIS technique and
LUMENs. The Egyptian Journal of Remote Sensing and Space Science, 19(2), 313-321.
Olin Albertson R. P. Bio. (2013). Stoney Creek Preliminary Watershed Health Assessment.
Nechako Environment & Water Stewardship Society (NEWSS). Avison Management
Services Ltd., Vanderhoof, BC
104

Oumer, H. A. (2009). Land use and land cover change, drivers and its impact: A comparative
study from Kuhar Michael and Lenche Dima of Blue Nile and Awash Basins of Ethiopia.
Ithaca, NY: Cornell University.
Paiboonvorachat, C. (2008). Using remote sensing and GIS techniques to assess land use/land
cover changes in the Nan Watershed, Thailand. Carbondale: Dept. of Geography and
Environmental Resources, Southern Illinois University.
Paiboonvorachat, C., & Oyana, T. J. (2011). Land-cover changes and potential impacts on soil
erosion in the Nan watershed, Thailand. International journal of remote sensing, 32(21),
6587-6609.
Parsa, V. A., Yavari, A., & Nejadi, A. (2016). Spatio-temporal analysis of land use/land cover
pattern changes in Arasbaran Biosphere Reserve: Iran. Modeling Earth Systems and
Environment, 2(4), 1-13.
Paul, S. S. (2013). Analysis of land use and land cover change in Kiskatinaw River watershed:
A remote sensing, GIS & modeling approach (Master Thesis, University of Northern
British Columbia, Prince George, BC, Canada).
Paul, S. S., Hasan, K., Akhter, S. H., & Li, J. (2012). Application of remote sensing and GIS
for investigating urbanization induced surface water body change in and around
Bangladesh. 12th International Environmental Specialty Conference. Edmonton:
Canadian society for Civil Engineering (CSCE).
Paul, S. S., Li, J., Wheate, R., & Li, Y. (2018). Application of Object Oriented Image
Classification and Markov Chain Modeling for Land Use and Land Cover Change
Analysis. Journal of Environmental Informatics, 31(1).
Pervaiz, W., Uddin, V., Khan, S. A., & Khan, J. A. (2016). Satellite-based land use mapping:
comparative analysis of Landsat-8, Advanced Land Imager, and big data Hyperion
imagery. Journal of Applied Remote Sensing, 10(2), 026004.
Pijanowski, B. C., Brown, D. G., Shellito, B. A., & Manik, G. A. (2002). Using neural networks
and GIS to forecast land use changes: a Land Transformation Model. Computer,
Environment and Urban Systems, 26 (6), 553-575.
Province of BC. (2016, September 16). Mountain pine beetle in BC. Retrieved from
https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forestresources/forest-health/forest-pests/bark-beetles/mountain-pine-beetle.
Qiang, Y., & Lam, N. S. (2015). Modeling land use and land cover changes in a vulnerable
coastal region using artificial neural networks and cellular automata. Environmental
monitoring and assessment, 187(3), 57.

105

Quinn.
J.
W.
(2001).
Band
combinations,
http://web.pdx.edu/~emch/ip1/bandcombinations.html.

Retrieved

from

website:

Rawat, J. S., and Kumar, M. (2015). Monitoring land use/cover change using remote sensing
and GIS techniques: A case study of Hawalbagh block, district Almora, Uttarakhand, India.
The Egyptian Journal of Remote Sensing and Space Sciences, 18, 77-84.
Regional District of Bulkley-Nechako (RDBN) (2014), Vanderhoof and surrounding area
profile. Retrieved from website: https://www.rdbn.bc.ca/economicdevelopment/regionalinformation/regional-profile.
Schulp, C., Nabuurs, G., & Verburg, P. (2008). Future carbon sequestration in Europe-effects
of land use change. Agriculture, Ecosystems & Environment, 127, 251-264.
Serneels, S., & Lambin, E. F. (2001). Proximate causes of land-use change in Narok District,
Kenya: a spatial statistical model. Agriculture, Ecosystems & Environment, 85 (1-3), 6581.
Shackelford, A. K., & Davis, C. H. (2003). A combined fuzzy pixel-based and object-based
approach for classification of high-resolution multispectral data over urban areas. IEEE
Transactions on GeoScience and Remote sensing, 41(10), 2354-2363.
Sörensen, R., Zinko, U., & Seibert, J. (2006). On the calculation of the topographic wetness
index: evaluation of different methods based on field observations. Hydrology and Earth
System Sciences Discussions, 10(1), 101-112.
Srivastava, P. K., Singh, S. K., Gupta, M., Thakur, J. K., & Mukherjee, S. (2013). Modeling
impact of land use change trajectories on groundwater quality using remote sensing and
GIS. Environmental Engineering & Management Journal (EEMJ), 12(12).
Stabile, M. (2012). Deconstructing the complexity of land use and cover classification and land
change modeling. New South Wales, Australia: Faculty of Agriculture, Food and Natural
Resources, The University of Sydney.
Steffen, W. L., Walker, B. H., Ingram, J. S., and Koch, G. W., 1992, Global change and
terrestrial ecosystems: the operational plan. IGBP Report No. 21, International
Geosphere–Biosphere Programme, Stockholm.
Steininger, M. K. (1996). Tropical secondary forest regrowth in the Amazon: age, area and
change estimation with Thematic Mapper data. International Journal of Remote
Sensing, 17(1), 9-27.
Sun, H., Forsythe, W., & Waters, N. (2007). Modeling urban land use change and urban sprawl:
Calgary, Alberta, Canada. Networks and Spatial Economics, 7, 353-376.
Tran, D. X., Pla, F., Latorre-Carmona, P., Myint, S. W., Caetano, M., & Kieu, H. V. (2017).
106

Characterizing the relationship between land use land cover change and land surface
temperature. ISPRS journal of photogrammetry and remote sensing, 124, 119-132.
Turner, B., Lambin, E., & Reenberg, A. (2007). The emergence of land change science for
global environmental change and sustainability. Proceedings of the National Academy of
Sciences. 104, pp. 20666-20671. PNAS.
Turner, M. G., & Ruscher, C. L. (1988). Changes in landscape patterns in Georgia,
USA. Landscape ecology, 1(4), 241-251.
U.S. Geological Survey. (2007, February 13). Global Geographic Information Systems.
Retrieved
May 29,
2013,
from
U.S.
Geological
Survey website:
http://webgis.wr.usgs.gov/globalgis/tutorials/what_is_gis.htm
U.S. Geological Survey. (2012). Landsat–A Global Land–Imaging Mission: U.S. Geological
Survey Fact Sheet: 2012-3072, 4 p., https://doi.org/10.3133/fs20123072.
U.S. Geological Survey. (2015, September 09). Frequently Asked Questions about the Landsat
Missions.
from
U.S.
Geological
Survey
website:
http://landsat.usgs.gov/best_spectral_bands_to_use.php.
U.S. Geological Survey. (2016). Landsat–Earth observation satellites (ver. 1.1, August 2016):
U.S.
Geological
Survey
Fact
Sheet
2015–3081,
4
p.,
http://dx.doi.org/10.3133/fs20153081.
U.S. Geological Survey. (2018). Earth Resources Observation and Science (EROS) Archive –
Landsat Archives: Landsat 4-5 Thematic Mapper (TM) Level-1 Data Products, Retrieved
from U.S. Geological Survey website: https://www.usgs.gov/centers/eros/science/usgseros-archive-landsat-archives-landsat-4-5-thematic-mapper-tm-level-1-data?qtscience_center_objects=0#qt-science_center_objects
U.S. Geological Survey. (2018). Earth Resources Observation and Science (EROS) Archive –
Landsat Archives: Landsat 8 OLI (Operational Land Imager) and TIRS (Thermal Infrared
Sensor),
Retrieved
from
U.S.
Geological
Survey
website:
https://www.usgs.gov/centers/eros/science/usgs-eros-archive-landsat-archives-landsat-45-thematic-mapper-tm-level-1-data?qt-science_center_objects=0#qtscience_center_objects
United Nations. (1986). Principles relating to remote sensing of the Earth from space. United
Nations, General Assembly, A/RES/41/65.
Urban, D. L., & Wallin, D. O. (2017). Introduction to Markov models. In Learning Landscape
Ecology (pp. 129-142). Springer, New York, NY.
Verburg, P. H., Eickhout, B., & van Meijl, H. (2008). A multi-scale, multi-model approach for
analyzing the future dynamics of European land use. The annals of regional science, 42(1),
107

57-77.
Verburg, P. H., Schot, P. P., Dijst, M. J., & Veldkamp, A. (2004). Land use modelling: current
practice and research priorities. GeoJournal , 61, 309-324.
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, Forthcoming.
Weng, Q. (2002). Land use change analysis in the Zhuajiang Delta of China using satellite
remote sensing, GIS and stochastic modeling. Journal of Environmental Management
(2002), 64, 273-284.
Yang, X. J., and Liu, Z. (2005). Using satellite imagery and GIS for land-use and land-cover
change mapping in an estuarine watershed. International Journal of Remote Sensing,
26(23), 5275-5296
Yeh, A. G. O. and Li, X. (1997). An integrated remote sensing-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.
Zhao, Y., Zhang, K., Fu, Y., & Zhang, H. (2012). Examining land-use/land-cover change in the
Lake Dianchi Watershed of the Yunnan-Guizhou Plateau of southwest China with remote
sensing and GIS techniques: 1974–2008. International Journal of Environmental
Research and Public Health, 9(11), 3843-3865.
Zhou, W., & Troy, A. (2008). An object-oriented approach for analysing and characterizing
urban landscape at the parcel level. International Journal of Remote Sensing, 29(11),
3119-3135.

108