UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Monitoring miombo woodlands of Southern Africa with multi-source information in a model-based framework Halperin, James J. 2017

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
24-ubc_2017_may_halperin_james.pdf [ 1.74MB ]
Metadata
JSON: 24-1.0342942.json
JSON-LD: 24-1.0342942-ld.json
RDF/XML (Pretty): 24-1.0342942-rdf.xml
RDF/JSON: 24-1.0342942-rdf.json
Turtle: 24-1.0342942-turtle.txt
N-Triples: 24-1.0342942-rdf-ntriples.txt
Original Record: 24-1.0342942-source.json
Full Text
24-1.0342942-fulltext.txt
Citation
24-1.0342942.ris

Full Text

MONITORING MIOMBO WOODLANDS OF SOUTHERN AFRICA WITH MULTI-SOURCE INFORMATION IN A MODEL-BASED FRAMEWORK by James J. Halperin B.A., Political Science, University of Southern Maine, 1993 A.A.S., Geographic Information Systems, Central Oregon Community College, 2000 M.S., Forestry, North Carolina State University, 2003 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (FORESTRY) THE UNIVERSITY OF BRITISH COLUMBIA (VANCOUVER)  February 2017 © James J. Halperin, 2017 ii  Abstract Acquiring forest resources information for tropical developing countries is challenging due to financial and logistical constraints, yet this information is critical for enhancing management capability and engaging in initiatives such as Reducing Emissions from Deforestation and forest Degradation (REDD+). In this dissertation, I investigated innovative approaches to monitoring forest resources and deepened understanding of multi-source information (i.e., remote-sensing, environmental, and disturbance data) needs by examining methods using a model-based framework for assessment of a jurisdictional landscape in the miombo ecoregion of Zambia. I focused on percent canopy cover (CC) and above-ground biomass (AGB) within a National Forest Inventory (NFI) context because both are important for land management and REDD+. First, I compared multi-source information and four modeling methods to estimate CC and total forest area. Landsat outperformed RapidEye and a generalized additive model was most precise. Available soil water content (AWC), slope, distance to district capital, and texture of remotely sensed data were crucial predictor variables in improving estimates. Second, I used multi-source information and compared methods with and without predicted CC in three modeling methods to estimate total AGB. A nonlinear sigmoidal model was most precise when using predicted CC, AWC, pH, occurrence of late season fire, and the Normalized Difference Moisture Index as predictor variables. Third, I developed an innovative monitoring framework using time series classification and a stock-difference approach to estimate change in land cover and AGB over a 16-year annual time series of Landsat data. Forest/nonforest change trajectories were used to develop change classes relevant to underlying biotic and abiotic factors and provided ecologically meaningful context for land cover change and AGB change. Overall, predictor iii  variables related to soil moisture, topography, shortwave infrared bands and texture of vegetation indices from remotely sensed data are vital to accurate models of CC and AGB. Genetic algorithms provide an opportune method for predictor variable selection across a diversity of modeling methods. Further, robust change estimates are feasible when using annual monitoring methods based on freely available, multi-source information. In conclusion, the model-based framework provides a precise, statistically sound, approach to estimating and mapping forest attributes within an NFI context.  iv  Preface The contributions of committee members to each research chapter are summarized in the tables below. Chapter 2. Categories James Halperin Valerie LeMay Nicholas Coops Peter Marshall Louis Verchot Total Identify research problem 86 8 2 2 2 100 Designing research 88 6 2 2 2 100 Analyzing data 94 6 0 0 0 100 Manuscript writing 85 6 3 3 3 100  Chapter 3. Categories James Halperin Valerie LeMay Nicholas Coops Peter Marshall Louis Verchot Total Identify research problem 94 4 2 0 0 100 Designing research 96 4 0 0 0 100 Analyzing data 98 2 0 0 0 100 Manuscript writing 95 2 1 1 1 100   v  Chapter 4. Categories James Halperin Valerie LeMay Nicholas Coops Peter Marshall Louis Verchot Total Identify research problem 95 2 1 1 1 100 Designing research 97 1 2 0 0 100 Analyzing data 100 0 0 0 0 100 Manuscript writing 92 2 2 2 2 100  A version of Chapter 3 was published as: Halperin, J., LeMay, V., Coops, N., Verchot, L., Marshall, P., & Lochhead, K. (2016). Canopy cover estimation in miombo woodlands of Zambia: Comparison of Landsat 8 OLI versus RapidEye imagery using parametric, nonparametric, and semiparametric methods. Remote Sensing of Environment, 179, 170-182. A version of Chapter 4 was published as: Halperin, J., LeMay, V., Chidumayo, E., Verchot, L., & Marshall, P. (2016). Model-based estimation of above-ground biomass in the miombo ecoregion of Zambia. Forest Ecosystems, 3(1), 14. A version of Chapter 5 has been submitted for publication as follows: Halperin, J., Chidumayo, E., LeMay, V., Coops, N., Verchot, L., & Marshall, P. A framework for estimating change in forest area and above-ground biomass in the miombo ecoregion.  vi  Table of Contents  Abstract .......................................................................................................................................... ii Preface ........................................................................................................................................... iv Table of Contents ......................................................................................................................... vi List of Tables ................................................................................................................................ ix List of Figures ............................................................................................................................... xi List of Abbreviations ................................................................................................................. xiii Acknowledgements ......................................................................................................................xv Dedication ................................................................................................................................. xviii Chapter 1: Introduction ................................................................................................................1 1.1 Background ..................................................................................................................... 1 1.2 Research goal and dissertation overview ........................................................................ 6 Chapter 2: Study area and ground data ....................................................................................10 2.1 Study area...................................................................................................................... 10 2.2 Ground-based plot measurements ................................................................................. 15 Chapter 3: Canopy cover and forest area modeling .................................................................17 3.1 Introduction ................................................................................................................... 17 3.2 Methods......................................................................................................................... 21 3.2.1 Canopy cover per ground plot................................................................................... 21 3.2.2 Satellite imagery pre-processing ............................................................................... 23 3.2.3 Predictor variables .................................................................................................... 28 vii  3.2.4 Canopy cover models ................................................................................................ 33 3.2.4.1 Effects of auxiliary variables ............................................................................ 33 3.2.4.2 Modeling methods ............................................................................................ 37 3.2.4.3 Model evaluations ............................................................................................. 41 3.3 Results ........................................................................................................................... 43 3.3.1 Effects of auxiliary variables .................................................................................... 43 3.3.2 Landsat 8 OLI versus RapidEye ............................................................................... 44 3.3.3 Comparison of methods ............................................................................................ 46 3.4 Discussion ..................................................................................................................... 50 Chapter 4: Above-ground biomass modeling ............................................................................57 4.1 Introduction ................................................................................................................... 57 4.2 Methods......................................................................................................................... 59 4.2.1 AGB per ground plot ................................................................................................ 59 4.2.2 Predictor variables .................................................................................................... 59 4.2.3 AGB models.............................................................................................................. 60 4.2.3.1 Models using CC only ...................................................................................... 60 4.2.3.2 Models using multiple predictor variables ........................................................ 62 4.2.3.3 Final variable selection ..................................................................................... 62 4.2.3.4 Uncertainty assessment and model selection .................................................... 63 4.3 Results ........................................................................................................................... 65 4.3.1 Models using CC only .............................................................................................. 65 4.3.2 Models using multiple predictor variables ................................................................ 67 4.3.3 Final variable and model selection ........................................................................... 68 viii  4.4 Discussion ..................................................................................................................... 72 Chapter 5: Estimating change in forest area and above-ground biomass ..............................80 5.1 Introduction ................................................................................................................... 80 5.2 Methods......................................................................................................................... 85 5.2.1 Ground plot data ....................................................................................................... 85 5.2.2 Landsat time series .................................................................................................... 86 5.2.3 Time-invariant models .............................................................................................. 88 5.2.4 Land cover change .................................................................................................... 89 5.2.5 AGB change .............................................................................................................. 93 5.2.6 Spatial context of change .......................................................................................... 94 5.3 Results ........................................................................................................................... 95 5.3.1 Time-invariant models .............................................................................................. 95 5.3.2 Land cover change .................................................................................................... 96 5.3.3 AGB change ............................................................................................................ 102 5.4 Discussion ................................................................................................................... 107 Chapter 6: Conclusion ...............................................................................................................112 6.1 Overall significance and contributions ....................................................................... 112 6.2 Limitations .................................................................................................................. 114 6.2.1 Data issues .............................................................................................................. 114 6.2.2 Modeling issues ...................................................................................................... 115 6.3 Future research directions ........................................................................................... 116 Bibliography ...............................................................................................................................118  ix  List of Tables Table 2.1. Characteristics of the four dominant vegetation types found in Nyimba District ....... 14 Table 3.1. Models to estimate tree crown radius (m). .................................................................. 22 Table 3.2. Remotely sensed data used for models to estimate canopy cover. .............................. 24 Table 3.3. Landsat 8 OLI predictor variables for estimating percent canopy cover ..................... 29 Table 3.4. RapidEye predictor variables for estimating percent canopy cover ............................ 30 Table 3.5. Disturbance predictor variables for estimating percent canopy cover ......................... 31 Table 3.6. Environmental variables for estimating percent canopy cover .................................... 33 Table 3.7. Comparison of model fit statistics from R2 versus GA variable selection................... 43 Table 3.8. Goodness-of-fit and validation statistics with predictor variables by estimation method for percent canopy cover. .............................................................................................................. 45 Table 4.1. Goodness-of-fit and validation statistics for two binomial generalized additive models to estimate percent canopy cover .................................................................................................. 67 Table 4.2. Above-ground biomass (AGB, t/ha) root mean square prediction error (RMSPE) using different predictor variable dataset combinations and models ..................................................... 67 Table 4.3. Fit and validation statistics for predicting miombo above-ground biomass (AGB, t/ha) using three models. ....................................................................................................................... 70 Table 4.4. Bootstrap results for mean and total above-ground biomass (AGB) with corresponding 95% confidence intervals .............................................................................................................. 71 Table 5.1. Land cover change status, classes and categories for trajectory category definitions representing forest/nonforest trajectories for the 2000 – 2015 time period .................................. 91 x  Table 5.2. Area of land cover categories characterizing annual forest/nonforest trajectories between 2000 and 2015, by distance to improved roads. ........................................................... 100 Table 5.3. Changes in above-ground biomass tonnes (ΔAGB t) by land cover category, class, and status. .......................................................................................................................................... 104 Table 5.4. Above-ground biomass change (ΔAGB t) by land cover class, status and distance to improved road class. ................................................................................................................... 105  xi  List of Figures Figure 1.1. Overview of dissertation structure................................................................................ 9 Figure 2.1. The miombo ecoregion in Southern and Eastern Africa with the location of Nyimba District, Eastern Province, Zambia. .............................................................................................. 11 Figure 2.2. Nyimba District, Eastern Province, Zambia. .............................................................. 11 Figure 2.3. Ground plots measured in Nyimba District. ............................................................... 16 Figure 3.1. Percent canopy cover per plot, as determined from tree point locations and crown radii models. .................................................................................................................................. 22 Figure 3.2. Coverage of RapidEye imagery over Nyimba District in 2013, with date of acquisition and percent cloud cover for each image. .................................................................... 25 Figure 3.3. Extents of Landsat 8 OLI imagery over Nyimba District with date of acquisition and Path/Row identifier. ...................................................................................................................... 26 Figure 3.4. 2013 Landsat 8 OLI mosaic covering Nyimba District.............................................. 27 Figure 3.5. 2013 RapidEye image mosaic covering Nyimba District. ......................................... 27 Figure 3.6. Observed vs predicted percent canopy cover (CC) using RapidEye imagery with environmental and disturbance variables ...................................................................................... 48 Figure 3.7. Observed vs predicted percent canopy cover (CC) using Landsat 8 OLI imagery with environmental and disturbance variables ...................................................................................... 49 Figure 3.8. Map of predicted canopy cover per 0.1 ha pixel and forest cover per 0.5 ha pixel .... 50 Figure 4.1 Scatterplot for canopy cover (CC) and above-ground biomass (AGB). ..................... 66 Figure 4.2. Observed versus estimated above-ground biomass (AGB) using predicted canopy cover, environmental, disturbance, and remotely-sensed predictor variables. ............................. 70 xii  Figure 4.3. Map of estimated above-ground biomass (AGB) per pixel (0.1 ha) and the Coefficient of Variation per pixel (CV). ....................................................................................... 72 Figure 5.1. Historical Landsat satellite image availability ........................................................... 87 Figure 5.2. Observed vs. predicted graphs from models for two forms of percent canopy cover (CCVERT, CCHEMI) and above-ground biomass (AGB). ................................................................ 96 Figure 5.3. Forest/nonforest maps for each year between 2000 and 2015 ................................... 98 Figure 5.4. Land cover class per 0.5 ha pixel in, line symbology as in Fig. 2.2. ........................ 101 Figure 5.5. (a) Close-up of deforestation by year. Landsat false color image composite (SWIR1, NIR, Red) for the year 2000 (b) and the year 2015 (c) ............................................................... 102 Figure 5.6. Above-ground biomass (AGB) change in tonnes (t) between 2000 and 2015 ......... 106 Figure 5.7. (a) Close-up of deforestation areas in relation to AGB change ................................ 107 xiii  List of Abbreviations AGB – above-ground biomass  ALS – airborne laser scanning A/R – afforestation/reforestation ARVI – Atmospherically Resistant Vegetation Index AWC – available soil water content CC – canopy cover  CI – confidence interval CO2 – carbon dioxide CTI – Compound Topographic Index CV – Coefficient of Variation  DBH – diameter at breast height DTW – Dynamic Time Warping  ETM – Enhanced Thematic Mapper EVI – Enhanced Vegetation Index FAO – Food and Agriculture Organization of the United Nations GA – genetic algorithm  GAM – generalized additive model  GLM – generalized linear model GOFC-GOLD – Global Observation of Forest and Land Cover Dynamics GPS – Global Positioning System GRZ – Government of the Republic of Zambia ha - hectare ILUA2 – Integrated Land Use Assessment 2 xiv  IPCC – Intergovernmental Panel on Climate Change IRMAD – Iteratively Reweighted Multivariate Alteration Detection kNN – k-Nearest Neighbors L1T – Level 1 Terrain Correction LEDAPS – Landsat Ecosystem Disturbance Adaptive Process MD – mean difference MODIS – Moderate Resolution Imaging Spectroradiometer MPD – mean prediction difference MSN – Most Similar Neighbor NDMI – Normalized Difference Moisture Index NDVI – Normalized Difference Vegetation Index NFI – National Forest Inventory OLI – Operational Land Imager OLS – ordinary least squares REDD+ – Reducing Emissions from Deforestation and forest Degradation  RF – Random Forests RMSE – root mean square error  RMSPE – root mean square prediction error  SAR – Synthetic Aperture Radar SWIR – Shortwave infrared t – metric tonne TM – Thematic Mapper USGS – United States Geological Survey yr – year xv  Acknowledgements The culmination of this dissertation relied on the assistance of many individuals. First, I would like to express an enormous amount of gratitude to my academic supervisor and mentor, Prof. Valerie LeMay. Her steadfast commitment to excellence in research and teaching gave me the opportunity to learn from a true leader in forest biometrics. For this, I am extremely grateful. My committee members at UBC, Prof. Peter Marshall and Prof. Nicholas Coops always provided insightful, nuanced and timely advice and support from which I profited greatly. My committee member from the Center for International Tropical Agriculture, Dr. Louis Verchot, originally introduced me to the opportunity for conducting research on forest monitoring in the miombo woodlands while we were both at the Center for International Forestry Research (CIFOR) in Indonesia. I have learned much from his insights in the areas of international tropical forest monitoring policy and tropical soils, as well as his exemplary leadership style.  I received funding for field work in Zambia and a one year stipend while at UBC from CIFOR, which was provided by the United States Agency for International Development under grant number BFS-G-11-00002, entitled the Nyimba Forest Project. Dr. Davison Gumbo of CIFOR-Zambia deserves special recognition for his foresight and perseverance in acquiring this funding to help promote improved forest monitoring tools and methods in the miombo ecoregion. Prof. Emmanuel Chidumayo helped me to gain insight into the complex ecology of miombo woodlands. His extensive research in this regard is unparalleled and I was very fortunate to work with him as a co-author. Mr. Kaala Moombe of CIFOR-Zambia helped me to understand the intricate linkages within the social aspects of forest issues in Zambia. I worked with a top notch crew while collecting data and learning about miombo woodland species identification, working xvi  in areas with dangerous wildlife, and local customs. Mr. Joel Lwambo, Mr. Martin Lyambai, and Mr. Andrew Goods Nkoma became my family for three months, sharing experiences we will never forget. Our fearless and tireless drivers, Mr. Geoffrey Tebuho and Mr. Susiku Muyapekwa delivered us safely to and from areas deep in the expanses of the miombo woodlands and I am deeply thankful for their efforts. Mrs. Bertha Kauseni, Ms. Rhoda Chiluba, and Mr. Shadreck Ngoma provided outstanding administrative and logistical support for our inventory work, without which the inventory could not have succeeded. In Nyimba District, we were warmly welcomed by Chief Nyalugwe, Chief Ndake, Chieftaness Mwape, and Chief Luembe, who all had the foresight to recognize that enhanced information on natural resources can only benefit improved management. Others in Zambia I would like to acknowledge for their contributions are Mr. Sylvester Siame, Mr. Smart Lungu, Mr. Saule Lungu, Mr. Moses Ngulube, Dr. Julian Fox, Mr. Abel Siampale, and the residents of Nyimba District, Zambia. While at The University of British Columbia, I received a Four-Year Fellowship to conduct analysis and writing. Without this support, it certainly would have taken more time to finalize this important work. While at the university, I have had the opportunity to meet many other professors, graduate students and staff and have always benefitted from the rich and stimulating environment on campus. However, Kyle Lochhead and Suborna Ahmed stand out as graduate student peers with whom I have engaged in countless discussions on forest biometrics topics. I have learned much from them and am sure I will continue to learn from them in the future.  All of the analyses in my research were conducted using the open source software platforms of R and Python. I utilized these software platforms as a way to demonstrate the utility of these powerful, yet freely available, tools for forest monitoring research and application. As such, I xvii  would like to express gratitude towards the open source statistical software development community for their ongoing efforts.  Even though I have traveled across North America, Africa, and Southeast Asia conducting forestry and natural resource work, my family has always been close, supportive, and a source of inspiration for helping to make the world a better place. Through these travels, I met the love of my life, Yulia, and she has provided never ending support, love, and encouragement; her dedication and outlook are amazing in every way and I look forward to many more adventures together. Lastly, I would like to acknowledge my son, Elijah, who already inspires me to continue the journey to help conserve our amazing global environment for future generations. xviii  Dedication  To those who seek a better world through improved knowledge and understanding of our environment, may this dissertation help in their journey.1  Chapter 1: Introduction 1.1 Background Information on forest resources in tropical developing countries is challenging to acquire due to financial and logistical constraints (Holmgren and Marklund 2007, Skutsch and Ba 2010, Stringer et al. 2012). Nevertheless, it is recognized that tropical forest resources are being exploited and converted to nonforest uses at rates that have been shown to outpace the capacity for regrowth (Shearman et al. 2012, Brink et al. 2014, Sawe et al. 2014, Suberu et al. 2014). This is especially true in the miombo ecoregion of Southern and Eastern Africa (Cabral et al. 2010, Mayes et al. 2015), where forest productivity is marginal (Frost 1996). Despite this low productivity, natural vegetation in this ecoregion is an important global biodiversity hotspot (Timberlake and Chidumayo 2011), contributes to regional watershed environmental services (Gumbo and Chidumayo 2010), and provides soil nutrient recycling functions important to agricultural production (Nyathi and Campbell 1993). The importance of forest resources in the miombo ecoregion to rural economies is also well documented (Chidumayo and Gumbo 2010, Kalaba et al. 2013, Ryan et al. 2016). Rural communities rely on forest products for both subsistence and trade (Shackleton and Gumbo 2010). Increasing population is placing greater demands on the productive capacity of the forests and on the land base in general (Hervey 2012). These pressures result in forest and land use changes that are occurring at an increasingly rapid pace (Hansen et al. 2013, Mayes et al. 2015). Understanding the distribution and amount of forest resources at local and national levels is a critical first step in being able to improve forest management (Stringer et al. 2012), for planning 2  effective policies and measures (Kleinn et al. 2010), and for engagement in emerging international initiatives such as Reducing Emissions from Deforestation and forest Degradation, or REDD+ (GOFC-GOLD 2015). Forest monitoring is a key component of developing and furthering this understanding. Forest monitoring programs can be defined as the collection and analysis of measurements over time for the purpose of estimating change in a resource (Schreuder et al. 1993). Monitoring programs can also be linked with progress towards achieving specific management goals (Elzinga et al. 2009). However, in either case there are two main informational components to monitoring programs: point-in-time estimates of forest attributes and estimates of change over time for those attributes.  Institutional efforts to improve forest monitoring in tropical countries often focus on development and/or enhancement of National Forest Inventories (NFI). Historically, information derived from NFIs has been used to support national forest policy related to timber harvesting (Köhl et al. 2006, Tomppo et al. 2010). More recently, NFIs have also been used to support forest resource decision making at sub-national levels, international reporting, and in providing information concerning the contribution of forests to the carbon cycle (Tomppo et al. 2010). In order to promote use of the information derived from NFIs by a broad range of stakeholders, NFI data collection and analysis needs to be methodologically and scientifically sound, cost-efficient, transparent, and harmonized with commonly used definitions (Köhl et al. 2006, Tomppo et al. 2010).  3  NFIs are used to describe forest populations, where the populations are all the elements over a given area. Two basic frameworks exist for making inference about a population: the design-based and the model-based framework (Gregoire 1998, Köhl et al. 2006, Kangas et al. 2006).    Within these two frameworks are approaches that utilize a mixture of sampling design and modeling, such as the design-based, model-assisted approach (Särndal et al. 1992) and the model-based, randomization-assisted approach (Särndal 2010),  also known as the hybrid approach (Ståhl et al. 2016).  The design-based framework assumes elements in a finite population have fixed quantities and the selection of elements (i.e., sampling units) in a sample is random, where each element has a positive inclusion probability (Gregoire 1998, Ståhl et al. 2016). Design-based inference is based on the variation of all possible samples of elements of a certain size that can be selected from the population when using a probability-based sample selection method (Kangas et al. 2006). Further, design-based inference makes no assumptions concerning the probability distribution or underlying structure of elements in the population (Gregoire 1998, Ståhl et al. 2016). Design-based estimators are generally unbiased, make few assumptions, and are well-studied (Schreuder et al. 1993, Köhl et al. 2006, McRoberts et al. 2010). NFIs often employ systematic selection of sampling units because of cost-efficiency and to ensure representative coverage over large areas (Köhl et al. 2006, Gregoire and Valentine 2007). Systematic sampling estimators are biased because they do not guarantee positive inclusion probabilities (Mandallaz 2007, Gregoire and Valentine 2007), and as such, simple random sampling estimators are often applied in their place (Kangas 2006). This application may be acceptable if the attributes of interest are distributed 4  randomly throughout the population (Köhl et al. 2006), in which case the variance tends to be over-estimated (Gregoire and Valentine 2007).  The model-based framework assumes that the elements in a population are random variables and that the population arises from a superpopulation that can be characterized using a model (Gregoire 1998, Ståhl et al. 2016). Model-based inference is based on the probability distribution specified by the model. A sample of elements from a population is used to estimate the model parameters; however, the method of selecting elements should not be based on the attribute of interest (Gregoire 1998). Further, the sample selection method is not a required feature of model-based inference, though it could be employed under a model-based, randomization-assisted approach (Särndal 2010). Theoretically, the estimate of a population parameter from a design-based framework is not equivalent to the estimate of a population parameter from model-based framework because of the assumption that the population is random under the latter framework (Schreuder et al. 1993).  Model-based estimators are considered model-unbiased if the true model that describes the underlying process can be identified (Köhl et al. 2006). If the true model is not identified, then precision may be optimistic (i.e., over-estimated) (Schreuder et al. 1993). When the model does hold, model-based estimators are often considered more efficient than design-based estimators in terms of precision, and useful inferences can be derived with small sample sizes (Schreuder et al. 1993). In the model-based framework, a wide definition of a model can be employed such that the model is “any assumption about the structure of the population” (Smith 1994, in Gregoire 1998).   5  Historically, NFI’s commonly employed design-based estimators (Lawrence et al. 2010, McRoberts et al. 2010); however, there has been an increasing trend towards used of model-based estimators because of efficiency advantages in terms of precision and cost (Magnussen 2015a, Ståhl et al. 2016). In the miombo ecoregion, a design-based framework has been employed for NFI sampling in several countries (e.g., GRZ 2014, Tomppo et al. 2014). However, the initial cost of implementing NFIs in tropical countries is generally high and funding is often dependent on donor support (Köhl et al. 2006, Arnold et al. 2014). Further, when change estimates are needed remeasurement costs can be prohibitive (FAO 2008). Recent work has examined a variety of design-based estimators for point-in-time characterization of forest attributes in the miombo ecoregion (Ene et al. 2016, Næsset et al. 2016). However, there is need for further research investigating the use of model-based estimators for providing both point-in-time and change estimates of forest attributes in this region. In both design-based and model-based frameworks, auxiliary information correlated with the forest attributes is commonly used in the design stage and/or the estimation stage (Schreuder et al. 1993, Köhl et al. 2006). In an NFI context, this is often referred to as multi-source information and includes information from remotely sensed data sources such as satellite imagery, aerial photography, or Airborne Laser Scanner (ALS) (Kangas et al. 2006, Tomppo et al. 2008). Environmental or disturbance data are also important components of multi-source information (Köhl et al. 2006), although these data have not been commonly employed in the miombo ecoregion. Multi-source information provides the ability to estimate forest attributes by sampling smaller areas and at a lower cost than estimation which does not use these data sources 6  (Tomppo et al. 2008, Tomppo and Tuomainen 2010). In a model-based framework, multi-source data are often used as predictor variables in models to estimate the forest attributes of interest (Köhl et al. 2006). There are a vast number of potential predictor variables that can be derived from multi-source datasets and a large number of ways in which to identify their usefulness in these models. Research is required before recommendations on data sources, predictor variables, and predictor variable selection methods can be provided for obtaining precise model-based estimates of forest attributes in the miombo ecoregion. 1.2 Research goal and dissertation overview  The main goal of this dissertation was to contribute to the development of methods and approaches to monitoring forest attributes. While the miombo ecoregion is studied in this dissertation, the methods and approaches were designed for broad application to other regions. To do this, I used a model-based framework incorporating multi-source information. The main research question I addressed was: What methods and data are most suited for estimating forest attributes at one point-in-time and change over time in the miombo ecoregion?  I focused on canopy cover (CC) and above-ground biomass (AGB) as two key forest attributes for the following reasons. CC is readily linked to forest area under commonly used forest definitions. When CC is estimated over time, forest area change can be described and this is related to deforestation and afforestation/ reforestation (A/R). Information on the amount and distribution of AGB is needed for fuelwood analyses, habitat assessments, and when estimated over time, is a main component of greenhouse gas accounting initiatives.  7  Specific objectives and research questions were identified and linked together in support of achieving the main goal (Fig. 1.1). The following overview provides an outline of the dissertation, and describes how the chapters are linked together. In Chapter 2, I describe characteristics of forest vegetation in the miombo ecoregion and the setting of the specific study area landscape. The ground-based plot measurements based on NFI protocol are also described.  In Chapter 3, I investigate models and data sources for estimating percent CC and forest area. Specifically, I compare the use of optical satellite imagery from two sensors and four modeling approaches. CC is defined as the percent of the ground covered by a vertical projection of the outermost perimeter of non-overlapping tree crowns on an NFI ground plot. The particular research questions asked were:  1. Does addition of environmental and/or disturbance predictor variables enhance model performance over using remotely-sensed data alone? 2. Does the use of high spatial resolution imagery improve model performance over medium resolution imagery? and, 3. Which modeling approach provides the best estimates?  In Chapter 4, I investigate the use of predicted CC in models to estimate AGB and compare these results with models that use multi-source information. I compared two forms of predicted CC as single predictor variables to estimate AGB. The first form of CC is the vertical projection view from Chapter 3, while the second form of CC is the result of a model to predict CC based 8  on ground plot measurements from a spherical densiometer (i.e., hemispherical projection view). The particular research questions asked were: 1. Which form of CC is more useful as a single predictor variable to estimate AGB, the hemispherical projection view or the vertical projection view? 2. Does AGB estimation accuracy improve if CC is combined with other predictor variables? 3. What AGB model performs best in terms of fit and validation statistics? and,  4. How do model-based methods compare to a simple design-based estimate?  In Chapter 5, I utilize time-invariant models, based on models from Chapters 3 and 4, to predict CC, forest area, and AGB at annual time steps using a time series of satellite images. I derived annual change in land cover from forest/nonforest trajectories and time series classification, and a stock-difference approach to estimate annual AGB change. The particular research questions asked were: 1. How much and what type of land cover change occurred?  2. When and where did land cover change occur?  3. What is the total change in AGB, the change in AGB within categories of land cover, and the annual change per hectare? and,  4. Where did AGB change occur? In Chapter 6, I provide an overall conclusion to integrate the results, highlighting contributions, strengths, limitations, and future research directions.   9   Figure 1.1. Overview of dissertation structure.   10  Chapter 2: Study area and ground data 2.1 Study area The study area landscape is Nyimba District (14°-15°S, 30°-31°E; 1,039,269 ha), Eastern Province, Zambia and located at the center of the miombo ecoregion that covers over 3.6 million km2 in Southern and Eastern Africa (Fig. 2.1, Frost 1996, Timberlake and Chidumayo 2011). Nyimba is bisected by the Luangwa River at the southwest extremity of the Great Rift Valley (Fig. 2.2, Dalal-Clayton et al. 1985). The miombo ecoregion is comprised of locally heterogeneous vegetation types, mainly characterized by tree species from the Caesalpinioideade sub-family of leguminous plants (Timberlake et al. 2010). The dominant vegetation type is miombo woodland, covering over 2.5 million km2 (Timberlake and Chidumayo 2011).  The name miombo is commonly used as a local name throughout the ecoregion and in reference to the main miombo woodland tree species (i.e., Brachystegia, Julbernardia and Isoberlinia, Timberlake et al. 2010). Other woodland types are present in the ecoregion (i.e., mopane woodland, munga woodland, Baikiaea woodland) (Timberlake and Chidumayo 2011). In general, woodlands of the miombo ecoregion are seasonally dry, deciduous or semi-deciduous forest types with highly variable CC (Frost 1996). However, semi-evergreen or evergreen forest types, such as riparian gallery forests, are also present, as well as grasslands, wetlands, and shrublands (Timberlake and Chidumayo 2011). The distribution of these vegetation types generally depends on rainfall, soil characteristics, and level of historic disturbance (Frost 1996). Soils in the miombo ecoregion are generally nutrient poor (Frost 1996).  11   Figure 2.1. The miombo ecoregion in Southern and Eastern Africa with the location of Nyimba District, Eastern Province, Zambia.  Figure 2.2. Nyimba District, Eastern Province, Zambia.  12  Miombo woodlands are generally divided into two distinct forest types: dry miombo woodland and wet miombo woodland (Timberlake and Chidumayo 2011). Wet miombo woodland generally receives more than 1,000 mm annual rainfall per year and has an average tree canopy height greater than 15 m, whereas dry miombo woodland generally receives less than 1,000 mm rainfall per year and has a lower average tree canopy height. Nyimba District receives an annual rainfall of 600-900 mm per year and most rain falls in the December to February period, with a distinct dry season from May to November (Timberlake et al. 2010). Dominant vegetation types found in Nyimba District include: dry miombo woodland, mopane woodland, munga woodland, and riparian forest (Table 2.1, GRZ 1976, FAO 2014a). However, other land cover types are also found, including agriculture, grassland, and settlements (GRZ 1976, Gumbo et al. 2016). The soils in Nyimba District are characterized as Lithosol-Cambisols on the hills and plateaus and Fluvisol-Vertisols in the valleys (Hengl et al. 2015). Elevation ranges from 450 m along the Luangwa River valley bottom to 1,000 m on the plateau, and reaches over 1,500 m on the mountain tops.  Most of the human population is located in the vicinity of the district capital (also called Nyimba), the sub-district capitals, and along the Great East Road, which is the national highway that connects Nyimba to the national capitol, Lusaka and to the Eastern Province capital, Chipata. West of the Luangwa River is the Petauke Game Management Area (GMA), where population density is very low. In 2010, the population of Nyimba District was 81,025, with a density of 8.1 people per km2 (CSO 2010). Population growth in the District between 2000 and 2010 was lower than provincial and national averages (1.9% versus 2.6% and 2.8%, respectively) (CSO 2010). 13  Explicit research on forest area change in Nyimba District has not been conducted; however, local evidence indicates that conversion of forest to nonforest is increasing due to agricultural expansion and tree harvesting (Gumbo et al. 2013, Watson et al. 2015, Gumbo et al. 2016). Agricultural expansion is associated with subsistence agriculture and cash cropping (Gumbo et al. 2016). Tree harvesting occurs for both fuel (charcoal) and timber. Charcoal use is balanced between consumption within the district and export to Lusaka and Chipata, the provincial capital (Gumbo et al. 2013). Tree harvesting for timber is only done on a small scale for producing planks to export out of the district (Gumbo et al. 2013).  Fire is a frequent disturbance, with an estimated return interval of one to three years (Frost 1996, Timberlake and Chidumayo 2011). Fires generally result from people engaged in agricultural land clearing, charcoal making, and hunting, with early dry-season fires exhibiting less intensity than late dry-season fires (Frost 1996).   14  Table 2.1. Characteristics of the four dominant vegetation types found in Nyimba District , Eastern Province, Zambia. Vegetation type Dominant species Structure Phenology Fire susceptibility of dominant species Soils and topography Source Dry miombo woodland Brachystegia spiciformis, B. boehmii and Julbernardia globiflora Open to closed two layer canopy, with mean height < 15 m; dense grass layer deciduous or semi-deciduous  fire-tolerant Nutrient-poor, plateaus and hills Timberlake and Chidumayo (2011), GRZ (2014) Mopane woodland Colophospermum mopane (generally monospecific) Open one layer canopy, with  height ranging from 6 – 18 m; sparse grass layer deciduous  fire-intolerant Clay dominated, on wide valley bottoms Timberlake 1995, Makhado et al. (2014), GRZ (2014) Munga woodland  Vechellia sp., Senegalia sp.,  Combretum sp., and trees associated with the Papilionoideae subfamily Open one or two layer canopy, emergents up to 18 m; dense grass layer deciduous  fire-tolerant Nutrient-rich and well-drained, plateaus Timberlake and Chidumayo (2011), GRZ (2014) Riparian forest Mixed Closed three layer canopy, up to 25 m; vines commonly occurring evergreen fire-intolerant Restricted to buffer zones around significant rivers GRZ (2014) 15  2.2 Ground-based plot measurements Field measurements were collected at 64, 0.1 ha permanent plots in Nyimba District from May to December, 2013 (Fig. 2.3). The plots are located on a systematic 10 km x 10 km grid corresponding to the Zambian Integrated Land Use Assessment 2 (ILUA2), a NFI program implemented by the Zambian Forest Department in collaboration with the FAO (GRZ 2014).  For ILUA2, a cluster of four 20 m x 50 m plots are located at approximately 20% of grid intersections across the country. Each cluster plot has four 20 m x 50 m rectangular subplots, and each subplot has a nested 20 m x 10 m microplot starting at the same origin as the subplot. In this study, ground plot data were collected according to Zambian NFI protocol only for the first subplot of each cluster to increase the spatial coverage, given the constraints of available field measurement time. As such, the 20 m x 50 m area where ground measurements were collected is referred to as the plot, and the 20 m x 10 m area as the subplot.  For each live tree with a diameter at breast height (DBH; outside bark diameter at 1.3 m above ground) ≥ 10 cm within the 20 x 50 m plot boundary, the following information was recorded: species, DBH (cm); total height (m); and the x, y location (m) at the center of the tree base relative to the plot starting point. At the plot level, the dominant vegetation type (i.e., miombo, mopane, munga, or riparian) was recorded; however, ‘nonforest’ was recorded if the ground plot was dominated by nonforest land use with none of the vegetation types present. CC measurements were collected using a spherical densiometer according to the NFI protocol. Ground plot coordinates were determined using a Trimble GeoXT Global Positioning System (GPS), which were differentially corrected (+/- 1 m).  Lastly, GPS line data was collected for roads, and point data for villages. Most of the ground plots (49) were visited between 5 May and 16  30 July 2013, while the remaining plots (15) were visited between 27 October and 6 December 2013.  Figure 2.3. Ground plots measured in Nyimba District.    17  Chapter 3: Canopy cover and forest area modeling 3.1  Introduction The world’s forests provide a number of socio-economic benefits to humans (FAO 2014b) and are critical to both biodiversity conservation and ecosystem services (Harrison et al. 2014). Consequently, there are many programs to monitor the state of forests (Hansen et al. 2010), with the simplest attribute being the forest area. However, having a single, world-wide definition of forest versus non-forest is difficult as it may not align with current standards (FAO 2010). If forest is considered as areas with only a few trees, then a number of cities would be considered as forest. Conversely, under other definitions, open, treed areas would be considered as non-forest, potentially leaving entire countries with no forested area. Therefore, global efforts which adopt a universal definition of forest area (e.g., Hansen et al. 2010) may not meet the needs of practitioners and policy makers at national levels. The miombo ecoregion of Southern Africa is a prime example where global definitions may not apply and where forest information is lacking (Mayes et al. 2015). This ecoregion is a diverse assemblage of tropical dry forest types and land uses (Frost 1996).  A number of environmental and disturbance factors interact contributing to this diversity including: human disturbances, wildlife impacts, fire, rainfall, and site characteristics such as geomorphology and soil (Frost 1996, Sankaran et al. 2008, Ahrends et al. 2010, Timberlake et al. 2010).  These woodlands are home to natural resource-dependent communities (Turpie et al. 2015) and contain highly diverse assemblages of megafauna (Timberlake et al. 2010).  Deforestation appears to be increasing at a rapid rate across the ecoregion (Cabral et al. 2010, Mayes et al. 2015).  Because of the 18  importance of these woodlands and the threats to resource values due to increased deforestation (Turpie et al. 2015), timely monitoring of forest area is essential for land managers and is key for implementing programs such as Reducing Emissions from Deforestation and forest Degradation (REDD+, FAO 2014b).   Since forest area definitions may differ over agencies and over time, monitoring percent CC is a better alternative as any definition of forest/non-forest can be later applied.  CC information can be used for a number of purposes including:  1) as an input variable to vegetation dynamics models to forecast productivity (Desanker and Prentice 1994) and ecosystem process models (Caylor and Shugart 2004); 2) to provide insights into habitat integrity needed for conservation planning (Pintea et al. 2003); 3) as a predictor in estimating forest biomass (Wu et al. 2013); and 4) to analyze the cost and benefits of ecosystem services  (Turpie et al. 2015).   Overall, percent CC as a forest attribute offers a number of benefits over identifying an area only as forest or non-forest.   However, measuring percent CC over a land area is difficult using forest inventory information.  In forest inventories, percent CC is commonly defined as the percent of ground covered by the vertical projection of tree crowns, and has also been labelled as canopy closure, crown cover, or crown closure (Jennings et al. 1999, FAO 2010). Percent CC has traditionally been measured on stereo-pairs of aerial photographs that cover a land area of interest. However, since aerial photographs are expensive to acquire, stands may be large, and the definition of stand boundaries can change over time, monitoring CC is often cost-prohibitive using this traditional method. One alternative is to process Airborne Laser Scanning (ALS) data (Urbazaev et al. 2015) to measure 19  percent CC at a pixel level. However, ALS may not provide a cost-efficient source of repeated measurements needed for long-term monitoring and data are not available for historical analysis. Another alternative for monitoring is to measure CC on field plots and then to use remotely sensed data and other information to estimate CC for non-sampled areas.  Although radar data has been used for this purpose (e.g. Naidoo et al. 2015, Urbazaev et al. 2015), radar has limitations which include lack of systematic availability, cost of acquisition, challenging processing requirements, and complex interactions with forest structure that are band-dependent (Sinha et al. 2015).   Passive optical satellite imagery will likely remain a principle data source for developing predictive models due to availability, repeatability, cost, and well-established processing methods (Kennedy et al. 2014). Assessments of CC in tropical dry forests of Africa have been conducted using medium resolution Landsat satellite imagery through a variety of methods. Traditional approaches delineated CC classes using visual interpretation (Mukosha and Siampale 2009) or supervised classification (Cabral et al. 2010). However, modeling methods such as parametric models (Wu et al. 2013) and regression trees (Karlson et al. 2015) can be used to estimate CC as a continuous variable rather than as a class. Using modeling approaches, Root Mean Square Errors (RMSE) ranged between 7 and 13 percent CC across several studies (Hansen et al. 2002, Adjorlolo and Mutanga 2013, Wu et al. 2013, Karlson et al. 2015). Use of high resolution satellite imagery has been suggested to reduce these error rates by providing spectral reflectance values at a finer spatial resolution (Hojas-Gascon et al. 2015). In this context, RapidEye satellite imagery has shown promise for CC estimation given the 5 m spatial resolution and the addition of a red-edge band. Ozdemir (2014) used RapidEye to develop different versions of parametric models for CC 20  in subtropical dry forests of Turkey, achieving RMSE of 7 – 10%. Similar research efforts to estimate percent CC in tropical, dry African forests are not often linked to plot-based measurements from permanent monitoring programs such as a National Forest Inventory (NFI) (Hansen et al. 2002, Bucini et al. 2010, Adjorlolo and Mutanga 2013, Wu et al. 2013, Karlson et al. 2015, Naidoo et al. 2015, Urbazaev et al. 2015). Further, few studies compare a range of methods for estimating CC, and forest area, in the miombo woodlands using optical remotely sensed data.  The main objective in this chapter was to recommend an accurate and cost-efficient method for estimating forest area, in the context of the existing Zambian NFI (GRZ 2013). This chapter compares the use of optical satellite imagery from two sensors and four modeling approaches for estimating percent CC in the miombo ecoregion of Zambia.  Percent CC was defined as the percent of the ground covered by a vertical projection of the outermost perimeter of non-overlapping tree crowns, aligning with the national definition of forest area in Zambia (FAO 2010, 2014a). RapidEye and Landsat 8 OLI imagery were compared using parametric (Generalized Linear Model), non-parametric (k-Nearest Neighbor imputation; Random Forests regression trees), and semiparametric (Generalized Additive Model) approaches.  Readily-available auxiliary variables that characterize disturbance and environmental factors were also examined to assess improvement in results, including: fire occurrence, distance to population, slope, and soil variables. The following specific questions were addressed: 1) Do auxiliary variables enhance model performance over using remotely sensed data alone?; 2) Does the use of high spatial resolution imagery (RapidEye) improve percent CC estimation over medium resolution imagery (Landsat 8 OLI)?; and 3) Which modeling approach provides the best 21  estimates? One approach was selected and used to estimate CC by pixel for the study area.  Forest/non-forest was then mapped based on the percent CC predictions and total forest area was estimated with a 95% confidence interval.  3.2 Methods 3.2.1 Canopy cover per ground plot Percent CC was calculated for each plot using the DBH measures and existing crown radius models following approaches used by Gill et al. (2000) and by Westfall and Morin (2013).  For this purpose, species-specific models were used when possible (Table 3.1).  Since trees < 10.0 cm DBH were measured only on a 20 m x 10 m subplot according to Zambian NFI protocol, crowns associated with these smaller trees were not included.  Also, the Burrows and Strang (1964) models are based on basal diameter instead of DBH. However, these authors noted that there is little taper in the three species with which they worked. Sixty of the 64 plots fell within the miombo or mopane woodland type, whereas the four remaining plots fell into munga or riparian woodland types (less than 5% of the tree data).  Given no available crown width or crown area models for these types, the miombo woodland type model was used. For each tree record, the tree location x,y points were buffered by the predicted radius from the tree species or woodland type crown radii model (Fig. 3.1). Overlapping buffer polygons were merged in the process. Finally, the canopy area polygons were clipped to the NFI ground plot boundaries to arrive at CC per plot.  22  Table 3.1. Models to estimate tree crown radius (m).  Tree species or woodland type Model1 Source Location Mean annual rainfall (mm) Brachystegia spiciformis   Burrows and Strang 1964 near Harare, Zimbabwe 825  Julbernardia globiflora     Parinari curatellifolia  Miombo      Isango 2007 Iringa District, Tanzania 565-900 Mopane  Mean crown radius: 1.58 m Fuller et al. 1997 Eastern Province, Zambia 750-950 1 cr: crown radius (m); bd: basal diameter (m); ht: height (m); DBH (m)   Figure 3.1. Percent canopy cover per plot, as determined from tree point locations and crown radii models. cr̂ = (-0.302  + 40.703 × bd        − 26.220 × bd2) / 2   cr̂   = (-0.258              + 35.266 × bd) / 2     cr̂= 0.073 + 0.113 ×       DBH + 0.136 × ht         cr̂= (0.324 + 22.931 × bd) / 2          23  3.2.2 Satellite imagery pre-processing Optical satellite imagery was acquired over the study area for June to July 2013 corresponding to the early dry season when most of the field measurements were collected (Table 3.2). This time period was identified as preferred for satellite image data acquisition in the miombo ecoregion (Fuller et al. 1997, Mayes et al. 2015).  Each RapidEye scene (Fig 3.2) was converted to surface reflectance using the Fast Line-of-sight Atmospheric Analysis of Hypercubes (FLAASH) algorithm in ENVI 5.1 (Adler-Golden et al. 1999). Some of the RapidEye scenes contained haze contamination; therefore, radiometric normalization to cloud- and haze-free scenes was implemented, utilizing the Iteratively Reweighted Multivariate Alteration Detection algorithm (IRMAD) as detailed in Canty (2014). Landsat 8 OLI data (Fig. 3.3) from the Landsat Ecosystem Disturbance Adaptive Process (LEDAPS) archive was downloaded. LEDAPS data is already processed to surface reflectance (Masek et al. 2006); therefore, radiometric normalization was not performed as both scenes were haze and cloud free. One image mosaic was created for each satellite image data source. After visual inspection of the Landsat 8 OLI mosaic (Fig. 3.4), it was found to be in agreement with the road network data collected during the field campaign and was not further processed. However, the RapidEye mosaic (Fig. 3.5) needed to be geometrically corrected and registered to the road network. Ground control points (112) were delineated on the RapidEye mosaic and a 4th order nearest neighbor polynomial transformation was implemented. The resulting Root Mean Square Error (RMSE) was 2.97 m. Any remaining clouds and cloud shadows were manually delineated from the RapidEye mosaic using on-screen digitizing, equal to 2.5% of Nyimba 24  District. These areas were removed from further analyses. The gaps did not cover the ground plot locations; therefore, the gaps were not filled in with other remotely sensed data. Table 3.2. Remotely sensed data used for models to estimate canopy cover. Satellite image  data source Data description Collection time frame Level 3A RapidEye (RapidEye 2012)   Pixel size: 5 m Image size: 25 km2 Mosaic of nine full scenes   and 23 partial scenes Final cloud coverage: 2.5% % coverage by date: June 21, 2013:    6% June 29, 2013:  25% July    3,  2013:    7% July   21, 2013:  61% Landsat 8 OLI  Pixel size: 30 m Image size: 185 km2 Mosaic of two partial scenes Path 170, Row 70 Path 171, Row 70 Final cloud coverage: < 1% % coverage by date: June    6, 2013: 80% June  29, 2013: 20%  25   Figure 3.2. Coverage of RapidEye imagery over Nyimba District in 2013, with date of acquisition and percent cloud cover for each image. Line symbology as in Fig. 2.2. 26   Figure 3.3. Extents of Landsat 8 OLI imagery over Nyimba District with date of acquisition and Path/Row identifier.  Line symbology as in Fig. 2.2.  27   Figure 3.4. 2013 Landsat 8 OLI mosaic covering Nyimba District.  Band combination is 6 (near-infrared), 5 (red), 4 (green).  Figure 3.5. 2013 RapidEye image mosaic covering Nyimba District.  Band combination is 5 (near-infrared), 3 (red), 2 (green).   28  3.2.3 Predictor variables The use of predictor variables developed from optical satellite imagery in models to estimate forest attributes has a long history and many combinations of band values and their transformations exist (Lu et al. 2014). RapidEye captures spectral information for five optical bands in the 0.440-0.850 μm range, including a unique red-edge band which was found useful for detecting stress in temperate dry forest trees (Eitel et al. 2011). Landsat 8 OLI includes multispectral information for seven optical bands in the 0.433 -2.300 µm range; however, the ‘ultra-blue’ band was not included as its primary use is for collecting information in marine and coastal environments (Irons et al. 2012).  For Landsat 8 OLI, a raster dataset was created for each of 24 predictor variables (Table 3.3). To account for Landsat 8 OLI’s spatial mismatch between the plot size and image pixel size, the pixel data was resampled to a 20 m x 50 m pixel size as recommended by Kattenborn et al. (2015) before calculating each variable. Resampling was done with a ‘nearest neighbor’ technique and not interpolation, so as to avoid changing the original pixel values when calculating indices. For RapidEye, a raster dataset was created for each of 33 predictor variables using a 4 pixel x 10 pixel moving window, calculating means and standard deviations (Table 3.4). In this process, each resultant RapidEye dataset was resampled to a 20 m x 50 m pixel size that corresponded with the plot size. With RapidEye-based and Landsat 8 OLI-based raster datasets now all at the same spatial resolution (20 m x 50 m pixel size), the variable Texture 2 was calculated as the standard deviation of a 3 x 3 window around each pixel (however, not on RapidEye Texture 1 datasets). Pixel values were extracted for each raster dataset at each plot, where the pixel center was closest to the plot center. There are more predictors derived from 29  RapidEye due to the use of within-plot textural measures (i.e., Texture 1 standard deviations) for each band and index.    Table 3.3. Landsat 8 OLI predictor variables for estimating percent canopy cover at the plot level.  Predictor variable Description Minimum Maximum Mean Citation Albedo Sum of all bands 0.4156 1.1433 0.7103 Lu et al. (2004)       Albedo.texture2*  0.0112 0.6469 0.0661 Hansen et al. (2002)       Atmospherically Resistant Vegetation Index (ARVI) (NIR – (2 x Red – Blue)) / (NIR + (2 x Red – Blue)) -0.1329 0.1974 0.0494 Kaufman and Tanre (1992)       ARVI.texture2  0.0039 0.1675 0.0283        Blue band  0.0115 0.0623 0.0324        Blue.texture2  0.0008 0.0556 0.0049        Enhanced Vegetation Index (EVI) 2.5 x ((NIR – Red) / (NIR + (6 x Red) –  (7.5 x Blue) + 1) 0.9530 1.6688 1.4284 Justice et al. (1998)       EVI.texture2  0.0079 0.1875 0.0545        Green band  0.0297 0.0966 0.0547        Green.texture2  0.0006 0.0786 0.0064        Normalized Burn Ratio (NBR) NIR – SWIR2 / NIR + SWIR2 -0.1119 0.5962 0.3084 Key and Benson (2006)       NBR.texture2  0.0054 0.3623 0.0561        Normalized Difference Moisture Index (NDMI) NIR – SWIR1 / NIR + SWIR1 -0.2136 0.3247 0.0475 Jin and Sader (2005)       NDMI.texture2  0.0047 0.2289 0.0419        Normalized Difference Vegetation Index (NDVI) NIR – Red / NIR + Red 0.3133 0.8169 0.5398 Rouse et al. (1974)       NDVI.texture2  0.0091 0.3330 0.0482        NIR band Near infrared 0.1187 0.2974 0.2263        NIR.texture2  0.0034 0.0455 0.0152        Red band  0.0247 0.1368 0.0685        Red.texture2  0.0016 0.1150 0.0106        SWIR1 band Shortwave infrared 1  0.1071 0.3333 0.2073        SWIR1.texture2  0.0028 0.1624 0.0205        SWIR2 band Shortwave infrared 2  0.0497 0.2527 0.1212        SWIR2.texture2   0.0021 0.1910 0.0156   * All variables denoted with .texture2 were calculated using the standard deviation for each band or index using a 3x3 pixel window.  30  Table 3.4. RapidEye predictor variables for estimating percent canopy cover at the plot level. Predictor variable Description Minimum Maximum Mean Citation Albedo Sum of all bands 0.2448 0.9403 0.5559 Lu et al. (2004) Albedo.texture2*  0.0117 0.6245 0.0682 Hansen et al. (2002) Albedo.texture1*  0.0236 0.5203 0.0873  Atmospherically Resistant Vegetation Index (ARVI) (NIR – (2 x Red – Blue)) /( NIR + (2 x Red – Blue)) -0.3396 0.5198 0.0390 Kaufman and Tanre (1992) ARVI.texture2  0.0067 0.5723 0.0881  ARVI.texture1  0.0322 0.4628 0.0999  Blue band  0.0026 0.0661 0.0207  Blue.texture2  0.0010 0.0735 0.0060  Blue.texture1  0.0037 0.0590 0.0091  Enhanced Vegetation Index (EVI) 2.5 x ((NIR – Red) / (NIR + (6 x Red) –  (7.5 x Blue) + 1) 0.0217 0.3834 0.1937 Justice et al. (1998) EVI.texture2  0.0054 0.2086 0.0308  EVI.texture1  0.0097 0.1287 0.0351  Green band  0.0185 0.1382 0.0605  Green.texture2  0.0018 0.1373 0.0099  Green.texture1  0.0043 0.1058 0.0135  NDRE Red edge – Red / Red edge + Red 0.0149 0.5356 0.2268 Eitel et al. (2011) NDRE.texture2  0.0039 0.2185 0.0350  NDRE.texture1  0.0172 0.1515 0.0460  NDREGI Red edge – Green / Red edge + Green 0.2382 0.5731 0.4193 Rana et al. (2014) NDREGI.texture2  0.0093 0.3259 0.0343  NDREGI.texture1  0.0208 0.2461 0.0617  Normalized Difference Vegetation Index (NDVI) NIR – Red / NIR + Red 0.0789 0.7239 0.3978 Rouse et al. (1974) NDVI.texture2  0.0048 0.3722 0.0615  NDVI.texture1  0.0256 0.3003 0.0697  NIR band Near infrared 0.0843 0.3374 0.2290  NIR.texture2  0.0052 0.0828 0.0220  NIR.texture1  0.0053 0.0930 0.0290  Red band  0.0282 0.1718 0.0998  Red.texture2  0.0040 0.1753 0.0172  Red.texture1  0.0061 0.1374 0.0214  Red edge band  0.0464 0.2502 0.1459  Red edge.texture2  0.0031 0.1604 0.0186  Red edge.texture1   0.0061 0.1287 0.0238    * All variables denoted with .texture1 were calculated using the standard deviation for each band or index within a 4x10 pixel window. All variables denoted with .texture2 were calculated using the standard deviation for each band or index using a 3x3 pixel window, where pixels were resampled to a 20m x 50m pixel size.  31  In addition to predictor variables derived from remotely sensed imagery, other auxiliary variables representing disturbance factors that potentially affect CC were included (Table 3.5).  Distances to roads and villages were determined using the field GPS data, since these distances have been shown to indirectly impact forest structure (Helmer et al. 2008, Ahrends et al. 2010). Missing road and village data was checked by overlaying the lines and points onto the RapidEye image mosaic, after spatial registration. Missing road or village locations were then digitized on-screen. Table 3.5. Disturbance predictor variables for estimating percent canopy cover at the plot level.  Predictor variable Description Minimum Maximum Mean Distance from ground plot to district capital (D2cap) Euclidean distance (m)  5,704 90,146 40,635 Distance from ground plot to nearest village or settlement (D2vill) Euclidean distance (m)  542 42,902 8,948 Distance from ground plot to nearest improved road (D2road) Euclidean distance (m)  192 59,363 12,042 Fire history – early season (Fire.early) Total number of fires per pixel, Apr. 1 – Jul. 31  0.0 6.0 0.6 Fire history – late season (Fire.late) Total number of fires per pixel, Aug. 1 – Nov. 30  0.0 8.0 1.8  Monthly fire history data was included at 500 m pixel size for 2003 – 2013 (Table 3.5), derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) Burned Area Product (Boschetti et al. 2013), since both the frequency and intensity of fires impact forest structure (Frost 1996). Low intensity burns are more characteristic of the early dry season when grass is still curing, while higher intensity burns are more likely to occur in the late dry season. The fire history data was processed into two raster datasets: 1) sum of fire occurrences per pixel between 32  April 1 and July 31 representing frequency of low intensity, early dry season fire; and, 2) sum of fire occurrences per pixel between August 1 and November 30 representing frequency of high intensity, late dry season fire. Lastly, environmental variables were included that may affect CC (Table 3.6). Miombo woodland attributes may be influenced by proximity to perennial water courses (Timberlake et al. 2010). The compound topographic index (CTI) can be used to characterize soil moisture given the effects of slope and upslope catchment area (Moore et al. 1993).  CTI and distance to perennial rivers/streams were calculated based on 30m Shuttle Radar Topography Mission Digital Elevation Model (USGS 2004). Soils information is considered to be a strong determinant influencing vegetation composition and structure in miombo woodland (Frost 1996, Timberlake et al. 2010). ISRIC – World Soil Information is a consistent and comprehensive data source for soils information covering the continent of Africa (Hengl et al. 2015). This raster soils data (250 m pixel size) is provided as estimates at six depths (2.5 cm, 10 cm, 22.5 cm, 45 cm, 80 cm, and 150 cm). These data were downloaded, clipped to the study area boundary, and means across all depths were calculated for each variable.  All disturbance and environmental raster datasets were also resampled to a 20 m x 50 m pixel size. The value of each predictor variable was assigned to a plot where the pixel center was closest to the plot center. Before selecting pixels to match with field plots for all predictor datasets, major rivers and paved roads were delineated by on-screen digitizing using a combination of the Landsat 8 OLI and RapidEye data. These areas were not subjected to any further analyses.  In a final step, each raster dataset was clipped to the Nyimba District boundary. 33  Preparation of spatial predictor variables was conducted using the R raster package (Hijmans et al. 2015) and ArcGIS v. 10.1 (ESRI 2012). Table 3.6. Environmental variables for estimating percent canopy cover at the plot level.  Predictor variable Description Minimum Maximum Mean Available soil water capacity (AWC) Volumetric fraction 12.1 15.2 14.0 Bulk density (Bulk) kg/m3 1,299.6 1,501.2 1,378.9 Carbon Organic soil carbon content (g/kg) 4.3 11.5 7.0 Cation Exchange Capacity (CEC) cmolc/kg 7.7 21.8 12.3 Clay Percent volume 21.5 35.1 28.7 Coarse fragments (Crs.frg) Percent volume 2.2 17.0 9.5 Compound Topographic Index (CTI) Unitless. High values indicate higher potential soil moisture. 5.8 39.6 8.2 Distance from ground plot to nearest perennial river/stream (D2hydro) Euclidean distance (m)  91 26,194 7,203 pH Unitless 59.0 67.8 62.6 Sand Percent volume 43.4 66.2 54.3 Silt Percent volume 9.9 22.8 17.0 Slope Percent  0.0 43.8 8.4  3.2.4 Canopy cover models 3.2.4.1 Effects of auxiliary variables To evaluate whether auxiliary variables improved model performance, models were built with remotely sensed variables only (Landsat 8 OLI or RapidEye) and then with auxiliary variables. Because of the large number of possible predictor variables, variable selection methods were used to find the best models. Stepwise variable selection methods have been widely employed for models using remote sensing data, particularly for linear models. However, issues with using 34  stepwise methods include:  i) the number and order of predictor variables entering the model affect the final choice; and ii) although the fit of a model will fail if predictor variables are multicollinear, the algorithms will continue if sets of variables are nearly multicollinear causing instability in parameter estimates (Pope and Webster 1972). This latter issue can be a problem for remotely sensed data particularly, where a number of variables are derived from the sensor band values.  Instead, two alternative variable selection methods were used: the R2 method that examines all possible combinations of predictor variables and the genetic algorithm (GA) that searches for the optimal set of predictor variables.   Each of these was applied to models with and without auxiliary variables.   For parsimony and to reduce the possibility of near multicollinearity, only models with up to six variables were considered.   R2 selection The ‘best subsets’ R2 selection method in the R leaps package (Lumley 2009) was used to calculate the coefficient of determination (R2) values for linear models fitted using Ordinary Least Squares (OLS) and the entire dataset.  The objective of the R2 selection method is to determine the variables that maximize R2 values given a fixed number of predictor variables (Morgan and Tatar 1972). For these models, percent CC was converted to a proportion (0, 1) and a logit transformation was applied.  Because the logit transformation is not defined for 0, a small adjustment factor of 0.005 was included for two plots which had 0 cover.  For each of the four variable sets (i.e., Landsat 8 OLI or RapidEye, with or without auxiliary variables), one model was selected with the highest R2.    35  Genetic Algorithm The GA can be applied to a wide variety of modeling methods and has shown promise for variable reduction in high-dimensional settings (Garcia-Gutierrez et al. 2014). It has near optimality in identifying predictor variables when using remotely sensed data to estimate forest attributes (Tomppo and Halme 2004, Tomppo 2006). Since using GAs in remote sensing applications is still relatively new, a brief summary is given here.   GAs were developed using the principle of survival of the fittest and borrow ideas from biological evolution (Scrucca 2013).   First, an individual is a particular model (i.e., a particular subset of predictor variables) out of a population of all possible models. The model is described by a binary string of 1’s (predictor variables in the model) and 0’s (predictor variables not in the model), described as a chromosome. The goal is to identify the best performing individual (i.e., the best chromosome) given a goodness-of-fit (called fitness) criterion. GAs use mechanisms of selection, crossover, and mutation through multiple generations. The initial selection of variables is random, whereas subsequent selections are based on results of “mating” the two best performing chromosomes in a given “generation”.  Cross-over of the two best chromosomes swaps genes between them, allowing subsequent chromosomes, known as children, to inherit the best traits (i.e., predictor variables) from their parents. Random noise created through mutation aids in avoiding local optima. This process is implemented until a user-defined stopping criterion is achieved.  For initial GA selection, a generalized linear model (GLM) fitted using maximum likelihood with a normal distribution, logit link on CC as a proportion, and the default variance function (i.e., equal variances) was used (McCullagh and Nelder 1989). This model is asymptotically equivalent to that used for the R2 selection (i.e., OLS with a logit transformation).  Within the 36  GA process, a key decision lies in the choice of an objective criterion, often termed fitness. The fitness was set to minimize the 20-fold cross-validated Root Mean Square Prediction Error (RMSPE) for each model under consideration, calculated as:  RMSPE = √∑(𝑦𝑖𝑣 −  ?̂?𝑖𝑣)2𝑛𝑛𝑖=1 ,   (3.1)  where 𝑦𝑖𝑣 is the percent cover of a plot not used for model building in a fold, ?̂?𝑖𝑣 is its corresponding prediction, and n = 64 field plots. When applied in the variable selection process, cross-validation has been reported to reduce the risk of overfitting (Hawkins 2004, Packalén et al. 2012). Leave-One-Out Cross-Validation (LOOCV) has been commonly applied in the remote sensing literature (Nelson et al. 2007, Bouvier et al. 2015, Vega et al. 2016). LOOCV is known to have a high computational burden (Hastie et al. 2009); therefore, 20 fold cross-validation was implemented as a compromise. After initial tests, suitable stopping criteria were established at 20 consecutive generations with no improvement in fitness, or a maximum of 100 generations. The single-point crossover probability between pairs of chromosomes was set to 0.8 and the uniform random mutation probability in a parent chromosome was set to 0.4 based on sensitivity tests.  Tests were also conducted to determine whether shifting these values impacted the GA implementation or run time; however no improvements were noticed.  The GA was implemented 100 unique times for models with remote sensing variables alone (Landsat 8 OLI or RapidEye) and then including the auxiliary variables. The GA package in R was used to conduct the GA predictor variable selection analysis (Scrucca 2013).  37  Comparisons Since the two selection methods use different fit criteria, a 20-fold RMSPE (Eq. 3.1) was also calculated for final models using the R2 selection method. Similarly, a Pseudo-R2 fit statistic (Schabenberger and Pierce 2002, p. 226) for percent CC was calculated for all models:   Pseudo-𝑅2  = 1 −∑ (𝑦𝑖 − ?̂?𝑖)2𝑛𝑖=1∑ (𝑦𝑖 −  ?̅?)2   𝑛𝑖=1 ,   (3.2)  where ?̂?𝑖 is the predicted percent cover based on the model using all observations. The RMSPE and Pseudo-R2 were compared for models using Landsat 8 OLI or RapidEye, with and without the auxiliary variables.   3.2.4.2 Modeling methods Parametric: Generalized linear models Linear models with one random effect (i.e., the error term) fitted using OLS have been used for many years to estimate forest attributes using remotely sensed data (Fuller et al. 1997).   For percent CC, a more suitable approach is to use a GLM (McCullagh and Nelder 1989), specified as:   𝑔(µ𝑖) = 𝑙𝑜𝑔𝑒 (µ𝑖1 − µ𝑖) =  𝛽0 + 𝛽1𝑥1𝑖 + 𝛽2𝑥2𝑖 … + 𝛽𝑘𝑥𝑘𝑖, (3.3)  where 𝑔(µ𝑖) is the link function used to transform the dependent variable.  For CC expressed as a proportion, the binomial distribution logit link, and default variance function (i.e., variances change with the predicted values) were used (GLM-B; Heiskanen and Kivinen, 2008). Since CC 38  as a proportion may be approximately normally distributed, a GLM with the normal distribution, logit link, and default variance function (i.e., equal variances) was also fitted (GLM-N). The GA was used to select variables for both GLM-N and GLM-B. All models were restricted to predictor subsets of six or fewer variables. Nonparametric: k-Nearest Neighbor imputation  k-Nearest Neighbor imputation (kNN) uses distances in variable space between records in target (i.e., dependent variable is not known) versus reference (i.e., dependent variable is measured) datasets to find the best matches.   The dependent variable for the selected reference records are then imputed to the target records (Crookston and Finley 2008). Since there is no probability distribution used in implementing kNN, this method is distribution-free and is commonly labelled as “nonparametric” in application literature (Eskelson et al. 2009). However, kNN can also be considered as a model-based method (Magnussen et al. 2009).  For this study, the Most Similar Neighbor (MSN) distance metric was used and calculated as:  𝑑𝑖𝑗2 =  (𝑥𝑖 −  𝑥𝑗)ΓΛ2Γ(𝑥𝑖 −  𝑥𝑗)′,    (3.4)  where xi is a (1 x p) vector of predictor x variables for the ith target record, xj is a (1 x p) vector of x variables for the jth reference record, 𝚪 is the matrix of canonical coefficients of predictor variables and 𝚲 is the diagonal matrix of the squared canonical correlations (Moeur and Stage 1995).  Logit-transformed CC as a proportion was the single dependent variable in this kNN imputation.  As a result, the canonical coefficients are the OLS coefficients between the logit of CC and the predictor variables, and the squared canonical correlation is the R2 from the OLS fit.   39  The distances were therefore based on the predicted logits of CC for the target record versus the reference records.   Using this distance metric, initial tests to assess the value of k using 20-fold cross-validation revealed that k = 3 generally resulted in the lowest RMSPE. kNN was implemented using the yaImpute package in R (Crookston and Finley 2008) and using the GA to select subsets with six or fewer variables. Nonparametric: Random Forests Regression trees have been used in many studies to impute forest attributes using remotely sensed data (Hansen et al. 2002, Bucini et al. 2010, Naidoo et al. 2015). Thorough descriptions of regression tree methodologies and applications are given in Liaw and Wiener (2002) and Moisen and Frescino (2002).   As with kNN, RF is distribution-free but has been labelled as nonparametric in application literature (Cutler et al. 2007, Su et al. 2016).  RF is a unique implementation of regression trees where the user can set the number of trees to be grown and the number of predictor variables to assess at each split in a tree (Liaw and Wiener 2002). Each implementation of RF is a bootstrapped aggregation (i.e., bagging) of all trees. Error estimates are based on the observations not included in the bootstrapped sample for a particular tree (i.e., out-of-bag) and are, therefore, considered similar to k-fold cross validation (Liaw and Wiener 2002).  To examine how results are affected by the number of trees, trials were performed using 20-fold cross-validation altering the number of trees from 1 to 500, by 50. The RF default settings (500 trees) performed well enough to warrant no change in subsequent RF runs. As with kNN, the dependent variable in RF was logit-transformed CC. RF models were implemented in the 40  Random Forests v4.6-10 package of R (Liaw and Wiener 2002) and using the GA to select subsets with six or fewer variables. Semiparametric: Generalized Additive Models Semiparametric Generalized Additive Models (GAMs) are well-suited to modeling complex relationships where an underlying nonlinear relationship is not known a priori (Wood 2006). The generalized form of the GAM is:  𝑔(𝜇𝑖) =  𝛽0 + ∑ 𝑓𝑗𝑝𝑗=1(𝑥𝑖𝑗),    (3.5)  where 𝛽0 is the intercept, and fj  is the smoothing function for each predictor variable.   Essentially, GAM’s are based on a linear combination of piecewise polynomials (Wood 2006). The polynomials characterizing the contribution of each predictor variable are combined into separate functions, where each function is defined as a smoothing spline (Wood 2006).  Because of their flexibility, GAMs have been used with promising results to build models for characterizing forest attributes by employing remotely sensed data (Moisen and Frescino 2002, Staver et al. 2011).  GAM’s can assist in providing insight into the curvilinear relationships between a set of remotely sensed predictor variables and a dependent variable for use in refining parametric models (Moisen and Edwards 1999), or be used as predictive models in their own right (Kattenborn et al. 2015). For example, Armston et al. (2009) employed OLS with natural splines (a basic form of GAM, Wood 2006) to model tree cover in Australia using Landsat 5 TM and 7 ETM+ data , and Moisen and Frescino (2002) modeled CC, among other dependent variables, with GAM’s using AVHRR data.  41  A GAM with 𝑔(µ𝑖) again as the logit link of CC as a proportion was employed. As with GLM, two alternative model specifications were investigated: i) a normal distribution with the default variance function; and ii) a binomial distribution with the default variance. Cubic regression splines were used with smoothing parameter estimation via Restricted Maximum Likelihood, which are computationally efficient and include parameters which are directly interpretable (Wood 2006).    Initially, the number of basis dimensions that characterize each predictor variable were allowed to be estimated up to the maximum number allowed in fitting. As such, these models consistently had very low RMSE (i.e., < 4%), yet exhibited very high RMSPE (i.e., > 25%). With GAM’s, this is one indication of overfitting (Wood 2006). The number of basis dimensions was restricted to increasingly smaller values while computing the cross-validated RSMPE in order to protect against overfitting in subsequent predictor variable selection runs. It was found that by restricting the number of basis dimensions per variable to be < 2, models exhibited reasonably good RMSE and RMSPE. All GAM models were fitted with the mgcv package in R (Wood 2014), using the GA to select subsets with six or fewer variables. 3.2.4.3 Model evaluations The main objective was to find the best method for predicting percent CC in miombo woodlands. To assist in evaluating this objective, fit statistics for Root Mean Square Error (RMSE) and mean difference were computed as: 42      RMSE = √∑(𝑦𝑖 −  ?̂?𝑖)2𝑛𝑛𝑖=1, (3.6)   Mean difference = ∑(𝑦𝑖  −  ?̂?𝑖 )𝑛𝑛𝑖=1 ,   (3.7)  as well as the Pearson correlation coefficient for observed (𝑦𝑖) and predicted (?̂?𝑖) values. Since fit statistics can provide an optimistic view of model results (Packalén et al. 2012), validation statistics using measures of predictive performance calculated through cross-validation were also computed. For each method, 20-fold cross-validation was performed and the following validation statistics were calculated: RMSPE (Eq. 1) and the mean prediction difference using  Mean prediction difference = ∑(𝑦𝑖𝑣  −  ?̂?𝑖𝑣 )𝑛𝑛𝑖=1 ,   (3.8)  where 𝑦𝑖𝑣 and ?̂?𝑖𝑣 are defined the same as in Eq. 1. Based on these fit and validation statistics, one model was chosen and used to estimate total forest area for the entire district along with a corresponding 95% confidence interval. This was done using a bootstrapping pairs approach, following recommendations from McRoberts (2011). For each of 500 bootstrap samples, the selected model was fitted, creating 500 raster datasets of predicted percent CC (i.e., predicted proportion expressed as a percent) for the entire study area.  Using the sum of cover estimates in a 0.5 ha moving window, each 0.5 ha pixel was classified as forest (sum > 10%) or non-forest based on the forest area definition in Zambia (FAO 2014a). Total forest area was then calculated as the average of the 500 forest area estimates, along with the 95% confidence interval. 43  3.3 Results 3.3.1 Effects of auxiliary variables Remotely-sensed variables were first assessed alone to ascertain if they would provide comparable performance to using the full set of possible predictors that includes environmental and disturbance variables. Considering remotely-sensed data only, the average Pseudo-R2 for each sensor was 0.35 for RapidEye and 0.34 for Landsat 8 OLI and the average RMSPE for RapidEye was 13.9% and for Landsat 8 OLI 14.1% (Table 3.7). Even though RapidEye alone may have a slight advantage over Landsat 8 OLI, on average, the latter sensor performs slightly better when using the GA selection method (Pseudo-R2: 0.51; RMSPE: 12.7%).  Table 3.7. Comparison of model fit statistics from R2 versus GA variable selection using remotely sensed, environmental, and disturbance data to predict miombo percent canopy cover. Sensor Method RMSPE (3.1) Pseudo-R2 (3.2) Remotely senseda Environmental Disturbance RapidEye R2 14.9% 0.25 Blue, Green, NIR, NIR.tex1, Albedo.tex1, NDVI.tex1 NA R2 13.7% 0.47 Blue, Green, EVI Bulk, Carbon, Slope None chosen GA 12.9% 0.45 NDVI,  EVI.tex1, NDREGI.tex1, Blue.tex2, NDVI.tex2, SR.tex2 NA   GA 9.5% 0.73 ARVI, EVI AWC, Crs.frg, pH, CTI None chosen Landsat 8 OLI   R2 15.6% 0.17 Blue, SWIR2, Albedo, ARVI, Albedo.tex2, SWIR1.tex2 NA R2  13.1% 0.48 Green, Albedo Bulk, Carbon, Slope Dist2cap GA 12.7% 0.51  Green, Red, SWIR2, ARVI.tex2, Green.tex2, NIR.tex2 NA GA  8.6% 0.70 Green, Red, SWIR2.tex2 Bulk, Carbon, Crs.frg None chosen aSee Tables 3.3, 3.4, 3.5, and 3.6 for descriptions of variables. 44  The addition of environmental and disturbance variables improved the Pseudo-R2 and RMSPE statistics for both sensors (Table 3.7). When using auxiliary data across variable selection methods, the average RMSPE decreased 3.3% and 2.7% for Landsat 8 OLI and RapidEye, respectively, while the average Pseudo-R2 increased to 0.25 for both sensors. However, average results were similar for both sensors (mean RMSPE: 10.8% vs. 11.2%; mean Pseudo-R2: 0.59 vs. 0.60; Landsat 8 OLI and RapidEye, respectively). Consistently, the addition of soil variables decreased RMSPE and increased the Pseudo-R2 for both sensors. 3.3.2 Landsat 8 OLI versus RapidEye Landsat 8 OLI outperformed RapidEye in four of the six estimation methods using remotely sensed and auxiliary data (Table 3.8). RapidEye exhibited better performance only for the two nonparametric methods. Overall, the average RMSE differed by only 0.3% and the average Pearson correlation coefficient was the same across the two sensors; however, the average RMSPE for Landsat 8 OLI was almost 1% lower than RapidEye across all estimation methods. RapidEye exhibited slightly less average mean bias and mean prediction bias; however, these statistics were highly dependent on the estimation method under consideration.       45  Table 3.8. Goodness-of-fit and validation statistics with predictor variables by estimation method for percent canopy cover.                Fit statistics Validation statistics  Sensor Method* RMSE (3.6) MD (3.7) Correlation RMSPE (3.1) MPD   (3.8) Predictor variables RapidEye   GLM-N 8.2% 1.5% 0.86 9.5% 1.4% ARVI, EVI, AWC, Crs.frg, pH, CTI GLM-B 8.9% 0.0% 0.82 10.1% 0.3% EVI, EVI.tex2, Green.tex2, AWC, Crs.frg, pH GAM-N 7.8% 0.9% 0.87 9.4% 1.0% ARVI, EVI, AWC, Crs.frg, pH, CTI GAM-B 8.0% 0.0% 0.86 9.9% -0.3% EVI, NDREGI, NIR.tex2, AWC, pH, Slope kNN 11.2% 1.7% 0.70 14.6% 1.7% EVI, ARVI.tex1, NDRE.tex1, NDVI.tex1, Red.tex2, CEC RF 12.1% 3.7% 0.69 12.3% 3.8% Green.tex1, NDRE.tex1, SR.tex1, NDREGI.tex2, Bulk, Dis2hydro     Landsat 8 OLI   GLM-N 7.4% 1.1% 0.89 8.6% 1.2% Green, Red, SWIR2.tex2, Bulk, Carbon Crs.frg GLM-B 8.5% 0.0% 0.84 9.6% 0.1% Green, Red, SWIR1.tex2, Bulk, Carbon, Crs.frg GAM-N 6.8% 1.4% 0.91 8.1% 1.5% Green, Red, SWIR2.tex2, Bulk, Clay, Dist2popcap GAM-B 6.6% 0.0% 0.91 8.4% -0.1% NDVI, ARVI.tex2, SWIR1.tex2, Dist2popcap, AWC, Slope kNN 13.1% 2.2% 0.59 14.2% 3.4% NDMI, EVI.tex2, SWIR2.tex2, Carbon, Slope RF 12.2% 3.9% 0.68 12.0% 4.0% Albedo.tex2, ARVI.tex2, Blue.tex2, AWC, Slope, Dist2hydro * GLM-N and GLM-B: generalized linear model with a normal and a binomial distribution, respectively; GAM-N and GAM-B: generalized additive model with a normal and a binomial distribution, respectively; kNN: k-Nearest Neighbor imputation; RF: Random Forest regression tree; RMSE: Root Mean Square Error; RMSPE: Root Mean Square Prediction Error from 20-fold cross-validation; MD: mean difference; Correlation: Pearson correlation coefficient between observed and predicted values; MPD: mean prediction difference. See Tables 3.3, 3.4, 3.5, and 3.6 for descriptions of variables. 46  3.3.3 Comparison of methods None of the nonparametric estimation methods performed well in this study. On average across sensors, the RMSE and RMSPE was almost 4% higher for nonparametric over parametric and semiparametric methods. Likewise, the average mean bias and mean prediction bias was between 2.5% and 3.0% higher, with a strong tendency to under-predict CC. GLM’s demonstrated acceptable predictive performance with average biases under 1%, and were simple to fit. However, three of the four GLM’s slightly under-predicted CC at the high end of the observed values (Fig. 3.6, Fig. 3.7). GAM’s demonstrated the best average RMSE and RMSPE (7.3% and 9.0%), though average mean bias and mean prediction bias were nearly identical to GLM’s. Overall, Landsat 8 OLI GAM’s outperformed all other sensor/fitting method combinations.  The predictive performance of Landsat 8 OLI GAM-N and GAM-B were similar, although the mean bias and mean prediction bias of GAM-B were clearly better than GAM-N. This same trend occurred for all model comparisons of normal versus binomial distributions. For Landsat 8 OLI GAM-N, five of the six variables were initially estimated with a basis dimension of one. This implies that these variables should enter the model as parametric terms (Wood 2006). When this model was re-fit as such, the fit and validation statistics remained unchanged, so this model was fit with five parametric variables and one nonparametric variable, distance to district capital. All Landsat 8 OLI GAM-B variables were estimated as nonparametric variables. The Landsat 8 OLI GAM-N demonstrated a slightly better RMSPE than the GAM-B, and predictions at the high end of the range of CC were more balanced (Fig. 3.7). However, the RMSE was slightly lower with GAM-B, while the mean bias and mean prediction bias were 47  approximately 1.5% less than GAM-N. Therefore, the following model for predicting CC was selected:    𝑐𝑎𝑛𝑜𝑝𝑦 𝑐𝑜𝑣𝑒𝑟% ̂    = −1.8747 + 𝑓(NDVI) + 𝑓(𝐴𝑅𝑉𝐼. 𝑡𝑒𝑥2) + 𝑓(𝑆𝑊𝐼𝑅1. 𝑡𝑒𝑥2) + 𝑓(𝑑𝑖𝑠𝑡2𝑝𝑜𝑝𝑐𝑎𝑝)+ 𝑓(𝐴𝑊𝐶) + 𝑓(𝑠𝑙𝑜𝑝𝑒).   (3.9)  Using this model and 500 bootstrap samples, the total forest area and associated 95% confidence interval was 758,100 ha +/- 3,956 ha for the 1,039,269 ha study area. Approximately 73% of the district was estimated as forest for the year 2013 (Fig. 3.8). Areas with higher CC were in more mountainous terrain on the western half of the study area. Although human disturbance predictor variables were not chosen in the final model, the indirect effect of distance to the district capital is noticeable in the southeast corner of the study area which was estimated as non-forest. Most of the population lives in this area, which is more densely settled and intensively farmed than other regions in the district. Other non-forest areas away from the capital and roads occurred at low elevation valley bottoms which receive less rainfall and have higher mean annual temperatures.    48   Figure 3.6. Observed vs predicted percent canopy cover (CC) using RapidEye imagery with environmental and disturbance variables  (n=64). The solid line indicates 1:1 correlation between observed and predicted values. Acronyms defined in Table 3.7 footnote.  0 10 20 30 40 50 600204060Obs CC (%)Pred CC (%)GLM-N0 10 20 30 40 50 600204060Obs CC (%)Pred CC (%)GLM-B0 10 20 30 40 50 600204060Obs CC (%)Pred CC (%)GAM-N0 10 20 30 40 50 600204060Obs CC (%)Pred CC (%)GAM-B0 10 20 30 40 50 600204060Obs CC (%)Pred CC (%)kNN0 10 20 30 40 50 600204060Obs CC (%)Pred CC (%)RF49   Figure 3.7. Observed vs predicted percent canopy cover (CC) using Landsat 8 OLI imagery with environmental and disturbance variables  (n=64). The solid line indicates 1:1 correlation between observed and predicted values. Acronyms defined in Table 3.7 footnote. 0 10 20 30 40 50 600204060Obs CC (%)Pred CC (%)GLM-N0 10 20 30 40 50 600204060Obs CC (%)Pred CC (%)GLM-B0 10 20 30 40 50 600204060Obs CC (%)Pred CC (%)GAM-N0 10 20 30 40 50 600204060Obs CC (%)Pred CC (%)GAM-B0 10 20 30 40 50 600204060Obs CC (%)Pred CC (%)kNN0 10 20 30 40 50 600204060Obs CC (%)Pred CC (%)RF50    Figure 3.8. Map of predicted canopy cover per 0.1 ha pixel and forest cover per 0.5 ha pixel for Nyimba District, Zambia. Map results from applying a binomial GAM using Landsat 8 OLI with soils and disturbance data.  3.4 Discussion Texture indices were consistently selected as important predictor variables contributing to model performance for either sensor. For the RapidEye imagery, both models without auxiliary data employed within-plot texture indices (Texture 1) to characterize the variability of spectral responses for a particular band or index within the plot area. However, texture around the plot figured more prominently in the better performing RapidEye models with no bias (GLM-B and GAM-B) and was employed in all of the Landsat 8 OLI models. With RapidEye, there was no 51  clear consensus on which texture index exhibited the most influence; however, with Landsat 8 OLI texture of either SWIR1 or SWIR2 appeared in five of the six models. Interestingly, the Landsat 8 OLI parametric and semiparametric models with SWIR texture outperformed each of their RapidEye counterparts. Miombo woodlands are known for exhibiting high heterogeneity in forest cover (Frost 1996), and these findings concerning the use of these variables are in line with those reached by Hansen et al. (2002), who also found texture to be an important variable when assessing CC in western Zambia.  While RapidEye texture variables were not systematically chosen, the vegetation index EVI was selected in five of the six RapidEye models. EVI has been shown to be well correlated with field observations of canopy phenology in miombo woodlands (Ryan et al. 2014a). For Landsat 8 OLI, the green and red bands, in combination with SWIR texture, were selected in three models, while only two vegetation indices were included as predictor variables (NDVI and NDMI, one method each). However, NDVI has been shown to be a good predictor of CC in dry forests of Africa (Wu et al. 2013) and Fuller et al. (1997) demonstrated the relationship between NDVI and CC in Zambian miombo woodlands using a radiative transfer model.  Including soils information contributed to the best performing model. The most frequently chosen soil variables were: available soil water capacity (AWC, six models), coarse fragments (five models), bulk density (four models), and pH (four models). These soil characteristics have been shown to be correlated with forest cover and land use types in the miombo ecoregion. Tinley (1982) demonstrated that soil moisture, in combination with soil texture, is a major driver of miombo woodland structural characteristics. Campbell et al. (1998) found that AWC was higher, while coarse fragments (i.e., particles > 2mm) and bulk density were lower, in forested 52  areas versus non-forested areas. Data in this research shows some similar relationships, where higher AWC is associated with higher CC while lower bulk density is related with higher CC (Pearson correlation coefficients 0.49 and -0.59, respectively).  For the RapidEye imagery, the interaction of EVI with AWC, coarse fragments and pH seemed to drive model performance; while for Landsat 8 OLI the most frequently occurring combination was the green and red bands with SWIR texture and bulk density. While the best performing model (Landsat 8 OLI GAM-B) is driven by a different combination of NDVI, texture, AWC, slope and distance to the district capital, this model is sensible given that disturbance from human activity plays a major role in shaping miombo woodland characteristics (Frost 1996, Ahrends et al. 2010). Linkages between variables such as soils, topography, and disturbance with African dry forest attributes, such as CC, have been investigated at local, regional and continental scales largely for inferential purposes (Mapaure and Campbell 2002, Sankaran et al. 2008, Staver et al. 2011). Separate work on predictive methods to estimate CC in these dry forests using optical remotely sensed data alone has shown utility (Hansen et al. 2002, Wu et al. 2013, Adjorlolo and Mutanga 2013), although few studies have investigated the combination of optical remotely sensed data with auxiliary variables for this purpose. As such, comparisons with other works are challenging due to differing definitions, data sources, and uncertainty measures. Nevertheless, these results align with other estimates, and are competitive with assessments using active sensor data, such as radar, in dry forests. RMSE estimates of CC models using medium or high resolution optical satellite imagery for tropical and sub-tropical dry forests are generally within the range of 7 to 13% (Hansen et al. 2002, Armston et al. 2009, Adjorlolo and Mutanga 2013, Ozdemir 2014, 53  Karlson et al. 2015). Recent studies on modeling CC in dry forest types of South Africa using radar data (Urbazaev et al. 2015, Naidoo et al. 2015), or radar data combined with Landsat 7 ETM+ (Bucini et al. 2010) achieved RMSE estimates of 7.9%, 8.8%, and 8.9%, respectively. The RMSPE of the best performing method for this study, assessed through 20-fold cross-validation, was 8.4%, while the corresponding RMSE was 6.6%.   Methods to estimate CC in dry forests of Africa frequently use nonparametric, distribution-free methods such as Random Forests and kNN imputation (Hansen et al. 2002, Bucini et al. 2010, Naidoo et al. 2015). While these methods require few or no theoretical assumptions and are relatively easy to implement (Temesgen and Ver Hoef 2015), they exhibited poor fit and validation statistics for both sensors in this study.  Moisen and Frescino (2002) and Temesgen and Ver Hoef (2015) also found that other methods outperformed Random Forest regression trees. Limitations leading to poor performance of distribution-free methods such as these may result from small sample sizes and imbalanced data (LeMay and Temesgen 2005; Temesgen and Ver Hoef 2015). Field measures used in this study were limited to a very small sampling intensity (64 plots distributed over a ~1,000,000 ha study area). The GAM’s employed in this study were able to accurately and precisely estimate percent CC. However, it was found that GAM’s have poor predictive capability when the number of basis dimensions is allowed to be estimated using default settings in the R package mgcv. In this case, the algorithm tended to maximize all available model degrees of freedom, which resulted in very high RMSPE. Poor predictive performance indicates that the model was over-fitting the data (Wood 2006). Assessing the appropriate number of basis dimensions requires an extra step to assess the maximum number of dimensions possible without overfitting. By performing variable 54  selection first with the highest number of basis dimensions possible then restricting this number to increasingly fewer dimensions and computing the cross-validated RMSPE, the number of basis dimensions was optimized so that subsequent variable selection runs were not overfit. Given the high heterogeneity in miombo woodlands and the relative importance of texture variables, models using the higher spatial resolution RapidEye imagery would generally be expected to perform better than the Landsat 8 OLI imagery. However, since the RapidEye pixel values were averaged over the plot area, the advantages of increased spatial resolution were limited to the Texture 1 variables (i.e., within-plot standard deviation). Models using RapidEye band values and vegetation indices were built at the plot scale, where model performance was mostly enhanced by using soils and disturbance predictor variables. Improvements in CC models might occur by using RapidEye at the pixel level, if the response variable is fractional tree cover per pixel. However, the focus of this research is on use of NFI data and there are challenges which would limit the effectiveness of high resolution pixel analyses in this context. First, the NFI protocol stipulates collection of x,y location data at the center of each tree stem to the nearest meter. Relating RapidEye pixels to tree stem, and tree cover, would likely require collecting tree stem location data to the nearest 0.1 m. Second, there is uncertainty in the RapidEye image mosaic geometric registration process. This issue introduces more overall spatial uncertainty for RapidEye, since many more images are required to cover the study area than for Landsat 8 OLI. While the average error of each of these images is within acceptable limits, their combination may introduce too much spatial uncertainty to effectively link tree cover polygons with RapidEye imagery at the pixel level.  55  The specific models identified in this research should be applied with caution when attempting implementation outside the range of conditions (forest types, soil, rainfall, disturbance patterns, etc.) in which the models were created (Foody et al. 2003).   However, application of the following suggestions may allow wider use of Zambian NFI data to CC and forest area modeling. First, the DBH-crown radius models used in this study were mostly generated outside of Zambia, and their predictive capacities are not well understood within this study area. These models could be updated to account for forest types in Zambia and applied to NFI data relatively easily. Second, overlap of CC on the plot which may have occurred from trees located outside the plot boundary was not taken into account during ground plot data collection. Similar issues have been noted with the use of ALS data in forest inventories, where the general recommendation is to use a larger plot size (McRoberts et al. 2014a). Edge effect correction methods are also well known (Gregoire and Valentine 2007), and their application in this context could warrant further investigation. Third, new technology in the form of unmanned aerial vehicles (i.e., drones) outfitted with cameras could be tested to determine if CC per plot can be estimated from aerial photography (Pajares 2015). CC is an important forest attribute that is useful for many purposes, such as ecosystem process models, conservation planning, and monitoring forest area. However, building an accurate and parsimonious model to estimate CC in miombo woodlands can be a difficult task. A wide variety of remotely sensed data sources and processing methods are available. In this chapter, Landsat 8 OLI and RapidEye remotely sensed data sources were compared, with and without environmental and disturbance data, across several modeling methods. Several insights can be gained through this process. First, RapidEye can be a competitive source of satellite imagery 56  when considered without auxiliary data; however, Landsat 8 OLI when combined with soils and disturbance information demonstrated consistently better performance over RapidEye and it is freely available. Second, CC can be predicted using a binomial Generalized Additive Model with a Root Mean Square Prediction Error of < 8.4%, which is competitive with models using active sensor remote sensing data (e.g., radar). Third, because of multiple spatial uncertainties, plot level analyses are likely the finest scale estimations which can be produced with NFI data. Lastly, updating DBH-crown radius models for miombo woodland species provides an opportunity to extend use of NFI data to forest area mapping across Zambia.  57  Chapter 4: Above-ground biomass modeling 4.1  Introduction AGB is a key variable of interest in forest monitoring programs, where estimates of AGB are necessary for assessing fuel wood and timber availability, as well as for monitoring forest carbon stocks. Integrating ground plot data from a forest inventory with environmental and/or remotely-sensed predictor variables to predict AGB has shown increasing utility for estimating AGB (Moisen et al. 2006, McRoberts et al. 2010, Lu et al. 2014, GOFC-GOLD 2015). Optical remotely-sensed data alone have been used to estimate miombo AGB with varying success (Samimi and Kraus 2004, Kashindye et al. 2013, Næsset et al. 2016). The use of these data is complicated by complex, nonlinear relationships between AGB and vegetation indices (Lu et al. 2014). Active sensor remotely-sensed data such as Airbone Laser Scanning (ALS) or radar may provide improved accuracy of AGB models (Ryan et al. 2012, Mitchard et al. 2013, Mauya et al. 2015, Næsset et al. 2016), although this is not guaranteed (Solberg et al. 2015). Another way to improve AGB estimation is to incorporate percent CC as a predictor variable (Tiwari and Singh 1984, Lefsky et al. 2002, Hall et al. 2006, González-Roglich et al. 2014, GOFC-GOLD 2015), particularly since percent CC is more readily estimated by optical remotely-sensed data (Lefsky and Cohen 2003, Lu et al. 2014).  González-Roglich and Swenson (2016) capitalized on the relationship between AGB and CC for modeling forest carbon in temperate savannahs of South America. They first developed a model to estimate CC using Landsat 5 Thematic Mapper along with topographic and climatic variables, and then used the estimated CC as the sole predictor to estimate forest carbon. Using this 58  process, the percent root mean square error (RMSE) for prediction of forest carbon was 35%.  Similar studies in arid woodlands of Australia (Suganama et al. 2006) and Sudan (Wu et al. 2013) also reported favorable results using estimated CC to estimate AGB. One complication is that the definition of CC varies (Jennings et al. 1999). In some cases, CC was defined as the percentage of the sky blocked by tree crowns over a hemispherical view (Jennings et al. 1999), often measured using a spherical densiometer. Alternatively, the vertical projection of tree crowns onto the ground is currently used in many countries as the basis for minimum forest area definition (FAO 2010). González-Roglich and Swenson (2016) used the former definition, while Suganama et al. (2006) and Wu et al. (2013) used the latter. Therefore, it is an open question as to which measure of CC may better relate to AGB.  Models to estimate AGB that also incorporate information on topography, soils, and disturbances may further improve the accuracy (Moisen and Frescino 2002, Powell et al. 2010, Pflugmacher et al. 2014, GOFC-GOLD 2015). In the miombo ecoregion, as in other ecosystems, ecological and physiological factors limit vegetation establishment and survival, while disturbances result in changes to vegetation structure and composition (Frost 1996, Chidumayo et al. 1996, Sankaran et al. 2008). Soil characteristics, fire regimes, and anthropogenic use have all been identified as principal determinants of vegetation structure (Timberlake and Chidumayo 2011, Ryan and Williams 2011). Investigations of these and other variables for estimating miombo AGB appear to be sparse, and further exploration is warranted.  The main objective of this chapter is to investigate methods for estimating total AGB for a study area within the miombo ecoregion by using ecologically relevant predictor variables. To achieve this objective, the following specific questions were addressed: 1) Which definition of CC is 59  more useful as a single predictor variable to estimate AGB, the hemispherical view or the vertical projection view?; 2) Does AGB estimation accuracy improve if CC is combined with other predictor variables?; 3) What model-based AGB estimation method performs best in terms of fit and validation statistics?; and, 4) How do model-based methods compare to a simple design-based estimate of the sample mean? This research is intended as a case study to help inform development of forest monitoring programs which apply NFI data to estimation of AGB at sub-national or national scales. 4.2 Methods 4.2.1 AGB per ground plot Tree measurements described in 2.2 were used to calculate tree-level AGB, using species-specific allometric models developed by Chidumayo (2012) and utilized by the Zambian Forest Department (GRZ 2014). In some cases, species-specific models were not available and woodland or forest type models were used (Chidumayo 2012). Tree-level AGB in plots and subplots was summarized to obtain plot-level AGB in tonnes per hectare (t/ha). 4.2.2 Predictor variables Two types of estimated CC were employed as a single predictor variable in separate models to estimate AGB. First, CC was estimated as the vertical projection of tree crowns onto the ground using a binomial GAM developed for the same study area (Halperin et al. 2016a), referred to as  𝐶?̂?VERT. Second, CC was estimated as the percentage of the sky blocked by tree crowns over a hemispherical projection view by using the spherical densiometer. For this purpose, a binomial GAM was built to estimate percent CC using the same predictor variable selection methods 60  described in Halperin et al. (2016a), and refer to its estimates as 𝐶?̂?HEMI. Results of estimating AGB using either  𝐶?̂?VERT or 𝐶?̂?HEMI were compared, and the better performing model was used as the base upon which to add other predictor variables. The other predictor variables investigated for use in AGB models include: remotely sensed predictor variables from the Landsat 8 OLI mosaic (Table 3.3), disturbance predictor variables (Table 3.5), and environmental predictor variables (Table 3.6). 4.2.3 AGB models 4.2.3.1 Models using CC only As noted, AGB was estimated first using either  𝐶?̂?VERT or  𝐶?̂?HEMI as a single predictor variable using three models. A semiparametric GAM was used since it is well-suited to modeling complex relationships where an underlying nonlinear relationship is not known a priori (Wood 2006). Because of their flexibility, GAMs have been used to estimate forest attributes from remotely-sensed data (Moisen and Frescino 2002, Kattenborn et al. 2015).  The form of the GAM is:  log(𝐴𝐺?̂?𝑖) = β0 +  ∑ 𝑓𝑗𝐽𝑗=1(x𝑖𝑗)  +  ε𝑖 ,     (4.1)  where 𝑖 is the ground plot, 𝐴𝐺𝐵𝑖 is in t/ha, β0 is the intercept, ε𝑖 is the residual error, and fj  is the smoothing function for a predictor variable x𝑗 created via a regression spline (Wood 2006).  The GAM was fit with a normal distribution, the log link, and the default variance function. Using the log link avoids negative AGB estimates. A value of 0.1 was added to the AGB for two ground plots with no trees. The maximum smoothing parameter was set at 2, in order to provide 61  sufficient flexibility to account for potential nonlinearity in the model (Halperin et al. 2016a), yet reduce the possibility of overfitting that may occur in GAM’s (Moisen et al. 2006). The GAM was fit with the mgcv package v. 1.8-9 in R v. 3.2.2 (Wood 2015, R Development Core Team 2015). Two nonlinear models that restricted the AGB estimate to be positive with an upper limit were also fit. First, a sigmoidal model (e.g., McRoberts et al. 2015) was used:  𝐴𝐺?̂?𝑖 =  α1 + exp (β0 +  ∑ β𝑗x𝑖𝑗𝐽𝑗=1 )+  ε𝑖 , (4.2)  where parameters are α and the β’s, and the other terms are as in Eq. 4.1. Second, an exponential model (e.g., Anaya et al. 2009) was used:  𝐴𝐺?̂?𝑖 =  α x exp (∑β𝑗x𝑖𝑗⁄𝐽𝑗=1) +  ε𝑖 . (4.3)  To fit these nonlinear models, the Gauss-Newton method in the nls function in R v.3.2.2 was used. Starting parameters were estimated using the linearized form of each nonlinear model (i.e., log(AGB𝑖)). After experimenting with a flexible or fixed asymptote α, it was found that setting α = 500 produced realistic estimates. This is in line with Næsset et al. (2016), who truncated miombo AGB estimates from linear ordinary least squares models at twice the maximum observed value in their sample data. For each fitted model (three models and  𝐶?̂?VERT or 𝐶?̂?HEMI), 20-fold cross-validation was used to calculate the root mean square prediction error (RMSPE) using Eq. 3.1 and RMSPE as a percent of the mean (RMSPE%), defined as: 62   RMSPE% = (RMSPE?̅?) ∗ 100, (4.4)  where ?̅? is the mean AGB (t/ha) calculated from the ground plots. To compare performance across the three models using  𝐶?̂?VERT or 𝐶?̂?HEMI, the mean RMSPE and mean RMSPE% were calculated. Based on these validation statistics, either  𝐶?̂?VERT or  𝐶?̂?HEMI was included in models with multiple predictor variables. 4.2.3.2 Models using multiple predictor variables Next, it was investigated whether expanding the models with increasingly diverse combinations of predictor variable sets could improve the accuracy of AGB estimates. For this purpose, models were compared using: 𝐶?̂?VERT or  𝐶?̂?HEMI alone; remotely-sensed (Landsat) variables alone; environmental and disturbance (Env/Dist) variables alone; and combinations of each predictor variable set (i.e., Landsat and Env/Dist, 𝐶?̂?VERT or  𝐶?̂?HEMI and Landsat, etc.). A genetic algorithm (GA) was used to perform predictor variable selection within each predictor variable set according to the same criteria in Chapter 3. For each predictor variable set (except  𝐶?̂?VERT or  𝐶?̂?HEMI alone), the GA was implemented one time and the lowest RMSPE was recorded. Mean RMSPE and mean RMSPE% for each predictor variable set were then calculated and compared, and the predictor variable set that minimized these two fit statistics was chosen.  4.2.3.3 Final variable selection Final predictor variable selection was performed using the predictor variable set that minimized RMSPE in the previous section. For this purpose, the GA was implemented 100 times for each model using the same fitness and stopping criteria as Chapter 3, and the subset of predictor 63  variables that minimized RMSPE across all 100 GA implementations was selected. For the sigmoidal and exponential models, commonly used transformations (i.e., square, square root, inverse) were also considered for each predictor variable, once the predictor variables were chosen through this process. Use of transformations for predictor variables has been examined in nonlinear models for improving model performance (Wang et al. 2007, Timilsina and Staudhammer 2012). To compare performance among the three models, the following fit statistics were computed: Pseudo-R2 (Eq. 3.2), RMSE (Eq. 3.6), MD (Eq. 3.7), and RMSE as a percent of the mean (RMSE%) was calculated as:  RMSE% = (RMSE?̅?) ∗ 100. (4.5)  20-fold validation statistics were also calculated, namely RMSPE (Eq. 3.1) and the mean prediction difference (MPD, Eq. 3.8).  4.2.3.4 Uncertainty assessment and model selection After choosing the best performing predictor variable subset for each model, the AGB total (tonnes), 95% confidence interval (CI) for the total, and AGB t/ha were estimated for the entire study area. This was accomplished using the bootstrap method described in Chapter 3. Because the GAM is not asymptotic on the upper end, the maximum estimated value for AGB (t/ha) was limited to the same value as the asymptote for the sigmoidal and exponential models, as is commonly done in non-asymptotic models (e.g., Næsset et al. 2016). For each model, total AGB (Ŷboot) was estimated using the following equation: 64   Ŷboot = ∑Ŷboot𝑘nbootnboot𝑘=1,    (4.6)  where Ŷbootk  is the estimated total AGB by summing all pixels (in tonnes) across the study area for bootstrap 𝑘, and nboot is 500. AGB/ha was calculated from dividing Ŷboot by the number of ha in the study area. The variance for Ŷboot was estimated using the following equation:  Var̂(Ŷboot) = ∑(Ŷbootk − Ŷboot)2nboot − 1nbootk=1. (4.7)  This variance estimate was used to calculate a 95% confidence interval (CI) for estimated total AGB based on the normal distribution (Efron and Tibshirani 1986), and reported the half-width of the CI for each model. To establish a baseline for comparing model performance, a null model was fit with no predictor variables. Using this null model, the AGB total, 95% CI for the total, and AGB (t/ha) were estimated (Bater et al. 2009). Then, the relative efficiency for each model was calculated by dividing the null model variance estimate by each model’s variance estimate, following Næsset et al. (2016). Finally, the four model-based estimates of total AGB, AGB/ha, and CI were compared to a simple design-based estimate of AGB total, 95% CI for the total, and AGB (t/ha) using the NFI systematic sample of ground plots, assuming equal probability sampling (Cochran 1977).  Based on these statistics, along with observed versus predicted AGB graphs, one model was chosen and used to estimate AGB for each 20m x 50m pixel in the study area. To assess the pixel-level variability of the AGB estimates, the Coefficient of Variation (CV) in percent was calculated (i.e., standard deviation per pixel across the 500 resamples divided by the corresponding mean pixel value, multiplied by 100, which was also used by Carreiras et al. 65  2013, Faßnacht et al. 2014, Kattenborn et al. 2015), a map of CV was presented, and the mean CV for the study area was reported.  4.3 Results 4.3.1 Models using CC only Observations of AGB (t/ha) and CCVERT (%) for the miombo woodland vegetation type follow a more regular pattern than the mopane woodland and riparian forest vegetation types, which exhibit a large range of observed AGB values at low observed values of CCVERT (Fig. 4.1). Observations of CCHEMI greater than 50% display a wide range of observed AGB values; however, observed AGB values corresponding to CCHEMI less than 50% follow a more regular pattern. As expected, the six nonforest ground plots which were measured have low observed AGB and low observed CC for both forms of CC. A binomial GAM was used to fit separate models for CCVERT and CCHEMI; the model to estimate  𝐶?̂?VERT had better fit and validation statistics than the model to estimate  𝐶?̂?HEMI (Table 4.1). However, in predicting AGB, clear differences are apparent in the usefulness of CC  as the percentage of the sky blocked by tree crowns over the hemispherical projection view (𝐶?̂?HEMI) versus CC as the vertical projection of tree crowns onto the ground (𝐶?̂?VERT). Across three model forms, 𝐶?̂?HEMI had an average RMSPE% that was 18% lower than for  𝐶?̂?VERT (Table 4.2). The exponential model performed poorly with RMSPE (72.1 t/ha) larger than the mean estimate of AGB t/ha from the sample data (64.1 t/ha). Considering only the GAM and sigmoidal models, the average difference in RMSPE between  𝐶?̂?VERT and  𝐶?̂?HEMI as a single predictor variable was 5 t/ha in favor of 𝐶?̂?HEMI. 66  Because of these statistics and relationships, 𝐶?̂?HEMI was chosen for estimating AGB in models with multiple predictor variables.  Figure 4.1 Scatterplot for canopy cover (CC) and above-ground biomass (AGB). Top row: Two forms of observed CC versus observed AGB, by vegetation type or land use. CCVERT is the vertical projection of tree crowns onto the ground, derived from crown area-DBH models; CCHEMI is percentage of the sky blocked by tree crowns measured using a spherical densiometer. Bottom row: Two forms of predicted CC versus observed AGB, by vegetation type or land use, with three models to estimate AGB based on 𝑪?̂?𝐕𝐄𝐑𝐓 or 𝑪?̂?𝐇𝐄𝐌𝐈.   𝐂𝐂𝐕𝐄𝐑𝐓 (%) 𝐂𝐂𝐇𝐄𝐌𝐈 (%) Obs AGB (t/ha) Obs AGB (t/ha) Obs AGB (t/ha) Obs AGB (t/ha) 𝑪?̂?𝐕𝐄𝐑𝐓 (%) 𝑪?̂?𝐇𝐄𝐌𝐈 (%) 67  Table 4.1. Goodness-of-fit and validation statistics for two binomial generalized additive models to estimate percent canopy cover based on crown radius-DBH models (𝐂?̂?𝐕𝐄𝐑𝐓) and measurements using a spherical densiometer (𝑪?̂?𝐇𝐄𝐌𝐈).  Fit statistics Validation statistics   Model RMSE (3.6) MD (3.7) RMSPE (3.1) MPD  (3.8) 𝐶?̂?VERT = −1.8747 + 𝑓(NDVI) + 𝑓(ARVI. texture2) +   𝑓(SWIR1. texture2)  +  𝑓(D2cap) + 𝑓(AWC) + 𝑓(slope)  8.5% 0.0% 8.4% 1.5%  𝐶?̂?HEMI =  0.6788 + Fire. late + D2cap + carbon+  𝑓(bulk + 𝑓(crs. frg)  13.4% -0.7% 17.2% -0.5% * See Tables 3.3, 3.5, and 3.6 for descriptions of predictor variables. RMSE: root mean square error; MD: mean difference; RMSPE: root mean square prediction error; MPD: mean prediction difference.  Table 4.2. Above-ground biomass (AGB, t/ha) root mean square prediction error (RMSPE) using different predictor variable dataset combinations and models , ordered by mean RMSPE and mean RMSPE%. Model 𝐶?̂?VERT 𝐶?̂?HEMI Landsat Env/ Dist Landsat + Env / Dist 𝐶?̂?HEMI + Landsat 𝐶?̂?HEMI + Env / Dist 𝐶?̂?HEMI + Env/Dist + Landsat GAM (4.1) 50.1 44.3 48.8 41.5 40.7 43.4 37.9 35.1 Sigmoidal (4.2) 49.8 45.2 47.8 47.0 45.6 44.2 39.8 37.5 Exponential (4.3) 72.1 46.3 45.9 46.6 44.3 42.8 40.5 39.4 Mean RMSPE 57.3 45.3 47.5 45 43.5 43.5 39.4 37.3 Mean RMSPE% 89% 71% 74% 70% 68% 68% 61% 58%   4.3.2 Models using multiple predictor variables A clear pattern in reduction of RMSPE can be seen by including more diverse combinations of predictor variable datasets (Table 4.2). Across models, 𝐶?̂?HEMI alone is a better predictor of AGB than models that use only Landsat 8 OLI predictor variables; however, the exponential 68  model that uses Landsat 8 OLI has a slight edge over  𝐶?̂?HEMI (RMSPE 45.9 vs. 46.6, respectively). The predictor variable combination of environmental and disturbance data alone and 𝐶?̂?HEMI alone are comparable with two of the three combinations that have two predictor variable datasets (Landsat + Env/Dist and  𝐶?̂?HEMI + Landsat), with a difference in mean RMSPE% of only three percent. Noticeable improvements in AGB RMSPE% are seen when combining 𝐶?̂?HEMI with environmental and disturbance data, as mean RMSPE% drops seven percent. However, the greatest improvement is when performing variable selection using all three predictor variable datasets combined, with a mean RMSPE% of 58%. Therefore, final variable selection was performed using this combination. 4.3.3 Final variable and model selection One influential observation, which likely resulted from haze in the remotely-sensed data, was evident across all models. As the influence could be considered uncommon given this dataset, the observation was discarded and continued with the analysis, following the recommendation of Kutner et al. (2005, p.438). Using the final selected variables, the use of predictor variable transformations for the sigmoidal and exponential models was explored. By incorporating squared terms in the sigmoidal model for two predictor variables (𝐶?̂?HEMI2 and pH2), an improvement of 4% in terms of RMSE% and an overall better fit in graphs of observed versus estimated values of AGB was observed. The exponential model did not improve by incorporating any transformed predictor variables. However, it was found that the fit and validation statistics of the exponential model remained the same after removing one variable (EVI.texture2); therefore, a subset of five predictor variables was used for this model. With the GAM, three variables were 69  estimated with a smoothing parameter of 1 (NBR, CTI, bulk density); therefore the model was fitted with these variables as parametric terms, as recommended by Wood (2006).  For each of the three models, average RMSE% was less than 60% and average RMSPE% was less than 65% (Table 4.3). The mean prediction difference (MPD) for the sigmoidal model was the highest at 3.2 t/ha; however, mean differences (MD) across all models were less than 1 t/ha. Despite the difference in MPD, the fit and validation statistics for the GAM and the sigmoidal models were otherwise similar, while the exponential model had the worst fit (Fig. 4.2). Each model had a tendency to under-estimate AGB at higher values, which was most pronounced in the exponential model; however, observations and estimations were relatively balanced along the 1:1 line up to approximately 100 t/ha.  As expected, 𝐶?̂?HEMI was selected as a predictor variable in all three models. Vegetation indices were chosen more consistently over raw band values or texture indices among the Landsat 8 OLI predictor variables. The Normalized Difference Moisture Index (NDMI) was chosen in both the GAM and the sigmoidal model, the Normalized Burn Ratio (NBR) was selected as another index in the GAM, while the common Normalized Difference Vegetation Index (NDVI) was selected as the only vegetation index in the exponential model. In the GAM and the sigmoidal model, NDMI was selected in conjunction with topographic or soil variables related to soil moisture (AWC, CTI, D2hydro). Available soil water capacity (AWC) was selected in both the sigmoidal and exponential models. Frequency of late season fire was the only disturbance variable selected, and only in the sigmoidal model.   70  Table 4.3. Fit and validation statistics for predicting miombo above-ground biomass (AGB, t/ha) using three models. Model RMSE (3.6) RMSE% (4.5) MD (3.7) Pseudo- R2 (3.2) RMSPE (3.1) RMSPE% (4.4) MPD (3.8) Variables* GAM    (4.1) 29.2 45% 0.8 0.65 36.8 57% 0 CĈHEMI, NBR, NDMI, Bulk, CTI, D2hydro Sigmoidal (4.2) 29.9 47% 0.9 0.63 36.9 58% 3.2 CĈHEMI, CĈHEMI2, NDMI, ARVI.texture2, Fire.late, AWC, pH, pH2 Exponential (4.3) 36.1 56% 0.6 0.47 39.7 62% 0.4 CĈHEMI, NDVI, AWC, Sand, Silt * See Tables 3.3, 3.5, and 3.6 for descriptions of remotely-sensed, environmental and disturbance variables.   Figure 4.2. Observed versus estimated above-ground biomass (AGB) using predicted canopy cover, environmental, disturbance, and remotely-sensed predictor variables. The solid line indicates 1:1 correlation between observed and estimated values. When each model was applied across the entire study area, AGB ranged from 61.9 to 67.3 t/ha and total AGB ranged from 64.3 million tonnes to 69.9 million tonnes (Table 4.4). For comparison, the null model estimate and the simple design-based estimate were both within this same range. All of the model-based methods, including the null model, exhibit confidence intervals that are more precise than the design-based estimate by an order of magnitude. Interestingly, even though the sigmoidal model exhibited the highest mean prediction bias, it had 71  the smallest confidence interval based on bootstrap resampling. Conversely, while the GAM demonstrated the best fit and validation statistics across the three models where variable selection was performed, it had a confidence interval wider than the other models. The sigmoidal model and the exponential model were both more efficient than the null model, while the GAM was less efficient. Since the sigmoidal model had fit and validation statistics similar to the GAM, yet had higher relative efficiency than the GAM, the sigmoidal model was selected to estimate AGB and the corresponding CV for each 20m x 50m pixel in the study area. Table 4.4. Bootstrap results for mean and total above-ground biomass (AGB) with corresponding 95% confidence intervals for four model-based methods, compared to a design-based method. Method Mean  AGB (t/ha) Total  AGB (t) 95% CI for total AGB (half-width) Relative efficiency GAM (1) 67.3 69,961,134 694,260 0.82 Sigmoidal (2) 62.1 64,526,209 477,730 1.20 Exponential (3) 61.9 64,316,094 500,070 1.15 Null model 62.8  65,244,018 573,365 - Design-based 63.0  65,472,297 13,055,331 -  AGB estimations per pixel using the sigmoidal model ranged from 0 to 390 t/ha, where lower values were more common around areas with higher population and higher values were more common in remote areas in the Luangwa River valley (Fig. 4.3). Estimates of CV per pixel generally were highest where AGB was the lowest, with a large area of high CV on the eastern boundary and two small valleys in the west half of the study area. The area of high CV on the eastern boundary is dominated by agricultural, nonforest land use; whereas the two small valleys appear to be dominated by grassland. The average CV per pixel across the whole study area was 31%, with a standard deviation of 14. 72   Figure 4.3. Map of estimated above-ground biomass (AGB) per pixel (0.1 ha) and the Coefficient of Variation per pixel (CV). AGB estimates derived from a sigmoidal model. 4.4 Discussion This research provides encouraging results for using a model-based approach to estimate miombo AGB with NFI data. The usefulness of CC as a predictor variable in estimating AGB was clear. Similar results were found by Hall et al. (2006) for boreal forests, by Wu et al. (2013) for tropical woodlands in Sudan, and by González-Roglich and Swenson (2016) for wooded savannas in Argentina. Hall et al. (2006) and González-Roglich and Swenson (2016) used CC as the percentage of the sky blocked by tree crowns over a hemispherical view, whereas Wu et al. used CC as the vertical projection of tree crowns onto the ground. In this research, CCHEMI (i.e., the hemispherical projection of CC) demonstrated a better relationship to AGB for two possible reasons. First, taller trees project more cover in a hemispherical measurement device such as the 73  spherical densiometer (Jennings et al. 1999), and tree height has been shown to be an important covariate in predicting tree-level biomass (Chave et al. 2014). Second, CCVERT (i.e., the vertical projection of CC) relies on DBH-crown radius models (Burrows and Strang 1964, Fuller et al. 1997, Isango 2007) to estimate CC per ground plot; these models are outdated and not localized in this research (Halperin et al. 2016a). Updating DBH-crown radius models may improve the usefulness of CCVERT as a predictor variable in AGB models for the miombo ecoregion and should be an area of active research. This is especially true for less common vegetation types such as mopane woodlands and riparian forests. Improvements in AGB estimates were observed when including a diverse range of ecologically meaningful covariates along with 𝐶?̂?HEMI. In this context, two aspects of moisture stand out. First, soil moisture was characterized by available soil water content (AWC, Hengl et al. 2015) in the sigmoidal and exponential models, or through landscape position via the compound topographic index (CTI) and distance to perennial water in the GAM. Soil moisture has been identified as an underlying factor contributing to miombo vegetation dynamics (Frost 1996, D'Odorico et al. 2007). For example, higher levels of soil moisture have been found under closed canopy miombo woodland compared with open canopy miombo woodland (Campbell et al. 1988). Second, vegetation moisture was captured through the NDMI, an important variable in both the GAM and sigmoidal models. NDMI has been shown to be correlated with tree canopy water content and is competitive with the more well-known NDVI for a variety of forest mapping purposes (Hunt and Rock 1989, Jin and Sader 2005). This study targeted collection of field and remotely-sensed data for the early dry season, which is considered a prime period for linking 74  these observations in the miombo ecoregion. In this timeframe, grass layers have generally senesced while the deciduous tree canopies have not yet fully shed their leaves (Fuller et al. 1997). As such, NDMI appears to be an important vegetation index for characterizing vegetation moisture in tree canopies, and higher canopy moisture could be related to greater levels of AGB (D'Odorico et al. 2007).  The dominant vegetation type in Nyimba District is dry miombo woodland where rainfall averages less than 1,000 mm/yr (Timberlake and Chidumayo 2011). Comparing these results with similar studies in the same vegetation type, Næsset et al. (2016) estimated mean AGB at 51.3 t/ha, Mauya et al. (2015) at 65.8 t/ha, Ryan et al. (2011) at 59.4 t/ha, and Ribeiro et al. (2008) at 70 t/ha, while Carreiras et al. (2013) estimated AGB at 25.7 t/ha. However, Carreiras et al. only studied areas with CC < 50%. The review by Frost (1996) provided an AGB estimate of 55 t/ha for dry miombo of Zambia and Zimbabwe, and in the most recent FAO Forest Resources Assessment for Zambia, the national estimate of AGB was 84 t/ha (FAO 2014a). Across models, these AGB estimates ranged from 61.9 to 67.3 t/ha, where the best performing model and the null model estimated 62.1 and 62.8 t/ha, respectively. Clearly, these estimates of AGB for a study area in the Zambian dry miombo woodland vegetation type are in line with other studies. Comparing accuracy statistics with other published results is challenging, given the number of studies using similar data sources, as well as the use of a wide variety of statistics. Therefore, a comparison of similar studies estimating AGB in other vegetation types of the miombo ecoregion, or other dry forest ecoregions, is necessary. Some researchers presented model fit statistics such as R2 (Wu et al. 2013, Kashindye et al. 2013), while others emphasized RMSE (Mitchard et al. 2013, Solberg et al. 2015). The latter two estimated RMSE at 12.8 t/ha (22% 75  CV) in Mozambique and 40.3 t/ha (78% CV) in Tanzania, respectively. However, fit statistics are known to be optimistic when concerning prediction (Guisan and Zimmermann 2000, Packalén et al. 2012). Using resampling methods to estimate prediction error for AGB in Mozambique, Ryan et al. (2012) and Carreiras et al. (2013) reported RMSPE of 19.6 t/ha and 5.0 t/ha, respectively, where the latter reported a mean CV of 25%. In a study site in Tanzania with vegetation and rainfall similar to Nyimba District, Mauya et al. (2015) and Næsset et al. (2016) reported RMSPE% of 47% and 62%, respectively. These latter studies all employed active sensor (radar or ALS) imagery as the single source for predictor variables.  Research on the use of optical sensor satellite imagery for estimating AGB in the miombo ecoregion is more limited. Næsset et al. (2016) used a linear ordinary least squares model and a stepwise predictor variable selection method. As potential predictor variables, these authors included CC predictions and CC gain or loss from Hansen et al. (2013), raw band values from Landsat 7 ETM+ and Landsat 8 OLI, NDVI calculated from the Landsat data, and square and square root transformations for all of these predictor variables. Based on their stepwise predictor variable selection procedure, squared CC was the only predictor variable chosen. While their finding helps to corroborate the use of CC and squared CC in the sigmoidal model, their model resulted in poor AGB prediction with an RMSPE% of 87%. It is possible that their model using Landsat data was disadvantaged by two aspects. First, they did not include other commonly used vegetation indices (e.g., ARVI, NDMI, NBR, etc.) or texture variables (i.e. standard deviation in a 3 x 3 window around the ground plot), which were found to be useful in this study, as have others (e.g., Lu et al. 2014). Second, the Landsat data that Næsset et al. used differed by at least a year from the time of their ground plot measurements.  76  Two studies that estimated miombo tree volume (m3/ha) using Landsat also provide a useful comparison. Employing cross-validation to compute RMSPE%, Pereira (2006) reported 48% in Mozambique using Landsat and k-Nearest Neighbor imputation, while Gara et al. (2015) averaged 62% using Landsat and nonlinear models. The chosen sigmoidal model has estimates of RMSE% at 47% and RMSPE% at 58%, and these fit and validation statistics are comparable to similar studies using large-scale inventories in the miombo ecoregion using optical remotely-sensed data (Pereira 2006) and active remote sensing data such as ALS or radar (Solberg et al. 2015, Mauya et al. 2015, Næsset et al. 2016). However, more research is needed to compare AGB confidence intervals to other studies because application of a model-based approach appears relatively uncommon in the miombo ecoregion. The use of GAMs and nonlinear models, such as the sigmoidal and exponential models in this study, for predicting AGB has shown promise in a number of studies in different biomes (e.g., Moisen and Frescino 2002, Hall et al. 2006, Anaya et al. 2009, McRoberts et al. 2015, Kattenborn et al. 2015). However, it appears that most approaches to estimate and map AGB in the miombo ecoregion used linear ordinary least squares (Ryan et al. 2012, Kashindye et al. 2013, Mitchard et al. 2013, Solberg et al. 2015, Næsset et al. 2016), despite the nonlinear relationships between forest attributes and remotely-sensed data (Sedano et al. 2008, Banskota et al. 2014, Lu et al. 2014, Gara et al. 2015). GAM’s were found to perform as well as nonlinear models with respect to fit and validation statistics; however, when applied to the population of predictor data the bootstrap confidence interval was less efficient than the null model.  The GAM was protected against overfitting by employing cross-validation in the variable selection process (Guisan and Zimmermann 2000, Hawkins 2004, Packalén et al. 2012); 77  however, the smoothing splines of GAM’s may still perform in an unpredictable manner when applied to population-level data which includes ranges of values not within the observation data (Guisan et al. 2002, Venables and Dichmont 2004). Because of this, nonlinear models may be preferred for several reasons. First, the sigmoidal and exponential models in this study are asymptotic, thereby preventing both negative and unreasonably high estimations (McRoberts et al. 2015). Second, nonlinear models are able to achieve more stable extrapolation at the ends of, and beyond, the ranges of the values of the predictor variables (Venables and Dichmont 2004). Third, judicious application of transformations may improve performance if there appears to be a lack of fit in nonlinear models.  An issue with the nonlinear models used in this study pertains to estimation of the asymptote. Initially a flexible asymptote () was investigated, with α estimated in the model fitting process. For this dataset, α was often estimated at unrealistically high values or the nonlinear models would not converge, both of which may indicate that the data does not support estimating a biologically meaningful asymptote. A flexible asymptote was further investigated by estimating α as a function of the predictor variables in the models (i.e., α0 +  α1𝑥𝑖); however, no improvements in model fit or biologically meaningful values of α were noted. In this context, Burkhart and Tomé (2012 p. 123) stated that if an asymptote cannot be realistically fit with the data at hand, then expert judgement may be employed to fix the asymptote at a reasonable value and continue with model fitting. AGB t/ha for the dry miombo woodland vegetation type tends to be less than 100 t/ha (Frost 1996), although upper limits generally appear to be understudied. Mauya et al. (2015) found an upper limit of 350 t/ha in dry miombo woodland for a study site in southeast Tanzania, while the riparian forest vegetation type may approach 500 t/ha (Chidumayo 78  personal communication 2016). Future modeling efforts should assess practical asymptotic limits for AGB in all miombo vegetation types. While this study made progress in terms of ecologically meaningful variable selection and highlighted model advantages and disadvantages, improvements may be realized in several aspects. First, it was found that CC as the percentage of the sky blocked by tree crowns over a hemispherical view (𝐶?̂?HEMI) was a better predictor for AGB than CC as the vertical projection of tree crowns onto the ground (𝐶?̂?VERT). However, values of CCVERT, used to predict 𝐶?̂?VERT, were derived from DBH-crown radius models (Burrows and Strang 1964, Fuller et al. 1997, Isango 2007) that should be updated for conditions in Zambia (Halperin et al. 2016a). Once updated, 𝐶?̂?VERT may prove to be as useful as 𝐶?̂?HEMI as a predictor of AGB. Second, the highest percent CV pixel values occurred where AGB was predicted to be the lowest, likely because there were very few ground plots with no trees. To accommodate such under-sampled areas, off-grid ground plots could be included along with the NFI ground plots in implementation of the model-based methods. Model-based methods are not restricted to one single sample design as long as the same variable of interest (i.e., AGB) is being sampled according to the same definitions; the main consideration is that sample selection is required to be uninformative in terms of Y (Gregoire 1998).  Third, active sensor remotely-sensed data allow tree heights to be estimated and height is known be correlated with tree-level (Chave et al. 2014) and plot-level biomass (Wulder et al. 2012a). ALS data is a proven source of height estimates, yet remains operationally out-of-reach for many countries (Ribeiro et al. 2012). Testing of this source of predictor variables for use in operational monitoring programs is just now underway in the miombo ecoregion (Mauya et al. 2015, Næsset 79  et al. 2016). Additionally, significant challenges remain in using more readily available radar data (Sinha et al. 2015, Solberg et al. 2015). However, space borne laser scanning systems are under development (Moussavi et al. 2014, GOFC-GOLD 2015) and may prove useful in future applications.  Fourth, a model-assisted approach may lend itself to use of the GAM, balancing the variability in population-level estimation with more conservative inference that takes advantage of the design-based properties of the sample (Opsomer et al. 2007, Næsset et al. 2016). Fifth, there is significant opportunity to explore model-assisted or model-based approaches in a small area estimation context for improving both estimates and inference at spatial scales where traditional design-based approaches are not possible (Magnussen et al. 2014). Lastly, if a map of AGB is a requirement, then use of georeferenced predictor variables is required. However, if a map of AGB is not required, then the null model would provide a viable alternative that could also be applied within a small area estimation framework.     80  Chapter 5: Estimating change in forest area and above-ground biomass 5.1 Introduction Globally, forest resources are important to the economic and cultural well-being of human society, as well as critical for maintaining biodiversity and environmental functions (FAO 2014b, Harrison et al. 2014). The permanent conversion of forest to nonforest and the associated loss of AGB are estimated to make large contributions to global carbon dioxide (CO2) emissions (Baccini et al. 2012, Houghton et al. 2012, Achard et al. 2014, Willcock et al. 2016), and these contributions are generally considered to exacerbate climate change (IPCC 2014). In some areas around the globe, however, research has demonstrated that AGB gains may outweigh losses and nonforest areas may be converting to forest, through either natural processes (e.g., increasing atmospheric CO2; Lewis et al. 2009a, Pan et al. 2011) or by human-induced causes (e.g., fire management, planting trees; Mitchard and Flintrop 2013, Payn et al. 2015). The knowledge of changes in forest area and AGB varies widely by continent and by biome (Houghton 2005, Köhl et al. 2015, Liu et al. 2015). The forests of Africa are of particular interest as recent findings demonstrate that they constitute a net sink of CO2 (Lewis et al. 2009b, Ciais et al. 2011, Valentini et al. 2014) despite losses in forest area (Achard et al. 2014).  In Africa, dry forests and woodlands cover more area than humid forest (Achard et al. 2014) and the miombo ecoregion of southern Africa encompasses the largest dry forest area on the continent. A large variety of vegetation formations are found in the miombo ecoregion, including several woodland types, wooded grasslands, bushlands and thickets (Chidumayo and Marunda 2010, Timberlake et al. 2010). The miombo ecoregion derives its name from the miombo 81  woodland, a type of dry forest that covers 2.4 million km2 and is the most widespread vegetation formation in the ecoregion (Timberlake et al. 2010). Miombo woodlands are typified by large heterogeneity in species composition and structure (i.e., tree height, CC), ranging from very sparse to dense CC (Frost 1996).  The spatiotemporal distribution of miombo woodlands arises from the complex interaction of both abiotic factors (i.e., climate, geomorphology, soil properties, and fire) and biotic factors (i.e., anthropogenic activities, tree-grass interactions, and grazing or browsing by large herbivores) (O’Brien 1993, Belsky 1995, Misana et al. 1996, Frost 1996, Scholes and Archer 1997, Chidumayo 2002, Sankaran et al. 2005, Trouet et al. 2010, Shekede et al. 2016).  Currently, anthropogenic activities are seen as the main biotic factor responsible for decreasing miombo woodland area and AGB (Grogan et al. 2013, Ryan et al. 2014b, Mayes et al. 2015, Watson et al. 2015, Jew et al. 2016). These activities are generally related to removal of tree cover for settlements, agriculture, and fuelwood production (Chidumayo and Marunda 2010). Agricultural land uses have been historically dominated by small-scale shifting agriculture (Timberlake and Chidumayo 2011), yet these agricultural land uses are transitioning to permanent cultivation (Timberlake et al. 2010). Even though fire is an abiotic factor, it is generally considered an anthropogenic activity in miombo woodlands which has been used as a tool for hunting, land clearing, and other socio-cultural purposes for hundreds, if not thousands, of years (Frost 1996). The frequency, timing, intensity, and severity of fire can be related to decreasing miombo woodland area and AGB by increasing tree mortality and promoting grass growth. Conversely, increasing miombo woodland area and AGB can be related to fire control, abandonment of agricultural lands, and land management (Frost 1996, Chidumayo 2002, 82  Timberlake et al. 2010, Mitchard and Flintrop 2013). Abiotic factors (i.e., rising atmospheric CO2 and increasing rainfall) have been linked with increasing miombo woodland area and AGB as well, in addition to the interaction of biotic and abiotic factors (Ciais et al. 2011, Sankaran et al. 2008, Timberlake et al. 2010). Improvements in monitoring are needed to quantify the change in miombo woodland area and AGB due to combined biotic and abiotic factors (Gumbo and Chidymayo 2010, Ribeiro et al. 2015), especially since research indicates the possibilities for both miombo woodland area loss (Mitchard and Flintrop 2013, Mayes et al. 2015, Jew et al. 2016) and increase in AGB (Lewis et al. 2009b, Ciais et al. 2011, Valentini et al. 2014). Monitoring programs are generally conducted under a stock-difference approach or a gain-loss approach (IPCC 2006). The stock-difference approach commonly relies on ground plot data collected at two or more points in time to make an estimate of change, often within the context of an NFI (GOFC-GOLD 2015). As an example of this approach in the miombo ecoregion, Zambia has implemented an NFI on a systematic grid using double sampling for stratification (GRZ 2016). When ground plot remeasurement data is available, the stock-difference approach will allow estimation of changes with well-known statistical properties (Schreuder et al. 1993).  The gain-loss approach focuses on estimating the balance between forest gains and forest losses, and generally relies on process-based inputs (IPCC 2006). For example, Canada implements a gain-loss approach through the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3) that utilizes information on land cover changes, the location and type of forest disturbance events, forest stand-level growth and yield models, and forest inventory data from at least one point in time (Kurz et al. 2009). However, some countries do not have permanent sample ground 83  plots necessary to develop growth and yield models, or forest harvest records detailing disturbance events. Without these data, a simpler version of the gain-loss approach may be applied that relies on a time series of forest/nonforest maps and a one-time forest inventory (GFOI 2013). Zambia has implemented a simple version of the gain-loss approach to assess change in forest area and AGB from 2000-2010, and 2010-2014 (GRZ 2016). This version of the gain-loss approach allows for future improvements, but currently does not quantify the area of land that changes from nonforest to forest nor the changes in AGB within areas that remain forest during the monitoring period. As an example of more refined monitoring in the miombo ecoregion, Ryan et al. (2012) demonstrated an approach to estimate land cover and AGB change over a three-year period for a 1,160 km2 landscape in central Mozambique. They measured different ground plots in each of the three years, and then combined ground plot data with annual L-band Synthetic Aperture Radar (SAR) satellite imagery. They then created annual regression models, using the SAR data as a source of predictor variables to estimate AGB per pixel per year and total change over the area for the three years. This research was limited by three factors. First, they employed a definition of forest based on the amount of AGB per pixel and not on CC per area as is more commonly used (FAO 2010, 2015). Second, a three-year time series is less likely to capture complex spatiotemporal patterns of change in miombo woodlands (Murwira and Skidmore 2006). Third, it is not common to have ground plot measurements for every year of a monitoring period.  Use of a framework that employs a longer time series with time-invariant models is another option (e.g., Kennedy et al. 2010, Powell et al. 2010). A time-invariant model is a model created 84  between the dependent data and a set of predictor variables for one year, and then applied to a time-series of predictor variables. Time-invariant models assume that there has been no change over time in the model parameters (i.e., regression coefficients) or the underlying properties of the predictor variables (Magnussen et al. 2015b) and have been successfully applied in wildlife population studies (Townsend et al. 2016), hydrological studies (Jones et al. 2016), and land use/land cover change (Andersson et al. 2016). Application of multiple, annual regression models (Ryan et al. 2012) or one time-invariant model (Powell et al. 2010) are both considered as stock-difference approaches because the quantification of change relies on estimates at specific points in time and not on process-based inputs as in a gain-loss approach.  Landsat satellite imagery is a standardized source of remotely sensed data (Claverie et al. 2015), making it suitable for use in time-invariant models. For example, Powell et al. (2010) implemented time-invariant models to generate annual AGB estimates for a 20-year period by combining NFI ground plot data for the United States with a Landsat time series. They then used time series segmentation to quantify AGB changes. However, this segmentation method may identify many small false positive changes (Kennedy et al. 2010), and does not lend itself well to an estimate of total change (McRoberts et al. 2014b). Other time series change analysis methods are also available, such as classification and clustering (Ratanamahatana et al. 2010). Classification is a well-known method for developing land cover maps; however, time series classification that applies to changes in forest area and AGB is becoming more common (Petitjean et al. 2012, Gómez et al. 2016, Lu et al. 2016).  In this context, there are two objectives for this chapter. First, the first objective aims to develop a framework to estimate annual land cover change from 2000 to 2015 using time series 85  classification for a study area within the miombo woodlands of Zambia. Land cover is explicitly defined in terms of forest or nonforest, along with biologically meaningful categories of change between forest and nonforest. The second objective aims to employ a stock-difference approach to estimate total AGB change within this time period, AGB change within categories of land cover, annual AGB change per hectare, and estimates of uncertainty. To achieve these objectives time-invariant models of CC and AGB were built and applied to a 16-year annual time series of Landsat imagery for a ~1 million ha landscape in eastern Zambia. Land cover changes and AGB changes within categories of land cover change were then assessed. With this information, the following specific questions were addressed: 1) How much and what type of land cover change occurred?; 2) When and where did land cover change occur?; 3) What is the total change in AGB, the change in AGB within categories of land cover, and the annual change per hectare?; and, 4) Where did AGB change occur? This framework was designed for broad technical application and to enhance understanding of change in woodland area and AGB as it applies to a monitoring context in the miombo ecoregion. 5.2 Methods 5.2.1 Ground plot data For estimating land cover and AGB change, three attributes per plot were required. First, percent CC measures that use a vertical projection of tree crowns onto the ground (CCVERT) were necessary to estimate forest cover consistent with the forest area definition in Zambia (GRZ 2016). CCVERT per plot was obtained by employing tree diameter-crown radius models described in Halperin et al. (2016a). Second, AGB per plot was obtained by calculating tree level AGB 86  using allometric models from Chidumayo (2012), summing AGB for all trees within each plot, and converting plot level AGB values into metric tonnes per hectare (t/ha). Third, another form of percent CC that uses a hemispherical projection of tree crowns onto the ground (CCHEMI) was necessary, as this has been shown to be an important predictor variable in miombo AGB models (Halperin et al. 2016b). CCHEMI per plot was obtained using the spherical densiometer measurements. 5.2.2 Landsat time series Time series of Landsat satellite imagery have been used to track changes in forest attributes in several studies (Powell et al. 2010, Kennedy et al. 2010). However, in the tropics creation of a Landsat time series faces additional challenges. The United States Geological Survey (USGS) Landsat archive carries fewer available images for areas outside North America, not all available images are geometrically corrected (i.e., Level 1 Terrain Corrected or L1T), cloud cover is more persistent, and phenology is more complex (Kovalskyy and Roy 2013, Ryan et al. 2014a, Hostert et al. 2015). Therefore, image compositing methods (Hermosilla et al. 2015) which utilize all available imagery that meet certain criteria are necessary.  All available Landsat imagery for the months of April/May (Julian dates 90 to 150) and for each year from 1990 – 2015 (Fig. 5.1) was downloaded from the United States Geological Survey archive covering the study area (Path 170/Row 70, Path 171/Row 70). Only imagery geometrically processed to L1T (Masek et al. 2006, Wulder et al. 2012b) was downloaded. Missing images were common during the 1990’s and the year 2000. Therefore, time series analysis focused on 2000 – 2015, utilizing imagery from 1998 and 1999 to augment imagery 87  from the base year of 2000. Each image was clipped to the study area; and major rivers and roads were masked out by on-screen digitizing. For each image, contaminated pixels (pixels dominated by clouds and cloud shadows identified in the Landsat Fmask data, Zhu and Woodcock 2012) were removed, as well as any pixel within a 50 pixel buffer around contaminated pixels (White et al. 2014). Radiometric normalization was then performed on each image utilizing the IRMAD algorithm (Canty 2014), using 2013 Landsat 8 OLI Path 170/Row 70 Julian date 125 as the base image.   Figure 5.1. Historical Landsat satellite image availability for the months of April and May from 1990 – 2015, covering Path 170/Row 70 (a) and Path 171/Row 70 (b). Landsat 5 Thematic Mapper ;       Landsat 7 Enhanced Thematic Mapper Plus, Scan Line Corrector on  ; Landsat 7 Enhanced Thematic Mapper Plus, Scan Line Corrector off ; Landsat 8 Operational Land Imager (OLI)  .     Contaminated pixels may also be caused by random noise, causing the spectral values to ‘spike’ (Kennedy et al. 2010). Therefore, contaminated pixels were identified and removed using the Kennedy et al. ‘spike’ detection algorithm. One composite image was created for each year by 88  taking the maximum value for each pixel and for each band within the two month period for every year (Holben 1986). Pixels with missing values were infilled by assigning the pixel value of the previous or the following year. Infilled pixels used the previous year if the standard deviation of the previous two years was lower than the standard deviation of the following two years, or vice versa, following recommendations by Hermosilla et al. (2015). All Landsat data processing steps were completed using Python 2.7.7 and the Geospatial Data Abstraction Library (GDAL 2016). 5.2.3 Time-invariant models Time-invariant models for CCVERT, CCHEMI, and AGB were developed using the methods of Halperin et al. 2016a and 2016b. The same predictor variables identified in previous models (Chapters 3 and 4) were not used because the previous models were based on June-July imagery and this time period in the miombo ecoregion can be impacted by burning (i.e., haze, burn scars; Archibald et al. 2010). As possible predictor variables, 24 remotely sensed predictor variables from the 2013 Landsat image composite (Table 3.3) and 12 environmental predictor variables (Table 3.6) were utilized. The environmental predictor variables are not expected to change within the time period of this study (McNicol et al. 2015). With these predictor variables, time-invariant models for three forest attributes were built: CCVERT (binomial GAM, Eq. 3.5), CCHEMI (binomial GAM, Eq. 3.5), and AGB (sigmoidal model, Eq. 4.2). For all models, predictor variable selection was conducted with a genetic algorithm (GA, Scrucca 2013), using the same settings and fitness criteria described in Chapter 3.  89  The chosen sigmoidal model in Chapter 4 was found to have a Mean Prediction Difference (MPD; i.e., the average difference between observed and predicted values from cross-fold validation) of 3.2 t/ha, indicating possible model heteroscedasticity. One way this can be addressed is by fitting a weighted regression that employs generalized nonlinear least squares (GNLS, Pinheiro and Bates 2000). After predictor variable selection the sigmoidal model was fit using an identity variance function, which uses a constant variance for all errors and is asymptotically equivalent to nonlinear least squares (Pinheiro et al. 2016).  Two other commonly used variance functions (power and exponential) were then investigated. The results of these model options were compared using fit and validation statistics (RMSE, Eq. 3.6; RMSE%, Eq. 4.4; MD, Eq. 3.7) and validation statistics (RMSPE, Eq. 3.1; RMSPE%, Eq. 4.4; MPD, Eq. 3.8) and the Akaike information criterion (AIC), where AIC = -2log(L) + 2 p. L is the estimate of the likelihood function for the model and p is the number of model parameters; smaller values of AIC indicate better performance (Wang et al. 2007). Tests revealed that the power function was not a good fit; therefore focus was placed on the identity and exponential functions. 5.2.4 Land cover change Monitoring programs to assess changes in forest area and AGB apply to all managed lands (IPCC 2006). All land in the study area is under some form of management regime (i.e., customary or government; Gumbo et al. 2016); therefore, the following definitions apply to the entirety of the land in the study area and the changes in forest area and AGB therein. The general definition of forest meets minimum CC, tree height, and area thresholds and also includes land exhibiting a temporary loss of CC (IPCC 2006, GOFC-GOLD 2015). As such, woodland and forest types in the miombo ecoregion meet these general criteria and are considered as forest. 90  From a satellite remote sensing perspective, the temporal aspect of this definition is challenging to apply (GFOI 2013, GOFC-GOLD 2015). In Zambia, the operational definition of forest does not include this aspect of temporary loss and focuses on identifying forest as tree CC > 10% over 0.5 ha, while the operational definition of nonforest includes all land cover types which specifically do not meet this definition of forest (i.e., agriculture, grassland, or settlement) (GRZ 2016). These operational definitions of forest and nonforest were also employed in this research. From these definitions, land cover was defined as a state of forest or nonforest and land cover change as the changes between forest and nonforest. Deforestation and afforestation/reforestation (A/R) are defined as human-induced land cover changes between forest and nonforest (IPCC 2006, GOFC-GOLD 2015). For the purposes of this research, deforestation and A/R were considered as a permanent change within the 16-year time series. To create a forest/nonforest time series for assessing land cover change using these definitions, the model for 𝐶?̂?VERT was applied to each year of the time series. Each pixel was resampled to 0.5 ha, using the sum of the 𝐶?̂?VERT estimates in a 0.5 ha moving window. Each pixel was then classified as forest (1) if the sum was > 10% or nonforest (0) if the sum was < 10%.  When using a 16-year time series of binary data, there are 216 possible transitions (i.e., trajectories) between forest and nonforest. Lambin (1997) proposed that generic trajectories can be defined in such a way as to meaningfully characterize land cover changes (e.g., Petit et al. 2001). Trajectory category definitions (Table 5.1) were created to characterize biologically meaningful change in land cover categories, based on forest-nonforest trajectories from 2000 to 2015. To characterize areas that exhibited no land cover change, two land cover categories were 91  established: forest remaining forest and nonforest remaining nonforest (Table 5.1). All other trajectories represent some type of land cover change which is either permanent or persistent.  Table 5.1. Land cover change status, classes and categories for trajectory category definitions representing forest/nonforest trajectories for the 2000 – 2015 time period. Each year in a trajectory category definition is represented by a 1 (forest), or a 0 (nonforest). Land cover change status Land cover class Land cover category Trajectory category definition No change Forest  Forest remaining forest 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1  Nonforest  Nonforest remaining nonforest 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0  Permanent change Afforestation/ Reforestation (A/R) 2003 0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1 2006 0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1 2009 0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1 2012 0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1 Deforestation 2003 1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0 2006 1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0 2009 1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0 2012 1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0 2014 1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0 Persistent change 2 & 3-year cyclical change Two-year cycle 1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0 Three-year cycle 1,1,1,0,0,0,1,1,1,0,0,0,1,1,1,0 5-year cyclical change Forest – nonforest – forest 1,1,1,1,1,1,0,0,0,0,0,1,1,1,1,1 Nonforest – forest – nonforest 0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0  Permanent change represents the conversion from nonforest to forest (A/R) or the conversion from forest to nonforest (deforestation) (Table 5.1). Permanent change patterns can be represented in annual time steps; however, this would create 30 separate land cover categories which would complicate interpretation and generalization of results. In a study of forest cover 92  changes in the United States, Jin and Sader (2005) recommended that five-year intervals were likely sufficient for assessing permanent change. Trajectory category definitions of permanent land cover change at three-year intervals were utilized as a compromise between annual change and five-year change. Continuity in land cover was ensured by establishing at least three years of the opposite land cover category at either end of the time series (i.e., three years of forest cover before deforestation is possible). Because recent forest cover loss due to agricultural expansion in the study area has been noted (Gumbo et al. 2016), a deforestation category for 2014 was established as an exception to capture these recent changes.  Persistent change represents a multiple stable-state, where biotic factors (i.e., anthropogenic activities, wildlife use) interact with abiotic factors (i.e., fire, soil properties, climate) to drive cyclical transitions between forest and nonforest (Belsky 1995, Frost 1996, Misana et al. 1996, Scholes and Archer 1997, Grogan et al. 2013). Two land cover classes were established representing two or three-year cyclical change and five-year cyclical change, with four trajectory category definitions characterizing cyclical change from forest to nonforest at two-year, three-year, and five-year cycles (Table 5.1).  With these trajectory category definitions, time series classification was employed to assign each pixel trajectory to a land cover category. Nearest Neighbor classification was implemented to determine the land cover category for each pixel, based on the Dynamic Time Warping (DTW) distance between each pixel’s trajectory and the trajectory category definitions. The DTW distance measure was utilized as it is a common choice for comparing differences between satellite image derived trajectories (Petitjean et al. 2012, Gómez et al. 2016, Lu et al. 2016). The DTW measure was computed between each pixel’s trajectory and each trajectory category and 93  classified the pixel’s land cover category based on the minimum DTW distance. In the event of a tie in minimum distance between two or more land cover categories, a category was assigned based a spatial weighting scheme using a 9 x 9 pixel neighborhood. After classification, calculations were made to describe the area of each land cover category, class, and status, and the net change in forest area as increase from A/R minus the decrease from deforestation (Coulston et al. 2013). Time series classification was implemented using Python 2.7.7, and the area calculations using the R raster package (Hijmans et al. 2015). 5.2.5 AGB change Change in AGB per land cover category was estimated by using a bootstrap resampling method with post-stratification (Efron and Tibshirani 1986, Tipton et al. 2013, Magnussen et al. 2015c). For this purpose, the sigmoidal model was refit using 500 bootstrap resamples, and each of the 500 resamples was applied to the entire study area for each year using the Landsat time series and environmental predictor variables. For each bootstrap resample, the pixels were rescaled to a 0.5 ha, where values of the pixels were calculated by using the sum of the 𝐴𝐺?̂? estimates in a 0.5 ha moving window and the units of the resampled pixels were in tonnes (t) per pixel. Total annual change in AGB per land cover category was calculated by:  ∆𝐴𝐺?̂?𝑗𝑏𝑜𝑜𝑡=  ∑ ∑ (( 𝐴𝐺?̂?1𝑖𝑗𝑘 −  𝐴𝐺?̂?2𝑖𝑗𝑘 ) + ⋯ +  ( 𝐴𝐺?̂?15𝑖𝑗𝑘 −  𝐴𝐺?̂?16𝑖𝑗𝑘 )) 500⁄500𝑘=1𝑛𝑗𝑖=1, (5.1)  where 𝐴𝐺?̂?1𝑖𝑗𝑘  is the predicted AGB using bootstrap resample 𝑘 for pixel 𝑖 in land cover category 𝑗 in 2000 and 𝐴𝐺?̂?16𝑖𝑗𝑘 in 2015 and 𝑛𝑗  is the number of pixels in land cover category 𝑗. Total and subtotal changes as the sum of ∆𝐴𝐺?̂?𝑗𝑏𝑜𝑜𝑡 across land cover categories by land cover class and 94  status were calculated. Change per hectare per year (t/ha/yr) was calculated by dividing ∆𝐴𝐺?̂?𝑗𝑏𝑜𝑜𝑡 for each category by its corresponding number of hectares, then the number of change periods (15). A model-based variance estimate for ∆𝐴𝐺?̂?𝑗𝑏𝑜𝑜𝑡 was calculated by:  Var̂ (∆𝐴𝐺?̂?𝑗𝑏𝑜𝑜𝑡) =  ∑ (∆𝐴𝐺?̂?𝑖𝑗𝑘−   ∆𝐴𝐺?̂?𝑖𝑗𝑏𝑜𝑜𝑡)2𝑛𝑗𝑖=1500 − 1⁄ . (5.2)  This variance estimate was used to calculate a 95% confidence interval (CI) for change in AGB per land cover category, class, and status, and report the half-width of the CI for each.  5.2.6 Spatial context of change Land cover change by distance from improved roads was quantified as previous research has indicated that change may be impacted by road proximity (Helmer et al. 2008, Ahrends et al. 2010). Some studies have found distances of 4 – 8 km to roads to be related to land cover changes or land use activities in miombo woodlands (Mertens and Lambin 1997, Sprague and Oyama 1999). However, these studies looked at bitemporal change in comparison to distance to improved and unimproved roads mapped for each year. As this research examines change over 16 years, it would be necessary to map the network of unimproved roads for the same time steps in order to assess detailed relationships between land cover, AGB and unimproved roads. As a compromise, five distance classes to improved roads were created which did not change during the time period (< 1.0 km, 1.0 to 4.9 km, 5.0 to 9.9 km, 10.0 to 14.9 km, and >= 15.0 km), similar to distance classes identified by a workshop of miombo researchers for a study of land cover change modeling in Tanzania (Swetnam et al. 2011). The area (ha) per land cover category was quantified for each distance-from-improved-road class and the change in AGB (t) for each 95  combination of distance-from-improved-road class and land cover class. For the latter, the full model estimator (Eq. 4.2) was applied, not the bootstrapped estimator (Eq. 5.1) because the objective of this spatial analysis is to provide context for the results and not to assess uncertainty of AGB change within distance classes. 5.3 Results 5.3.1 Time-invariant models Both GAM’s exhibited fit and validation statistics with low MD between observed and predicted values (Table 5.3, Fig. 5.2). In the GAM for CCVERT, coefficients for slope and carbon were estimated as fixed, while the other four predictor variables were all estimated as smoothed functions. In the GAM for CCHEMI, all six predictor variables were estimated as smoothed functions. Slope and the NIR band were chosen for both CC models in the GA, and each CC model had one soils variable (carbon for 𝐶?̂?VERT and crs.frg for 𝐶?̂?HEMI).  Table 5.3. Goodness-of-fit and validation statistics for two binomial GAM’s to estimate percent canopy cover based on crown radius-DBH models (𝑪?̂?𝐕𝐄𝐑𝐓) and measurements using a spherical densiometer (𝑪?̂?𝐇𝐄𝐌𝐈), and two sigmoidal models using different variance functions (identity and exponential) to estimate above-ground biomass (𝐀𝐆?̂? t/ha).  Fit statistics Validation statistics  Attribute RMSE (3.6) RMSE% (4.5) MD (3.7) RMSPE (3.1) RMSPE% (4.4) MPD (3.8) Predictor variables* 𝐶?̂?VERT  7.6 44% 0.0 9.4 54% 0.0 Carbon, Slope, NIR.texture2, NIR, SWIR1.texture2, SWIR2 𝐶?̂?HEMI   11.3 21% 0.0 16.3 31% 0.1 Crs.frg, Slope, ARVI, Blue, NBR, NIR AGB̂ ident.  30.9 48% 1.7 35.7 56% 1.5 𝐶?̂?HEMI, Crs.frg, pH, Slope, NIR.tex, SWIR2 AGB̂ expo.  31.4 49% 0.8 36.4 57% 0.2 * See Table 3.3 and 3.5 for descriptions of predictor variables.  96   Figure 5.2. Observed vs. predicted graphs from models for two forms of percent canopy cover (CCVERT, CCHEMI) and above-ground biomass (AGB). GA predictor variable selection for the AGB sigmoidal model included 𝐶?̂?HEMI while soils variables were also important to the estimation. The sigmoidal model fit using GNLS with an identity variance function resulted in positive MD and MPD and an AIC of 627. However by using an exponential function based on 𝐶?̂?HEMI, the MD was reduced by one half, MPD was close to zero, and AIC was 625. Given these results, the AGB model was chosen with an exponential variance function based on 𝐶?̂?HEMI. 5.3.2 Land cover change Forest/nonforest maps were created for each year across the time series, based on the time-invariant model for CCVERT (Fig. 5.3). Using these data, 72.3% of the land exhibited no change followed by persistent change (17.1%) then permanent change (10.6%) (Table 5.2). The largest land cover class is forest (58.2%), four times greater than the next largest class (nonforest, 14.1%). Considering persistent change, the area of 2 & 3-year cyclical change is more than four times greater than 5-year cyclical change (13.9% vs. 3.1%). Within permanent change, deforestation accounted for twice as much area as A/R (7.2% vs. 3.4%). Overall, there was a net 97  loss in forest area of 38,847 ha, equivalent to 2,590 ha/yr. Deforestation covered between 0.8% and 2.3% of the study area over time. Most deforestation occurred in the 2003 and 2012 time periods, averaging approximately 23,000 ha which were twice as much as any other time period. A/R ranged from 0.7% to 1.1% per year, with a slight decrease over time. Considering both persistent and permanent change, area in the 2 & 3-year cyclical change categories was at least twice as large as any other change category.    98   Figure 5.3. Forest/nonforest maps for each year between 2000 and 2015 for Nyimba District, Eastern Province, Zambia.     99  Most of the study area is > 15.0 km from an improved road (42.6%), within which the largest land cover category was forest remaining forest (30.5%) (Table 5.2). Within nonforest remaining nonforest, there was more area in the 1.0 to 4.9 km distance class (4.3%) than the other distance classes, though there was also a sizeable amount in the > 15.0 km distance class (3.5%). Area of A/R categories in the 1.0 to 4.9 and > 15.0 km distance classes (1.1% and 0.9%, respectively) was greater than the other distance classes. Percent total area of deforestation by distance class occurred fairly evenly among distance classes > 1.0 km and more than half of deforestation occurred in the 2003 and 2012 periods. The largest amount of deforestation occurred in the 5.0 to 9.9 km distance class. More than one half of all persistent change occurred > 10 km from improved roads, dominated by 2 & 3-year cyclical change in the > 15.0 km distance class (5.2%) (Table 5.2) Much of this persistent change occurred in remote valleys on the west side of the study area, the northern half of the Luangwa River valley and a concentrated area on the northern half of the east side of the study area (Fig. 5.4). Deforestation generally occurred at forest edges closer to population centers or improved roads. Further, deforestation over time and space is apparent in Fig. 5.5a, where it progressed further from the sub-district capital over time. Visual interpretation of Landsat image composites (Fig. 5.5b and 5.5c) from 2000 and 2015 demonstrates this change from forest to nonforest.      100  Table 5.2. Area of land cover categories characterizing annual forest/nonforest trajectories between 2000 and 2015, by distance to improved roads. Land cover change status Land cover  class Land cover  category           Area (ha)              Percent of total area   Distance to improved road class (km)   Distance to improved road class (km)  < 1.0 1.0 to 4.9 5.0 to 9.9 10.0 to 14.9 > 15.0 Total < 1.0 1.0 to 4.9 5.0 to 9.9 10.0 to 14.9 > 15.0 Total No change Nonforest  NRN 17,381 44,465 28,899 18,240 36,551 145,535 1.7 4.3 2.8 1.8 3.5 14.1 Forest FRF 22,849 98,749 95,719 69,397 314,657 601,371 2.2 9.6 9.3 6.7 30.5 58.2 Subtotal   40,230 143,214 124,618 87,636 351,208 746,906 3.9 13.9 12.1 8.5 34.0 72.3 Permanent change Afforestation/ Reforestation  (A/R) 2003 1,118 3,305 1,959 1,203 3,362 10,946 0.1 0.3 0.2 0.1 0.3 1.1 2006 696 1,993 1,245 747 3,254 7,934 0.1 0.2 0.1 0.1 0.3 0.8 2009 1,314 3,442 1,871 800 1,258 8,684 0.1 0.3 0.2 0.1 0.1 0.8 2012 1,068 2,857 1,476 656 1,656 7,712 0.1 0.3 0.1 0.1 0.2 0.7 Subtotal 4,195 11,596 6,550 3,405 9,530 35,275 0.4 1.1 0.6 0.3 0.9 3.4 Deforestation 2003 2,102 6,043 6,148 4,089 3,928 22,308 0.2 0.6 0.6 0.4 0.4 2.2 2006 950 2,697 2,498 1,143 1,295 8,583 0.1 0.3 0.2 0.1 0.1 0.8 2009 923 2,645 2,503 1,233 1,452 8,755 0.1 0.3 0.2 0.1 0.1 0.8 2012 1,826 4,727 6,411 4,954 5,992 23,909 0.2 0.5 0.6 0.5 0.6 2.3 2014 677 1,840 2,059 1,457 4,536 10,568 0.1 0.2 0.2 0.1 0.4 1.0 Subtotal 6,478 17,950 19,618 12,875 17,202 74,122 0.6 1.7 1.9 1.2 1.7 7.2 Subtotal   10,672 29,546 26,168 16,280 26,731 109,397 1.0 2.9 2.5 1.6 2.6 10.6 Persistent change 2 & 3-year cyclical change 2YC 5,844 18,760 15,907 14,562 37,908 92,981 0.6 1.8 1.5 1.4 3.7 9.0 3YC 3,390 10,976 11,628 9,010 15,855 50,857 0.3 1.1 1.1 0.9 1.5 4.9 Subtotal 9,234 29,736 27,535 23,571 53,763 143,838 0.9 2.9 2.7 2.3 5.2 13.9 5-year cyclical change FNF 1,554 4,408 3,627 2,117 3,361 15,066 0.2 0.4 0.4 0.2 0.3 1.5 NFN 1,222 4,155 3,755 2,971 5,250 17,351 0.1 0.4 0.4 0.3 0.5 1.7 Subtotal 2,775 8,562 7,381 5,088 8,611 32,417 0.3 0.8 0.7 0.5 0.8 3.1 Subtotal   12,009 38,298 34,916 28,659 62,373 176,254 1.2 3.7 3.4 2.8 6.0 17.1 Total     62,911 211,057 185,702 132,575 440,312 1,032,556 6.1 20.4 18.0 12.8 42.6 100.0 * NRN: Nonforest remaining nonforest; FRF: Forest remaining forest; 2YC: Two-year cyclical change; 3YC: Three-year cyclical change; FNF: Forest-nonforest-forest; NFN: Nonforest-forest-nonforest.  101   Figure 5.4. Land cover class per 0.5 ha pixel in, line symbology as in Fig. 2.2.    102   Figure 5.5. (a) Close-up of deforestation by year. Landsat false color image composite (SWIR1, NIR, Red) for the year 2000 (b) and the year 2015 (c). Line symbology as in Fig. 2.2. 5.3.3 AGB change From 2000 to 2015, total AGB increased by ~6 million t (bootstrapped confidence interval of +/- 215,049 t) (Table 5.3), representing an annual change of 0.4 t/ha/yr. Within the land cover status of no change, there was an increase of more than 7.1 million t and the vast majority of AGB increase was associated with the forest remaining forest land cover category (~7 million t, 0.8 t/ha/yr). There was a relatively small increase in AGB associated with the nonforest remaining nonforest land cover category (98,400 t). Within the land cover status of persistent change, 5-103  year cyclical change demonstrated an increase (104,811 t); however, this was more than offset by AGB decrease associated with 2 & 3-year cyclical change (-270,579 t). There was a decrease in AGB in both 2 & 3-year cyclical change; however, the decrease associated with the 3-year cycle was 3.5 times more than the 2-year cycle. Within the land cover status of permanent change, the increases in AGB associated with A/R (646,970 t) were more than offset by decreases associated with deforestation (-1,614,518 t). The greatest increases associated with A/R occurred early in the time period (< 2003). All deforestation categories exhibited decreases in AGB, where the largest decreases were in the 2003 and 2012 time periods. Across distance-to-improved-road classes, AGB change by land cover class generally demonstrated either a consistent increase (i.e., forest) or a consistent decrease (i.e., deforestation) (Table 5.4). However, this spatial pattern varied in three ways. First, nonforest demonstrated AGB decrease < 10 km from an improved road and AGB increase > 10 km from an improved road. Second, 2 & 3-year cyclical change demonstrated consistent AGB decrease < 15 km from an improved road and AGB increase > 15 km from an improved road. Third, the 5.0 to 9.9 km distance-to-improved-road class was the only distance class to demonstrate a decrease in total AGB. This can be seen in Fig. 5.6, where AGB decreases are generally further away from improved roads east of the Luangwa River. Patterns of decrease appear to follow small rivers on either side of the Luangwa and the improved road that runs due north of the district capital. AGB decreases were highest within areas identified as deforestation (Fig. 5.7a). Most AGB increases occurred west of the Luangwa River, where the greatest increases were seen in remote mountainous areas (Fig. 5.7a). Visual interpretation of the Landsat image composite from 2015 shows that large AGB decreases occurred in areas currently mapped as nonforest (Fig. 5.7b).  104  Table 5.3. Changes in above-ground biomass tonnes (ΔAGB t) by land cover category, class, and status. Land cover change status Land cover  class Land cover  category ΔAGB (t) 95% CI for change AGB (t) (half-width) ΔAGB       (t/ha/yr) No change Nonforest  Nonforest remaining nonforest 98,400 20,999 0.1 Forest Forest remaining forest 7,022,676 110,132 0.8 Subtotal   7,121,076 111,389 0.8 Permanent change Afforestation/ Reforestation  (A/R) 2003 210,033 2,797 1.3 2006 144,744 1,977 1.2 2009 149,912 2,052 1.2 2012 142,282 1,955 1.2 Subtotal 646,970 8,726 1.2 Deforestation 2003 -442,866 10,474 -1.3 2006 -227,146 4,565 -1.8 2009 -264,384 5,594 -2.0 2012 -458,660 11,873 -1.3 2014 -221,462 5,571 -1.4 Subtotal -1,614,518 37,877 -1.5 Subtotal   -967,548 29,928 -0.6 Persistent change 2 & 3-year cyclical change Two-year cyclical change -57,447 22,374 -0.1 Three-year cyclical change -213,131 10,308 -0.3 Subtotal -270,579 32,154 -0.1 5-year cyclical change Forest-nonforest-forest 72,364 1,880 0.3 Nonforest-forest-nonforest 32,446 2,493 0.1 Subtotal 104,811 4,125 0.2 Subtotal   -165,768 36,027 -0.1 Total AGB change  5,987,760 215,049 0.4     105   Table 5.4. Above-ground biomass change (ΔAGB t) by land cover class, status and distance to improved road class. Land cover change status Land cover  class ΔAGB (t) Distance to improved road (km) class < 1.0 1.0 to 4.9 5.0 to 9.9 10.0 to 14.9 > 15.0 No change Nonforest  -17,201 -120,203 -109,865 58,761 260,026 Forest 209,148 782,167 583,674 419,541 5,123,852 Subtotal 191,948 661,964 473,809 478,302 5,383,877 Permanent change Afforestation/ Forestation 85,791 225,775 104,128 69,692 186,208 Deforestation -175,464 -493,321 -533,463 -233,396 -248,302 Subtotal -89,674 -267,546 -429,335 -163,705 -62,095 Persistent change 2 & 3-year cyclical change -36,899 -176,860 -225,094 -40,141 168,220 5-year cyclical change 9,491 6,398 1,533 21,186 66,400 Subtotal -27,407 -170,461 -223,560 -18,956 234,620 Total AGB change 74,866 223,957 -179,086 295,642 5,556,402  106   Figure 5.6. Above-ground biomass (AGB) change in tonnes (t) between 2000 and 2015. Line symbology as in Fig. 2.2.  107   Figure 5.7. (a) Close-up of deforestation areas in relation to AGB change. Deforestation area polygons derived from the deforestation land cover class data in Fig. 5.4. (b) 2015 Landsat false color image composite (SWIR1, NIR, Red). Line symbology as in Fig. 2.2.  5.4 Discussion These results demonstrate the utility of estimating annual land cover and AGB changes in a framework that employs time series classification and a stock-difference approach. Permanent land cover changes associated with deforestation and A/R both occurred in the study area, and that the area of deforestation is twice as much as the area of forestation. Similar conclusions regarding land cover change in the miombo ecoregion have been reached by others (Prins and 108  Kikula 1996, Petit et al. 2001, Mitchard and Flintrop 2013, Bodart et al. 2013, Achard et al. 2014, Mayes et al. 2015). However, these findings highlight the opportunity for annual time series data and time-invariant models to quantify changes at spatiotemporal scales more relevant to biotic and abiotic factors underlying the distribution of miombo woodland area and AGB. For example, it is not possible to estimate area of persistent change, and the associated changes in AGB, with a bi-temporal time series. Deforestation is generally driven by multiple, interacting biotic factors which are commonly attributed to anthropogenic activities (Geist and Lambin 2002, Ahrends et al. 2010, Hosonuma et al. 2012). In general, deforestation at any time period closer to improved roads is more likely attributable to direct anthropogenic activity (Cabral et al. 2011, Ryan et al. 2014b). Deforestation at any time period further from improved roads could be attributable to direct anthropogenic activity, land management policies that affect wildlife populations, or an interaction of multiple biotic and abiotic factors (Frost 1996, Baxter et al. 2005, Watson et al. 2015). For example, intense browsing of trees by elephants in remote protected areas opens the woodland canopy, influencing grass growth that leads to more frequent, intense fires and a transition from woodland to wooded grassland (Timberlake et al. 2010).  The largest land cover class exhibiting change was associated with 2 & 3-year cyclical change. Persistent change closer to improved roads is more likely to be aligned with anthropogenic activities related to agriculture (Misana et al. 1996, Grogan et al. 2013). However, more than one-third of persistent change in the study area occurred > 15 km from improved roads. Persistent changes at these distances could also be related to agriculture in the eastern half of the study area, given its higher population density. The western half of the study area was more 109  populated in the past, but since becoming a Game Management Area now has much lower population density (Watson et al. 2015, Gumbo et al. 2016). Persistent change in this region is concentrated in lowland river valleys where a combination of biotic and abiotic factors likely interact, such as successional changes from historic land use patterns, elephant use, fire, seasonal flooding along major rivers, and/or soil nutrient status. These results indicated that increases in AGB within forest remaining forest far outweighed AGB decreases associated with permanent or persistent land cover change. Other broad research in Africa supports this finding (Lewis et al. 2009b, Ciais et al. 2011, Valentini et al. 2014) and comparisons with other studies provide several useful insights. In a review of eight studies from dry forests across Africa, Ciais et al. (2011) found that AGB increased by approximately 0.5 t/ha/yr. This corresponds well with my annual estimate of AGB change (0.4 t/ha/yr) in this research. The increase in AGB associated with A/R was estimated at 1.2 t/ha/yr. In Mozambique, Williams et al. (2008) and McNicol et al. (2015) found the change in AGB on land that had been previously cleared and used for agricultural purposes increased at 1.4 and 1.6 t/ha/yr, respectively. In Zambia, Chidumayo (1990, 2002) analyzed data from re-measured ground plots on sites previously cleared and abandoned, and estimated AGB increases at 1.7 and 1.0 t/ha/yr, respectively. Within forest remaining forest, an increase in AGB was estimated at 0.8 t/ha/yr, whereas Ek (1994) estimated growth as 1.6 t/ha/yr using re-measured ground plots in Tanzania. The differences in these estimates from these others could easily be accounted for by differences in site conditions (Chidumayo and Frost 1996, Chidumayo et al. 1996).  Comparing AGB decrease from deforestation with other studies is challenging since decreases depend on AGB values at the beginning of the time period, local environmental conditions, and 110  the types of local disturbances. As a general comparison, Ahrends et al. (2010) and Ryan et al. (2012) estimated AGB decreases associated with deforestation at -1.6 and -0.9 t/ha/yr, respectively. AGB decrease associated with deforestation was estimated at -1.5 t/ha/yr, which is biologically reasonable given land cover change from forest to nonforest. The decrease in AGB associated with deforestation was greatest in the 5.0 – 9.9 km distance from improved road class, and this is the only distance class with a total decrease in AGB.  This could indicate a frontier of deforestation, which is critical to monitor as a hotspot of opportunities for improving ecosystem services, land management, and biodiversity conservation (Salerno 2016, Ryan et al. 2016).   AGB changes within land cover classes of nonforest and 2 & 3-year cyclical change exhibited a complex spatial pattern, shifting from decreasing AGB to increasing AGB as the distance from roads increases. The AGB decrease at distances closer to roads is possibly attributable to anthropogenic activities related to agricultural land use, while the AGB increase at distances further from roads is possibly attributable to a similar interaction of biotic and abiotic factors related to the presence of persistent change. However, recently declining elephant populations in the region (Nyirendra et al. 2015), and the associated decrease in elephant browsing, could also lead to increasing AGB (Morrison et al. 2016).  AGB decrease within forest remaining forest is often highlighted as a source of GHG emissions (Houghton 2005, Putz et al. 2008, Ahrends et al. 2010, Federici et al. 2015) and has been shown to occur within miombo woodlands (Jew et al. 2016). However, this research indicates that forest remaining forest in miombo woodlands can demonstrate large AGB gains. All sources and sinks of GHG emissions on managed lands that make up > 5% of an overall GHG emissions profile should be monitored (IPCC 2006). However, it is common for monitoring programs to focus on 111  loss only (Mitchard and Flintrop 2013, GOFC-GOLD 2015) because coarse default methods assume no AGB change on land remaining in the same land cover category during a monitoring period (IPCC 2006). Estimates of both decreases and increases of AGB on all land cover categories of managed lands should be included in GHG monitoring programs.  By using an annual time series of satellite imagery, annual changes can be operationally monitored using land cover categories relevant to biotic and abiotic factors that underlie the spatiotemporal distribution of miombo woodland area and AGB. To further enhance operational GHG monitoring within the miombo ecoregion and to promote understanding of the complex human-environment interactions needed to improve policies, land management, and sustainable use of woodlands in rural communities that depend on natural resources for their livelihoods, the following research areas should be prioritized: 1) refining time-invariant models and validating change with repeat ground plot measurements; 2) developing statistical relationships to identify and attribute specific abiotic and biotic factors related to land cover categories and AGB change; 3) where factors can be identified as anthropogenic, identifying the direct and indirect anthropogenic activities and underlying socio-economic factors behind land cover change and AGB change; and, 4) developing scenarios of forecasted changes in land cover and AGB based on biotic factors, abiotic factors, and the combination of these factors.    112  Chapter 6: Conclusion 6.1 Overall significance and contributions Estimating the amount and distribution of forest area and AGB is crucial for improving land management, designing effective policies, and engaging in international initiatives such as REDD+. The main goal of the research presented in this dissertation was to contribute to the development of methods and approaches to monitor forest attributes, with a focus on AGB and CC. This was accomplished through model-based approaches that incorporated NFI data and multi-source information. While a jurisdictional landscape located in the center of the miombo ecoregion was the focus in this research, overall this research demonstrated that estimating the amount and distribution of forest attributes is enhanced by including ecologically meaningful predictor variables in a model-based framework for estimation. My research has contributed to the knowledge-base of NFI programs and forest monitoring in general by developing and refining approaches to: 1. estimate forest area by using CC models (Chapter 3); 2. estimate total AGB and uncertainty at large and small scales (Chapter 4); and,  3. monitor land cover and AGB using ecologically relevant post-stratification (Chapter 5). Further, the following insights for improving and enhancing estimation of CC, forest area, and AGB for monitoring purposes were generated: 1. freely available optical satellite imagery of medium spatial resolution (Landsat 8 OLI) is consistently better than costly satellite imagery of high spatial resolution (RapidEye); 113  2. a genetic algorithm provides an efficient way to perform predictor variable selection, especially when using cross-validated accuracy statistics as model selection criterion; 3. predictor variables related to soil moisture and texture of remotely sensed vegetation indices are important in models for both CC and AGB; 4. semi-parametric generalized additive models improve upon more well-known parametric or nonparametric methods for CC estimation when used judiciously; 5. accuracy of models that employ freely available multi-source information is competitive with models that use expensive active sensor remotely sensed data; 6. the biophysical relationship between AGB and CC can be exploited to improve AGB estimation; 7. estimating the model-based confidence interval using bootstrap resampling is a viable alternative to the traditional design-based confidence interval; 8. the first two months of the early dry season should be prioritized when using optical remotely sensed data as a source of predictor variables in forest attribute models in the miombo ecoregion; 9. annual changes in land cover can be operationally monitored using times series classification and land cover categories relevant to biotic and abiotic factors that underlie the spatiotemporal distribution of forest attributes; and 10. an annual time series allows estimation of area under persistent change, which is not possible when change analysis employs two or three time periods.  114  6.2 Limitations Procedures to obtain estimates of forest attributes have become increasingly reliant on model-based or design-based, model-assisted frameworks rather than the traditional design-based frameworks used in the past.  Since the traditional design-based framework has a long history, this approach has been well-studied.  However, the trend towards replacing the traditional design-based framework with model-based approaches is largely because resulting estimators can be more precise, they make use of multi-source information that has become increasingly accessible, and the computational requirements are more feasible; therefore, making model-based approaches more cost-efficient than in the past.  The specific methods evaluated in this research were constructed with multi-source data specific to the study area landscape. Although the recommended approaches can be applied in other regions, applying these specific models outside the range of conditions (i.e., vegetation types, soil characteristics, rainfall, disturbance patterns, etc.) represented in the study region should be avoided.  As noted, multi-source data were used in this research. New ground data were collected using the NFI protocols for Zambia and other data used were obtained from a variety of sources. As a result, improvements are possible. These improvements concern both data and modeling issues. 6.2.1 Data issues Predictor variable data issues (i.e., spatial and/or temporal registration, uncertainties from measurements and/or predictions) may contribute to uncertainty and/or decreased precision in an estimate. Overall data limitations of this research include the following issues. 115  1. Measures of CC (i.e., CCVERT) that specifically relate to a common forest area definition were not part of the NFI measurements protocols. Crown radius from DBH models (DBH-crown models) can be used to fill this data gap, and existing DBH-crown radius models from neighboring regions were used in this research.  Since crowns are affected by a large number of factors that may not all be reflected in the DBH-crown relationship, local data on measured crowns should be collected and used to produce local DBH-crown radius models.     2. The soils data used to obtain the soils predictor variables in this research were extracted from existing map databases based on continental-scale soils models. While these soils data proved vital in forest attribute models over large land areas, improvements in the soils data map accuracy and spatial resolution would likely improve forest attribute modeling results.  3. Since there were no repeated measurements of the NFI ground plots, validating change estimates for Chapter 5 was not possible.  6.2.2 Modeling issues A wide variety of modeling methods were investigated in this research. Specific advantages and disadvantages of the models were discussed in the research chapters. Overall modeling limitations of this research include the following issues.  1. Variances needed to construct model-based confidence intervals were approximated by bootstrap resampling. For some modeling methods (e.g., nonlinear models), there are analytical variance estimators available; however, for other methods (e.g., kNN, GAM) these estimators are not available.  116  2.  Nonlinear models proved useful in estimating AGB, although enhancements are possible. More effort is needed to estimate or choose biologically relevant asymptotes, which is currently difficult given the nature of multi-source datasets and large scale forest inventories. Such enhancements may require improving knowledge about the growth process of specific vegetation types.  6.3 Future research directions The research process is a dynamic enterprise, in continuous development for the purpose of improving our understanding of the world and our role in shaping the future. Forest inventory and monitoring plays a key role in this enterprise and requires essential advancements to provide cost-effective solutions. In order to continue improvements in forest monitoring, future research should prioritize the following areas.  1. New technology in the form of unmanned aerial vehicles outfitted with cameras should be investigated to examine if CC per plot can be reliably and consistently estimated from digital aerial photography.  2. When repeated ground plot measurements are available, examination of both design-based and model-based frameworks should be conducted to assess whether these provide improvements over the time-invariant model approach to change estimation that was used in this research. 3. Information from an NFI is often considered in decision-making processes for national and/or regional forest policy. Support for NFIs in developing countries often comes from donor support, yet the value of this information for its intended purpose is rarely assessed. There is a critical need to develop the methodological tools for assessing the 117  value of NFI information and to investigate this need in the context of a cost-plus-loss framework. 4. Forest policies and management planning require information concerning why past changes in land cover and AGB have occurred.  In order to enhance decision-making processes, research is needed in three key areas. First, the specific abiotic and biotic factors related to land cover and AGB change need to be identified. Second, where biotic factors can be identified as anthropogenic, research is needed to improve understanding of the direct and indirect anthropogenic activities and the underlying socio-economic factors behind land cover change and AGB change. Third, with this information scenarios of possible changes in land cover and AGB should be examined in relation to climate change, policy regimes, and socio-economic status of rural communities. 5. Given the estimation framework in Chapter 5, spatially and temporally precise estimation of specific changes in AGB within forest remaining forest is possible. The direction of such change may indicate growth (positive change) or decline (negative change). Natural and anthropogenically-induced decline of AGB within ecosystems is common (e.g., insect outbreaks, fire, harvesting, etc.); however, anthropogenically-induced AGB decline that exceeds a threshold over long time periods could be defined as degradation. In order to answer the question of whether degradation is occurring, research is needed in two key areas. First, the reasons for AGB decline within forest remaining forest need to be identified. Second, where these reasons can be attributed to human activities, the issue of whether these activities result in long-term AGB reduction needs to be addressed.  118  Bibliography Achard, F., Beuchle, R., Mayaux, P., Stibig, H.J., Bodart, C., Brink, A., Carboni, S., Desclée, B., Donnay, F., Eva, H.D. & Lupi, A. (2014) Determination of tropical deforestation rates and related carbon losses from 1990 to 2010. Global Change Biology, 20(8), 2540-2554. Adjorlolo, C., & Mutanga, O. (2013). Integrating remote sensing and geostatistics to estimate woody vegetation in an African savanna. Journal of Spatial Science, 58(2), 305-322. Adler-Golden, S. M., Matthew, M. W., Bernstein, L. S., Levine, R. Y., Berk, A., Richtsmeier, S. C., Acharya, P.K., Anderson G.P., Felde, J.W., Gardner, J.A., Hoke, M.L., Jeong, L.S., Pukall, B., Ratkowski, A.J., & Hsiao-hua, K. B. (1999). Atmospheric correction for shortwave spectral imagery based on MODTRAN4. In SPIE's International Symposium on Optical Science, Engineering, and Instrumentation (pp. 61-69). International Society for Optics and Photonics. Ahrends, A., Burgess, N. D., Milledge, S. A., Bulling, M. T., Fisher, B., Smart, J. C., Clarke, G., Mhorok, B.E., and Lewis, S. L. (2010). Predictable waves of sequential forest degradation and biodiversity loss spreading from an African city. Proceedings of the National Academy of Sciences, 107(33), 14556-14561. Anaya, J. A., Chuvieco, E., & Palacios-Orueta, A. (2009). Above-ground biomass assessment in Colombia: a remote sensing approach. Forest Ecology and Management, 257(4). 119  Andersson K., Lawrence D., Zavaleta J., Guariguata M.R. (2016). More Trees, More Poverty? The Socioeconomic Effects of Tree Plantations in Chile, 2001–2011. Environmental Management, 57(1), 123-136.  Archibald S., Scholes R. J., Roy D.P., Roberts G., Boschetti L. (2010). Southern African fire regimes as revealed by remote sensing. International Journal of Wildland Fire, 19(7), 861-878. Armston, J. D., Denham, R. J., Danaher, T. J., Scarth, P. F., & Moffiet, T. N. (2009). Prediction and validation of foliage projective cover from Landsat-5 TM and Landsat-7 ETM+ imagery. Journal of Applied Remote Sensing, 3(1), 033540-033540. Arnold, F. E., Rametsteiner, E., & Kleinn, C. (2014). User-oriented national forest monitoring planning: a contribution to more policy relevant forest information provision. International Forestry Review, 16(4), 389-404. Baccini, A. G. S. J., Goetz, S. J., Walker, W. S., Laporte, N. T., Sun, M., Sulla-Menashe, D., Hackler, J., Beck, P.S.A., Dubayah, R., Friedl, M.A. & Samanta, S. (2012). Estimated carbon dioxide emissions from tropical deforestation improved by carbon-density maps. Nature Climate Change, 2(3), 182-185. Banskota, A., Kayastha, N., Falkowski, M. J., Wulder, M. A., Froese, R. E., & White, J. C. (2014). Forest monitoring using Landsat time series data: a review. Canadian Journal of Remote Sensing, 40(5), 362-384. Bater, C. W., Coops, N. C., Gergel, S. E., LeMay, V., & Collins, D. (2009). Estimation of standing dead tree class distributions in northwest coastal forests using lidar remote sensing. Canadian Journal of Forest Research, 39(6), 1080-1091. 120  Baxter, P.W. & Getz, W.M. (2005). A model-framed evaluation of elephant effects on tree and fire dynamics in African savannas. Ecological Applications, 15(4), 1331-1341. Belsky, A.J. (1995). Spatial and temporal landscape patterns in arid and semi-arid African savannas. In Hansson, L., Fahrig, L. & Merriam, G. (Eds.), Mosaic landscapes and ecological processes (pp. 31-56). Springer Netherlands. Bodart, C., Brink, A. B., Donnay, F., Lupi, A., Mayaux, P., & Achard, F. (2013). Continental estimates of forest cover and forest cover changes in the dry ecosystems of Africa between 1990 and 2000. Journal of Biogeography, 40(6), 1036-1047. Boschetti, L., Roy, D., & Hoffmann, A. A. (2013). MODIS Collection 5.1 Burned Area Product-MCD45. User’s Guide. Ver, 3.0.1. University of Maryland.  Brink, A. B., Bodart, C., Brodsky, L., Defourney, P., Ernst, C., Donney, F., Lupi, A. & Tuckova, K. (2014). Anthropogenic pressure in East Africa—Monitoring 20 years of land cover changes by means of medium resolution satellite data. International Journal of Applied Earth Observation and Geoinformation, 28, 60-69. Bucini, G., Hanan, N. P., Boone, R. B., Smit, I. P. J., Saatchi, S., Lefsky, M. A., & Asner, G. P. (2010). Woody fractional cover in Kruger National Park, South Africa: remote-sensing-based maps and ecological insights. In Hill, M. J.  & Hanan, N. P. (Eds.), Ecosystem function in savannas: measurement and modeling at landscape to global scales (pp. 219-238). CRC Press. Burkhart, H. E., & Tomé, M. (2012). Modeling forest trees and stands. Springer Science & Business Media.  Burrows, P. M., & Strang, R. M. (1964). The relation of crown and basal diameters in Rhodesian Brachystegia woodland. The Commonwealth Forestry Review, 331-333. 121  Cabral, A. I. R., Vasconcelos, M. J., Oom, D., & Sardinha, R. (2010). Spatial dynamics and quantification of deforestation in the central-plateau woodlands of Angola (1990–2009). Applied Geography, 31(3), 1185-1193. Campbell, B. M., Swift, M. J., Hatton, J., & Frost, P. G. H. (1988). Small-scale vegetation pattern and nutrient cycling in miombo woodland. In: Verhoeven, J.T.A., Heil, G.W., & Werger, M. (Eds.), Vegetation structure in relation to carbon and nutrient economy (pp. 69-85.), SPB Academic Publishing, The Hague.  Campbell, B. M., Swift, M. J., Frost, P. G. H., & Kirchmann, H. (1998). Comparative ecosystem characteristics of a miombo woodland and an adjacent agricultural field (Zimbabwe). In Bergstrom, L., & Kirchmann, H. (Eds.), Carbon and nutrient dynamics in natural and agricultural tropical ecosystems (pp. 201-226). Cab International. Canty, M. J. (2014). Image analysis, classification and change detection in remote sensing: with algorithms for ENVI/IDL and Python, Third Edition. CRC Press.  Carreiras, J., Melo, J. B., & Vasconcelos, M. J. (2013). Estimating the above-ground biomass in miombo savanna woodlands (Mozambique, East Africa) using L-band synthetic aperture radar data. Remote Sensing, 5(4), 1524-1548. Caylor, K. K., & Shugart, H. H. (2004). Simulated productivity of heterogeneous patches in Southern African savanna landscapes using a canopy productivity model. Landscape Ecology, 19(4), 401-415. Chave, J., Réjou-Méchain, M., Búrquez, A., Chidumayo, E., Colgan, M.S., Delitti, W.B., Duque, A., Eid, T., Fearnside, P.M., Goodman, R.C. & Henry, M. (2014). Improved 122  allometric models to estimate the above-ground biomass of tropical trees. Global Change Biology, 20(10), 3177-3190. Chidumayo, E. N. (1988). A re-assessment of effects of fire on miombo regeneration in the Zambian Copperbelt. Journal of Tropical Ecology, 4(4), 361-372. Chidumayo, E. N. (1990). Above-ground woody biomass structure and productivity in a Zambezian woodland. Forest Ecology and Management, 36(1), 33-46. Chidumayo, E. N. (2002). Changes in miombo woodland structure under different land tenure and use systems in central Zambia. Journal of Biogeography, 29(12), 1619-1626. Chidumayo, E.N. (2012). Assessment of existing models for biomass and volume calculations for Zambia. Report prepared for FAO-Zambia and the Integrated Land Use Assessment (ILUA) Phase 2 Project. Lusaka. Chidumayo, E. N. & Frost, P. (1996). Population biology of miombo trees. In Campbell, B. (Ed.), The Miombo in transition: woodlands and welfare in Africa (pp. 59-71). Center for International Forestry Research, Bogor, Indonesia.  Chidumayo, E., Gambiza, J., & Grundy, I. (1996). Managing miombo woodlands. In Campbell, B. (Ed.), The Miombo in transition: woodlands and welfare in Africa (pp. 175-193). Center for International Forestry Research, Bogor, Indonesia.  Chidumayo, E. N. & Marunda C. (2010). Dry forests and woodlands in sub-Saharan Africa: Context and challenges. In Chidumayo, E.N. & Gumbo, D.J. (Eds.), The dry forests and woodlands of Africa: Managing for products and services (pp. 1-10). Earthscan, London. 123  Ciais, P., Bombelli, A., Williams, M., Piao, S.L., Chave, J., Ryan, C.M., Henry, M., Brender, P. & Valentini, R. (2011). The carbon balance of Africa: synthesis of recent research studies. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 369(1943), 2038-2057. Claverie, M., Vermote, E.F., Franch, B. & Masek J.G. (2015). Evaluation of the Landsat-5 TM and Landsat-7 ETM+ surface reflectance products. Remote Sensing of Environment, 169, 390-403. Cochran, W. G. (1977). Sampling techniques, 3rd edition. New York: Wiley.  Coulston, J.W., Reams, G.A., Wear, D.N. & Brewer C.K. (2013). An analysis of forest land use, forest land cover and change at policy-relevant scales. Forestry, 0, 1–10, doi:10.1093/forestry/cpt056. Crookston, N. L., & Finley, A. O. (2008). yaImpute: An R package for kNN imputation. Journal of Statistical Software, 23(10), 1-16. CSO (Central Statistics Office). (2012). Zambia 2010 Census of Population and Housing: Volume 11, National Descriptive Tables. Republic of Zambia, Central Statistics Office. Lusaka, Zambia. Cutler, D. R., Edwards Jr, T. C., Beard, K. H., Cutler, A., Hess, K. T., Gibson, J., & Lawler, J. J. (2007). Random forests for classification in ecology. Ecology, 88(11), 2783-2792. Dalal-Clayton, D.B., English, C., Williams, G.J. & Spaargarden, O.C. (1985). A geomorphic legend for Zambia, Technical Guide No. 15. Soil Survey Unit, Research Branch, Department of Agriculture. Lusaka, Zambia. Desanker, P. V., & Prentice, I. C. (1994). MIOMBO—A vegetation dynamics model for the miombo woodlands on Zambezian Africa. Forest Ecology and Management, 69(1), 124  87-95.ESRI. (2012). ArcGIS Desktop: Release 10.1. Redlands, CA: Environmental Systems Research Institute. D'Odorico, P., Caylor, K., Okin, G. S., & Scanlon, T. M. (2007). On soil moisture–vegetation feedbacks and their possible effects on the dynamics of dryland ecosystems. Journal of Geophysical Research: Biogeosciences (2005–2012), 112(G4). Efron, B., & Tibshirani, R. (1986). Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy. Statistical Science, 1(1), 54-75.  Eitel, J. U., Vierling, L. A., Litvak, M. E., Long, D. S., Schulthess, U., Ager, A. A., Krofcheck, D. J. & Stoscheck, L. (2011). Broadband, red-edge information from satellites improves early stress detection in a New Mexico conifer woodland. Remote Sensing of Environment, 115(12), 3640-3646. Ek, T. (1994). Biomass structure in Miombo woodland and semi-evergreen forest: development in 22 permanent plots in Morogoro, Tanzania. A thesis submitted in partial fulfillment of the requirement for the degree of Cand. Agric., Department of Forestry, Agricultural University of Norway.  Elzinga, C. L., Salzer, D. W., Willoughby, J. W., & Gibbs, J. P. (2009). Monitoring plant and animal populations: a handbook for field biologists. John Wiley & Sons. Ene, L.T., Næsset, E., Gobakken, T., Mauya, E.W., Bollandsås, O.M., Gregoire, T.G., Ståhl, G. and Zahabu, E. (2016). Large-scale estimation of aboveground biomass in miombo woodlands using airborne laser scanning and national forest inventory data. Remote Sensing of Environment, 186, 626-636.   Eskelson, B. N., Temesgen, H., LeMay, V., Barrett, T. M., Crookston, N. L., & Hudak, A. T. (2009). The roles of nearest neighbor methods in imputing missing data in forest 125  inventory and monitoring databases. Scandinavian Journal of Forest Research, 24(3), 235-246. ESRI. (2012). ArcGIS Desktop: Release 10.1. Redlands, CA: Environmental Systems Research Institute. FAO. (2008). NFMA approach and process: an analysis of cost and time. National Forest Monitoring and Assessment Working Paper 39. Rome. FAO. (2010). Global forest resources assessment, 2010. Terms and definitions. Working paper 144/E. Rome. FAO. (2014a). Zambia - Global Forest Resources Assessment 2015 – Country Report. Rome. 76p. FAO. (2014b). State of the world’s forests: enhancing the socioeconomic benefits from forests. Food and Agriculture Organization of the United Nations, Rome.  FAO. (2015). Global Forest Resources Assessment 2015: How have the world’s forests changed? Food and Agriculture Organization of the United Nations. Rome, Italy. Faßnacht, F. E., Hartig, F., Latifi, H., Berger, C., Hernández, J., Corvalán, P., & Koch, B. (2014). Importance of sample size, data type and prediction method for remote sensing-based estimations of aboveground forest biomass. Remote Sensing of Environment, 154, 102-114. Federici, S., Tubiello, F.N., Salvatore, M., Jacobs, H. & Schmidhuber J. (2015). New estimates of CO2 forest emissions and removals: 1990–2015. Forest Ecology and Management, 352, 89-98. 126  Foody, G. M., Boyd, D. S., & Cutler, M. E. (2003). Predictive relations of tropical forest biomass from Landsat TM data and their transferability between regions. Remote Sensing of Environment, 85(4), 463-474. Frost, P. (1996). The ecology of miombo woodlands. In Campbell, B. (Ed.), The Miombo in transition: woodlands and welfare in Africa (pp. 11-57). Center for International Forestry Research, Bogor, Indonesia.  Fuller, D. O., Prince, S. D., & Astle, W. L. (1997). The influence of canopy strata on remotely sensed observations of savanna-woodlands. International Journal of Remote Sensing, 18(14), 2985-3009. Gara, T. W., Murwira, A., Ndaimani, H., Chivhenge, E., & Hatendi, C. M. (2015). Indigenous forest wood volume estimation in a dry savanna, Zimbabwe: exploring the performance of high-and-medium spatial resolution multispectral sensors. Transactions of the Royal Society of South Africa, 1-9. Garcia-Gutierrez, J., Gonzalez-Ferreiro, E., Riquelme-Santos, J. C., Miranda, D., Dieguez-Aranda, U., & Navarro-Cerrillo, R. M. (2014). Evolutionary feature selection to estimate forest stand variables using LiDAR. International Journal of Applied Earth Observation and Geoinformation, 26, 119-131. GDAL. (2016). GDAL - Geospatial Data Abstraction Library: Version 2.0.2, Open Source Geospatial Foundation, http://gdal.osgeo.org. Geist, H.J. & Lambin, E.F. (2002). Proximate causes and underlying driving forces of tropical deforestation. BioScience, 52(2), 143-150. GFOI. (2013). Integrating remote-sensing and ground-based observations for estimation of emissions and removals of greenhouse gases in forests: methods and guidance from 127  the Global Forest Observations Initiative. Group on Earth Observations, Geneva, Switzerland. 163p. Gill, S. J., Biging, G. S., & Murphy, E. C. (2000). Modeling conifer tree crown radius and estimating canopy cover. Forest Ecology and Management, 126(3), 405-416. GOFC-GOLD. (2015). A sourcebook of methods and procedures for monitoring and reporting anthropogenic greenhouse gas emissions and removals associated with deforestation, gains and losses of carbon stocks in forests remaining forests, and forestation. GOFC-GOLD Report version COP21-1. GOFC-GOLD Land Cover Project Office, Wageningen University, The Netherlands. Gómez, C., White, J.C. & Wulder, M.A. (2016). Optical remotely sensed time series data for land cover classification: A review. ISPRS Journal of Photogrammetry and Remote Sensing, 116, 55-72. González-Roglich, M., Swenson, J. J., Jobbágy, E. G., & Jackson, R. B. (2014). Shifting carbon pools along a plant cover gradient in woody encroached savannas of central Argentina. Forest Ecology and Management, 331, 71-78. González-Roglich, M., & Swenson, J. J. (2016). Tree cover and carbon mapping of Argentine savannas: scaling from field to region. Remote Sensing of Environment, 172, 139-147.  Gregoire, T. G. (1993). Estimation of forest growth from successive surveys. Forest Ecology and Management, 56(1), 267-278. Gregoire, T. G. (1998). Design-based and model-based inference in survey sampling: appreciating the difference. Canadian Journal of Forest Research, 28(10), 1429-1447. 128  Gregoire, T. G., & Valentine, H. T. (2007). Sampling strategies for natural resources and the environment. CRC Press.  Grogan, K., Birch-Thomsen, T. & Lyimo, J. (2013). Transition of shifting cultivation and its impact on people’s livelihoods in the Miombo woodlands of northern Zambia and south-western Tanzania. Human Ecology, 41(1), 77-92.  GRZ (Government of the Republic of Zambia). (1968). The Republic of Zambia – National rainfall map, 1:3,000,000. Government Printer, Lusaka. Accessed online at http://eusoils.jrc.ec.europa.eu/esdb_archive/eudasm/africa/lists/czm.htm on 11 April, 2014. GRZ (Government of the Republic of Zambia). (1976). Vegetation map, 1:500,000. Sheets 1 – 9. Survey Department. Lusaka, Zambia. GRZ (Government of the Republic of Zambia). (1986). The Republic of Zambia – National soil map, 1:3,000,000. Government Printer, Lusaka. Accessed online at http://eusoils.jrc.ec.europa.eu/esdb_archive/eudasm/africa/lists/czm.htm on 11 April, 2014. GRZ (Government of the Republic of Zambia). (2014). Integrated Land Use Assessment Phase 2 - Zambia Biophysical Field Manual, version March 25, 2014. Forestry Department, Ministry of Lands, Natural Resources and Environmental Protection in cooperation with Food and Agriculture Organization (FAO). Lusaka, Zambia.  GRZ (Government of the Republic of Zambia). (2016). Submission on proposed forest reference level (4 January 2016). Accessed online at http://redd.unfccc.int/submissions.html on 20 July 2016.  129  Guisan, A., & Zimmermann, N. E. (2000). Predictive habitat distribution models in ecology. Ecological Modelling, 135(2), 147-186. Guisan, A., Edwards, T. C., & Hastie, T. (2002). Generalized linear and generalized additive models in studies of species distributions: setting the scene. Ecological modelling, 157(2), 89-100.  Gumbo, D.J. & Chidumayo, E.N. (2010). Managing dry forests and woodlands for products and services: a prognostic synthesis. In Chidumayo, E.N. & Gumbo, D.J. (Eds.), The dry forests and woodlands of Africa: Managing for products and services (pp. 261-280). Earthscan, London. Gumbo, D. J., Moombe, K.B., Kandulu, M.M., Kabwe, G., Ojanen, M., Ndhlovu, E. & Sunderland T.C. (2013). Dynamics of the charcoal and indigenous timber trade in Zambia: A scoping study in Eastern, Northern and Northwestern provinces. Occasional Paper 86, Center for International Forestry Research. Bogor, Indonesia. Gumbo, D. J., Mumba, K. Y., Kaliwile, M. M., Moombe, K. B., & Mfuni, T. I. (2016). Agrarian changes in the Nyimba District of Zambia. In: Deakin, L., Kshatriya, M., & Sunderland, T. (Eds.), Agrarian change in tropical landscapes (pp. 234-268.) Center for International Forestry Research, Bogor. Hall, R. J., Skakun, R. S., Arsenault, E. J., & Case, B. S. (2006). Modeling forest stand structure attributes using Landsat ETM+ data: Application to mapping of aboveground biomass and stand volume. Forest Ecology and Management, 225(1), 378-390.  Halperin, J., LeMay, V., Coops, N., Verchot, L., Marshall, P., & Lochhead, K. (2016a). Canopy cover estimation in miombo woodlands of Zambia: Comparison of Landsat 8 130  OLI versus RapidEye imagery using parametric, nonparametric, and semiparametric methods. Remote Sensing of Environment, 179, 170-182. Halperin, J., LeMay, V., Chidumayo, E. N., Verchot, L. & Marshall, P. (2016b). Model-based estimation of above-ground biomass in the miombo ecoregion of Zambia. Forest Ecosystems, 3:14. Hansen, M. C., DeFries, R. S., Townshend, J. R. G., Marufu, L., & Sohlberg, R. (2002). Development of a MODIS tree cover validation data set for Western Province, Zambia. Remote Sensing of Environment, 83(1), 320-335. Hansen, M. C., Stehman, S. V., & Potapov, P. V. (2010). Quantification of global gross forest cover loss. Proceedings of the National Academy of Sciences, 107(19), 8650-8655. Hansen, M.C., Potapov, P.V., Moore, R., Hancher, M., Turubanova, S.A., Tyukavina, A., Thau, D., Stehman, S.V., Goetz, S.J., Loveland, T.R. & Kommareddy, A. (2013). High-resolution global maps of 21st-century forest cover change. Science, 342(6160), 850-853. Harrison, P. A., Berry, P. M., Simpson, G., Haslett, J. R., Blicharska, M., Bucur, M., Dunford, R., Egohe, B., Garcia-Llorente, M., Geamănă, N., Geertsemah, W., Lommeleni, E., Meiresonnei, L., & Turkelboom, F. (2014). Linkages between biodiversity attributes and ecosystem services: A systematic review. Ecosystem Services, 9, 191-203. Hawkins, D. M. (2004). The problem of overfitting. Journal of Chemical Information and Computer Sciences, 44(1), 1-12. 131  Heiskanen, J., & Kivinen, S. (2008). Assessment of multispectral,-temporal and-angular MODIS data for tree cover mapping in the tundra–taiga transition zone. Remote Sensing of Environment, 112(5), 2367-2380. Helmer, E. H., Brandeis, T. J., Lugo, A. E., & Kennaway, T. (2008). Factors influencing spatial pattern in tropical forest clearance and stand age: implications for carbon storage and species diversity. Journal of Geophysical Research, 113, G02S04, DOI:10.1029/2007JG000568. Hengl, T., Heuvelink, G. B., Kempen, B., Leenaars, J. G., Walsh, M. G., Shepherd, K. D., Sila, A., MacMillan, R.A., de Jesus, J.M., Tamene, L. & Tondoh, J. E. (2015). Mapping soil properties of Africa at 250 m resolution: Random forests significantly improve current predictions. PloSone, 10(6), e0125814. Hermosilla, T., Wulder, M.A., White, J.C., Coops, N.C. & Hobart, G.W. (2015). An integrated Landsat time series protocol for change detection and generation of annual gap-free surface reflectance composites. Remote Sensing of Environment, 158, 220-234. Hijmans, R. J., van Etten, J., Mattiuzzi, M., Sumner, M., Greenberg, J., Lamigueiro, O., Bevan, A., Racine, E., & Shortridge, A. (2015). Raster: geographic analysis and modeling with raster data. R package version 2.3-41. http://cran.r-project.org/web/packages/raster/index.html. Hojas-Gascon, L., Cerutti, P. O., Eva, H., Nasi, R., & Martius, C. (2015). Monitoring deforestation and forest degradation in the context of REDD+: lessons from Tanzania. Infobrief 124. Center for International Forestry Research (CIFOR), Bogor, Indonesia. 132  Holben, B.N. (1986). Characteristics of maximum-value composite images from temporal AVHRR data. International Journal of Remote Sensing, 7(11), 1417-1434. Holmgren, P., and L. G. Marklund. F. (2007). National forest monitoring systems: purposes, options and status. IN: Freer-Smith, P. H., Broadmeadow, M. S., & Lynch, J. M. (Eds.), Forestry and climate change (pp. 163-173). CAB International, UK. Hosonuma, N., Herold, M., De Sy, V., De Fries, R.S., Brockhaus, M., Verchot, L., Angelsen, A. & Romijn, E. (2012). An assessment of deforestation and forest degradation drivers in developing countries. Environmental Research Letters, 7(4), 044009. Hostert, P., Griffiths, P., van der Linden, S. & Pflugmacher, D. (2015). Time series analyses in a new era of optical satellite data. In Kuenzer, C., Dech, S. & Wagner, W. (Eds.), Remote Sensing Time Series: Revealing Land Surface Dynamics, Vol. 22 (pp. 25-41). Springer. Houghton, R.A. (2005). Aboveground forest biomass and the global carbon balance. Global Change Biology, 11(6), 945-958. Houghton, R.A., House, J.I., Pongratz, J., Van der Werf, G.R., DeFries, R.S., Hansen, M.C., Quéré, C.L. & Ramankutty, N. (2012). Carbon emissions from land use and land-cover change. Biogeosciences, 9(12), 5125-5142. Hunt, E. R., & Rock, B. N. (1989). Detection of changes in leaf water content using near-and middle-infrared reflectances. Remote Sensing of Environment, 30(1), 43-54. IPCC. (2006). IPCC Guidelines for National Greenhouse Gas Inventories, prepared by the National Greenhouse Gas Inventories Programme (eds Eggleston HS, et al.), Institute for Global Environmental Strategies, Hayama. 133  IPCC. (2014). Climate Change 2014: Mitigation of Climate Change. Contribution of Working Group III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (eds Edenhofer OR, Pichs-Madruga Y, Sokona E, et al.), Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. Irons, J. R., Dwyer, J. L., & Barsi, J. A. (2012). The next Landsat satellite: The Landsat data continuity mission. Remote Sensing of Environment, 122, 11-21. Isango. J.A. (2007). Stand structure and tree species composition of Tanzania miombo woodlands: a case study from miombo woodlands of community based forest management in Iringa District. Working Papers of the Finnish Forest Research Institute 50: 43–56. Jennings, S. B., Brown, N. D., & Sheil, D. (1999). Assessing forest canopies and understory illumination: canopy closure, canopy cover and other measures. Forestry, 72(1), 59-74. Jew, E.K., Dougill, A.J., Sallu, S.M., O’Connell, J. & Benton T.G. (2016). Miombo woodland under threat: Consequences for tree diversity and carbon storage. Forest Ecology and Management, 361, 144-153. Jin, S., & Sader, S. A. (2005). Comparison of time series tasseled cap wetness and the normalized difference moisture index in detecting forest disturbances. Remote Sensing of Environment, 94(3), 364-372. Jones, O.D., Lane, P.N., Sheridan, G.J. (2016). The stochastic runoff-runon process: Extending its analysis to a finite hillslope. Journal of Hydrology, http://dx.doi.org/10.1016/j.jhydrol.2016.06.056  Justice, C. O., Vermote, E., Townshend, J.R.G., DeFries, R., Roy, D.P., Hall, D.K., Salomonson, V.V., Privette, J.L., Riggs, G., Strahler, A., Lucht, W., Myneni, R.B., 134  Knyazikhin, Y., Running, S.W., Nemani, R.R., Zhengming W., Huete, A.R., Van Leeuwen, W., Wolfe, R.E., Giglio, L., Muller, J.-P., Lewis, P., Barnsley, M.J. (1998). The Moderate Resolution Imaging Spectroradiometer (MODIS): land remote sensing for global change research. Geoscience and Remote Sensing, IEEE Transactions on, 36(4), 1228-1249. Kalaba, F. K., Quinn, C. H., & Dougill, A. J. (2013). The role of forest provisioning ecosystem services in coping with household stresses and shocks in Miombo woodlands, Zambia. Ecosystem Services, 5, 143-148. Kangas, A. (2006). Design-based sampling and inference. In: Kangas, A., & Maltamo, M. (Eds.), Forest inventory: methodology and applications (pp. 13-38). Springer Science & Business Media.  Kangas, A., Gove, J., & Scott, C. (2006). Introduction. In: Kangas, A., & Maltamo, M. (Eds.), Forest inventory: methodology and applications (pp. 3-8). Springer Science & Business Media.  Karlson, M., Ostwald, M., Reese, H., Sanou, J., Tankoano, B., & Mattsson, E. (2015). Mapping tree canopy cover and aboveground biomass in Sudano-Sahelian woodlands using Landsat 8 and Random Forest. Remote Sensing, 7(8), 10017-10041. Kashindye, A., Mtalo, E., Mpanda, M. M., Liwa, E., & Giliba, R. (2013). Multi-temporal assessment of forest cover, stocking parameters and above-ground tree biomass dynamics in miombo woodlands of Tanzania. African Journal of Environmental Science and Technology, 7(7), 611-623. Kattenborn, T., Maack, J., Faßnacht, F., Enßle, F., Ermert, J., & Koch, B. (2015). Mapping forest biomass from space – fusion of hyperspectral EO1-hyperion data and Tandem-X 135  and WorldView-2 canopy height models. International Journal of Applied Earth Observation and Geoinformation, 35, 359-367. Kaufman, Y. J., & Tanre, D. (1992). Atmospherically resistant vegetation index (ARVI) for EOS-MODIS. Geoscience and Remote Sensing, IEEE Transactions on, 30(2), 261-270. Kennedy, R., Yang, Z. & Cohen W.B. (2010). Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr—Temporal segmentation algorithms. Remote Sensing of Environment, 114(12), 2897-2910. Kennedy, R.E., Andréfouët, S., Cohen, W.B., Gómez, C., Griffiths, P., Hais, M., Healey, S.P., Helmer, E.H., Hostert, P., Lyons, M.B. & Meigs, G.W. (2014). Bringing an ecological view of change to Landsat-based remote sensing. Frontiers in Ecology and the Environment, 12(6), 339-346. Key, C. H., & Benson, N. C. (2006). Landscape assessment (LA). FIREMON: fire effects monitoring and inventory system. Gen. Tech. Rep. RMRS-GTR-164-CD, Fort Collins, CO: US Department of Agriculture, Forest Service, Rocky Mountain Research Station. Kleinn, C., Corrales, L., & Morales, D. (2002). Forest area in Costa Rica: a comparative study of tropical forest cover estimates over time. Environmental Monitoring and Assessment, 73(1), 17-40. Kleinn, C., Ståhl, G., Fehrmann, L., & Kangas, A. (2010). Issues in forest inventories as an input to planning and decision processes. Allgemeine Forst- und Jagdzeitung, 181(11/12):205-210. Köhl, M., Magnussen, S. S., & Marchetti, M. (2006). Sampling methods, remote sensing and GIS multiresource forest inventory. Springer Science & Business Media.  136  Köhl, M., Lasco, R., Cifuentes, M., Jonsson, Ö., Korhonen, K.T., Mundhenk, P., de Jesus Navar, J. & Stinson, G. (2015). Changes in forest production, biomass and carbon: Results from the 2015 UN FAO Global Forest Resource Assessment. Forest Ecology and Management, 352, 21-34. Kovalskyy, V. & Roy, D. (2013). The global availability of Landsat 5 TM and Landsat 7 ETM+ land surface observations and implications for global 30m Landsat data product generation. Remote Sensing of Environment, 130, 280-293. Kurz, W.A., Dymond, C.C., White, T.M., Stinson, G., Shaw, C.H., Rampley, G.J., Smyth, C., Simpson, B.N., Neilson, E.T., Trofymow, J.A. & Metsaranta, J. (2009). CBM-CFS3: a model of carbon-dynamics in forestry and land-use change implementing IPCC standards. Ecological Modelling, 220(4), 480-504. Kutner, M. H., Nachtsheim, C. J., Neter, J., & Li, W. (2005). Applied linear statistical models, 5th ed. McGraw Hill.  Lambin, E.F. (1997). Modelling and monitoring land-cover change processes in tropical regions. Progress in Physical Geography, 21(3), 375-393. Lawrence, M., McRoberts, R. E., Tomppo, E., Gschwantner, T., & Gabler, K. (2010). Comparisons of national forest inventories. In: Tomppo, E., Gschwantner, T., Lawrence, M., & McRoberts, R. E. (Eds.), National forest inventories: Pathways for Common Reporting (pp. 19-32). Springer.  Lefsky, M. A., Cohen, W. B., Harding, D. J., Parker, G. G., Acker, S. A., & Gower, S. T. (2002). Lidar remote sensing of above-ground biomass in three biomes. Global ecology and biogeography, 11(5), 393-399. 137  Lefsky, M. A., & Cohen, W. B. (2003). Selection of remotely sensed data. In: Wulder, M., & Franklin, S. E. (Eds.), Remote sensing of forest environments: concepts and case studies (pp. 13-46). Springer Science & Business Media. Springer US.  LeMay, V., & Temesgen, H. (2005). Comparison of nearest neighbor methods for estimating basal area and stems per hectare using aerial auxiliary variables. Forest Science, 51(2), 109-119. Lewis, S., Lloyd, J., Sitch, S., Mitchard, E. & Laurance, W. (2009a). Changing ecology of tropical forests: evidence and drivers. Annual Review of Ecology, Evolution, and Systematics, 40, 529-549.  Lewis, S.L., Lopez-Gonzalez, G., Sonké, B., Affum-Baffoe, K., Baker, T.R., Ojo, L.O., Phillips, O.L., Reitsma, J.M., White, L., Comiskey, J.A. & Ewango, C.E. (2009b). Increasing carbon storage in intact African tropical forests. Nature, 457(7232), 1003-1006. Liaw, A., & Wiener, M. (2002). Classification and regression by randomForest. R News 2(3), 18-22. Liu, Y.Y., Van Dijk, A.I., De Jeu, R.A., Canadell, J.G., McCabe, M.F., Evans, J.P. & Wang, G. (2015). Recent reversal in loss of global terrestrial biomass. Nature Climate Change, 5(5), 470-474. Lohr, S., & Rao, J. K. (2006). Estimation in multiple-frame surveys. Journal of the American Statistical Association, 101(475), 1019-1030. Lu, D., Chen, Q., Wang, G., Liu, L., Li, G., & Moran, E. (2014). A survey of remote sensing-based aboveground biomass estimation methods in forest ecosystems. International Journal of Digital Earth, 1-43. DOI:10.1080/17538947.2014.990526. 138  Lu, Y., Coops, N. & Hermosilla, T. (2016). Regional assessment of pan-Pacific urban environments over 25 years using annual gap free Landsat data. International Journal of Applied Earth Observation and Geoinformation, 50, 198-210. Lumley, T. (2009). R package ‘leaps’, version 2.9, July 2, 2014.  http://cran.r-project.org/web/packages/leaps/index.html Magnussen, S., McRoberts, R. E., & Tomppo, E. O. (2009). Model-based mean square error estimators for k-nearest neighbour predictions and applications using remotely sensed data for forest inventories. Remote Sensing of Environment, 113(3), 476-488. Magnussen, S., Mandallaz, D., Breidenbach, J., Lanz, A., & Ginzler, C. (2014). National forest inventories in the service of small area estimation of stem volume. Canadian Journal of Forest Research, 44(9), 1079-1090. Magnussen, S. (2015a). Arguments for a model-dependent inference? Forestry 88(3): 317–325, http://dx.doi.org/10.1093/forestry/cpv002 Magnussen, S., Næsset & E., Gobakken, T. (2015b). LiDAR-supported estimation of change in forest biomass with time-invariant regression models. Canadian Journal of Forest Research, 45(11), 1514-1523. Magnussen, S., Andersen, H.E. & Mundhenk, P. (2015c). A second look at endogenous poststratification. Forest Science, 61(4), 624-634. Makhado, R. A., Mapaure, I., Potgieter, M. J., Luus-Powell, W. J., & Saidi, A. T. (2014). Factors influencing the adaptation and distribution of Colophospermum mopane in southern Africa’s mopane savannas – a review. Bothalia, 44(1), 9p. Mandallaz, D. (2007). Sampling techniques for forest inventories. CRC Press. 139  Mapaure, I. N., & Campbell, B. M. (2002). Changes in miombo woodland cover in and around Sengwa Wildlife Research Area, Zimbabwe, in relation to elephants and fire. African Journal of Ecology, 40(3), 212-219. Masek, J. G., Vermote, E. F., Saleous, N. E., Wolfe, R., Hall, F. G., Huemmrich, K. F., Gao, F., Kutler, J. and Lim, T. K. (2006). A Landsat surface reflectance dataset for North America, 1990-2000. Geoscience and Remote Sensing Letters, IEEE, 3(1), 68-72. Mauya, E. W., Ene, L. T., Bollandsås, O. M., Gobakken, T., Næsset, E., Malimbwi, R. E., & Zahabu, E. (2015). Modelling aboveground forest biomass using airborne laser scanner data in the miombo woodlands of Tanzania. Carbon Balance and Management, 10(1), 1-16. Mayes, M. T., Mustard, J. F., & Melillo, J. M. (2015). Forest cover change in miombo woodlands: modeling land cover of African dry tropical forests with linear spectral mixture analysis. Remote Sensing of Environment, 165, 203-215. McCall, M. K., Chutz, N., & Skutsch, M. (2016). Moving from Measuring, Reporting, Verification (MRV) of forest carbon to community Mapping, Measuring, Monitoring (MMM): Perspectives from Mexico. PloS one, 11(6), e0146038. McCullagh, P., & Nelder, J. A. (1989). Generalized Linear Models. Chapman and Hall, London. McNicol, I.M., Ryan, C.M. & Williams, M. (2015). How resilient are African woodlands to disturbance from shifting cultivation? Ecological Applications, 25(8), 2320-2336. McRoberts, R. E. (2011). Satellite image-based maps: scientific inference or pretty pictures? Remote Sensing of Environment, 115(2), 715-724. 140  McRoberts, R. E., Tomppo, E. O., & Næsset, E. (2010). Advances and emerging issues in national forest inventories. Scandinavian Journal of Forest Research, 25(4), 368-381. McRoberts, R. E., Magnussen, S., Tomppo, E. O., & Chirici, G. (2011). Parametric, bootstrap, and jackknife variance estimators for the k-Nearest Neighbors technique with illustrations using forest inventory and satellite image data. Remote Sensing of Environment, 115(12), 3165-3174. McRoberts, R. E., Andersen, H. E., & Næsset, E. (2014a). Using airborne laser scanning data to support forest sample surveys. In Maltamo, M. Næsset, E. & Vauhkonen, J. (Eds.), Forestry applications of airborne laser scanning (pp. 269-292). Springer Netherlands. McRoberts, R.E., Bollandsås, O.M., & Næsset, E. (2014b). Modeling and estimating change. In Maltamo, M. Næsset, E. & Vauhkonen, J. (Eds.), Forestry applications of airborne laser scanning (pp. 293-313). Springer Netherlands. McRoberts, R. E., Næsset, E., & Gobakken, T. (2015). The effects of temporal differences between map and ground data on map-assisted estimates of forest area and biomass. Annals of Forest Science, 1-9.  Mertens, B. & Lambin, E. (1997). Spatial modelling of deforestation in southern Cameroon: spatial disaggregation of diverse deforestation processes. Applied Geography, 17(2), 143-162. Misana, S., Mung'ong'o, C. & EMukamuri, B. (1996). Miombo woodlands in the wider context: macro-economic and inter-sectoral influences. In: Campbell, B. (Ed.). The Miombo in transition: woodlands and welfare in Africa (pp. 73-99). Center for International Forestry Research, Bogor, Indonesia.  141  Mitchard, E.T. & Flintrop, C.M. (2013). Woody encroachment and forest degradation in sub-Saharan Africa's woodlands and savannas 1982–2006. Phil. Trans. R. Soc. B, 368(1625), 20120406. Mitchard, E. T., Meir, P., Ryan, C. M., Woollen, E. S., Williams, M., Goodman, L. E., Mucavele, J., Watts, P., Woodhouse, I., & Saatchi, S. S. (2013). A novel application of satellite radar data: measuring carbon sequestration and detecting degradation in a community forestry project in Mozambique. Plant Ecology & Diversity, 6(1), 159-170. Moeur, M., & Stage, A. R. (1995). Most similar neighbor: an improved sampling inference procedure for natural resource planning. Forest Science, 41(2), 337-359. Moisen, G. G., & Edwards Jr, T. C. (1999). Use of Generalized Linear Models and digital data in a forest inventory of northern Utah. Journal of Agricultural, Biological, and Environmental Statistics, 372-390. Moisen, G. G., & Frescino, T. S. (2002). Comparing five modelling techniques for predicting forest characteristics. Ecological Modelling, 157(2), 209-225. Moisen, G. G., Freeman, E. A., Blackard, J. A., Frescino, T. S., Zimmermann, N. E., & Edwards Jr, T. C. (2006). Predicting tree species presence and basal area in Utah: a comparison of stochastic gradient boosting, generalized additive models, and tree-based methods. Ecological Modelling, 199(2), 176-187. Moore, ID., P.E. Gessler, G.A. Nielsen, & G.A. Petersen. (1993). Terrain attributes: estimation methods and scale effects. In Jakeman, A.J., Beck, M.B. & McAleer, M.  (Eds.), Modeling Change in Environmental Systems (pp. 189 – 214). Wiley, London.  142  Morgan, J. A., & Tatar, J. F. (1972). Calculation of the residual sum of squares for all possible regressions. Technometrics, 14(2), 317-325. Morrison, T.A., Holdo, R.M. & Anderson, T.M. (2016). Elephant damage, not fire or rainfall, explains mortality of overstorey trees in Serengeti. Journal of Ecology, 104(2), 409-418. Moussavi, M. S., Abdalati, W., Scambos, T., & Neuenschwander, A. (2014). Applicability of an automatic surface detection approach to micro-pulse photon-counting LiDAR altimetry data: implications for canopy height retrieval from future ICESat-2 data. International Journal of Remote Sensing, 35(13), 5263-5279. Mukosha, J. & Siampale, A. (2009). Integrated Land Use Assessment Zambia, 2005–2008. Ministry of Tourism, Environment and Natural Resources and Forestry Department, Government of the Republic of Zambia and the Food and Agriculture Organization of the United Nations. Lusaka, Zambia. Murwira, A. & Skidmore, A.K. (2006). Monitoring change in the spatial heterogeneity of vegetation cover in an African savanna. International Journal of Remote Sensing, 27(11), 2255-2269.  Næsset, E., Ørka, H.O., Solberg, S., Bollandsås, O.M., Hansen, E.H., Mauya, E., Zahabu, E., Malimbwi, R., Chamuya, N., Olsson, H. & Gobakken, T. (2016). Mapping and estimating forest area and aboveground biomass in miombo woodlands in Tanzania using data from airborne laser scanning, TanDEM-X, RapidEye, and global forest maps: A comparison of estimated precision. Remote Sensing of Environment, 175, 282-300. 143  Naidoo, L., Mathieu, R., Main, R., Kleynhans, W., Wessels, K., Asner, G., & Leblon, B. (2015). Savannah woody structure modeling and mapping using multi-frequency (X-, C-and L-band) Synthetic Aperture Radar data. ISPRS Journal of Photogrammetry and Remote Sensing, 105, 234-250. O'Brien, E.M. (1993). Climatic gradients in woody plant species richness: towards an explanation based on an analysis of southern Africa's woody flora. Journal of Biogeography, 181-198. Opsomer, J. D., Breidt, F. J., Moisen, G. G., & Kauermann, G. (2007). Model-assisted estimation of forest resources with generalized additive models. Journal of the American Statistical Association, 102(478), 400-409. Ozdemir, I. (2014). Linear transformation to minimize the effects of variability in understory to estimate percent tree canopy cover using RapidEye data. GIScience & Remote Sensing, 51(3), 288-300. Packalén, P., Temesgen, H., & Maltamo, M. (2012). Variable selection strategies for nearest neighbor imputation methods used in remote sensing based forest inventory. Canadian Journal of Remote Sensing, 38(5), 557-569. Pajares, G. (2015). Overview and current status of remote sensing applications based on Unmanned Aerial Vehicles (UAVs). Photogrammetric Engineering & Remote Sensing, 81(4), 281-330. Pan, Y., Birdsey, R.A., Fang, J., Houghton, R., Kauppi, P.E., Kurz, W.A., Phillips, O.L., Shvidenko, A., Lewis, S.L., Canadell, J.G. & Ciais, P. (2011). A large and persistent carbon sink in the world’s forests. Science, 333(6045), 988-993. 144  Payn, T., Carnus, J.M., Freer-Smith, P., Kimberley, M., Kollert, W., Liu, S., Orazio, C., Rodriguez, L., Silva, L.N. & Wingfield, M.J. (2015). Changes in planted forests and future global implications. Forest Ecology and Management, 352, 57-67.  Pereira, C. R. (2006). Estimating and mapping forest inventory variables using the k-NN method: Mocuba District case study – Mozambique. PhD dissertation, Universita Tuscia.  Petit, C., Scudder, T. & Lambin, E. (2001). Quantifying processes of land-cover change by remote sensing: resettlement and rapid land-cover changes in south-eastern Zambia. International Journal of Remote Sensing, 22(17), 3435-3456. Petitjean, F., Inglada, J. & Gançarski, P. (2012). Satellite image time series analysis under time warping. IEEE transactions on geoscience and remote sensing, 50(8), 3081-3095. Pflugmacher, D., Cohen, W. B., Kennedy, R. E., & Yang, Z. (2014). Using Landsat-derived disturbance and recovery history and lidar to map forest biomass dynamics. Remote Sensing of Environment, 151, 124-137.  Pinheiro, J. & Bates, D. (2000). Mixed-effects models in S and S-PLUS. Springer, NY.  Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D. & R Core Team (2016). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-128. Pintea, L., Bauer, M. E., Bolstad, P. V., & Pusey, A. (2003). Matching multiscale remote sensing data to inter-disciplinary conservation needs: the case of chimpanzees in western Tanzania. In S. Morain & A. Budge (Eds.), ISPRS Comission I Mid-Term Symposium in conjunction with Pecora 15/Land Satellite Information IV Conference, Integrated Remote Sensing at the Global, Regional and Local Scale, November 10-15, 2002, Denver, CO USA. 145  Pope, P.T. & Webster, J.T. (1972). Use of an F-statistic in stepwise regression procedures. Technometrics, 14, 327–340. Powell, S. L., Cohen, W. B., Healey, S. P., Kennedy, R. E., Moisen, G. G., Pierce, K. B., & Ohmann, J. L. (2010). Quantification of live aboveground forest biomass dynamics with Landsat time-series and field inventory data: A comparison of empirical modeling approaches. Remote Sensing of Environment, 114(5), 1053-1068. Prins, E. & Kikula, I.S. (1996). Deforestation and regrowth phenology in miombo woodland-assessed by Landsat Multispectral Scanner System data. Forest Ecology and Management, 84(1), 263-266. Putz, F.E., Zuidema, P.A., Pinard, M.A., Boot, R.G., Sayer, J.A., Sheil, D., Sist, P. & Vanclay, J.K. (2008). Improved tropical forest management for carbon retention. PLoS Biol, 6(7), e166. R Development Core Team. (2015). R: A Language and Environment for Statistical Computing, version 3.2.2 (August 2015). R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org Rana, P., Tokola, T., Korhonen, L., Xu, Q., Kumpula, T., Vihervaara, P., & Mononen, L. (2013). Training area concept in a two-phase biomass inventory using Airborne Laser Scanning and RapidEye satellite data. Remote Sensing, 6(1), 285-309. RapidEye. (2012). Satellite Imagery Product Specifications, Version 4.1, September 2012. https://earth.esa.int/documents/10174/896711/RE_Product_Specifications_ENG.pdf Ratanamahatana, C., Lin, J., Gunopulos, D., Keogh, E., Vlachos, M. & Das G. (2010). Mining time series data. In Maimon, O. & Rokach, L. (Eds.), Data Mining and Knowledge Discovery Handbook, 2nd edition (pp. 1049-1077). Springer US.  146  Ribeiro, N. S., Saatchi, S. S., Shugart, H. H., & Washington-Allen, R. A. (2008). Aboveground biomass and leaf area index (LAI) mapping for Niassa Reserve, northern Mozambique. Journal of Geophysical Research: Biogeosciences (2005–2012), 113(G3). Ribeiro, N., Cumbana, M.,  Mamugy, F., & Chaúque, A. (2012). Remote sensing of biomass in the miombo woodlands of southern Africa: opportunities and limitations for research. In Fatoyinbo (Ed.), Remote Sensing of Biomass - Principles and Applications (pp. 77-98). ISBN: 978-953-51-0313-4, InTech.  Ribeiro, N., Syampungani, S., Matakala, N., Nangoma, D. & Ribeiro-Barros, A.I. (2015). Chapter 19: Miombo woodlands research towards the sustainable use of Ecosystem services in Southern Africa. In Lo, Y-H, Blanco, J. & Roy, S. (Eds.), Biodiversity in Ecosystems-Linking Structure and Function (pp. 475-491). InTech.  Rouse, J. W., Jr., Haas, R. H., Schell, J. A. & Deering, D. W. (1974). Monitoring vegetation systems in the Great Plains with ERTS. In Freden, S.C., Mercanti, E.P. & Becker, M.A. (Eds.),  Proceedings of the third earth resources tech. satellite-1 symposium 10–14 Dec. 1973, Vol. 1: Technical Presentations Sec. A (pp. 309–317). NASA Science and Technology Information Office, Washington, DC. Ryan, C. M., & Williams, M. (2011). How does fire intensity and frequency affect miombo woodland tree populations and biomass?  Ecological Applications, 21(1), 48-60. Ryan, C. M., Williams, M., & Grace, J. (2011). Above-and belowground carbon stocks in a miombo woodland landscape of Mozambique. Biotropica, 43(4), 423-432. Ryan, C. M., Hill, T., Woollen, E., Ghee, C., Mitchard, E., Cassells, G., Grace, J., Woodhouse, I.H., & Williams, M. (2012). Quantifying small-scale deforestation and 147  forest degradation in African woodlands using radar imagery. Global Change Biology, 18(1), 243-257. Ryan, C.M., Williams, M., Hill, T.C., Grace, J. & Woodhouse, I.H. (2014a). Assessing the phenology of southern tropical Africa: a comparison of hemispherical photography, scatterometry, and optical/NIR remote sensing. IEEE Transactions on Geoscience and Remote Sensing, 52(1), 519-528, 10.1109/TGRS.2013.2242081. Ryan, C.M., Berry, N.J., Joshi, N. (2014b). Quantifying the causes of deforestation and degradation and creating transparent REDD+ baselines: a method and case study from central Mozambique. Applied Geography, 53, 45-54. Ryan, C.M., Pritchard, R., McNicol, I., Owen, M., Fisher, J.A. & Lehmann, C. (2016). Ecosystem services from southern African woodlands and their future under global change. Phil. Trans. R. Soc. B, 371(1703), 20150312. Salerno, J. (2016). Migrant decision-making in a frontier landscape. Environmental Research Letters, 11(4), 044019. Samimi, C., & Kraus, T. (2004). Biomass estimation using Landsat-TM and-ETM+. Towards a regional model for Southern Africa? GeoJournal, 59(3), 177-187. Sankaran, M., Hanan, N.P., Scholes, R.J., Ratnam, J., Augustine, D.J., Cade, B.S., Gignoux, J., Higgins, S.I., Le Roux, X., Ludwig, F. & Ardo, J. (2005). Determinants of woody cover in African savannas. Nature, 438(7069), 846-849. Sankaran, M., Ratnam, J., & Hanan, N. (2008). Woody cover in African savannas: the role of resources, fire and herbivory. Global Ecology and Biogeography, 17(2), 236-245. Särndal, C. E., Swensson, B., & Wretman, J. (1992). Model-assisted survey sampling Springer. New York. 148  Särndal, C. E. (2010). Models in survey sampling. In Carlson, M., Nyquist, H., & Villani, M. (Eds), Official Statistics – Methodology and Applications in Honour of Daniel Thorburn (pp. 15-28). Department of Statistics, Stockholm University. Sawe, T. C., Munishi, P. K., & Maliondo, S. M. (2014). Woodlands degradation in the Southern Highlands miombo of Tanzania: implications on conservation and carbon stocks. International Journal of Biodiversity and Conservation, 6(3), 230-237. Schabenberger, O., & Pierce, F. J. (2002). Contemporary statistical models for the plant and soil sciences. CRC Press.  Scholes, R.J. & Archer, S.R. (1997). Tree-grass interactions in savannas. Annual Review of Ecology and Systematics, 517-544.  Scrucca, L. (2013). GA: a package for Genetic Algorithms in R. Journal of Statistical Software, 53(4), 1-37. Sedano, F., Gómez, D., Gong, P., & Biging, G. S. (2008). Tree density estimation in a tropical woodland ecosystem with multiangular MISR and MODIS data. Remote Sensing of Environment, 112(5), 2523-2537. Shearman, P., Bryan, J., & Laurance, W. F. (2012). Are we approaching ‘peak timber’in the tropics? Biological Conservation, 151(1), 17-21. Shekede, M., Murwira, A., Masocha, M. & Zengeya, F. (2016). Decadal changes in mean annual rainfall drive long-term changes in bush-encroached southern African savannas. Austral Ecology. Sinha, S., Jeganathan, C., Sharma, L. K., & Nathawat, M. S. (2015). A review of radar remote sensing for biomass estimation. International Journal of Environmental Science and Technology, 12(5), 1779-1792. 149  Skutsch, M. M., & Ba, L. (2010). Crediting carbon in dry forests: The potential for community forest management in West Africa. Forest Policy and Economics, 12(4), 264-270. Smith, T. M. F. (1994). Sample surveys 1975-1990; an age of reconciliation? International Statistical Review, 62(1), 5-34. Solberg, S., Gizachew, B., Næsset, E., Gobakken, T., Bollandsås, O. M., Mauya, E. W., Olsson, H., Malimbwi, R. & Zahabu, E. (2015). Monitoring forest carbon in a Tanzanian woodland using interferometric SAR: a novel methodology for REDD+. Carbon Balance and Management, 10(1), 1-14. Sprague, D.S. & Oyama, S. (1999). Density and distribution of citemene fields in a Miombo woodland environment in Zambia. Environmental Management, 24(2), 273-280.  Ståhl, G., Saarela, S., Schnell, S., Holm, S., Breidenbach, J., Healey, S. P., Patterson, P. L., Magnussen, S., Næsset, E., McRoberts, R. E. & Gregoire, T. G. (2016). Use of models in large-area forest surveys: comparing model-assisted, model-based and hybrid estimation. Forest Ecosystems, 3(1), 1. Staver, A. C., Archibald, S., & Levin, S. (2011). Tree cover in sub-Saharan Africa: rainfall and fire constrain forest and savanna as alternative stable states. Ecology, 92(5), 1063-1072. Stringer, L. C., Dougill, A. J., Thomas, A. D., Spracklen, D. V., Chesterman, S., Speranza, C. I.,  Rueff, H., Riddell, M.,  Williams, M., Beedy, T.,  Abson, D.J., Klintenberg, P., Syampungani, S., Powell, P.,  Palmer, A. R., Seely, M. K., Mkwambisi, D. D., Falcao, M., Sitoe, A., Ross, S., & Kopolo, G. (2012). Challenges and opportunities in linking 150  carbon sequestration, livelihoods and ecosystem service provision in drylands. Environmental Science & Policy, 19, 121-135. Su, Y., Guo, Q., Xue, B., Hu, T., Alvarez, O., Tao, S., & Fang, J. (2016). Spatial distribution of forest aboveground biomass in China: Estimation through combination of spaceborne lidar, optical imagery, and forest inventory data. Remote Sensing of Environment, 173, 187-199. Suberu, M. Y., Bashir, N., & Mustafa, M. W. (2014). Over use of wood-based bioenergy in selected sub-Saharan Africa countries: review of unconstructive challenges and suggestions. Journal of Cleaner Production. doi:10.1016/j.jclepro.2014.04.014 Suganuma, H., Abe, Y., Taniguchi, M., Tanouchi, H., Utsugi, H., Kojima, T., & Yamada, K. (2006). Stand biomass estimation method by canopy coverage for application to remote sensing in an arid area of Western Australia. Forest Ecology and Management, 222(1), 75-87. Swetnam, R.D., Fisher, B., Mbilinyi, B.P., Munishi, P.K.T., Willcock, S., Ricketts, T., Mwakalila, S., Balmford, A., Burgess, N.D., Marshall, A.R. & Lewis, S.L. (2011). Mapping socio-economic scenarios of land cover change: A GIS method to enable ecosystem service modelling. Journal of Environmental Management, 92(3), 563-574. Temesgen, H., & Ver Hoef, J. M. (2015). Evaluation of the spatial linear model, Random Forest and Gradient Nearest-Neighbour methods for imputing potential productivity and biomass of the Pacific Northwest forests. Forestry, 88(1), 131-142. Timberlake, J. R. (1995). Colophospermum mopane: annotated bibliography and review. The Zimbabwe Bulletin of Forestry Research, No. 11. ISBN: 0-7974-1420-7.  151  Timberlake, J., Chidumayo, E.N., & Sawadogo, L.  (2010). Distribution and characteristics of African dry forests. In Chidumayo, E. N.  & Gumbo, D.J. (Eds.), The dry forests and woodlands of Africa: managing for products and services (pp. 11-42). Earthscan.  Timberlake, J., & Chidumayo, E.N. (2011). Miombo ecoregion vision report. Occasional Publications in Biodiversity No. 20 Biodiversity Foundation for Africa, Bulawayo, Zimbabwe. Timilsina, N., & Staudhammer, C. L. (2012). Individual tree mortality model for slash pine in Florida: a mixed modeling approach. Southern Journal of Applied Forestry, 36(4), 211-219. Tinley, K. L. (1982). The influence of soil moisture balance on ecosystem patterns in southern Africa. In Walker, B. H. (Ed.), Ecology of tropical savannas (pp. 175-192). Springer Berlin Heidelberg. Tipton, J., Opsomer, J. & Moisen, G. (2013). Properties of endogenous post-stratified estimation using remote sensing data. Remote Sensing of Environment, 139, 130-137. Tiwari, A. K., & Singh, J. S. (1984). Mapping forest biomass in India through aerial photographs and nondestructive field sampling. Applied Geography, 4(2), 151-165. Tomppo, E., & Halme, M. (2004). Using coarse scale forest variables as ancillary information and weighting of variables in k-NN estimation: a genetic algorithm approach. Remote Sensing of Environment, 92(1), 1-20. Tomppo, E. (2006). The Finnish multi-source National Forest Inventory – small area estimation and map production. In: Kangas, A., & Maltamo, M. (Eds.), Forest inventory: methodology and applications (pp. 195-224). Springer Science & Business Media. 152  Tomppo, E., Haakana, M., Katila, M., & Peräsaari, J. (2008). Multi-source national forest inventory: methods and applications. Springer Science & Business Media.  Tomppo, E. & Tuomainen, T. (2010). Development of Finland’s National Forest Inventory. In: Tomppo, E., Gschwantner, T., Lawrence, M., & McRoberts, R. E. (Eds.), National forest inventories: Pathways for Common Reporting (pp. 185-206). Springer.  Tomppo, E., Schadauer, K., McRoberts, R., Gschwantner, T., Gabler, K., & Ståhl, G. (2010). Introduction. In: Tomppo, E., Gschwantner, T., Lawrence, M., & McRoberts, R. E. (Eds.), National forest inventories: Pathways for Common Reporting (pp. 1-18). Springer.  Tomppo, E., Malimbwi, R., Katila, M., Mäkisara, K., Henttonen, H.M., Chamuya, N., Zahabu, E. & Otieno, J. (2014). A sampling design for a large area forest inventory: case Tanzania. Canadian Journal of Forest Research, 44(8), 931-948. Townsend, A.K., Cooch, E.G., Sillett, T.S., Rodenhouse, N.L., Holmes, R.T. & Webster M.S. (2016). The interacting effects of food, spring temperature, and global climate cycles on population dynamics of a migratory songbird. Global Change Biology, 22(2), 544-555. Trouet, V., Esper, J. & Beeckman, H. (2010). Climate/growth relationships of Brachystegia spiciformis from the miombo woodland in south central Africa. Dendrochronologia, 28(3), 161-171. Turpie, J., Warr, B., & Carter, J. (2015). Benefits of forest ecosystems in Zambia and the role of REDD+ in a green economy transformation. UNEP, Nairobi.  Urbazaev, M., Thiel, C., Mathieu, R., Naidoo, L., Levick, S. R., Smit, I. P., Asner, G.P. & Schmullius, C. (2015). Assessment of the mapping of fractional woody cover in 153  southern African savannas using multi-temporal and polarimetric ALOS PALSAR L-band images. Remote Sensing of Environment, 166, 138-153. USGS (United States Geological Survey). (2004). Shuttle Radar Topography Mission, 1 Arc Second SRTM s14e29, s14e30, s14e31, s15e29, s15e30, s15e31, s16e29, s16e30, s16e31, Unfilled Unfinished 2.0, Global Land Cover Facility, University of Maryland, College Park, Maryland, February 2000. Valentini, R., Arneth, A., Bombelli, A., Castaldi, S., Cazzolla Gatti, R., Chevallier, F., Ciais, P., Grieco, E., Hartmann, J., Henry, M. & Houghton, R.A. (2014). A full greenhouse gases budget of Africa: synthesis, uncertainties, and vulnerabilities. Biogeosciences, 11, 381-407. Venables, W. N., & Dichmont, C. M. (2004). GLMs, GAMs and GLMMs: an overview of theory for applications in fisheries research. Fisheries research, 70(2), 319-337. Wang, Y., LeMay, V. M., & Baker, T. G. (2007). Modelling and prediction of dominant height and site index of Eucalyptus globulus plantations using a nonlinear mixed-effects model approach. Canadian Journal of Forest Research, 37(8), 1390-1403. Watson, F.G., Becker, M.S., Milanzi, J. & Nyirenda, M. (2015). Human encroachment into protected area networks in Zambia: implications for large carnivore conservation. Regional Environmental Change, 15(2), 415-429. Westfall, J. A., & Morin, R. S. (2013). A cover-based method to assess forest characteristics using inventory data and GIS. Forest Ecology and Management, 298, 93-100. White, J.C., Wulder, M.A., Hobart, G.W., Luther, J.E., Hermosilla, T., Griffiths, P., Coops, N.C., Hall, R.J., Hostert, P., Dyk, A. & Guindon, L. (2014). Pixel-based image 154  compositing for large-area dense time series applications and science. Canadian Journal of Remote Sensing, 40(3), 192-212. Willcock, S., Phillips, O.L., Platts, P.J., Swetnam, R.D., Balmford, A., Burgess, N.D., Ahrends, A., Bayliss, J., Doggart, N., Doody, K. & Fanning, E. (2016). Land cover change and carbon emissions over 100 years in an African biodiversity hotspot. Global Change Biology, 22, 2787–2800. Williams, M., Ryan, C., Rees, R., Sambane, E., Fernando, J. & Grace J. (2008). Carbon sequestration and biodiversity of re-growing miombo woodlands in Mozambique. Forest Ecology and Management, 254(2), 145-155. Wood, S.N. (2006). Generalized Additive Models: an introduction with R. Chapman & Hall/CRC Press. Chicago.  Wood, S.N. (2015). Reference manual for mgcv: Mixed GAM Computation Vehicle with GCV/AIC/REML smoothness estimation, v 1.8-9 (30 Oct 2015). Available at: http://cran.r-project.org/web/packages/mgcv/index.html Wu, W., De Pauw, E., & Helldén, U. (2013). Assessing woody biomass in African tropical savannahs by multiscale remote sensing. International Journal of Remote Sensing, 34(13), 4525-4549. Wulder, M. A., White, J. C., Nelson, R. F., Næsset, E., Ørka, H. O., Coops, N. C., Hilker, T., Bater, C. W. & Gobakken, T. (2012a). Lidar sampling for large-area forest characterization: a review. Remote Sensing of Environment, 121, 196-209. Wulder, M.A., Masek, J.G., Cohen, W.B., Loveland, T.R. & Woodcock C.E. (2012b). Opening the archive: How free data has enabled the science and monitoring promise of Landsat. Remote Sensing of Environment, 122, 2-10. 155  Zhu, Z. & Woodcock, C.E. (2012). Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sensing of Environment, 118, 83-94.  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0342942/manifest

Comment

Related Items