CREATING SPATIAL PROBABILITY DISTRIBUTIONS FOR LONGLEAF PINE ECOSYSTEMS ACROSS EAST MISSISSIPPI, ALABAMA, THE PANHANDLE OF FLORIDA, AND WEST GEORGIA Except where reference is made to the work of others, the work described in this thesis is my own or was done in collaboration with my advisory committee. This thesis does not include proprietary or classified information. ________________________________________ John Scott Hogland Certificate of Approval: _________________________ _________________________ Ralph S. Meldahl Mark D. MacKenzie, Chair Associate Professor Assistant Professor Forest Biometrics Forest Biology and Ecology _________________________ _________________________ James B. Grand Nedret Billor Associate Professor Associate Professor Wildlife Sciences Mathematics and Statistics _________________________ Stephen L. McFarland Acting Dean Graduate School CREATING SPATIAL PROBABILITY DISTRIBUTIONS FOR LONGLEAF PINE ECOSYSTEMS ACROSS EAST MISSISSIPPI, ALABAMA, THE PANHANDLE OF FLORIDA, AND WEST GEORGIA John Scott Hogland A Thesis Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Master of Science Auburn, Alabama December 16, 2005 iii CREATING SPATIAL PROBABILITY DISTRIBUTIONS FOR LONGLEAF PINE ECOSYSTEMS ACROSS EAST MISSISSIPPI, ALABAMA, THE PANHANDLE OF FLORIDA, AND WEST GEORGIA John Scott Hogland Permission is granted to Auburn University to make copies of this thesis at its discretion, upon request of individuals or institutions and at their expense. The author reserves all publication rights. ______________________________ Signature of Author ______________________________ Date of Graduation iv THESIS ABSTRACT CREATING SPATIAL PROBABILITY DISTRIBUTIONS FOR LONGLEAF PINE ECOSYSTEMS ACROSS EAST MISSISSIPPI, ALABAMA, THE PANHANDLE OF FLORIDA, AND WEST GEORGIA John Scott Hogland Master of Science, December 16, 2005 (B.S., Auburn University, 2001) 139 Typed Pages Directed by Mark D. MacKenzie Longleaf ecosystems have severely decreased in total area since pre-European settlement. These dramatic losses are the principle reason for the listing of many plants and animals as endangered, and have been the driving factor for recent longleaf ecosystem restoration efforts. While studies have documented the regional decline of longleaf ecosystems, they provide little information on fine scale fragmentation patterns and current locations. This lack of information often limits the efficacy of longleaf ecosystem management, monitoring, and restoration. To aid longleaf restoration efforts we developed a series of fine grain (30 m) ecosystem probability distributions using multitemporal Landsat enhanced thematic mapper plus imagery, digital elevation models, field data, ancillary data sets, polytomous logistic regression, and a hierarchical classification scheme. Using our ecosystem v probability distributions, resource managers can identify the most probable locations for longleaf ecosystems, locate potential restoration sites, prioritize restoration efforts, and estimate ecosystem area. vi ACKNOWLEDGEMENTS The author would like to thank Dr. Mark D. MacKenzie, Dr. James B. Grand, Dr. Ralph S. Meldahl, and Dr. Nedret Billor, for their guidance and help, Melissa Reynolds, Kevin Kleiner, and John Gilbert for their insights and advice, and fianc?e Melissa and family members Butch, Judy, and Dan for their support throughout this project. vii Style manual or journal used Photogrammetric Engineering and Remote Sensing Computer software used SAS? version 8.2 ArcGIS? version 8.3 ESRI?s Spatial Analyst? extension ERDAS IMAGINE? version 7 Garmin eTrex Legend? GPS Microsoft; Excel? Word? PowerPoint? Access? viii TABLE OF CONTENTS LIST OF TABLES ix LIST OF FIGURES xii Introduction 1 Chapter 1: Comparing polytomous logistic regression and discriminant analysis: a remote sensing perspective John S. Hogland, Mark D. MacKenzie, and Nedret Billor 4 Chapter 2: Bringing images to a common radiometric scale using Aggregate No Change Regression: a comparison between radiometric normalization techniques John S. Hogland and Mark D. MacKenzie 36 Chapter 3: Creating longleaf ecosystem probability distributions across the Southeast John S. Hogland, Mark D. MacKenzie, Ralph S. Meldahl, James B. Grand, and Nedret Billor 66 Cumulative Bibliography 108 Appendix 1 120 Appendix 2 122 Appendix 3 124 ix LIST OF TABLES CHAPTER 1: Table 1. Cross walk between our stage 1 land cover classes and NLCD classes. 22 Table 2. Landsat ETM + Scene, date, and number of classes used to compare maximum likelihood classification and PLR classifications. 23 Table 3. Number of samples per class type for each Landsat ETM+ scene. Approximately 70% of the samples were used to generate each classification model, while the remaining 30% of the data were used to test each classification methodology. 24 Table 4. Model fit statistics for the PLR classification method using the Wetland class as the baseline response variable. Model one represents our top model. 25 Table 5. Maximum likelihood beta estimates, standard errors, chi-square values, p-values, and odds ratio for the top ranked probabilistic classification model. Band 3 for the deciduous and evergreen classes and band 6 for the evergreen class have p-values > 0.05. While these variables are not significant at ? = 0.05, for their respective class, they are significant for the other class, thus making them significant in the overall model. 26 CHAPTER 2: Table 1. Method, source, and normalization issues addressed 56 x Table 2. Initial mean regression estimates for relative normalization procedures. AS-NCR intercept and slope estimates correspond to changes in sensor reflectance. Values inside () indicate digital number (DN) equivalents. 57 Table 3. Elevation Distance and digital number values of selected dark objects used in ASR and DOS1 normalization techniques. 58 Table 4. Summarized validation statistics. Values represent averaged estimates across 11 scene comparisons. AS-NCR and ASR intercept estimates correspond to changes in sensor reflectance (RF). DOS1 intercept estimates corresponds to changes in surface RF. AG-NCR and NCR intercept estimates correspond to changes in digital number (DN) values. Values inside () indicate the DN equivalents for validation intercept and slope estimates for models that use RF. 59 CHAPTER 3: Table 1. Cross walk between our level 1 land cover classes and NLCD + classes. 88 Table 2. PLR model statistics for each level and stage of our hierarchical classification scheme (HCS). W and S/F identify winter and summer/fall seasonality of the imagery. 89 Table 3. Slope and intercept estimates for each stage of level 1 in our hierarchical classification scheme (HCS). The land cover class type ?Water? was used as xi the baseline category. W and S/F identify winter leaf-off and summer/fall seasons, respectively. 90 Table 4. Slope and intercept estimates for each stage of level 2 in our hierarchical classification (HCS). The forested ecosystem Coastal Plain Longleaf was used as the baseline class. W and S/F identify winter leaf-off and summer/fall seasons, respectively. 91 Table 5. Predicted 95% class count confidence intervals (CI) vs. observed class counts for level 1 validation data set. 92 Table 6. Predicted 95% confidence intervals (CI) for class counts vs. weighted observed class counts for level 2 validation data sets. The observed number of samples (actual) were weighted to remove the effect of prior probabilities. 93 xii LIST OF FIGURES CHAPTER 1 Figure 1. Cross sectional view of a hypothetical, multidimensional, 5 class probability distribution. The vertical dashed lines identify the location of each class? maximum likelihood allocation rule while holding variables X i-1 constant and allowing variable X 1 to vary. 27 Figure 2. Seasonality of each Landsat ETM+ scene and the spatial location of each sample point (total of 15 different scenes). 28 Figure 3. Estimated mean kappa values (measure of agreement among classes) and 95% confidence intervals for each Landsat ETM+ scene (Path Row) for the PLR and MLC methods. 29 Figure 4. Comparison between an observed hard classification accuracy assessment and a predicted hard classification accuracy assessment using the top ranked PLR probabilistic classification model and a MLAR. 30 CHAPTER 2: Figure 1. Scenes (path and row) and overlaps used in the comparison between normalization methods. 60 xiii Figure 2. Scatter plots for simulated, image overlapping regions where a slave image is a linear function of a master image and has a constant geometric rectification error of 1, 2, and 4 pixels. Expected intercept and slope values are 5 and 0.875 respectively. 61 CHAPTER 3 Figure 1. An overlay of the historical range of longleaf ecosystems (after Little, 1971) and our study area. The extent of our study area was defined as the counties within the state of Alabama and the counties that intersected United States Fish and Wildlife Service econumber 29 and 30 (USFWS, 2000). 94 Figure 2. An overlay of our study area (outlined in gray), and Landsat scenes by path row (PR; outlined in red) showing the dates of ETM+ image acquisition for a winter leaf-off (W), a summer/fall leaf-on (F), and a spring leaf-on (S) season. Dates followed by an M or X indicate master images and images that could not be brought to a common radiometric scale, respectively. N/A identifies scenes for which we were unable to obtain a suitable image representing a spring leaf-on season. 95 Figure 3. Hierarchical classification scheme for two seasons. W and S/F identify land use change (LUC) and temporal feature (TF) classes found in the normalized ETM+ imagery representing a winter leaf-off and summer/fall leaf-on season, respectively. 96 xiv Figure 4. Three maps illustrating the spatial configuration of forested ecosystem prior abundance estimates. Maps A, B, and C correspond to the prior abundance estimates, relative to a maximum abundance estimate of 100, for Mountain Longleaf, Coastal Plain Longleaf, and Slash ecosystems, respectively. Prior abundance estimates for Hardwood, Mixed Hardwood/Pine, and Loblolly ecosystems were set at 100 across our study area. 97 Figure 5. Level 1 Land cover and land use change probabilistic classifications. 98 Figure 6. Level 2 Forested ecosystem probabilistic classifications. 99 Figure 7. Average area (m 2 ) per pixel for Coastal Plain Longleaf Ecosystems within Blackwater State Forest. Blackwater State Forest is located in the northwest portion of the Florida panhandle and is depicted as the red area in the location map. 100 Figure 8. Probability distribution of the Coastal Plain Longleaf ecosystem within and around Conecuh National Forest located in southern Alabama (identified as the red area within the location map). 101 1 INTRODUCTION Currently, longleaf ecosystems are estimated to occupy 1.2 X 10 6 ha across the Southeastern United States, a mere 5% of the 24.3 X 10 6 ha pre-European settlement estimate (Outcalt and Sheffield, 1996). This dramatic loss of habitat has had a substantial impact on numerous plants and animals, and is the primary reason that many Southeastern species have been listed as threatened or endangered (Tuldge, 1999). Due to this loss of habitat, there is a strong need for the conservation and restoration of these critically endangered ecosystems (Noss et al., 1995). While conservation and restoration efforts have begun, they have been limited, in part, by the lack of information depicting the current location of these ecosystems. Long-term studies such as the Forest Inventory Analysis have been useful in identifying trends in longleaf ecosystem decline (Kelly and Bechtold, 1990; Outcalt and Sheffield, 1996), but are ill-suited to provide meaningful information at fine spatial scales. Due to the coarse nature of these data sets (e.g., 20 km grain), organizations have had to take a broad based approach towards longleaf ecosystem management, monitoring, and restoration, often limiting the efficacy of their efforts. To become more effective, these organizations need accurate, fine scale data sets that identify forested ecosystem types and depict the current location and distribution of longleaf ecosystems. Remotely sensed data (e.g., satellite imagery and digital elevation models) provide a unique opportunity to generate such a data set by linking fine grain (30 m) 2 spectral information with spatially explicit examples of forested ecosystem types. Due to the amount of spectral overlap among coniferous ecosystems in the Southeast, though, few analysts have successfully differentiated longleaf ecosystems from other coniferous ecosystems using common classification (e.g., maximum likelihood classifiers, clustering, classification trees, and artificial neural networks) and radiometric normalization (e.g., at-sensor reflectance, dark object subtraction, and ridge regression) techniques. This suggests either one of two scenarios: there are no differences (spectrally or in elevation) between longleaf and other coniferous ecosystems, or current methodologies may be too restrictive and/or inappropriate to depict the differences between coniferous ecosystems in the Southeast. Given the visual and structural differences of longleaf ecosystems (e.g., relatively sparse overstories, and diverse understories, composed of many grasses and forbs), we believed the latter scenario to be the case. Therefore, we adopted a probabilistic classification procedure (polytomous logistic regression), created a new radiometric normalization procedure, and implemented these procedures with a hierarchical classification scheme to differentiate between forested ecosystems types. While these procedures theoretically provided the flexibility and accuracy needed to distinguish among forested ecosystems, they have not been thoroughly tested against commonly used classification and normalization techniques in remote sensing. Therefore, we compared each procedure against other commonly used classification and normalization procedures. To present our findings in a clear and concise manner we separated each comparison into distinct chapters. Chapter 1 directly compares polytomous logistic regression (PLR) with a maximum likelihood classifier, also known 3 as discriminant analysis (i.e., linear and quadratic discriminant functions), to determine if there are any benefits to adopting the PLR process. Chapter 2 directly compares our radiometric normalization procedure, Aggregate No Change Regression, with 4 other radiometric normalization procedures in terms of bringing Landsat enhanced thematic mapper plus imagery to a common radiometric scale. Finally, in Chapter 3, we introduce our hierarchical classification scheme and describe how the procedures introduced in Chapters 1 and 2 were used to predict the distributions of longleaf ecosystems. Each chapter can be read independently and is styled in publication format according to the editor?s guidelines of Photogrammetric Engineering and Remote Sensing. Our findings indicate that, on average, longleaf ecosystems are spectrally different than other forested ecosystems. Furthermore, by incorporating our newly developed radiometric normalization technique and adopting PLR in a hierarchical framework, these differences can be used to accurately and precisely predict the probability distribution of longleaf and other forested ecosystems, thereby providing resource managers with the information needed to begin addressing fine scale questions pertaining to longleaf ecosystem management, monitoring, and restoration. 4 CHAPTER 1 COMPARING POLYTOMOUS LOGISTIC REGRESSION AND DISCRIMINANT ANALYSIS: A REMOTE SENSING PERSPECTIVE ABSTRACT Maximum likelihood classification, also known as discriminant analysis, is a popular supervised technique used by remote sensing analysts. While this classification procedure has been embraced by the remote sensing community, it has some distinct drawbacks such as being limited to continuous data, the assumptions of multivariate normality and equal covariance, limited modeling diagnostics, and few techniques that address model parsimony. An alternative classification approach, which is less restrictive and has been successfully used in other fields to distinguish among class types, is called polytomous logistic regression (PLR). To assess the utility of PLR in image classification, we compared and contrasted PLR with the traditional maximum likelihood procedure (using linear and quadratic discriminant functions). Our findings indicate PLR is a flexible alternative to the traditional maximum likelihood classification. 5 INTRODUCTION Remote sensing analysts commonly use a maximum likelihood classification technique to perform supervised classifications (Lillesand and Kiefer, 1994; Erdas, 1997; Jensen, 2000). This technique, also referred to as dicriminant analysis in the statistical literature (Press and Wilson, 1978; Johnson and Wichern, 2002), distinguishes among classes by estimating the distance between each class, using either a linear or a quadratic discriminant function, for a given set of explanatory variables (Erdas, 1997; Johnson and Wichern, 2002; Metternicht, 2003). From these distance measures, analysts either create a hard classification (i.e., maximum likelihood classification) by generating rules that allocate class types to new observations based on minimizing class distance or maximizing posterior probability (Jensen, 1986; Johnson and Wichern, 2002), or by incorporating posterior probabilities into a fuzzy classification (Foody, 1996; Benz et al., 2004; Metternicht, 2003). Discriminant analysis, though, assumes multivariate normality and equal covariance for linear dicriminant functions, and multivariate normality for quadratic functions. As these assumptions are often difficult to satisfy (Press and Wilson, 1971; Foody, 1996), discriminant analysis tends to overestimate the magnitude of association among classes (Halpern et al., 1971; Press and Wilson, 1978; Hosmer et al, 1983; Hosmer and Lemeshow, 1989) and produces misleading posterior probabilities (Press and Wilson, 1978; Hosmer and Lemeshow, 2000; Johnson and Wichern, 2002). To circumvent these issues, analysts have employed a number of different techniques ranging from ignoring the statistical assumptions of discriminant analysis and performing a maximum likelihood classification (Ram?rez-Garc?a et al., 1998; Keuchel et al., 2003; Galv?o et al, 2005) to using alternative procedures that do not rely on the 6 assumption of multivariate normality (Yoshida and Omatu, 1994; Gopal and Woodcock, 1996; Brown de Colstoun et al., 2003; Pal and Mather, 2003). The accuracy of their classification is then assessed using a different analysis (e.g., class proportions and/or a kappa statistic), which compares the number of times classes were correctly identified to the number of times they were misclassified (Foody, 2002). Assessing class accuracy using class proportions and kappa statistics (Foody, 2002; Agresti, 2002) assumes a multinomial distribution and large sample normality, respectively, (Agresti, 2002) which are typically easier to satisfy. An alternative classification technique, which assumes a multinomial distribution and large sample normality, is polytomous logistic regression (PLR). This technique has been successfully used in numerous other fields including geography (Wrigley, 1985), engineering (Hasegawa and Kurita, 2002), biological and molecular sciences (Bailey et al., 2003), education (Peng and Nichols, 2003), and environmental sciences (Mahopatra and Kant, 2005). Some of the benefits of using a PLR classification include the ability to use probabilistic classifications, it has relatively few statistical assumptions, the ability to use both continuous and categorical data as explanatory variables, and it focuses on directly modeling class probabilities (Agresti, 2002). Using PLR, hard classifications can be generated by setting probability thresholds ranging from identifying minimum or maximum probabilities for each class to using a maximum likelihood allocation rule (MLAR) for each observation (i.e., the same allocation rule used in discriminant analysis). PLR hard classifications that use a MLAR, however, may not always provide good separation among classes when there is substantial overlap among explanatory 7 variables and can bias class area estimates. For example, consider the hypothetical probability distribution of 5 class types illustrated in Figure 1. In this scenario, using a MLAR (represented by the dashed lines) to classify each observation results in substantial misclassification (misclassification rate > 25%) between the following classes: Water Deciduous, Deciduous and Evergreen, Evergreen and Field, and Field and City. Moreover, class area estimates, which are calculated by summing the multiple of the number of observations allocated to each class (using a MLAR) by the area of each observation, are substantially biased because misclassification errors are far from symmetric for the majority of class comparisons. Using class probabilities, however, increases the ability to distinguish among class types at different ranges of the explanatory variable. In addition, unbiased estimates of class area can be calculated by weighting the area of each observation by the probability of each class and then summing the weighted area estimates by class for all observations. Many have made the argument that hard classifications do not always adequately describe class transitions (Foody, 1996; Benz et al., 2004; Metternicht, 2003). To address this issue they have incorporated posterior probabilities, calculated from the discriminant analysis procedure, into their classification. These probabilities, however, are biased when statistical assumptions are not met (Hosmer and Lemeshow, 2000; Johnson and Wichern, 2002). In contrast, the PLR method does not have the same restrictive assumptions as discriminant analysis and directly models class probabilities, potentially providing analysts with a more accurate estimate of class probabilities. To illustrate the difference between posterior probabilities of discriminant analysis and PLR probabilities, it is useful to look at how each method derives these 8 estimates. In discriminant analysis, class probabilities are indirectly calculated based on class distance measurements and the assumption that each class has a multivariate normal distribution. Maholanobis distances (MD), quadratic discriminant score (QDS), and posterior probabilities (PP) are calculated as follows: ()( ) p1,..., jn 1,...,i )(MD i ==? ? ?== ? xxSxx 1 ij iD (1) () ()() ,..., lnlnQDS class jprior for p nj piD j jjijjijj = = +? ? ???== ? 1 2 1 xxSxx 2 1 S 1 (2) () )( 2 1 exp )( 2 1 exp | 1 ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? = j i j j iD iD ijPP (3) where x i and x are the observed vector and mean vector, respectively, and S and S j are the pooled sampled covariance matrix and the sample covariances, respectively (Johnson and Wichern, 2002). The PLR method, however, directly estimates class probabilities and assumes that multi-category responses have a multinomial distribution with asymptotic errors around the linear form of the natural log transformation of class odds (logits). Estimated class (response) probabilities {? j (x)} are determined by manipulating baseline category logits as follows: given that ( ) , , 111 ,...,J-j j j == ? x? 9 then () ().-1 hJ xx ? ? = = 1 1 J h ?? Using the linear form of the logit for p covariates and a constant term, denoted by the vector, x, of length p + 1 ( ) () []1,...pi , ln =?=??= ? ? ? ? ? ? ? ? ij J j xx?x x x 1 ? ? and substituting () Jh for -1 ?? x ? ? = 1 1 J h ? j and ? J are solved as follows: () ( ) () ? ? = ?+ ? = 1 1 1 J h h j j ?x ?x x exp exp ? (4) () () exp ? ? = ?+ = 1 1 1 1 J h h J ?x x? (5) Where ? j (x) is the mean probability of group j, with class means and variances equal to n? j and n? j (1-? j ), respectively (Agresti, 2002). Maximum likelihood estimates of beta values are determined in an iterative fashion, typically using the Newton-Raphson or Fisher scoring methods, and the standard errors for each beta are based on profile likelihood functions or asymptotic normality (Agresti, 1990). While the PLR method has a number of benefits over the classical supervised approach, there are some drawbacks. The first drawback deals with the efficiency of the PLR method when the multivariate normality assumption holds (Bull and Donner, 1987). Multivariate normality, however, rarely holds in remote sensing (Foody, 1996). 10 Secondly, PLR cannot obtain a maximum likelihood estimate when there is no overlap in class explanatory values. While an unsolvable maximum likelihood estimate may be troubling in terms of mathematic complexity and model fit estimates (Agresti, 2002), viewed from a classification perspective this situation means that some of the class types can be separated from the rest of the class types with 100% accuracy given a set of rules. In this situation, a probabilistic classification is not required. Instead, class types can be assigned using means and/or Mahalanobis distance. For classes that do have overlap in the explanatory values, a maximum likelihood estimate can be solved and a probabilistic classification can be generated. From a theoretical standpoint, PLR is a very robust classification technique that should provide a better depiction of class distributions when compared with discrimant analysis. However, few have directly compared these techniques (Bull and Donner, 1987; Hossain et al., 2002) and to our knowledge, no one has applied PLR in a remote sensing framework. Therefore, we compared and contrasted the two techniques from the standpoint of a hard classification (using a MLAR) and a probabilistic classification. We then further demonstrate the utility and flexibility of PLR using an example of 3 land cover types that share a significant portion of spectral space. METHODS To contrast these two classification techniques, we compared hard classifications generated by the PLR method using a MLAR, with hard classifications generated by both linear and quadratic discrimant functions, for 15 different Landsat enhanced thematic mapper plus (ETM+) scenes. Landsat ETM+ scenes were rectified using Multi-Resource Land Characteristics Consortium data processing level 1t (NASA, 2005) procedures. 11 Class types consisted of generalized National Land-Cover Database categories (Table 1; after Homer et. al., 2004) and temporal features (i.e. clouds, burnt areas, shadows, and smoke) and were visually interpreted at the spatial scale of one ETM+ pixel using ETM+ imagery and digital ortho quarter quads. To maintain consistency across the study area, one image interpreter identified all class types. Classification errors potentially caused by image acquisition dates representing different seasons were accounted for by randomly selecting 1 of 3 phonologies for each scene comparison; leaf-off winter season, leaf-on spring growing season, leaf-on fall season (Figure 2; Table 2). Digital number (DN) values occurring at the same spatial location as image interpreted samples were extracted on a nearest pixel basis, by band, using the sample command in Environmental Systems Research Institute?s (ESRI) Spatial Analyst extension (ESRI, 2005). Samples were then randomly partitioned into a training (~ 70% of the data for each scene) and validation (~ 30% of the data for each scene) data set and imported into Statistical Analysis Software (SAS) version 8.2 to perform all analyses (Discriminant and Logistic procedures were used to perform the maximum likelihood classification and PLR classifications, respectively). Total sample size, by class and scene, are listed in Table 3. In studies that have contrasted PLR and discrimant analysis (Bull and Donner, 1987; Hossain et al., 2002), only 3 classes were used to compare the two methods. In remote sensing, however, there are often many more than 3 classes. To determine whether the number of classes had an appreciable effect on the accuracy of the hard classification method, as many as 10 class types were sampled from some scenes, while 12 other scenes had as few as 4 class types sampled. In all, 10,350 image-interpreted samples were collected across 15 different scenes. Validation accuracies were used to estimate the level of agreement between observed and predicted classes using kappa statistics (Foody, 2002; Erdas, 1997; Agresti, 1990). Each classification method?s mean kappa estimate and corresponding lower and upper kappa confidence limits were compared on a scene-by-scene basis to determine if the number of classes had an appreciable effect on the accuracy of a particular classification method. To test general accuracy trends, we compared mean estimates of kappa for all scenes among all classification methods using a one-way nonparametric analysis of variance. In situations where there was complete spectral separation among some of the classes, the logistic procedure was allowed to continue using the last maximum likelihood iteration to determine fit statistics. While fit statistics in these situations may be misleading due to an unsolvable maximum likelihood estimate, predicted class types, based on probabilities and a MLAR, can still be used to compare the class accuracy of the two methods. After comparing the classification methods using a hard classification, the PLR procedure was performed again for a randomly chosen scene (path/row 21/37) and the classes Evergreen, Deciduous, and Wetland to develop and interpret a parsimonious probabilistic classification. Training data for this scene were used to generate a suite of classification models, from which the best fitting and most parsimonious model was selected using Akaike?s Information Criterion (AIC; Akaike, 1973). The top ranked model was then used to perform a probabilistic classification for that scene. To validate 13 our probabilistic classification, we compared the summed class frequency estimates (lower and upper 95% confidence limits) against observed class frequencies using the validation data set for that scene. Class frequency confidence intervals were calculated for each observation using the delta method (Agresti, 1990; Agresti, 2002). Under this scenario, a PLR model that contains a predefined proportion (i.e., 95%) of observed class frequencies within its predefined confidence interval (i.e, 95% confidence interval), suggests a good model fit and that the model can be generalized to the rest of the population. RESULTS Hard Classification Kappa estimates for the PLR hard classifications and Quadratic maximum likelihood classification were significantly larger than kappa estimates for linear maximum likelihood classification for scenes 19/38 and 20/38 (Figure 3). Additionally, the PLR hard classifications had a significantly larger kappa value when compared to the linear maximum likelihood classification method for scenes 20/37 and 20/39. There were no significant differences in mean kappa estimates, however, between PLR hard classifications and quadratic maximum likelihood classification. These results suggest that the number of classes and the sample size of each comparison did not have an appreciable effect on the accuracy of any particular classification method (Table 2; Figure 3). In addition, there were no significant differences in hard classification accuracy trends among the different classification methods (all method comparisons: Kruskal-Wallace X 2 df = 2 = 1.6935; p-value = 0.43, PLR vs Linear maximum likelihood 14 classification: Kruskal-Wallis X 2 df =1 = 0.7230; p-value = 0.40, PLR vs Quadratic maximum likelihood classification Kruskal-Wallis X 2 df = 1 = 0.1897; p-value = 0.66). Probabilistic Classification Reapplying the PLR classification method for the classes Evergreen, Deciduous, and Wetland in scene 21/37, we found only 6 spectral bands were needed to sufficiently describe the probability transition among classes (Table 4). Using Landsat ETM+ bands 2 through 7, we were able to generate a statistically significant model (X 2 df = 12 = 227.59; p-value < 0.0001) that explained the majority of information within our training data (max rescaled 8211.0 ~ 2 =R ; SAS 2005). Interpreting the beta estimates, in terms of odds ratio (calculated as follows; ratio, odds= ? ij x e where B j = the slope estimate for band j and x i = the DN value for band j ), we identified the effect each Landsat ETM+ band had on class probabilities (Table 5). For example, with every incremental increase in DN values for band 4, while holding bands 2, 3, 5, 6, and 7 constant, the average odds (chances) of the classes Deciduous and Evergreen increased by multiples of 1.455 and 1.202, respectively, when compared with the Wetland Class. Comparing the Deciduous and the Evergreen class, however, was less straightforward. Odds ratio between the Deciduous and Evergreen classes were calculated as follows: Given that the log odds of, ; ln,ln j Wetland Deciduous Wetland Evergreen and ?x?= ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? where i = DN values for bands 2 ? 7 and J= classes Evergreen and Deciduous, and that, ()() EvergreenDeciduous Wetland Deciduous Wetland Evergreen Deciduous Evergreen ?x?x ???= ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? lnlnln (7) 15 Using equation 7 and exponentiating the results (Agresti, 2002), the estimated average odds of the class Deciduous increased by a multiple of 1.211 for every incremental increase in band 4 DN values (Table 4). Using the delta method (Agresti, 1990; Agresti, 2002) to generate 95% confidence intervals for each observation in the validation data set and summing the lower and upper probability estimates for all observations, we predicted (within the lower and upper confidence limits) the total number of observations for each class (Figure 4). In addition, we were able to predict (within the lower and upper bounds of a 95% confidence interval) the cell counts of a hard classification accuracy assessment by using a MLAR and summing the lower and upper probability confidence interval for each class (Figure 4). User accuracies (row proportions; Figure 4), were calculated by dividing the area under a specified class distribution (the sum of class probabilities for a given class type), within the bounds of a given MLAR, by the total area of that class distribution (the sum of class probabilities for a given class type across all MLARs). Producer accuracies (column proportions; Figure 4) were calculated by dividing the area of each specified class distribution, within the bounds of each MLAR, by the area of all class distributions within the bounds of that MLAR. DISCUSSION For all 15 scene comparisons, hard classification accuracies were similar among PLR, linear discriminant analysis, and quadratic discriminant analysis. Similar results have been documented in other studies (binary case: Press and Wilson, 1978; multinomial case: Hossain et al., 2002). Under different circumstances, such as when categorical variables provide important class separating information, the PLR 16 classification would have a significant advantage over the linear or quadratic maximum likelihood classification method given that categorical variables can be directly incorporated into a PLR model (Hosmer et al., 1983). In cases where maximum likelihood classification assumptions are met, both the PLR and maximum likelihood classification methods will give equivalent results (Hosmer and Lemeshow, 2000), the primary difference among the methods is the efficacy of the algorithms (binary case: Efron, 1975; multinomial case: Bull and Donner, 1987). Historically, the PLR classification method may have been less appealing due to the amount of processing time required to calculate a maximum likelihood estimate. With today?s computer processors, however, this is no longer an issue. When performing a maximum likelihood classification, sometimes analysts ignore the statistical assumptions of the classification technique and rely on an independent accuracy assessment to estimate classification accuracy. While an accuracy assessment does assess the hard classification rule, it does not consider issues such as over fitting the data, modeling assumptions, multicollinearity, or model parsimony in the initial statistical model. Failure to address these issues can mask the relationship between explanatory variables and response variables, reducing the overall accuracy of the classification, and limiting our understanding of the driving components of the model (Press and Wilson, 1978; Foody, 1996; Hosmer and Lemeshow, 2000; Johnson and Wichern, 2002). The PLR method differs from the discriminant analysis method in its modeling assumptions, its ability to incorporate both categorical and continuous variables, its focus on directly modeling class probabilities, and its ability to estimate model error (i.e., error 17 in the estimated probabilities). While the discriminant method can be used to derive posterior probabilities, these probabilities are based on the assumption of multivariate normality, which seldom holds for remote sensing data (Foody, 1996). Moreover, posterior probabilities only represent the mean probability given a set of explanatory variables values. By definition, they do not include measures of probability error. In contrast, the PLR method estimates probability error, which provides analysts and end users with a much finer level of detail. In essence, a PLR classification is not limited to one map. Instead, a probabilistic map, with corresponding lower and upper confidence limits can be produced for each respective class and observation. For example, each observation in scene 21/37 represents a Landsat ETM+ pixel. Class probabilities for each pixel in that scene are interpreted as the mean proportion of times that one would expect to find each class at a pixel with given spectral values. Assume there were 100 pixels with DN values of 47, 36, 65, 69, 124, and 33 for bands 2 through 7, respectively. Using equations 4 and 5 and the delta method, we would expect on average 28 (12-45, 95% confidence intervals), 7 (0-17, 95% confidence intervals), and 65 (47-83, 95% confidence intervals) pixels out of the original 100 pixels to be Evergreen, Deciduous, and Wetland, respectively. These types of maps not only depict the probability of each class for each pixel, but also maintain the modeling errors associated with class probabilities, which can be directly incorporated into other predictive models. A number of techniques exist to validate PLR models (Neter et al., 1996). For example, a classical model validation approach is to generate models for both training and validation data sets and compare beta estimates and significant model parameters. 18 The difficulty with this approach is that it requires a substantial amount of data. We demonstrated (Figure 4) how conditional probabilities and probability confidence interval can be used to assess the accuracy of a PLR model when sample size is limited. One can then infer whether PLR probabilities can be generalized to the rest of the population (i.e., the rest of the pixels inside a Landsat ETM+ image) by comparing these predicted estimates with observed independent data. Evaluating PLR fit statistics and taking into consideration model parsimony, Landsat ETM+ bands 2 through 7 were useful in separating the classes Evergreen, Deciduous, and Wetland in scene 21/37. ETM+ band 1 did not provide any additional information and was subsequently removed from the analysis. Interpreting the results for scene 21/37 in terms of scaled energy (i.e., DN values) and odds ratios, we see that as the amount of energy being reflected and transmitted from the earth?s surface in ETM+ bands 3 and 6 increases (wavelength ranges of 0.630 ?m - 0.690 ?m and 10.40 ?m ? 12.50 ?m, respectively), while the energy being reflected and transmitted in band 2 remains low (wavelength range 0.525 ?m ? 0.605 ?m) the odds of the Evergreen class increase (Table 4). Conversely, as the amount of energy being reflected and transmitted in bands 3 and 6 decrease while maintaining band 2 at low levels the odds of the Deciduous class increase. These results have biological implications. For example, given the date of image acquisition and that certain objects use different portions of the electromagnetic spectrum, these results may suggest that in the spring, Deciduous class types absorb a relatively large portion of visible light (bands 2 and 3) when performing photosynthesis, 19 and that these class types maintain cooler daily temperatures (band 6) when compared with Evergreen and Wetland classes. Whether remote sensing analysts are primarily focused on developing a hard classification, building a series of probability maps, or understanding the driving components of a classification, it is critical that modeling assumptions are checked. Failing to do so may result in biased classification models. If the assumptions of a classification method are met, such as in the PLR example for scene 21/37, then class probabilities can be used as an alternative to a hard classification. Moreover, these class probabilities (and estimates of error) can provide a more accurate depiction of class area (and area confidence intervals). This is not to imply that every pixel has a certain proportion of its area allocated to each class. Instead, it indicates that, on average, we expect a certain number of pixels to be allocated to each class type, given a set of explanatory values. This means that as the number of pixels used to calculate class area for a specified region increases, the accuracy of our class area estimate should also increase. Identifying class location, however, is less straightforward. Given that class probabilities do not represent the amount of class area within each pixel, a class allocation rule must be specified to identify class location. Arguably, one of the easiest rules to apply is a MLAR. While this type of rule will ensure that every pixel is classified into one of the predefined categories, in some instances this rule may allocate class types to pixels where the probability of being a specific class could be as low as 100% divided by the number of classes. Alternatively, with the PLR method users can select class probability thresholds that identify class locations. While this approach may 20 leave some pixels unclassified, it provides users of the data set with a level of precision (in terms of probabilities) for each class location. By maintaining class probabilities for each pixel, data set users have the flexibility to address numerous different scenarios. For example, if a user was interested in identifying locations in scene 21/37 where there was a high probability of the class Evergreen (x > 66%) and medium to high probability of the class Wetland (33% < x < 66%), they could perform a GIS query to find such locations. To determine the accuracy of each potential class, a user could apply GIS zonal functions (ESRI, 2005) to average the summed probabilities of each class combination. CONCLUSION The PLR technique has desirable qualities for classifying remotely sensed data which include relatively unrestrictive model assumptions, the ability to incorporate both continuous and categorical variables directly into the classification scheme, relatively easy techniques for model comparison (e.g., AIC), an intuitive relationship between class types and explanatory variables (odds ratios), and a focus on directly modeling class probabilities (model error). The outputs from the PLR method are class probability distributions. Using class probability distributions and estimates of model error, users can accurately calculate the amount of area for each class type within a predefined level of confidence. In addition, they can query each class probability distribution to generate a wide variety of maps identifying class location that are tailored to specific questions. These maps can include a MLAR map, but are not limited to this type of map. Overall, the PLR classification method provides an extremely flexible alternative to the classical supervised approach. 21 ACKNOWLEDGEMENTS We would like to acknowledge the following people who have provided valuable comments and insight in this paper; Ralph Meldahl, Barry Grand, Melissa Reynolds, Kevin Kleiner, and John Gilbert. Funding for this study was provided by Alabama Gap Analysis Project. 22 Table 1. Cross walk between our stage 1 land cover classes and NLCD + classes. Land Cover Classes NLCD Classes Water 11 Open Water 21 Developed, Open Space 22 Developed, Low Intensity 23 Developed, Medium Intensity 24 Developed, High Intensity 31 Barren Land (Rock/Sand/Clay) Urban / Transportation / Bare Ground 32 Unconsolidated Shore 41 Deciduous Forest Forested * 42 Evergreen Forest 43 Mixed Forest Deciduous Evergreen 52 Shrub/Scrub 71 Grassland/Herbaceous 81 Pasture/Hay Field 82 Cultivated Crops 90 Woody Wetlands Wet Vegetated Area 95 Emergent Herbaceous Wetlands + NLCD classes 12, 51, 72, 73, and 74 were not applicable to our study *Forested classes where split into 2 groups Evergreen and Deciduous based on the dominance of "Evergreen" and "Deciduous" trees. 23 Table 2. Landsat ETM + Scene, date, and number of classes used to compare maximum likelihood classification and PLR classifications. Scene (Path/Row) Season Date Number of Samples Number of Classes 19/37 Spring 4/5/2000 251 5 19/38 Spring 4/5/2000 2155 8 19/39 Winter 12/20/2001 345 5 20/36 Spring 6/18/2001 579 5 20/37 Fall 10/8/2001 1371 7 20/38 Fall 10/8/2001 1421 10 20/39 Winter 1/25/2001 1005 6 21/36 Winter 3/5/2001 406 4 21/37 Spring 4/19/2000 897 9 21/38 Winter 2/15/2000 633 8 21/39 Spring 5/24/2001 318 7 22/36 Fall 10/3/2000 88 4 22/37 Spring 5/15/2001 405 8 22/38 Spring 5/15/2001 282 7 22/39 Fall 11/7/2001 194 6 Table 3. Number of samples per class type for each Landsat ETM+ scene. Approximately 70% of the samples were used to generate each classification model, while the remaining 30% of the data were used to test each classification methodology. Scene (Path Row) Burn City Clouds Deciduous Evergreen Field Shadow Smoke Water Wetland 19/37 44 46 0 0 48 0 0 47 66 0 19/38 0 310 262 249 313 268 267 218 268 0 19/39 59 0 0 0 71 73 0 0 71 71 20/36 0 0 110 111 109 139 0 0 110 0 20/37 192 195 191 197 204 196 0 0 196 0 20/38 128 128 128 128 240 136 133 140 128 132 20/39 170 157 0 169 172 0 0 0 168 169 21/36 0 0 0 88 87 146 0 0 85 0 21/37 86 240 80 77 81 78 0 74 101 80 21/38 65 70 0 63 63 63 0 65 64 180 21/39 44 44 0 50 45 44 0 0 47 44 22/36 0 0 22 0 0 22 0 0 22 22 22/37 69 44 78 44 34 45 0 0 37 44 22/38 41 40 0 40 40 0 0 41 40 40 22/39 33 0 29 41 33 0 29 0 0 29 24 Table 4. Model fit statistics for the PLR classification method using the Wetland class as the baseline response variable. Model one represents our top model. Model # Landsat ETM+ Bands Used in the Model DF AIC ?AIC AIC Model Weight X 2 Nested Comparison P-value 1 bands 2-7 12 182.277 0.000 0.500 2 vs 1 0.14266 2 bands 1-7* 14 182.382 0.105 0.474 NA NA 3 bands 2, 4-7 10 188.195 5.918 0.026 3 vs 2 3 vs 1 0.07921 0.07048 4 bands 1-6 12 204.660 22.383 0.000 4 vs 2 0.00000 5 bands 2-6 10 207.000 24.723 0.000 5 vs 2 0.00000 6 bands 3-7 10 219.876 37.599 0.000 6 vs 2 0.00000 * Model with most parameters (i.e., all Landsat ETM+ Bands) 25 Table 5. Maximum likelihood beta estimates, standard errors, chi-square values, p-values, and odds ratio for the top ranked probabilistic classification model. Band 3 for the deciduous and evergreen classes and band 6 for the evergreen class have p-values > 0.05. While these variables are not significant at ? = 0.05 for their respective class, they are significant for the other class, thus making them significant in the overall model. Variables Class DF Beta Estimates Standard Error X 2 P-value Odds ratio (Class/Wetland Odds ratio (Deciduous/Evergreen) Intercept deciduous 1 263.2 59.218 19.7474 <.0001 Exp(263.2) * Intercept evergreen 1 66.4517 23.3085 8.128 0.0044 Exp(66.45) * Exp(196.75) * Band 2 deciduous 1 -1.2422 0.3525 12.4202 0.0004 0.289 * Band 2 evergreen 1 -1.0443 0.2176 23.0397 <.0001 0.352 * 0.820 Band 3 deciduous 1 -0.5096 0.2806 3.2996 0.0693 0.601 Band 3 evergreen 1 0.2172 0.1315 2.7259 0.0987 1.243 0.483 * Band 4 deciduous 1 0.3752 0.1062 12.4725 0.0004 1.455 * Band 4 evergreen 1 0.1839 0.0593 9.6166 0.0019 1.202 * 1.211 Band 5 deciduous 1 -0.3569 0.1599 4.9818 0.0256 0.7 * Band 5 evergreen 1 -0.4247 0.1066 15.8722 <.0001 0.654 * 1.070 Band 6 deciduous 1 -1.8291 0.4746 14.8544 0.0001 0.161 * Band 6 evergreen 1 -0.2898 0.1794 2.6094 0.1062 0.748 0.215 * Band 7 deciduous 1 1.1635 0.3278 12.5954 0.0004 3.201 * Band 7 evergreen 1 0.8261 0.2043 16.3519 <.0001 2.284 * 1.401 * Odds ratio that is significantly different than 1 at ? = 0.05 26 Figure 1. Cross sectional view of a hypothetical, multidimensional, 5 class probability distribution. The vertical dashed lines identify the location of each class? maximum likelihood allocation rule while holding variables X i-1 constant and allowing variable X 1 to vary. 0% 100% Water Deciduous Evergreen Field City Class Probability Variable X 1 27 Figure 2. Seasonality of each Landsat ETM+ scene and the spatial location of each sample point (total of 15 different scenes). 28 017034085 Kilometers ? Scene Seasonality Winter Fall Spring Samples 20/36 22/36 21/36 20/37 22/37 19/3721/37 21/3822/38 20/38 19/38 19/39 20/39 22/39 21/39 0.4 0.5 0.6 0.7 0.8 0.9 1 1937 1938 1939 2036 2037 2038 2039 2136 2137 2138 2139 2236 2237 2238 2239 Landsat ETM+ Scene (Path Row) Ka p p a e s t i m a t e PLR Method Linear MLC Method Quadratic MLC Method Figure 3. Estimated mean kappa values (measure of agreement among classes) and 95% confidence intervals for each Landsat ETM+ scene (Path Row) for the PLR and MLC methods. 29 30 Figure 4. Comparison between an observed hard classification accuracy assessment and a predicted hard classification accuracy assessment using the top ranked PLR probabilistic classification model and a MLAR. Predicted Classes Deciduous Evergreen Wetland Totals Counts Row Proportion Column Proportion Lower predicted (95% CI) Observed (predicted) Upper predicted (95% CI) Lower predicted (95% CI) Observed (predicted) Upper predicted (95% CI) Lower predicted (95% CI) Observed (predicted) Upper predicted (95% CI) Lower predicted (95% CI) Observed (predicted) Upper predicted (95% CI) 15.49 18.00 (18.48) 20.98 0.17 0.00 (2.03) 4.43 0.00 4.00 (0.91) 2.32 15.66 22.00 (21.42) 27.73 0.72 0.82 (0.86) 0.98 0.01 0.00 (0.09) 0.21 0.00 0.18 (0.04) 0.11 0.73 1.00 1.30 Deciduous 0.74 0.86 (0.88) 1.00 0.01 0.00 (0.11) 0.23 0.00 0.17 (0.04) 0.10 0.00 2.00 (1.56) 3.72 11.36 18.00 (14.60) 17.42 0.57 4.00 (3.08) 5.81 11.93 24.00 (19.24) 26.95 0.00 0.08 (0.08) 0.19 0.59 0.75 (0.76) 0.91 0.03 0.17 (0.16) 0.30 0.62 1.00 1.40 Evergreen 0.00 0.10 (0.07) 0.18 0.60 0.95 (0.77) 0.92 0.02 0.17 (0.13) 0.24 0.00 1.00 (0.96) 2.36 0.52 1.00 (2.37) 4.55 16.43 16.00 (20.01) 23.08 16.95 18.00 (23.34) 29.99 0.00 0.06 (0.04) 0.10 0.02 0.06 (0.10) 0.19 0.70 0.89 (0.86) 0.99 0.72 1.00 1.28 Obse rve d Cla sses Wet land 0.00 0.05 (0.05) 0.11 0.03 0.05 (0.12) 0.24 0.68 0.67 (0.83) 0.96 15.49 21.00 27.06 12.05 19.00 26.4 17 24.00 31.21 44.54 64.00 84.67 Totals 0.74 1.00 1.29 0.64 1.00 1.39 0.70 1.00 1.30 31 REFERENCES Agresti, A. 1990. Categorical data analysis. Wiley-Interscience. New York. 558 pp. Agresti, A. 2002. Categorical data analysis. Wiley-Interscience. Hoboken, New Jersey. 710 pp. Akaike H. 1973. Information theory and an extension of the maximum likelihood principle. In: B.Petrov and F. Cazakil editers, Proceedings of the 2 nd international Symposium on information Theory, Aakademiai Kidao Budapest. Bailey, N., T. Clemets, J. T. Lee, and S. Thompson, 2003. Modeling soil series data to facilitate targeted habitat restoration: a polytomous logistic regression approach. Journal of Environmental Management. 67: 395-407. Benz, U.C., P. Hofmann, G. Willhauck, I. Lingenfelder, and M. Heynen. 2004. Multi- resolution, object-oriented fuzzy analysis of remote sensing data for GIS-ready information. Journal of Photogrammetry and Remote Sensing. 58: 239-258. Brown de Colstoun, E. C., M. H. Story, C. Thompson, K. Commisso, T. G. Smith, J. R. Irons. 2003. National park vegetation mapping using multitemporal landsat 7 data and a decision tree classifier. Remote Sensing of Environment. 85: 316-327. Bull, S. B. and A. Donner. 1987. The efficiency of multinomial logistic regression compared with multiple group discriminant analysis. Journal of the American Statistical Association, 82(400): 1118-1122. Efron, B. 1975. The efficiency of logistic regression compared to normal discriminant analysis. Journal of the American Statistical Association. 70(352): 892-898. ERDAS IMAGINE. 1997. ERDAS field guide; fourth edition. ERDAS Inc. Atlanta, Georgia. 656 pp. 32 ESRI 2005. ArcGIS Desktop Help. URL: http://webhelp.esri.com/arcgisdesktop/9.1/index.cfm?TopicName=welcome (last date accessed: 7 June 2005). Foody, G. M. 1996. Fuzzy modeling of vegetation from remotely sensed imagery. Ecological Modeling. 85: 3-12. Foody, G. M. 2002. Status of land cover classification accuracy assessment. Remote Sensing of Environment. 80: 185-201. Galv?o, L. S., A. R. Formaggio, and D. A. Tistot. 2005. Discrimination of sugarcane varieties in Southeastern Brazil with EO-1 Hyperion data. Remote Sensing of Environment. 94: 523-534. Gopal, S. and C. Woodcock. 1996. Remote sensing of forest change using artificial neural networks. IEEE Transactions on Geoscience and Remote Sensing. 34(2): 398-404. Halpern, M., W. C. Blackwelder, and J. I. Verter. 1971. Estimation of the multivariate logistic risk function: a comparison of the discriminant function and maximum likelihood approaches. Journal of Chronic Disease. 24: 125-158. Hasegawa, O. and T. Kurita. 2002. Face and non-face classification by multinomial logit models and kernel feature compound vectors. Proceedings of the 9 th International Conference on Neural Information Processing. 2: 996-1000. Homer, C. C. Huang, L. Yang, B. Wylie and M. Coan. 2004. Development of a 2001 National Landcover Database for the United States. Photogrammetric Engineering and Remote Sensing. 70(7): 829-840. 33 Hosmer, T.D. W. Hosmer, and L. L. Fisher, 1983. A comparison of the maximum likelihood and discriminant function estimators of the coefficients of the logistic regression model for mixed continuous and discrete variables. Communications in Statistics. B12: 577-593. Hosmer, D. W. and S. Lemeshow.1989. Applied logistic regression. Wiley-Interscience New York. 307 pp. Hosmer, D. W., S. Lemeshow. 2000. Applied logistic regression, second edition. Wiley- Interscience New York. 375 pp. Hossain, M., S.Wright, and L. A. Petersen. 2002. Comparing performance of multinomial logistic regression and discriminant analysis for monitoring access to care for acute myocardial infarction. Journal of Clinical Epidemiology. 55: 400-406. Jensen, J.R. 1986. Introductory digital image processing: a remote sensing perspective. Prentice-Hall; Englewood Cliffs, N.J. 378 pp. Jensen, J.R. 2000. Remote sensing of the environment: an earth resource perspective. Prentice Hall; 2nd edition, Upper Saddle River, N.J. 544 pp. Johnson, R.A. and D.W. Wichern. 2002. Applied multivariate statistical data analysis. Prentice Hall; 5 th edition, Upper Saddle River, N.J. 767 pp. Keuchel, J., S. Naumann, M. Heiler, and A. Siegmund. 2003. Automatic land cover analysis for Tenerife by supervised classification using remotely sensed data. Remote Sensing of Environment. 86: 530-541. Lillesand, T. M. and R. W. Kiefer. 1994. Remote sensing and image interpretation. Wiley and Sons; 3 rd edition, New York. 750 pp. 34 Mahopatra, K. and S. Kant. 2005. Tropical deforestation: a multinomial logistic model and some country-specific policy prescriptions. Forest Policy and Economics. 7: 1-24. Metternicht, G.I. 2003. Categorical fuzziness: a comparison between crisp and fuzzy class boundary modeling for mapping salt-affected soils using Landsat TM data and a classification based on anion ratios. Ecological Modeling. 168: 371-389. [NASA] National Aeronautics and Space Adiministration. ?Landsat 7 Science Data Users Handbook? [internet online] (accessed 5 January 2005); available from http://ltpwww.gsfc.nasa.gov/IAS/handbook/handbook_toc.html. Neter, J., M. Kutner, C.J. Nachtsheim, and W.Wasserman. 1996. Applied linear statistical models. McGraw-Hill. Boston, Massachusetts. 1408 pp. Pal, M. and P. M. Mather. 2003. An assessment of the effectiveness of decision tree methods for land cover classification. Remote Sensing of Environment. 86: 554- 556. Peng, C. J. and R. N. Nichols. 2003. Using multinomial logistic models to predict adolescent behavioral risk. Journal of Modern Applied Statistical Methods. 2(1): 1-13. Press, S. J. and S. Wilson. 1978. Choosing between logistic regression and discriminant analysis. Journal of the American Statistical Association. 73(364): 699-705. Ram?rez-Garc?a, P., J. L?pez-Blanco, and D. Oca?a. 1998. Mangrove vegetation assessment in the Santiago River Mouth, Mexico, by means of supervised classification using Landsat TM imagery. Forest Ecology and Management. 105: 217-229. 35 SAS/STAT Software: Changes and Enhancements, Release 8.2, Cary, NC: SAS Institute Inc., 2001. SAS/STAT 2005 SAS OnlineDoc version eight URL: http://v8doc.sas.com/sashtml/ (last date accessed: 20 June 2005). Wrigley, N., 1985.Categorical data analysis for geographers and environmental scientists. Longman New York. 385 pp. Yoshida, T. and S. Omatu. 1994. Neural network approach to land cover mapping. IEEE Transactions on Geoscience and Remote Sensing. 32(5): 1103-1109. 36 CHAPTER 2 BRINGING IMAGES TO A COMMON RADIOMETRIC SCALE USING AGGREGATE NO CHANGE REGRESSION: A COMPARISON BETWEEN RADIOMETRIC NORMALIZATION TECHNIQUES ABSTRACT Aggregate no change regression (AG-NCR) is a relative normalization procedure that addresses changes in the reflectance between images by aggregating the digital number values of each image, which have not experienced land use change or temporal features such as clouds, and regressing those values against one another. What separates this normalization procedure from other common relative normalization techniques is that it explicitly addresses spectral differences between images caused by geo-rectification errors. Compared with 2 absolute normalization models, a top ranking relative normalization model, and a combination of absolute and relative normalization model, AG-NCR significantly reduced the variability between 11 Landsat ETM+ scene comparisons, accounting for approximately 96% of the total variation between images. Using this technique, multiple images representing the same season, can be brought to a 37 common radiometric scale for the purpose of generating fine grain regional classifications. INTRODUCTION In recent years, the ability to produce fine grain (30 m) classifications that span across large extents has dramatically increased (Vogelmann et al., 2001; Skole et al., 1997). Many of these regional classifications have been successfully performed by mosaicing adjacent satellite images acquired from sensors such as Landsat 7 (Homer et al., 1997, 2004). Adjacent images, though, often have significantly different spectral reflectance due to changes in the environment and changes in the sensor as a function of image acquisition dates (Hall et al., 1991; Du et al., 2002). If not accounted for, these spectral differences can increase class variability, causing the precision and accuracy of a classification to decrease (Song et al., 2001; Hall et al., 1991). To compensate for spectral differences between images, analysts have employed a number of different techniques ranging from classifying on a scene-by-scene basis (Vogelmann et al., 2001), to normalizing images prior to classification (Pax-Lenney et al., 2001). While the first approach addresses spectral differences between scenes, it does so in an indirect fashion, leading to issues such as edge matching, increased sampling, and multiple classification accuracies (many classification models), making the final map less appealing, more costly, and difficult to assess. Transforming imagery to a common radiometric scale, however, addresses spectral differences directly, potentially removing the issues associated with scene-by-scene classifications. For these reasons, radiometric normalization is often preferred over the time intensive scene-by-scene classification technique (Pax-Lenney et al., 2001; Song et al., 2001; Olthof et al., 2005a, b). Currently 38 though there is no commonly accepted radiometric normalization technique, and existing procedures do not fully address the sources of spectral differences between images. The purpose of this paper is to compare and contrast four commonly used radiometric normalization techniques with our newly developed technique that directly considers geo- rectification errors. Background Radiometric normalization can be separated into two broad categories, absolute and relative. Absolute radiometric normalization attempts to transform spectral reflectance to surface reflectance using empirically derived models that account for atmospheric conditions and changes in the sensor given the temporal nature of the imagery. These models, referred to as radiometric transfer codes (RTC), have been shown to be effective at converting at-sensor reflectance to surface reflectance (Holm et al., 1989; Moran et al., 1992), but typically require accurate atmospheric information at the time of image acquisition, which is seldom available for historical imagery and can be costly to collect for multiple scenes (Chavez, 1996; Du et al., 2002; Velloso and De Souza, 2002; Canty et al., 2004). Therefore, alternative absolute normalization routines have been developed, such as dark object subtraction (Chavez, 1988, 1989, 1996; Song et al., 2001), which assumes the reflectance of certain objects (e.g., water and shadows) are near or at zero and that reflectance values deviating from zero, for these objects, are caused by atmospheric interference. In this methodology, dark objects are often manually selected through image interpretation (Chavez, 1989). In instances where manual selection is too time consuming and/or expensive, automated routines have been developed (Teillet and Fedosejevs, 1995; McDonald et al., 1998; Song et al., 2001). 39 Relative radiometric normalization attempts to transform spectral reflectance of one image (slave image) to match the spectral reflectance of another image (master image). The focus of this approach is not to transform image spectral reflectance to surface reflectance but to bring images to a common radiometric scale. This approach not only accounts for atmospheric and sensor conditions, but also accounts for changes in the biotic and abiotic environment (e.g., differences in soil moisture and growth in vegetation). Relative normalization procedures include algorithms such as histogram equalization and matching (Erdas, 1997), linear regression (Jensen, 1983), multiple regression (Olsson, 1993), non-parametric smoothing (Velloso and De Souza, 2002), and piecewise regression (Erdas, 1997). This approach typically uses features that remain constant through time to capture the spectral relationship, by band, between two images. Using these relationships, relative normalization then transforms the spectral reflectance of the slave image to match the spectral reflectance of the master image. Identifying spectrally constant features in the imagery is integral to the relative normalization approach. Numerous ways have been proposed to achieve this task ranging from manual interpretation (Schott et al., 1988) to automated techniques (Hall et al., 1991). Others have suggested pixels occupying the same geometric location between images be used to transform the digital number (DN) values of one image to match the DN values of another image (Jensen, 1983; Song et al., 2001). In this approach, pixels that have been affected by land use change, given different image acquisition dates, are thought to have a minimal effect on the image transformation process. Yuan and Elvidge (1996), however, have demonstrated that a no change regression (NC) procedure reduces the effects of land use change and/or temporal features, such as clouds between image 40 acquisition dates and improves the normalization process. They also found that pixels that have undergone land use change or that have temporal features could have a significant, negative impact on the normalization process. Given that the objectives of absolute and relative normalization differ, few have directly compared the 2 methodologies. Song et al., (2001) evaluated 4 dark object subtraction (DOS) methods (absolute), 2 dense dark vegetation (DDV) methods (absolute), a path radiance (PARA) method (absolute), and a regression method (relative) in terms of overall classification accuracies and found that 2 DOS methods and the regression technique were the top performing normalization models. However, no statistical tests were performed to determine if overall differences in classification accuracy could be attributed to a particular normalization method (e.g., compare variability in overall accuracies, perform ANOVA, or perform non-parametric ANOVA). Furthermore, errors in the classification may have been an artifact of the classification methodology and or classification rule(s) as opposed to the normalization technique, potentially masking the effect of the normalization process. Alternatively, one could evaluate both approaches by comparing how well each method transforms one image to match the spectral values of another image. Comparing absolute and relative normalization techniques in this manner, Yuan and Elvidge (1996) found that NC regression outperformed Chavez?s (1988) haze correction algorithm (early DOS), Jensen?s simple regression (1983), Hall et al.?s (1991) dark and bright regression, Schott et al.?s (1988) pseudoinvariant features normalization, mean-standard deviation normalization, and a minimum maximum normalization procedure. 41 While the NC regression showed a remarkable improvement in bringing two MSS images (27 years between scene acquisition dates) to a common radiometric scale, it makes assumptions that may limit the utility of the technique. For example, only MSS bands 3 and 4 were used to identify no change pixels. Different bands, however, accentuate different features in satellite imagery (Jensen, 2000). Pixels identified as having constant reflectance or changed reflectance in bands 3 and 4 may not represent pixels that have constant reflectance or changed reflectance in other spectral bands. Furthermore, initial estimates that are determined by locating the centers of land and water ellipsoids from band 3 and 4 scattergrams can potentially be influenced by extreme land and water values. Finally, NC regression assumes the geometric rectification between images is accurate (i.e., the spatial location of a pixel in one scene represents the same spatial location in another scene). Geometric accuracies ranging from + 1 to + 8 pixels for ETM+ imagery, depending on whether the source data have been terrain corrected or systematically corrected, are commonly encountered (NASA, 2005), potentially limiting the validity of this assumption. Given the assumptions of the NC regression, recent improvements in DOS methods (Chavez, 1996; Song et al., 2001), and that the Multi-Resource Land Characteristics Consortium (MRLC) has adopted at-satellite reflectance transformation of Landsat imagery (Huang et al., 2002), we compared and contrasted the ability of five normalization methods to bring images, acquired over different dates with the same representative season, to a common radiometric scale (Table 1). 42 METHODOLOGY Overview This section provides an overview of the study design and methods. Sections following the overview describe different aspects of our analyses in detail. Eleven overlapping areas of 9 Landsat enhanced thematic mapper plus (ETM+) images (Figure 1), acquired over a 2-year period representing a leaf-off winter season, were used to compare and contrast 5 radiometric normalization procedures (Table 1). Each Landsat image was preprocessed by the MLRC using processing level 1T procedures (NASA, 2005). Because extremely bright objects (e.g., sand, rocks, and concrete) have the potential to saturate individual Landsat ETM+ bands in one image without necessarily saturating the same band in another image (given changes in sun elevation angle, solar illumination, and atmospheric condition) we removed pixels that had DN values > to 240 for image pairs from the analysis. In addition, areas that had substantial cloud cover were removed from the analysis. For the remaining pixels in each of the overlapping regions, we developed a non-parametric change detection method and used it to identify no change pixels. Each no change pixel was then randomly assigned into a training (~ 50%) and validation (~ 50%) data set. Spectral values were extracted from each image, by band, on a pixel basis, using the no change subset, and imported into Statistical Analysis Software version 8.3 (SAS). The appropriate transformation coefficients for each normalization method were determined and applied to the validation data set. Master image spectral values for the validation data sets were then compared against transformed slave images spectral values for each method, on a pixel basis, using 43 ordinary least squares regression. Techniques that have brought images to a common radiometric scale in this analysis should predict intercept (?) and slope (?) values approximately equal to 0 and 1, respectively, and should explain a large amount of the variation between spectral values in image overlaps [i.e., have a large coefficient of determination (R 2 )]. To test whether there was a significant improvement among normalization methods, we compared predicted intercept and slope estimates with expected values for each method and performed a non-parametric one-way analysis of variance using R 2 estimates. R 2 estimates were compared by band as opposed to root mean squared error estimates because they are independent of scale, thereby allowing for direct comparison among the different normalization techniques (Olsson, 1993). Identifying No Change Pixels After removing potentially saturated DN values (i.e., > 240 for both images) from the overlap between two images, we subtracted the DN values of the slave image from the DN values of the master image by band on a pixel basis. Using an equal area slicing algorithm we ranked differenced values from smallest to largest and then grouped rankings into 100 classes, each with approximately equal area (ESRI, 2005; ERDAS, 1997). Assuming that 90% of the area in the overlap between scenes had not experienced land use change or temporal features we then selected no change pixels by querying sliced categories > 5 and < 95. Ordinary Least Squares Regression Ordinary least squares regression was used to determine how well the spectral values of each transformed slave image matched the spectral values of each master image. Assuming that transformed slave images have been brought to the same 44 radiometric scale as master images, one would expect a perfect 1 to 1 linear relationship that crosses through the origin (i.e., intercept of 0, slope of 1, and R 2 = 1). Using ordinary least squares regression and the validation data set, we estimated intercept (?), slope (?), and R 2 parameters by minimizing the deviation (Q) between transformed slave image and master image spectral values as follows (Neter et al., 1996): () ? = ??= n i ikkkik xy 1 2 ?? k Q (1) where y and x represent the master and slave spectral values for Landsat ETM+ band k respectively, and ? and ? are solved as follows: ()() () kkkk n i kik n i kikkik k xy xx yyxx ??? ?= ? ?? = ? ? = = 1 2 1 (2) To determine the amount of variation explained between master and transformed slave image we use the following equation: ()[] () ? ? = = ? ?+ = n i kik n i kikkk k yy yx R 1 2 1 2 2 ?? (3) The top performing normalization method was determined by comparing estimated ?, ?, and R 2 values against expected values. Normalization Methods Dark Object Subtraction (DOS) ? DOS assumes that there are objects in the imagery that are completely dark. However, due to atmospheric scattering these objects will not appear absolutely dark (Chavez, 1988). Song et al. (2001) indicated that the top 2 45 DOS models, in terms of classification accuracy, were a method that assumed no atmospheric transmittance loss and no diffuse downward radiation (DOS1) and a method that assumes transmittance loss, diffuse downward radiation, and Rayleigh scattering (DOS3). Given that both methods gave comparable results in their study, we selected DOS1, the simpler of the two methods for our comparison. To transform DN values (a measure of at-satellite radiance) to surface reflectance, using DOS1, we needed to identify the location of dark objects and account for the additive effect of atmospheric interference. Dark objects (DNmin) were located using the lowest DN value within an image, by band, which had > 1000 pixels (Teillet and Fedosejevs, 1995; Song et al., 2001). Assuming a 1 % surface reflectance for these objects (Chavez, 1989; Moran et al., 1992; Song et al., 2001) path radiance (Lp k ) was calculated as follows; ( ){ } 7 5,-1 Bands ETM k ?sinESUN. BDNminGLp kkkk +=?+?= ? 010 (4) where G , B, ESUN, and ? are band specific measures of gain, bias, solar exoatmospheric irradiance, and sun elevation angle, respectively. Using equation 4, surface reflectance (sr k ) was derived using the following equation (after Song et al., 2001): ()?sinESUN )Lp?(Lsat sr kk k ? ? = (5) where Lsat k is at-satellite radiance and is calculated as follows (Markham and Barker; 1986; Huang et al., 2002): kkkk BDNGLsat +?= (6) For this method, the overlapping areas in slave and master scenes were transformed to surface reflectance and then regressed against one another using the validation data set. 46 At-Sensor Reflectance (ASR) ? ASR is a first order normalization routine adopted by the MRLC (Huang et al., 2002). This routine directly addresses issues of sensor change and illumination geometry. ASR is calculated using equation 6 to convert DN values to at-satellite radiance and then converting at-satellite radiance to sensor reflectance as follows (Huang et al., 2002): ()?sinESUN dLsat? ASR k k ? ?? = 2 (7) where d is the earth sun distance in astronomical units and can be calculated according to Iqbal (1983). This approach is similar to DOS1 except that Lp k term is removed and the square of earth sun distance is added to the reflectance formula. In this method, the overlapping areas in slave and master scenes were transformed to sensor reflectance and then regressed against one another using the validation data set. No Change Regression (NCR) ? NCR is a relative normalization procedure that attempts to bring adjacent overlapping images to a common radiometric scale. Unlike the previously described procedures, this technique requires images to share a portion of geographic space. After removing pixels that represent land use change or temporal features, NCR regresses, on a pixel basis, the DN values of a slave image onto a master image using ordinary least squares regression (Yuan and Elvidge, 1996). Estimated intercept and slope coefficients are then used to transform DN values of the slave image to match DN values of the master image. Elvidge et al. (1995) suggest identifying no change pixels by first locating the centers of land and water pixels in scattergrams of Landsat multispectral scanner (MSS) bands 3 and 4 (approximate ETM+ spectral equivalent of band 4) and deriving initial estimates for each linear transformation. Using 47 these estimates and a half-perpendicular width of 10 DN (approximate ETM+ equivalent, 20 DN), they generated a series of Boolean equations separating change verses no change pixels. Extreme water and land values (i.e. land use change or temporal features), however, can skew the center of water and land pixels, adversely affecting the initial intercept and slope estimates of the linear transformation. Furthermore, pixels identified as no change in certain bands may not represent no change pixels in other bands given that each band accentuates different environmental characteristics for both MSS and ETM+ sensors (Jensen, 2000). Therefore, we adopted our non-parametric approach of identifying no change pixels described early in this paper. Using the training data set and equations 1 and 2 we estimated linear transformation coefficients by band for each overlapping area. Applying these estimates to the DN values of the slave image validation data set, we then regressed transformed slaved image DN values against master image DN values, by band, for all image comparisons. At-Sensor No Change Regression (AS-NCR) ? AS-NCR is a merger of the absolute ASR and the relative NCR normalization methods. The primary difference between this method and ASR and NCR is that before applying NCR we transformed DN values to at-sensor reflectance. After converting DN values to at-sensor reflectance, we regressed slave image reflectance onto the master image reflectance using the training data set. Estimated slope and intercept coefficients were then used to transform the slave image?s sensor reflectance values, for the validation data set by band for all 11 image comparisons. Finally, these transformed sensor reflectance were regressed against master sensor reflectance by band for all image comparisons. 48 Aggregate No Change Regression (AG-NCR) ? AG-NCR is a relative normalization procedure that, in addition to accounting for changes in the biotic and abiotic environment, also addresses errors in the geographic rectification process. This technique is identical to the NCR procedure with one exception; prior to extracting the spectral values from the overlapping region of two images, AG-NCR performs a mean aggregating function (ESRI, 2005; ERDAS, 1997). The aggregating function helps to limit the effect of geographic rectification inaccuracies between two images by minimizing the importance of the spatial location of any one individual pixel. The assumption is that factors affecting differences between scene DN values occur at coarser scales than an individual ETM+ pixel (grain size 30 meters), and that geo-rectification errors are such that a pixel occupying the same spatial location in two scenes may not represent the spectral reflectance of the same objects. To determine an appropriate aggregation size, a series of empirical tests were performed, simulating geo-rectification errors between the overlapping areas of two images. DN values (0-255) were randomly assigned to pixels within a master scene overlap, assuming a worst-case scenario between adjoining pixels (i.e., no spatial correlation between spectral reflectance). Slave scene pixel values were then calculated using an arbitrary linear transformation of the master scene. ( ) 875.0 5? = DN DN master slave (8) Next, constant spatial shifts of 1, 2, and 4 pixels were applied to the slave image to mimic geo-rectification errors. Finally, mean aggregates of 0, 25, 100, 400, 625, 2500, 10000, and 40000 pixels were compared, using ordinary least squares regression analysis. From 49 the regression analyses mean estimated intercept, slope, and R 2 were obtained for each aggregation level and were plotted against the size of the aggregating function for each corresponding pixel shift (Figure 2). From these plots a conservative aggregation size of 50 (2500 pixels) was selected based on identifying a point where all estimated values (intercept, slope, and R 2 , respectively) begin to reach a local asymptote (5, 0.875, and 1, respectively), for the three simulated pixel shifts. After aggregating the bands of each image, we randomly assigned 50% of the aggregated pixels to a training and validation data set. From the training data set, we generated intercept and slope coefficients by band for each image comparison, using equations 1 and 2. Using these estimates and the validation data sets, we transformed the DN values of the slave image and regressed those values against the master image DN values by band for each image comparisons. RESULTS Transformation Coefficients Transformation intercepts and slope estimates for NCR and AS-NCR (DN equivalent) were very similar (Table 2). Comparing mean R 2 values for these methods revealed that both techniques explain equal amounts of variation between slave and master image spectral reflectance (Table 2). Given at-sensor reflectance is a linear transformation of DN values, these results were expected. On average, AG-NCR had smaller intercepts than NCR and AS-NCR (DN equivalents) and had larger slopes than NCR (? = 0.05). In addition, AG-NCR had initial R 2 values that were significantly larger than the other two relative methods (Kruskal- Wallace X 2 df = 2 = 35.66; p-value < 0.0001; Table 2). 50 For ASR and DOS1 methods, transformation coefficients were determined by band using image header files and the lowest DN value with at least 1000 pixels within each image (Teillet and Fedosejevs, 1995; Song et al., 2001). Transformation values for each parameter within ASR and DOS1 can be found in Table 3. Validation Coefficients The amount of variation explained between master spectral values and transformed slave values for DOS1, ASR, NCR, and AS-NCR methods were equal for all bands, suggesting that given a linear transformation, these 4 methods would equivalently bring each image to common radiometric scale. The appropriate linear transformation for NCR and AS-NCR would be, on average, an intercept and slope value equal to 0 and 1, respectively, implying that these two methods are at a common radiometric scale (Table 4). DOS1 and ASR, however, had intercept and slope values different than 0 and 1, indicating that an additional additive and multiplicative transformation is necessary to bring each image to a common radiometric scale (Table 4). In addition to having average intercept and slope values approximately equal to 0 and 1, respectively, the AG-NCR method showed a significant improvement (Kruskal- Wallace X 2 df = 4 = 44.41; p-value < 0.0001) in the amount of variation explained between master and transformed slave spectral reflectance (Table 4). On average, AG-NCR explained an additional 6% of the variation (averaged across all ETM+ bands) between master spectral values and transformed slave values. DISCUSSION Comparing transformed slave spectral values with master spectral values, we determined that AG-NCR out-performed the rest of the normalization procedures in 51 terms of bringing each scene to a common radiometric scale. NCR and AS-NCR also brought each Landsat ETM+ image to common radiometric scale but did not explain as much of the variability between each image comparison, thereby suggesting geo- rectification accuracy is an additional source of spectral variability between images. Furthermore, when comparing AS-NCR and NCR it was apparent both methods explained the same amount of variation between master and transformed slave spectral values, indicating the additional step of converting DN to at-sensor reflectance prior to performing regression is unnecessary. The DOS1 and ASR methods did not bring adjacent scenes to a common radiometric scale. This finding implies that these models may not be appropriate, by themselves, to mosaic adjacent images for the purposes of classification. When transformed slave values were regressed against master slave values, DOS1 and ASR explained equal amounts of the variation in the overlap, again suggesting that converting to surface and satellite reflectance is unnecessary. While DOS1 and ASR performed poorly in our tests, it should be noted the objective of these methods are to estimate surface and satellite reflectance, not to bring images to common radiometric scale. Studies interested in quantifying the affects of biological growth between two images acquired at different time periods may be better suited by a normalization approach such as DOS1, ASR, dark object bright object regression (Hall et al., 1991), pseudo invariant features (Schott et al., 1988), or a combination of these techniques with AG-NCR. Prior to comparing normalization methods, we identified and removed pixels that had experienced land use change and/or temporal features using a non-parametric slicing algorithm. This algorithm was developed to address issues of automation, outlier 52 influence, and individual band differences. Using our non-parametric differencing technique, one can conservatively estimate the amount of land use change and/or temporal features within the overlapping areas of two scenes and implement those estimates directly into the model. In our study, we chose to mask the majority of temporal features before running our non-parametric differencing technique. This additional step, however, could have been directly addressed using our non-parametric differencing procedure. For example, given a short time period between ETM+ image acquisition dates (1 year) one might hypothetically estimate that 6% of the overlapping area between two images has experienced land use change. In addition, for the visible and near infrared portions of the electromagnetic spectrum (ETM+ bands 1-4), temporal features, such as clouds, appear to be covering 7% of the slave scene in the overlapping area. Given that the slave scene is subtracted from the master scene in our technique, we can set slicing thresholds as follows: valueslice5,7 bands valueslice valueslice4-1 bands valueslice pixels NC 3 pixels NC 10 97 97 <> <> The assumption in this example is that land use change is equally distributed (3% of the area in the overlap) to upper and lower extreme differenced values, while differenced values associated with temporal features are allocated to only the lower extreme differences in ETM+ bands 1-4. Conservative (i.e., inflated) estimates of land use change and temporal features should be used when applying this approach and special care should be given when allocating thresholds for temporal features given that clouds in a slave scene would manifest as negative differenced values where as clouds in a master scene would manifest as positive differenced values in ETM+ bands 1-4. By inflating 53 land use change and temporal features estimates, analysts can be assured of removing pixels that have experienced extreme spectral change from the regression analysis. Interestingly, the amount of average variation unexplained in the AG-NCR method was very close to the amount of variation left unexplained for an aggregation size of 50 and a one pixel shift in our simulation study (Figure 2). This similarity suggests further improvements to bring images to a common radiometric scale can be achieved through increased geo-rectification accuracies. Alternatively, one could increase the aggregation size, but this would more than likely have a minimal impact on the amount of variance explained given that R 2 began to reach a local asymptote at an aggregate size of 50 in our simulation (Figure 2). In a recent study (Velloso and De Souza, 2002), there has been some question as to the validity of linear regression assumptions (i.e., normally distributed errors and a linear relationship between two images acquired at different dates). Checking regression diagnostics for our study, these assumptions appeared to be valid. Furthermore, almost all of the variation in DN values for each master scene was explained by the DN values of each transformed slave scene (mean R 2 = 0.96). Deviations from a linear relationship and or low R 2 values (< 0.8) may suggest that slave and master images represent different seasons and/or pixels representing land use change or temporal features may still be affecting the regression model. In the former case, our approach would be inappropriate, given that different vegetation types use light differently at various times of the year. In the latter case, non-parametric differencing thresholds could be adjusted to alleviate the problem. 54 Given numerous reasons for variation between two images, we suggest that R 2 values be used to determine if two images have been brought to a common radiometric scale. Conservatively speaking, when using AG-NCR, R 2 values less than 0.90 for any band should be checked to see if the assumption of linearity is being met and values less than 0.80 more than likely suggest that images have not been brought to a common radiometric scale. CONCLUSION Compared with 3 popular radiometric normalization techniques and a combination of at-sensor reflectance and no change regression, Aggregate No Change Regression (AG-NCR) showed a significant improvement in bringing 11 adjacent scenes to a common radiometric scale. This method can easily be automated and provides a quantitative assessment of how well two images have been normalized (R 2 ). The advantages of AG-NCR include: 1) it directly addresses geo-rectification errors, 2) it uses a non-parametric differencing technique to identify changed pixels, 3) and it works on an individual band basis. The intent of AG-NCR is to bring images acquired under similar seasons and solar illumination geometries, to a common radiometric scale. Often, however, subtle changes in illumination and seasonality can go undetected with visual inspection of the imagery alone. Using AG-NCR and estimates of model fit, remote sensing analysts can identify these situations and determine which images have been brought to a common radiometric scale. Properly normalized imagery can then be mosaiced together and used to perform accurate fine grain, regional classifications, in a cost and time efficient manner. 55 AKNOWLEDGEMENTS We would like to acknowledge the following people who have provided valuable comments and insight in this paper; Ralph Meldahl, Barry Grand, Nedret Billor, Melissa Reynolds, Kevin Kleiner, and John Gilbert. Funding for this study was provided by Alabama Gap Analysis Project. Table 1. Method, source, and normalization issues addressed Method Acronym Source Issues Addressed Dark object subtraction 1 DOS1 Song et al., 2001 sensor and atmosphere At-Sensor Reflectance ASR NASA, 2001 sensor No Change Regression NCR after Yuan and Elvidge, 1996 sensor, atmosphere, changes in biotic and abiotic condition At-Sensor No Change Regression AS-NCR NASA, 2001 after Yuan and Elvidge, 1996 sensor, atmosphere, changes in biotic and abiotic condition Aggregate No Change Regression AG-NCR This Paper sensor, atmosphere, changes in biotic and abiotic condition, and geometric rectification 56 Table 2. Initial mean regression estimates for relative normalization procedures. AS-NCR intercept and slope estimates correspond to changes in sensor reflectance. Values inside () indicate digital number (DN) equivalents*. * DN equivalents were estimated by back transforming sensor reflectance to DN values. Average sun elevation and solar distance of 34.96? and 0.9864, respectively, were used to calculate DN intercept and slope values. Normalization Technique Parameter Band 1 Band 2 Band 3 Band 4 Band 5 Band 7 intercept 3.400 2.278 2.309 2.862 2.457 1.544 slope 0.952 0.957 0.945 0.961 0.955 0.957 AG-NCR R 2 0.946 0.960 0.967 0.941 0.968 0.963 intercept 5.892 3.192 3.099 4.654 4.082 2.725 slope 0.907 0.935 0.927 0.931 0.932 0.930 NCR R 2 0.861 0.907 0.910 0.889 0.919 0.911 intercept 0.012 (5.95) 0.006 (3.01) 0.005 (2.74) 0.013 (4.04) 0.010 (3.48) 0.006 (2.23) slope 0.944 0.972 0.965 0.976 0.972 0.970 AS-NCR R 2 0.861 0.907 0.910 0.889 0.919 0.911 57 Table 3. Elevation Distance and digital number values of selected dark objects used in ASR and DOS1 normalization techniques. Path Row Date Elevation Distance Dark Object (Band 1) Dark Object (Band 2) Dark Object (Band 3) Dark Object (Band 4) Dark Object (Band 5) Dark Object (Band 7) 1937 12/20/2001 29.1600 0.9840 41 25 18 11 3 4 1938 12/20/2001 30.3500 0.9840 42 27 19 12 3 4 2036 2/26/2001 40.0700 0.9902 48 29 21 12 2 3 2037 12/27/2001 28.9700 0.9836 41 25 18 10 2 3 2038 1/25/2001 33.1700 0.9846 44 28 21 12 2 3 2039 1/25/2001 34.2800 0.9846 46 28 18 10 2 2 2136 3/5/2001 42.5600 0.9919 53 33 25 14 3 4 2137 2/15/2000 37.5200 0.9878 47 29 21 13 2 3 2138 2/15/2000 38.5700 0.9878 49 32 24 14 3 4 58 Table 4. Summarized validation statistics. Values represent averaged estimates across 11 scene comparisons. AS-NCR and ASR intercept estimates correspond to changes in sensor reflectance (RF). DOS1 intercept estimates corresponds to changes in surface RF. AG-NCR and NCR intercept estimates correspond to changes in digital number (DN) values. Values inside () indicate the DN equivalents* for validation intercept and slope estimates for models that use RF. *DN equivalents were estimated by back transforming RF to DN values. Average sun elevation (34.96?), solar distance (0.9864), and dark objects values for ETM+ bands 1-5 and 7 (46, 28, 21, 12, 2, and 3, respectively) were used to calculate the intercept and slope values for each appropriate normalization model. Normalization Technique Parameter Band 1 Band 2 Band 3 Band 4 Band 5 Band 7 intercept -0.02 0.015 0.181 -0.02 0.066 0.119 slope 1.001 1.000 0.997 1.001 1.000 0.998 AG-NCR R 2 0.947 0.962 0.968 0.942 0.968 0.963 intercept 0.006 0.00 -0.00 0.006 0.00 0.001 slope 1.000 1.000 1.000 1.000 1.000 1.000 NCR R 2 0.861 0.907 0.909 0.888 0.919 0.910 intercept 0.000 (0.01) 0.00 (0.00) 0.00 (0.00) 0.000 (0.01) 0.00 (0.00) 0.000 (0.00) slope 1.000 1.000 1.000 1.000 1.000 1.000 AS-NCR R 2 0.861 0.907 0.909 0.888 0.919 0.910 intercept 0.011 (5.95) 0.006 (3.01) 0.005 (2.74) 0.012 (4.04) 0.009 (3.48) 0.005 (2.23) slope 0.944 0.972 0.965 0.976 0.972 0.970 ASR R 2 0.861 0.907 0.910 0.888 0.919 0.910 intercept 0.006 (5.37) 0.006 (-0.63) 0.005 (-1.31) 0.013 (1.47) 0.011 (0.56) 0.006 (-0.96) slope 0.946 (0.95) 0.973 (0.96) 0.967 (0.96) 0.978 (0.96) 0.974 (0.97) 0.971 (0.96) DOS1 R 2 0.861 0.907 0.910 0.888 0.919 0.910 59 60 2137/2138 2136/2137 2038/2138 2037/2137 2136/2036 2036/2037 2038/2039 2037/2038 1938/2038 1937/2037 1937/1938 Scene Overlaps Used in Comparisons 2039 2138 2038 1938 2037 1937 2137 2036 2136 Figure 1. Scenes (path and row) and overlaps used in the comparison between normalization methods. 61 Figure 2. Scatter plots for simulated, image overlapping regions where a slave image is a linear function of a master image and has a constant geometric rectification error of 1, 2, and 4 pixels. Expected intercept and slope values are 5 and 0.875 respectively. Aggregate Size 2 Pixel Shift 1 Pixel Shift 4 Pixel Shift 0 10 20 30 40 50 60 70 80 90 100 110 120 130 Intercept 0.000 0.100 0.200 0.300 0.400 0.500 0.600 0.700 0.800 0.9000.9 R 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 20 40 60 80 100 120 140 160 180 200 Slope 1.000 62 REFERENCES Canty, M.J., A.A. Nielsen, and M. Schmidt. 2004. Automatic radiometric normalization of multitemporal satellite imagery. Remote Sensing of Environment. 91: 441-451. Chavez, P.S. Jr. 1988. An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data. Remote Sensing of Environment. 24: 459-479. Chavez, P.S. Jr. 1989. Radiometric calibration of Landsat Thematic mapper multispectral images. Photogrammetric Engineering and Remote Sensing. 55(9):1285-1294. Chavez, P.S. Jr. 1996. Image-based atmospheric corrections revisited and improved. Photogrammetric Engineering and Remote Sensing. 66(9): 1025-1036. Du, Y., P.M. Teillet, and J. Cihlar. 2002. Radiometric normalization of multitemporal high-resolution satellite images with quality control for land cover change detection. Remote sensing of Environment. 82: 123-134. Elvidge C.D., D. Yuan, R.D. Weerackoon, and R.S. Lunneta. 1995. Relative radiometric normalization of Landsat Multispectral Scanner (MSS) data using an automatic scattergram-controlled regression. Photogrammetric Engineering and Remote Sensing. 61(10):1255-1260. ERDAS. 1997. ERDAS field guide; fourth edition. ERDAS Inc. Atlanta, Georgia. 656 pp. ESRI. 2005. ArcGIS Desktop Help. URL: http://webhelp.esri.com/arcgisdesktop/9.1/index.cfm?TopicName=welcome (last date accessed: 7 June 2005). 63 Hall, F.G., D.E. Strebel, J.E. Nickeson, and S. J. Goetz. 1991. Radiometric rectification: toward a common radiometric response among multidate, multisensor images. Remote Sensing of Environment. 35: 11-27. Holm, R.G., R.D. Jackson, B. Yuan, M.S. Moran, P.N. Slater, and S.F. Bigger. 1989. Surface reflectance factor retrieval from Thematic Mapper data. Remote Sensing of Environment. 27: 47-57. Homer, C.G., R.D. Ramsey, T.C. Edwards, Jr., and A. Falconer, 1997. Land cover-type modeling using a multi-scene Thematic Mapper mosaic, Photogrammetric Engineering and Remote Sensing. 63: 59?67. Homer, C., C. Huang, L. Yang, B. Wylie and M. Coan. 2004. Development of a 2001 National Landcover Database for the United States. Photogrammetric Engineering and Remote Sensing. 70(7): 829-840. Huang, C., Z. Zhang, L. Yang, B. Wylie, and C. Homer. 2002. MRLC 2001 Image Preprocessing Procedure, USGS White Paper. Iqbal, M. 1983. An introduction to solar radiation. Academic Press; Toronto. 390 pp. Jensen, J.R. (Editor) 1983. Urban/suburban land use analysis. In: Manual of Remote Sensing, 2 nd ed. American Society of Photogrammetry. Vol. 2, pp.1571-1666. Jensen, J.R. 2000. Remote sensing of the environment: an earth resource perspective. Prentice Hall; 2nd edition, Upper Saddle River, N.J. 544 pp. Markham, B.L., and J.L. Barker. 1986. Landsat MSS and TM post-calibration dynamic ranges, exoatmospheric reflectances and at-satellite temperatures. EOSAT Landsat Technical Notes, 1: 3-8. 64 McDonald, A.J., F.M. Gemmell, and P.E. Lewis. 1998. Investigation of the utility of spectral vegetation indices for determining information on coniferous forests. Remote Sensing of Environment. 66: 250-272. Moran, M.S., R.D. Jackson, P.N. Slater, and P.M. Teillet. 1992. Evaluation of simplified procedures for retrieval of land surface reflectance factors from satellite sensor output. Remote Sensing of Environment. 41: 169-184. [NASA] National Aeronautics and Space Administration. Landsat 7 Science data Users Handbook, URL: http://ltpwww.gsfc.nasa.gov/IAS/handbook/handbook_toc.html (last date accessed: 5 January 2005). Neter, J., M. Kutner, C.J. Nachtsheim, and W. Wasserman. 1996. Applied linear statistical models. McGraw-Hill. Boston, Massachusetts. 1408 pp. Olsson, H. 1993. Regression functions for multitemporal relative calibration of Thematic Mapper data over boreal forest. Remote Sensing of Environment. 46: 89-102. Olthof, I., C. Butson, and R. Fraser. 2005a. Signature extension through space for northern landcover classification: a comparison of radiometric correction methods. Remote Sensing of Environment. 95: 290-302. Olthof, I., D. Pouliot, R. Fernandes, and R. Latifovic. 2005b. Landsat-7 ETM+ radiometric normalization comparison for northern mapping applications. Remote Sensing of Environment. 95: 388-398. Pax-Lenney, M., C.E. Woodcock, S.A. Macomber, S. Gopal, and C. Song. 2001. Forest mapping with a generalized classifier and Landsat TM data. Remote Sensing of Environment. 77: 241-250. 65 Schott, J.R., C. Salvaggio, and W.J. Volchok. 1988. Radiometric Scene Normalization using Pseudoinvarinat Features. Remote Sensing of Environment. 26: 1-16. Skole, D.L., C.O. Justice, J.R.G. Townshend, and A.C. Janetos. 1997. A land cover change monitoring program: strategy for an international effort. Mitigation Adaptation Strategies Global Change. 2:157-175. Song, C., C.E. Woodcock, K.C. Seto, M.P. Lenney, and S.A. Macomber. 2001. Classification and change detection using Landsat TM data: when and how to correct atmospheric effects. Remote Sensing of Environment. 75: 230-244. Teillet, P.M., and G. Fedosejevs. 1995. On the dark target approach to atmospheric correction of remotely sensed data. Canadian Journal of Remote Sensing. 21: 373-387. Velloso, M.L., and F.J De Souza. 2002. Non-Parametric smoothing for relative radiometric correction on remotely sensed data. Proceedings of the XV Brazilian Symposium on Computer graphics and Image Processing. IEEE International 2002. Vogelmann, J.E., S.H. Howard, L. Yang, C.R. Larson, B.K. Wylie, and N. Van Driel. 2001. Completion of the 1990s national land cover data set for the conterminous United States from Landsat Thematic Mapper data and ancillary data sources. Photogrammetric Engineering and Remote Sensing. 67, 650-662. Yuan, D., and C.D. Elvidge. 1996. Comparison of relative radiometric normalization techniques.ISPRS Journal of Photogrammetry and Remote Sensing. 51: 117-126 66 CHAPTER 3 CREATING LONGLEAF ECOSYSTEM PROBABILITY DISTRIBUTIONS ACROSS THE SOUTHEAST ABSTRACT Longleaf ecosystems have declined to a mere 5% of their original range since European settlement. These dramatic losses, in what was once the dominant pine ecosystem across the southeastern U.S., are the principle reasons for the listing of many plants and animals as threatened and endangered and have been the driving factor for recent longleaf ecosystem restoration efforts. While studies have documented the regional decline of longleaf ecosystems, they provide little information on fine scale fragmentation patterns and current ecosystem locations. This lack of information often limits the efficacy of restoration efforts. To aid longleaf restoration efforts we developed a series of fine grain (30 m) ecosystem probability distributions using multitemporal Landsat enhanced thematic mapper plus imagery, digital elevation models, field data, ancillary data sets, polytomous logistic regression, and a hierarchical classification scheme. Using our ecosystem probability distributions, resource managers can identify the most probable locations for 67 longleaf ecosystems, locate potential restoration sites, prioritize restoration efforts, and estimate ecosystem area. INTRODUCTION Longleaf ecosystems, some of the most species rich plant ecosystems outside of the tropics (Outcalt and Sheffield, 1996), have experienced increased attention over the last several decades. Recent studies on longleaf ecosystems have included topics ranging from structural characteristics (Noel et al., 1998), ecological classification (Peet and Allard, 1993; Carter et al., 1999), vegetation composition (Hedman et al., 1999), to current trends and condition (Outcalt and Sheffield, 1996), with the latter example becoming a primary focus for numerous organizations due to the ecological value of these ecosystems. Findings, based on the USDA Forest Service?s, Forest Inventory Analysis (FIA, Gillespie, 1996), have documented dramatic area reductions in longleaf ecosystems across the southeastern U.S. (henceforth Southeast) since 1935 and suggest that the ecosystem has been eliminated from northeastern North Carolina and southeastern Virginia (Outcalt and Sheffield, 1996). Historically, longleaf ecosystems ranged from Virginia to Texas (Frost, 1993; Figure 1), covering between 24-37 million hectares across the landscape (Outcalt and Sheffield, 1996; Frost, 1993), making it formally the dominant coniferous ecosystem in the Southeast (Frost, 1993; Peet and Allard, 1993). However, by the early 1930?s, longleaf ecosystems had been decimated to 8 million hectares (Outcalt and Sheffield, 1996) and by the 1980?s to less than 1.6 million hectares across the Southeast (Kelly and Bechtold, 1990) due to over harvesting, fire suppression, and land conversion. 68 Currently, longleaf ecosystems are estimated to occupy 1.2 million hectares across the Southeast, a mere 5% of their pre-European settlement estimate (Outcalt and Sheffield, 1996). This tremendous loss of habitat has had a dramatic impact on numerous plants and animals and is the primary reason that many species have become listed as threatened or endangered (Tuldge, 1999). Since the 1980?s, the amount of longleaf ecosystem area has remained fairly stable on public lands, while private land holdings have shown a declining area trend (Outcalt and Sheffield, 1996). This is especially alarming for states such as Alabama where longleaf ecosystems are primarily found on private lands (Outcalt and Sheffield, 1996). These findings indicate a strong need for the conservation of this critically endangered ecosystem (Noss et al., 1995). While conservation and restoration efforts have begun, they have been limited, in part, by the lack of information depicting the current location of these ecosystems. Long term studies such as the FIA have identified declining longleaf ecosystem trends, but are ill-suited to provide meaningful information at fine spatial scales. Due to the coarse nature of these data sets, organizations have had to take a broad based approach towards longleaf ecosystem management, monitoring, and restoration, often limiting the efficacy of their efforts. To become more effective, these organizations need accurate, fine scale data sets that identify coniferous ecosystem types and depict the current location and distribution of longleaf ecosystems. Many recent studies have focused on gradient analysis and ecosystem classification (Peet, and Allard, 1993; Carter et al., 1999; Weakley et al., 2000; Goebel et al., 2001) to better understand the abiotic and vegetated composition of longleaf ecosystems. While these types of studies have been used to quantify differences in 69 longleaf ecosystems, they typically do not address other coniferous ecosystems, and often fall short in identifying the location of longleaf ecosystems, particularly across large extents. Furthermore, patterns emerging at fine scales may not extend to coarser scales (Perry et al., 2002), thereby causing confusion in regional ecosystem classifications. Finally, these classification techniques are typically very labor intensive, can be challenging to interpret, and require a substantial amount of field work when trying to map fine scale class distributions, making these approaches logistically unattainable and too expensive to apply at regional scales. Alternatively, remote sensing offers a unique opportunity to link fine grain information (i.e., field samples) with large extents (i.e., the Southeast). Remote sensing satellites, such as Landsat 7, continuously quantify portions of the electro-magnetic spectrum reflected and transmitted from the earth?s surface at a grain size of 15, 30, or 60 m (225, 900, and 3600 m 2 , respectively) depending on the wavelength of the electromagnetic radiation, at an extent of 183 by 170 km (Jensen, 2000; NASA, 2005). Using spectral information obtained from remote sensing satellites and field sample locations, analysts have successfully identified land use (Yang et al., 2001), forest types (Jakubauskas, 1996; Schmidt and Skidmore, 2003; Foody and Hill, 1996; Brown de Colstoun et al., 2003; Yu et al., 1999), forest structure (Hall et al., 1995; Chen, 1996; McDonald et al., 1998; Makela et al., 2001), and many other variables across numerous landscapes. Several techniques and algorithms have been developed to analyze remotely sensed data. These techniques range from isoclustering algorithms with post identification to regression analysis, multivariate analysis, classification trees, artificial 70 neural networks, and polytomous logistic regression (PLR). Using PLR, Hogland et al. (in progress; see chapter 1) demonstrated the flexibility and utility of probabilistic classifiers when there was substantial spectral overlap between land cover types. Given the structural similarities between longleaf and other coniferous ecosystems in the Southeast, a probabilistic classification would be well suited to distinguish ecosystem types. To improve classification accuracy for projects with an extent that spans across multiple satellite images (Song et al., 2001; Vogelmann et al., 2001; Du et al., 2002), remote sensing analysts have developed techniques that radiometrically normalize sensor reflectance (Elvidge, 1995; Song et al., 2001; Hogland and MacKenzie, in progress; see chapter 2). Comparing normalization techniques, Hogland and MacKenzie (in progress; see chapter 2) found that an Aggregate No Change Regression (AG-NCR) procedure significantly improved the ability to bring images to a common radiometric scale. Once images have been brought to a common radiometric scale, reflectance in one image can be generalized to reflectance in another image, implying that relationships between vegetated ecosystems and the spectral values in one image can be used to identify the same relationships in other images. While remote sensing offers one of the most economically viable and regionally consistent means to identify the location of ecosystems across the Southeast, few analysts have focused on the critically endangered longleaf ecosystem. To aid in management, monitoring, and restoration of longleaf ecosystems, we generated a series of fine scale (grain size 30 m) forested ecosystem data sets depicting the probability distributions of Hardwood, Mixed Hardwood/Pine, Coastal Plain Longleaf, Mountain Longleaf, Slash, and Loblolly ecosystems across eastern Mississippi, Alabama, western Georgia, and the 71 panhandle of Florida (Figure 1). These data sets can be used with a geographic information system (GIS) to identify ecosystem locations, potential longleaf restoration areas, and corridors between longleaf ecosystems. They can also be used to manage existing longleaf ecosystems, determine the amount of longleaf ecosystem area, and/or be directly incorporated into existing ecosystem models. METHODS Forested ecosystem probability distributions were mapped across portions of the Southeast using the spatial locations of image and field interpreted samples, multitemporal Landsat enhanced thematic mapper plus (ETM+) imagery, ancillary data sets, PLR, and a hierarchal classification scheme. Our study area (Figure 1) is known for its diverse vegetation found on a wide variety of soil types ranging from clayey hills occurring in the Piedmont and Ridge and Valley portions of Alabama to xeric sandhill sites comprised of coarse sandy soils (Peet and Allard, 1993). Longleaf ecosystems across our study area are known to vary from relatively sparse overstory and understory vegetation in the more xeric sites to floristically rich savannas and flatwoods found along mesic sites (Peet and Allard, 1993). This pattern of species scarcity and abundance depending on soil moisture led Peet and Allard (1993) to hypothesize that longleaf composition would vary in an interpretable manner between physiographic regions. From their study, we hypothesized that spectral differences in forested ecosystems can be used to identified longleaf ecosystems. Our primary explanatory variables include Landsat ETM+ imagery (spectral bands 1-5 and 7) acquired from 1999 to 2002 (Figure 2). Landsat ETM+ imagery was chosen specifically because of its availability, cost, and resolution (spectral, spatial, and 72 temporal). Recent studies have found that multitemporal imagery improved classification accuracy (Reese et al., 2002; Homer et al., 2004). Therefore, image dates were selected to approximate three seasons: 1) leaf-on spring growing season, 2) leaf-on late summer/fall season, and 3) leaf-off winter season (Figure 2). In all, 64 ETM+ images, representing 22 Landsat scenes, were used to distinguish among forested ecosystem types. Ancillary information used in this study include FIA findings, digital elevation models (DEMs), and vector layers of known publicly and privately owned protected lands within our study area, tree distribution maps (after Little, 1971), and Level III Omernik ecoregions (Omernik, 1987). Due to the inherent spectral variability of multiple Landsat ETM+ images (Elvidge et al., 1995; Yuan and Elvidge, 1996; Song et al., 2001; Hogland and MacKenzie, in progress; see chapter 2), AG-NCR was used to bring ETM+ scenes of the same season, to a common radiometric scale (Hogland and MacKenzie, in progress; see chapter 2). AG-NCR is a relative normalization procedure that determines the optimal linear transformations needed to bring spectral bands of 2 images to common radiometric scales by regressing the overlapping aggregated digital number (DN) values of a slave and master image. The inputs to this normalization technique include a conservative percent area estimate of land use change and temporal features (e.g., clouds, smoke, burnt areas, and shadows) for the overlapping region of a slave and master image and an appropriate aggregation size. For this study, we estimated land use change and temporal features to account for 20% of the area in the overlapping region and used an aggregation size of 50 (2500 pixels). Large bodies of water and areas that had significant cloud cover were removed from this analysis prior to performing the AG-NCR procedure. Images that 73 could not be brought to a common radiometric scale, identified using measures of model fit (R 2 < 0.85 for any ETM+ bands), were removed from the classification. Once images representing the same season were brought to a common radiometric scale they were then merged together by band to create three large images representing each season. Bands from each season were then used as explanatory variables for our hierarchical classification scheme. Hierarchical Classification Scheme Our hierarchical classification scheme is a 2 level, multi-stage classification that uses the conditional probabilities of a preceding PLR model to constrain the conditional probabilities of subsequent PLR models. Using this approach, we were able to reduce the impact of temporal features in multitemporal imagery and were able to hierarchically organize our classification without losing spatially explicit class information. The benefits of our hierarchical classification scheme include fewer field samples, preserving modeling and classification errors, and the ability to account for confounding temporal features. The first stage in level 1 of our hierarchical classification scheme generated a series of land cover, land use change, and temporal features probability distributions at the resolution of one ETM+ pixel, using extracted spectral values visually interpreted from our normalized multitemporal ETM+ imagery and digital ortho quarter quads. Temporal feature probability distributions occurring in stage 1, which identify the season and type of temporal feature, were then allocated to land cover and land use change probability distributions by restricting the explanatory variables of each consecutive stage to seasonal imagery that did not have confounding temporal feature class types. For 74 example, a probabilistic land cover classification that uses multitemporal ETM+ imagery representing a leaf-on and leaf-off season generates 3 PLR models (i.e., stages; Figure 3). The first PLR model identifies the probability distribution of land cover, land use change, and temporal feature class types using the ETM+ bands from both seasons. The second and third PLR models (i.e., second and third stage) generate a series of land cover and land use change probability distributions using only leaf-on or leaf-off imagery as explanatory variables, respectively. Land cover and land use change probability distributions generated in the second and third PLR models are then multiplied by the temporal feature probability distributions of the first PLR model. The product of these multiplications are then added to the land cover and land use change probability distributions of the first PLR model, thereby constraining the land cover and land use change probability distributions of the second and third PLR model to the temporal features probability distributions of the first PLR model. The first stage in the second level of our hierarchical classification scheme generated a series of forested ecosystem probability distributions using extracted spectral values from ETM+ imagery interpreted from field samples and PLR (Figure 3). Similar to level 1 of our hierarchical classification, temporal features can have confounding effects on forested ecosystem probability distributions. To account for these effects, the same staging scheme as described in level 1 was used to generate a series of ecosystem PLR models. The probability distributions of each forested ecosystem stage were then multiplied by the probability distributions of Deciduous and Evergreen land cover types, for each corresponding level 1 stage. The product of these multiplications were then 75 summed across all ecosystem stages, thereby constraining forested ecosystem probabilities by Deciduous and Evergreen land cover types. In the second level of our hierarchical classification, field samples potentially could occur at the same geographic location as temporal features in the imagery. To account for this scenario, an initial hard classification (sensue Foody and Hill, 1996) of the field sample points was performed using the level 1, stage 1 PLR model and a maximum likelihood allocation rule. Field samples identified as a temporal features were then subset into groups, based on the seasonality of the confounding feature, and used in the appropriate forested ecosystem PLR stage. Image Sampling Scheme Land cover, land use change, and temporal features class types were visually interpreted at the spatial scale of one ETM+ pixel across the study area using the normalized ETM+ imagery and digital ortho quarter quads. Class types include generalized National Land-Cover Database categories (Table 1; after Homer et al., 2004), land use change (e.g., cropland converted to bare ground, clear cuts, and developing areas) and temporal features (e.g., clouds, burnt areas, shadows, and smoke). To maintain consistency across the study area, one image interpreter identified all class samples. Using ArcGIS version 8.3, Environmental Systems Research Institute?s Spatial Analyst extension and our normalized multitemporal ETM+ imagery, DN values occurring at the same spatial location as image interpreted sample plots were extracted on a nearest pixel basis by band using the sample command (ESRI, 2005). In all, 10,941 spectral samples, representing the DN values of land cover, land use change, and temporal features for each 76 band of our normalized multitemporal ETM+ imagery, were collected across our study area. Field Sampling Scheme Due to the presence of large, contiguous, coniferous stands on public lands (Outcalt and Sheffield, 1996) and issues of accessibility, field samples of forested ecosystems were primarily collected in national forests, national wildlife refuge areas, state forests, and military installations. Prior to field visits, potential sample sites were identified using aerial photography, stand data, and other ancillary data sets (e.g., spatial layers identifying forest boundaries, roads, and stream networks). After navigating to prioritized areas, sample locations were selected based on field observations. Locations that appeared to be composed of predominantly one overstory species ( > 75% of the overstory component) were sampled based on the following criteria; 1) sampling plots were at least 60 m from any road, 2) sampling plots represented at least a 3 by 3 pixel area of homogeneous overstory, and 3) sample locations had a minimum distance of 60 m from one another. In all, 1772 field samples were collected across our study area. At each sample location, a 0.04 ha plot (~ half an ETM+ pixel) was used to categorize ecosystem types by quantifying the relative density of tree species > 4 m in height. Plots that had > 75% relative density for a single pine tree species, were categorized as one of the following forested ecosystem; Coastal Plain Longleaf, Mountain Longleaf, Loblolly, or Slash. Plots that had > 75% relative density for hardwood tree species were grouped into one class called Hardwoods. Plots that had a mixture of hardwood and pine species (relative densities of hardwood and pine species < 75% and > 25% for a given group) were categorized as Mixed Hardwood Pine. Plots that 77 had mixture of pine species with relative densities < 75% and > 25% for given groups of pines were removed from our analysis because they were found to occupy the same spectral space as low class probabilities of the two most prevalent pine species within those groups, thus generating redundant information. By setting probabilistic thresholds for each pine dominant ecosystem (> 75% relative density for a given pine species), users can identify the probable locations of mixed pine ecosystems. Dominant pine ecosystems not represented in our classification scheme (e.g., shortleaf, sand, spruce and Virginia pine) accounted for a small portion of our study area and thus were not directly modeled. In addition to categorizing relative densities, geographic coordinates (WGS 84) were collected at each sample location using a Garmin eTrex Legend global positioning systems (GPS) unit. To maintain a close proximity to the true location of our sample plot, GPS coordinates were only recorded when location errors were less than 10 meters. Sample plots were then projected into Albers equal area projection using ArcGIS version 8.3 and were used to extract elevation and spectral data from our DEMs and our normalized multitemporal ETM+ imagery, on a nearest pixel basis, using the same process described in the previous section. PLR Models Prior to extracting spectral information for image and field interpreted samples, we randomly partitioned each data set into 2 groups, a training (~70% of the data) and a validation (~ 30% of the data) data set. Training data were used to develop PLR classifications for each stage of our hierarchical classification scheme by first removing redundant or insignificant explanatory variables via a stepwise procedure (SAS, 2005). Thresholds for variables entering and staying in each PLR model were set at a 78 significance level of 0.15 and 0.10, respectively. Next, using combinations of the remaining explanatory variables, a suite of potential PLR classifications were generated, from which the best fitting, most parsimonious models were selected based on Akaike's Information Criterion (AIC; Akaike, 1973; Table 2). The validation data set was then used to independently assess the predictive capabilities of our classification models using methods described in Hogland et al. (in progress; see chapter 1). To avoid extrapolating our PLR models when generating our probabilistic classification, we restricted each predictive model to the upper and lower values of sampled explanatory variables. In level 1 of our hierarchical classification scheme, all pixels within our study area fell within the spectral range of image interpreted samples. In level 2 of our hierarchical classification scheme, many pixels within our study area fell outside the range of our DN and elevation values sampled in the field. Subsequently, ecosystem probabilities for those pixels were not calculated and class probabilities for those pixels reverted to level 1 Evergreen and Deciduous class probabilities. Because the spatial location of each forested ecosystem type was not known prior to the beginning of our study, we approached each stage of level 2 in our hierarchical classification as a retrospective study (Agresti, 2002). Therefore, forested ecosystems were given an equal sampling weight for each PLR model. To account for forested ecosystem prior probabilities across our study area, we incorporated Bayesian statistics into our analyses. Bayesian statistics were generated in the form of prior abundance estimates, relative to a maximum abundance of 100, for each forested ecosystem using FIA findings (Outcalt and Sheffield, 1996) and the spatial layers of tree species range maps (after 79 Little, 1971), known publicly owned and privately owned protected lands, and Level III Omernik Ecoregions (Omernik, 1987). Prior abundance estimates can be converted to prior probabilities as follows; AbundancePrior AbundancePrior iesProbabilitPrior n 1i i i i ? = = (1) where i identifies each ecosystem type and n is the total number of ecosystems. Figure 4 illustrates the prior abundance estimates used at the intersection and extent of each spatial layer. RESULTS All winter leaf-off and summer/fall leaf-on ETM+ images were brought to a common radiometric scale for each respective season (R 2 > 0.85 for every ETM+ band; Figure 2). On average, the AG-NCR normalization procedure accounted for 96.6 % and 96.8% of the variation between ETM+ images representing a winter and summer/fall season, respectively. Six ETM+ images representing a spring leaf-on season could not be brought to a common radiometric scale (R 2 < 0.85 for at least one ETM+ band; Figure 2), and were subsequently removed from our hierarchical classification. Furthermore, 2 representative images, identified as path/row 18/39 and 19/36 for the spring leaf-on season (Figure 2), could not be obtained and consequently were not used in our hierarchical classification. On average, the AG-NCR normalization procedure accounted for 96.5% of the variation between the remaining 14 spring leaf-on ETM+ images. Models selected for each stage in our hierarchical classification can be found in Table 2. While the best fitting most parsimonious models were achieved using a 80 combination of spectral bands from all seasons, these models could not be generalized to our entire study area because of the problems with 8 spring leaf-on images previously mentioned. Therefore, we limited our final probabilistic classifications to PLR models which used winter and summer/fall normalized ETM+ imagery. Reducing the explanatory variables in our hierarchical classification had a minimal effect on model fit ( 2 R ~ ? ~ 0.01; SAS, 2005) and decreased the total number of land use change and temporal features classes (i.e., land use change and temporal features occurring in the spring imagery), making for a less complex probabilistic classification that could be generalized across our study area. PLR stages Slope estimates for each explanatory variable in our hierarchical classification (Tables 3 and 4) identify changes in the natural log odds ratio (logits) of each class comparison (Agresti, 2002; Hogland et al., in progress; see chapter 1), given a change in that explanatory variable. Odds ratios represent the number of times one class was found versus a baseline class (i.e., a class with which all other classes are compared against in the PLR classification), for a given set of explanatory variables. Alternatively, one can view these class counts as class probabilities (i.e., the counts of one class ?A? divided by the counts of another class ?B? is equivalent to the ratio of {A/{A+B}} / {B/{A+B}}). Therefore, the exponentiation of the slope estimates for each logit provides an intuitive description of the change in probabilities between one class versus the baseline class. This description is properly interpreted as a multiplicative increase/decrease in the ratio of 2 class probabilities for every incremental increase in a given explanatory variable. To determine the change in probabilities between one class and another class other than the 81 baseline class, one can exponentiate the difference between the beta estimates of a given explanatory variable. For example, in stage 1 of level 1 in our hierarchical classification (Table 3) the odds of the Deciduous class (i.e., the probability of the Deciduous class / the probability of the Water class) increase on average by a multiple of 1.59 (e (0.464) ) for every unit increase in ETM+ spectral band 4 during a winter leaf-off season. However, when comparing the Deciduous class with the Evergreen class the odds of the Deciduous class (i.e., the probability of the Deciduous class / the probability of the Evergreen class) decrease on average by a multiple of 0.94 (e (0.464 - 0.526) ) for every unit increase in spectral band 4 during a leaf-off season. Using odds ratios, analysts can quickly determine how changes in the spectral values affect the probability of one class given the probability of another class, they can objectively compare the slope estimate of explanatory variable when values have been standardized (Appendix 1), and they can generate class probability through a series of mathematical manipulations (Appendix 2A). Level 1 of the hierarchical classification scheme ? The top performing model in stage 1, level 1 used normalized spectral bands 1-5 and 7 from both seasons to calculate class probabilities of land cover, land use change, and temporal features (Table 2). Subsequent stages also used ETM+ bands 1-5 and 7 to calculate class probabilities, but were restricted to only one season for each PLR model (Figure 3). Each stage in level 1 generated statistically significant models (p-value < 0.0001) that explained the majority of the information within the data ( 2 R ~ values > 0.96). Slope estimates for each stage of level 1 can be found in table 3. 82 Assessing each stage?s PLR model using class 95% confidence intervals, a maximum likelihood allocation rule, and the validation data set (Hogland et al, in progress; see chapter 1), we found that each PLR model accurately predicted the number of observations occurring for each class (Table 5). This result suggests that our PLR models for level 1 were unbiased and could be applied across our study area. Applying our level 1 PLR models, we effectively removed the confounding temporal features in our multitemporal imagery while maintaining the spatially explicit nature of our PLR model and classification errors (Figure 5). Level 2 of the hierarchical classification scheme ? The top performing model in level 2, stage 1 used ETM+ spectral bands 1-5 from both the winter leaf-off and summer/fall leaf-on season and DEMs to differentiate between forested ecosystems types (Tables 2 and 4). Stage 2 and 3 PLR models used DEMs and ETM+ bands 1-5 and 7 for summer/fall leaf-on and winter leaf-off seasons, respectively, to distinguish forested ecosystems. Each stage of level 2 generated statistically significant models (p-value < 0.0001) that explained large portions of the information within the data ( 2 R ~ values > 0.93). Comparatively, both the Slash and Mountain Longleaf ecosystem classes are spectrally closest to the Coastal Plain Longleaf ecosystem class in the winter leaf-off season (i.e., slopes closest to zero in stage 1; Table 4 and appendix 1). However, these forested ecosystems quickly become distinguishable in band 1 of the summer/fall leaf-on season and low and high elevations, respectively (Table 4; appendix 1). Furthermore, the spectral differences between the Slash and Coastal Plain Longleaf ecosystem 83 dramatically increases in ETM+ bands 2, 3, and 4 representing a summer/fall leaf-on season. Assessing the accuracy of our forested ecosystem classification using the same methodology as in stage one, we were able to predict the number of observations occurring for each PLR iteration (Table 6). Constraining our ecosystem models by the conditional probabilities of stage 1 Evergreen and Deciduous classes, we were able to identify the probability distribution of each forested ecosystem type while maintaining stage 1 and stage 2 modeling and classification errors (Figure 6). DISCUSSION Using the spatial locations of image interpreted and field interpreted samples, multitemporal ETM+ imagery, ancillary data sets, PLR, and our hierarchical classification scheme, we have accurately (i.e., good model fit) depicted the probability distributions of longleaf and other pine ecosystems across a large portion of the Southeast. To account for differences among forested ecosystem prior probabilities, given spectral reflectance and the elevation of each pixel, we estimated prior abundance of each ecosystem, relative to a maximum abundance of 100, on a pixel basis, using recent FIA findings (Outcalt and Sheffield, 1996) and the spatial distribution of our ancillary data sets (Figure 4). While we believe these estimates to be an accurate depiction of the chance of finding each forested ecosystem, given the spectral reflectance and elevation of each pixel, they may be insufficient to account for site-specific cases. In instances where the prior abundance of a given ecosystem is greater or less than our assumed values, absolute ecosystem probabilities will be skewed to either the lower or higher end of the probability distribution, respectively. Relative values (i.e., the magnitude of change in the 84 logits between any two pixels within the same prior abundance estimates), however, remain constant across the spatial distribution of our prior abundance estimates. Therefore, individuals that have further knowledge of prior ecosystem abundance for a site specific location (in relation to the spectral reflectance and elevation of that location) can directly incorporate that information into the probabilistic ecosystem classification (Appendix 2B). Potential biological implications While the goal of this study was to generate a spatially explicit data set that accurately depicts the distribution of longleaf ecosystems across a large portion of the Southeast, many other biologically important questions, pertaining to ecosystem structure and the effects on spectral reflectance can be addressed. For example, given the dependence of longleaf ecosystems on fire (Lear et al., 2005) and that longleaf ecosystems typically have sparse overstories with a diverse understory (Peet and Allard, 1993) when compared to loblolly ecosystems, one would expect Coastal Plain Longleaf ecosystems to reflect and transmit more spectral energy than Loblolly ecosystems in the photosynthetically active portion of the electromagnetic spectrum (ETM+ bands 1 and 3). To test this prediction we compared the probability of finding a Loblolly ecosystem with the probability of finding a Coastal Plain Longleaf ecosystem (odds ratio; Table 4) using exponentiated slope estimates for ETM+ bands 1 and 3, acquired during a winter leaf-off and summer/fall leaf-on season using level 1, stage 1 in our hierarchical classification scheme. Surprisingly, we found a mixed effect when comparing these values. In ETM+ band 1, for both the winter and summer/fall periods, the mean odds of finding a loblolly ecosystem increased by multiples of 1.40 and 1.36, respectively, for every increase of 1 85 DN value. In ETM+ band 3, however, the odds of finding a loblolly ecosystem declined by multiples of 0.73 and 0.62, respectively, for every DN value increase. Interpreting these results, we see that as reflectance increase in the spectral range of ETM+ band 1, coastal plain longleaf ecosystems are less likely to occur when compared with loblolly ecosystems and that as reflectance increase in the spectral range of ETM+ band 3, longleaf ecosystems are more likely to occur when compared with loblolly ecosystems. This could suggest a relative trade off in the way each ecosystem absorbs electromagnetic energy for the purposes of photosynthesis (i.e., longleaf ecosystems prefer higher energy, shorter wavelength, electromagnetic energy, and loblolly ecosystems prefer lower energy, longer wavelength, electromagnetic energy). To determine if these results are due to chance alone, one could incorporate standard error estimates (Appendix 3) with mean slope estimates (Tables 4) to calculate the probability that each slope estimate is equivalent to zero. In this example, the slope estimates for both seasons and bands were significantly different from zero (p-value < 0.0001), meaning the relationship found between the probability of finding a loblolly vs. longleaf ecosystem is very unlikely to be caused by chance alone. Comparing the effect of changes in other spectral bands, such as the near and middle infrared bands (ETM+ bands 4, 5, and 7), on the odds of each forested ecosystem, analysts can begin to characterize other biological aspects of each ecosystem. It is important to remember, however, that our relationships are correlational not causal, and that further tests, designed as an experiment, should be performed to substantiate biological interpretations. 86 Management implications Using our forested ecosystem probability distributions and a GIS, resource managers can generate probability distribution maps for a specific location (Figure 8), set probability thresholds to identify the most probable locations of longleaf ecosystems, identify potential longleaf ecosystem restoration sites, and incorporate ancillary data sets to prioritize restoration locations. By weighting the area of each pixel by the probability of each ecosystem, managers can obtain an estimate of the amount of ecosystem area for a predefined location. For example, if a resource manager wanted to estimate the average amount of area within the boundary of the Blackwater State Forest in Florida that was classified as a coastal plain longleaf ecosystem type, they could use zonal statistics to sum the multiples of pixel area by ecosystem probability (Figure 7). In so doing, managers would estimate there to be approximately 19,500 hectares of Coastal Plain Longleaf ecosystem within Blackwater State Forest. This area estimate, however, would more than likely be inflated given that mixed pine ecosystems were not modeled. To account for mixed pine ecosystems we would suggest setting a minimum probability threshold of 25% when performing area calculations. Furthermore, this estimate represents only the mean estimate of ecosystem area. Using modeling error and the delta method (Agresti, 1990) managers could generate probability confidence intervals for each pixel within the study area, which could be used to generate a range of ecosystem area estimates representing a predefined level of ecosystem area accuracy. CONCLUSION We successfully mapped longleaf ecosystems at a fine spatial resolution (30 m grain), across a large portion of the Southeast. These probabilistic ecosystem 87 classifications provide resource managers with a level of detail that is statistically accurate and precise and flexible enough to begin addressing fine scale longleaf ecosystem restoration questions. In addition, model and classification errors have been maintained in a spatially explicit manner across our study area, thereby allowing other researchers to incorporate our model errors into their work. Future studies could potentially improve upon our results by incorporating ETM+ imagery from a spring leaf- on season, adding a textual component to the analysis, and/or directly incorporating the spatial locations of FIA data. Finally, probabilistic ecosystem classification will be made available through the Alabama Gap Analysis by the end of August 2005. AKNOWLEDGEMENTS We would like to acknowledge the following people who have provided valuable comments and insight in this paper; Melissa Reynolds, Kevin Kleiner, and John Gilbert. Funding for this study was provided by Alabama Gap Analysis Project. 88 Table 1. Cross walk between our level 1 land cover classes and NLCD + classes. Land Cover Classes NLCD Classes Water 11 Open Water 21 Developed, Open Space 22 Developed, Low Intensity 23 Developed, Medium Intensity 24 Developed, High Intensity 31 Barren Land (Rock/Sand/Clay) Urban / Transportation / Bare Ground 32 Unconsolidated Shore 41 Deciduous Forest Evergreen 42 Evergreen Forest 43 Mixed Forest Fore sted * Deciduous 52 Shrub/Scrub 71 Grassland/Herbaceous 81 Pasture/Hay Field 82 Cultivated Crops 90 Woody Wetlands Wet Vegetated Area 95 Emergent Herbaceous Wetlands + National Land-Cover Dataset (NLCD) classes 12, 51, 72, 73, and 74 were not applicable to our study. *Forested classes were split into 2 groups Evergreen and Deciduous based on the dominance (> 50% percent cover) of Evergreen and Deciduous trees. Table 2. PLR model statistics for each level and stage of our hierarchical classification scheme (HCS). W and S/F identify winter and summer/fall seasonality of the imagery. + All explanatory variables listed had values significantly different than 0 (? < 0.0001) * Model selected for each stage of our HCS Explanatory Variables + HCS Level HCS stage Top 3 models W ETM+ Bands S/F ETM+ Bands DEM Degrees of Freedom AIC ?AIC Model Weight X 2 Nested Comparison p-value 1* 1-5 and 7 1-5 and 7 no 168 7443.25 0.000 1.0000 - - 2 1, 3-5, and 7 1-5 and 7 no 154 7496.27 53.020 0.0000 1 vs.2 0.0000 1 3 1-5, and 7 1-5 no 154 7551.562 108.312 0.0000 1 vs.3 0.0000 1* - 1-5 and 7 no 36 3952.42 0.000 0.9994 - - 2 - 1-5 no 30 4180.83 228.410 0.0000 1 vs.2 0.0000 2 3 - 2-7 no 30 3967.41 14.990 0.0006 1 vs.3 0.0001 1* 1-5 and 7 - no 36 4839.211 0.000 1.000 - - 2 2-7 - no 30 5045.916 206.705 0.000 1 vs.2 0.0000 1 3 3 1-5 - no 30 5188.836 349.625 0.000 1 vs.3 0.0000 1 1-5 1-5 and 7 yes 60 3271.926 0.000 0.457 - - 2* 1-5 1-5 yes 55 3272.197 0.271 0.399 1 vs.2 0.0679 1 3 1 and 3-5 1-5 and 7 yes 55 3274.232 2.306 0.144 1 vs.3 0.0308 1* - 1-5 and 7 yes 35 4027.868 0.000 0.8921 - - 2 - 1-5 yes 30 4032.092 4.224 0.1079 1 vs.2 0.0142 2 3 - 1, 3-5, and 7 yes 30 4194.528 166.660 0.0000 1 vs.3 0.0000 1* 1-5 and 7 - yes 35 4104.671 0.000 0.8242 - - 2 1, 3-5, and 7 - yes 30 4107.761 3.090 0.1758 1 vs.2 0.0225 2 3 3 1-2, 4-5, and 7 - yes 30 4219.027 114.356 0.0000 1 vs.3 0.0000 89 Table 3. Slope and intercept estimates for each stage of level 1 in our hierarchical classification scheme (HCS). The land cover class type ?Water? was used as the baseline category. W and S/F identify winter leaf-off and summer/fall seasons, respectively. HCS stage Logits intercept W Band 1 W Band 2 W Band 3 W Band 4 W Band 5 W Band 7 S/F band 1 S/F band 2 S/F band 3 S/F band 4 S/F band 5 S/F band 7 LUC W, S/F -3.527 -0.493 -0.040 -0.170 0.278 -0.118 0.218 0.341 -0.202 -0.059 0.261 0.117 0.134 LUC W -5.503 -0.078 -0.315 0.133 0.292 -0.013 -0.010 0.204 -0.287 -0.099 0.306 -0.154 0.322 LUC S/F 0.356 -0.403 -0.069 -0.408 0.334 0.027 -0.110 0.477 -0.319 -0.098 0.252 0.106 0.224 Burnt W -9.352 -0.198 -0.219 -0.391 0.423 -0.410 0.536 0.423 -0.086 -0.184 0.331 0.098 0.068 Burnt S/F -5.316 -0.186 0.003 -0.405 0.434 -0.094 0.056 0.390 -0.126 -0.196 0.143 -0.024 0.253 Smoke W -33.550 0.892 -0.170 -0.526 0.545 -0.367 0.374 -0.279 -0.280 -0.348 0.338 0.077 0.325 Smoke S/F -21.892 -0.704 0.028 -0.213 0.474 -0.145 0.226 0.946 -0.083 -0.520 0.319 -0.005 0.240 Clouds S/F -29.553 -0.430 0.197 -0.248 0.385 -0.125 0.090 0.627 -0.060 -0.205 0.418 -0.057 0.230 Shadow S/F 13.476 -0.941 0.588 -0.669 0.834 0.138 0.130 0.755 -1.041 -0.096 -0.139 -0.348 -0.251 Fields -22.224 -0.432 0.245 -0.303 0.506 -0.292 0.333 0.249 -0.266 -0.065 0.467 0.099 0.175 Urban / Trans. / Bare -29.764 -0.194 0.146 -0.183 0.440 -0.294 0.268 0.419 -0.162 -0.009 0.388 -0.159 0.293 Deciduous -4.725 -0.162 -0.183 -0.209 0.464 -0.135 0.104 0.298 -0.233 -0.491 0.351 0.213 0.024 Evergreen 5.317 -0.354 0.143 -0.160 0.526 -0.186 -0.075 0.066 -0.168 -0.208 0.327 0.009 0.138 1 Wet Vegetated -21.056 0.005 -0.026 -0.082 0.236 -0.090 -0.038 0.273 0.195 -0.262 0.210 0.165 -0.101 LUC S/F -18.146 - - - - - - 0.291 -0.670 0.136 0.293 0.226 0.152 Fields -26.316 - - - - - - 0.192 -0.580 0.137 0.526 0.126 0.256 Urban / Trans. / Bare -30.296 - - - - - - 0.345 -0.268 0.126 0.400 -0.041 0.249 Deciduous -5.990 - - - - - - 0.240 -0.530 -0.214 0.338 0.388 -0.262 Evergreen -4.310 - - - - - - 0.187 -0.316 -0.058 0.338 0.068 -0.107 2 Wet Vegetated -23.247 - - - - - - 0.299 0.100 -0.154 0.149 0.336 -0.294 LUC W -0.435 0.073 -0.654 0.100 0.299 0.107 0.081 - - - - - - Fields -15.736 -0.085 -0.110 -0.169 0.555 -0.042 0.253 - - - - - - Urban / Trans. / Bare -20.152 0.085 -0.008 -0.083 0.490 -0.219 0.289 - - - - - - Deciduous 2.170 0.205 -0.800 -0.177 0.466 0.160 -0.001 - - - - - - Evergreen 14.310 -0.080 -0.639 -0.036 0.506 0.110 -0.249 - - - - - - 3 Wet Vegetated -12.737 0.400 -0.399 -0.073 0.266 0.161 -0.188 - - - - - - 90 Table 4. Slope and intercept estimates for each stage of level 2 in our hierarchical classification (HCS). The forested ecosystem Coastal Plain Longleaf was used as the baseline class. W and S/F identify winter leaf-off and summer/fall seasons, respectively. HCS stage Logits intercept W band 1 W band 2 W band 3 W band 4 W band 5 W band 7 S/F band 1 S/F band 2 S/F band 3 S/F band 4 S/F band 5 S/F band 7 Elevation Slash 0.923 -0.060 -0.051 -0.033 -0.043 0.020 - -0.160 0.374 -0.163 0.150 -0.026 - -0.009 Hardwood 1.626 0.225 0.320 -0.297 -0.254 0.250 - 0.332 -0.741 -0.617 0.206 0.106 - -0.002 Mixed Hardwood/Pine -6.648 0.267 0.206 -0.327 -0.146 0.236 - 0.115 -0.149 -0.305 0.057 -0.044 - 0.012 Mountain Longleaf 1.969 0.149 0.083 0.107 -0.021 -0.027 - -0.349 -0.103 -0.007 0.066 0.006 - 0.025 1 Loblolly -7.415 0.340 0.182 -0.312 -0.075 0.069 - 0.304 -0.173 -0.485 -0.162 0.106 - 0.024 Slash -0.508 - - - - - - -0.212 0.321 -0.151 0.105 -0.034 0.032 -0.012 Hardwood -4.733 - - - - - - 0.642 -0.875 -0.820 0.188 0.435 -0.169 0.012 Mixed Hardwood/Pine -10.663 - - - - - - 0.420 -0.186 -0.508 0.039 0.155 -0.009 0.019 Mountain Longleaf -0.650 - - - - - - -0.168 0.013 -0.013 0.055 0.087 -0.146 0.034 2 Loblolly -10.427 - - - - - - 0.580 -0.067 -0.555 -0.198 0.209 -0.177 0.030 Slash -2.315 -0.008 0.047 -0.041 0.055 -0.091 0.132 - - - - - - -0.010 Hardwood -16.141 0.459 0.205 -0.487 -0.253 0.346 -0.135 - - - - - - 0.012 Mixed Hardwood/Pine -13.610 0.361 0.197 -0.470 -0.127 0.235 -0.091 - - - - - - 0.015 Mountain Longleaf -13.584 0.102 0.004 0.090 0.025 0.016 -0.111 - - - - - - 0.031 3 Loblolly -12.831 0.423 0.211 -0.377 -0.150 0.178 -0.254 - - - - - - 0.024 91 Table 5. Predicted 95% class count confidence intervals (CI) vs. observed class counts for level 1 validation data set. * Land use change + Temporal features Stage 1 Stage 2 Stage 3 Land Cover, LUC * , and TF + classes Lower 95% CI Observed Upper 95% CI Lower 95% CI Observed Upper 95% CI Lower 95% CI Observed Upper 95% CI LUC W, S/F 86.9 107 147.4 - - - - - - Burnt W 90.1 118 151.4 - - - - - - LUC W 71.8 98 129.6 - - - 179.0 205 247.5 Smoke W 15.7 19 24.8 - - - - - - Burnt S/F 95.4 143 189.1 - - - - - - LUC S/F 89.7 115 131.4 202.7 222 254.1 - - - Clouds S/F 108.9 149 179.9 - - - - - - Shadow S/F 69.0 69 71.7 - - - - - - Smoke S/F 83.4 108 147.8 - - - - - - Fields 237.5 286 322.0 245.7 286 308.9 252.0 286 322.3 Urban / Trans. / Bare 228.9 246 287.2 228.2 246 266.9 234.2 246 283.0 Deciduous 272.9 315 389.7 294.2 315 372.6 270.0 315 368.4 Evergreen 292.8 356 379.0 307.8 356 380.9 308.5 356 368.5 Water 340.0 356 363.5 347.3 356 361.6 341.5 356 372.5 Wet Vegetated 129.7 182 219.1 144.4 182 211.8 135.5 182 209.5 92 Table 6. Predicted 95% confidence intervals (CI) for class counts vs. weighted observed class counts for level 2 validation data sets. The observed number of samples (actual) were weighted to remove the effect of prior probabilities. Stage 1 Stage 2 Stage 3 Forested Ecosystem classes Lower 95% CI Observed (actual) Upper 95% CI Lower 95% CI Observed (actual) Upper 95% CI Lower 95% CI Observed (actual) Upper 95% CI Slash 80.7 118 (55) 150.1 98.5 129 (58) 156.3 93.3 127 (58) 152.6 Hardwoods 95.9 123 (83) 142.7 109.2 136 (85) 153.0 105.0 129 (83) 159.9 Mixed Hardwood/Pine 81.1 111(72) 173.7 98.8 123 (78) 177.4 95.7 117 (74) 175.2 Mountain Longleaf 63.0 79 (19) 105.7 79.4 86 (19) 113.7 74.5 84 (19) 109.6 Coastal Plain Longleaf 78.3 126 (126) 154.5 94.9 143 (143) 159.5 84.5 136 (136) 148.3 Loblolly 79.9 120 (51) 150.8 101.7 136 (61) 166.7 88.9 128 (52) 155.3 93 Figure 1. An overlay of the historical range of longleaf ecosystems (after Little, 1971) and our study area. The extent of our study area was defined as the counties within the state of Alabama and the counties that intersected United States Fish and Wildlife Service econumber 29 and 30 (USFWS, 2000). Study Area Longleaf Range Map 94 Figure 2. An overlay of our study area (outlined in gray), and Landsat scenes by path row (PR; outlined in red) showing the dates of ETM+ image acquisition for a winter leaf-off (W), a summer/fall leaf-on (F), and a spring leaf-on (S) season. Dates followed by an M or X indicate master images and images that could not be brought to a common radiometric scale, respectively. N/A identifies scenes for which we were unable to obtain a suitable image representing a spring leaf-on season. PR 1838 W 12/24/1999 F 10/23/2000 S 04/30/2000 PR 1839 W 02/02/2003 F 10/23/2000 S N/A X PR 2336 W 02/18/2002 F 10/29/2001 S 05/14/2001 PR 2238 W 02/27/2002 F 11/07/2001 S 05/15/2001 PR 1936 W 03/07/2001 F 09/28/2000 S N/A X PR 2237 W 02/06/2000 F 10/03/2000 S 05/15/2001 PR 1937 W 12/20/2001 F 09/28/2000 S 04/05/2000 PR 1938 W 12/20/2001 F 10/01/2001 S 04/05/2000 PR 1939 W 12/20/2001 F 11/29/1999 S 04/05/2000 PR 2236 W 01/26/2002 F 10/03/2000 S 04/29/2001 PR 2039 W 01/25/2001 F 11/09/2001 S 04/18/2002 PR 2038 W 01/25/2001 F 10/08/2001 S 04/18/2002 PR 2036 W 02/26/2001 F 11/09/2001 S 06/18/2001 PR 2338 W 02/21/2000 F 10/24/1999 S 05/14/2001 X PR 2337 W 12/27/1999 F 10/24/1999 S 06/18/2002 X PR 1837 W 01/27/2001 F 10/26/2001 S 04/30/2000 X PR 2239 W 01/21/2000 F 11/07/2001 S 04/26/2000 X PR 2037 W 12/27/2001 F 10/08/2001 S 05/14/2000 X PR 2136 W 03/05/2001 F 10/15/2001 S 06/25/2001 X PR 2139 W 01/02/2002 F 10/15/2001 M S 05/24/2001 PR 2138 W 02/15/2000 M F 10/15/2001 M S 04/19/2000 M PR 2137 W 02/15/2000 M F 09/29/2001 S 04/19/2000 M 95 Figure 3. Hierarchical classification scheme for two seasons. W and S/F identify land use change (LUC) and temporal feature (TF) classes found in the normalized ETM+ imagery representing a winter leaf-off and summer/fall leaf-on season, respectively. 96 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 Water Wet Vegetated Areas Evergreen Deciduous Urban / Trans. / Bare Ground Fields Shadow S/F Clouds S/F Smoke S/F Burnt Areas S/F Smoke W Burnt Areas W LUC W LUC S/F LUC W, S/F Tem p o r al F e at u r e s 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 10 11 12 13 14 15310 11 12 13 14 15 2 Stage 1; Both W and S/F ETM+ imagery Stage 2; S/F ETM+ imagery Stage 3; W ETM+ imagery 1 2 3 Level 1 Probabilistic Land Cover and LUC Classification Stage 1 Stage 2 2 10 11 12 13 14 15 10 11 12 13 14 15 10 11 12 13 14 153Stage 3 1 2 3 10 11 12 13 14 15 + + Coastal Plain Longleaf21 Mountain Longleaf20 Loblolly19 Mixed Hardwood/Pine18 Hardwoods17 Slash16 16 17 18 19 20 21Stage 2; S/F ETM+ imagery 16 17 18 19 20 21Stage 3; W ETM+ imagery 16 17 18 19 20 21 Stage 1; Both W and S/F ETM+ imagery + + 16 17 18 19 20 21 Level 2 Probabilistic Ecosystem Classification Level 1 Level 2 Figure 4. Three maps illustrating the spatial configuration of forested ecosystem prior abundance estimates. Maps A, B, and C correspond to the prior abundance estimates, relative to a maximum abundance estimate of 100, for Mountain Longleaf, Coastal Plain Longleaf, and Slash ecosystems, respectively. Prior abundance estimates for Hardwood, Mixed Hardwood/Pine, and Loblolly ecosystems were set at 100 across our study area. Prior Abundance Estimates A B C 1 10 25 50 100 97 Figure 5. Level 1 Land cover and land use change probabilistic classifications. Land Cover and LUC Class Probabilities 0 - 10 10 - 20 20 - 30 30 - 40 40 - 50 50 - 60 60 - 70 70 - 80 80 - 90 90 - 100 Winter and Fall LUCWinter LUCFall LUC Urban / Transportation / Bare Ground Fields Deciduous Evergreen Wet Vegetated AreasWater 98 Figure 6. Level 2 Forested ecosystem probabilistic classifications. Evergreen Deciduous HardwoodMixedLoblolly SlashMountain Longleaf Coastal Plain Longleaf Forested Ecoystem Probabilities 0 - 10 10 - 20 20 - 30 30 - 40 40 - 50 50 - 60 60 - 70 70 - 80 80 - 90 90 - 100 99 Figure 7. Average area (m 2 ) per pixel for Coastal Plain Longleaf Ecosystems within Blackwater State Forest. Blackwater State Forest is located in the northwest portion of the Florida panhandle and is depicted as the red area in the location map. Coastal Plain Longleaf Area (sq m) 0 - 100 100 - 200 200 - 300 300 - 400 400 - 500 500 - 600 600 - 700 700 - 800 800 - 900 Location Map 100 Figure 8. Probability distribution of the Coastal Plain Longleaf ecosystem within and around Conecuh National Forest located in southern Alabama (identified as the red area within the location map). 0255012.5 Kilometers Coastal Plain Longleaf Probabilities 0 - 10 10 - 20 20 - 30 30 - 40 40 - 50 50 - 60 60 - 70 70 - 80 80 - 90 90 - 100 Location Map 101 102 REFERENCES Agresti, A. 1990. Categorical data analysis. Wiley-Interscience. New York. 558 pp. Agresti, A. 2002. Categorical data analysis. Wiley-Interscience. Hoboken, New Jersey. 710 pp. Akaike H. 1973. Information theory and an extension of the maximum likelihood principle In: B. Petrov and F. Cazakil editors, Proceedings of the 2 nd international Symposium on information Theory, Aakademiai Kidao Budapest. Brown de Colstoun, E.C., M.H. Story, C. Thompson, K. Commisso, T.G. Smith, J.R. Irons. 2003. National park vegetation mapping using multitemporal landsat 7 data and a decision tree classifier. Remote sensing of the Environment. 85: 316-327. Carter, R. E., M. D. MacKenzie, and D. H. Gjerstad. 1999. Ecological land classification in the Southern Loam hills of south Alabama. Forest Ecology and Management 114: 395-404. Chen, J.M. 1996. Canopy architecture and remote sensing of the fraction of photosynthetically active radiation absorbed by boreal conifer forests. Ieee transactions of geoscience and remote sensing. 34(6): 1353-1368. Du, Y., P. M. Teillet, and J. Cihlar. 2002. Radiometric normalization of multitemporal high-resolution satellite images with quality control for land cover change detection. Remote sensing of Environment. 82: 123-134. Elvidge C.D., D. Yuan, R.D. Weerackoon, and R.S. Lunneta. 1995. Relative radiometric normalization of Landsat Multispectral Scanner (MSS) data using an automatic scattergram-controlled regression. Photogrammetric Engineering and Remote Sensing. 61(10):1255-1260. 103 ESRI 2005. ArcGIS Desktop Help. URL: http://webhelp.esri.com/arcgisdesktop/9.1/index.cfm?TopicName=welcome (last date accessed: 7 June 2005). Foody, G.M., and R.A. Hill. 1996. Classification of tropical forest classes from Landsat TM data. International journal of remote sensing. 17(12): 2353-2367. Frost, C.C. 1993. Four centuries of changing landscape patterns in the longleaf pine ecosystem. In: The longleaf pine ecosystem: ecology, restoration and management. Proc. Tall Timbers Fire Ecol. Conf. No.18. 17-43. Gillespie, A.J.R. 1999. Rationale for a national annual forest inventory program. Journal of Forestry. 97(12):16-20. Goebel, P.C., B.J. Palik, K. Kirkman, and M.B. Drew. 2001. Forest ecosystems of a lower gulf coastal plain landscape: multifactor classification and analysis. Journal of the Torrey Botanical Society. 128(1): 47-75. Hall, F.G., Y.E. Shimabukuro, and K.F. Huemmrich. 1995. Remote sensing of forest biophysical structure using mixture decomposition and geometric reflectance models. Ecological applications. 5(4): 993-1013. Hedman, C.W., S.L. Grace, and S.E. King. 1999. Vegetation composition and structure of southern coastal plain pine forests: an ecological comparison. Forest Ecology and Management. 134: 233-247. Hogland, J.S., N. Billor, and M.D. MacKenzie. in progress. Comparing polytomous logistic regression and discriminant analysis: a remote sensing perspective. 104 Hogland, J.S., and M.D. MacKenzie. in progress. Bringing images to a common radiometric scale using Aggregate No Change Regression: a comparison between radiometric normalization techniques. Homer, C., C. Huang, L. Yang, B. Wylie and M. Coan. 2004. Development of a 2001 National Landcover Database for the United States. Photogrammetric Engineering and Remote Sensing. 70(7) 829-840. Jakubauskas, M.E. 1996. Thematic mapper characterization of lodgepole pine seral stages in Yellowstone National Park, USA. Remote Sensing of Environment 56: 118- 132. Jensen, J.R. 2000. Remote sensing of the environment: an earth resource perspective. Prentice Hall; 2nd edition, Upper Saddle River, NJ. 181-541. Kelly, J.F., and W.A. Bechtold. 1990. The longleaf pine resource. In: Proceedings of symposium on the management of longleaf pine, 1989. April 4-6; Long Beach, MS. Department of Agriculture, Forest Service, Southern Forest Experiment Station, 11-22. Lear, D. H., W.D. Carroll, P.R. Kapeluck, and R. Johnson. 2005. History and restoration of the longleaf pine-grassland ecosystem: implications for species at risk. Forest Ecology and Management. 211: 150-165. Little, E.L., Jr. 1971. Atlas of United States trees. Volume 1, Conifers and Important Hardwoods: U.S. Department of Agriculture Miscellaneous Publication 1146, 9 pp., 200 maps. 105 Makela, H., and A. Pekkarinen. 2001. Estimation of timber volume at the sample plot level by means of image segmentation and Landsat TM imagery. Remote Sensing of Environment. 77: 66-75. McDonald, A.J., F.M. Gemmell, and P.E. Lewis. 1998. Investigation of the utility of spectral vegetation indicies for determining information on coniferous forests. Remote Sensing of Environment. 66: 250-272. [NASA] National Aeronautics and Space Administration. Landsat 7 Science data Users Handbook, URL: http://ltpwww.gsfc.nasa.gov/IAS/handbook/handbook_toc.html (last date accessed: 5 January 2005). Noel, J.M., W.J. Platt, and E.B. Moser. 1998. Structural characteristics of old-and second-growth stands of longleaf pine (Pinus palustris) in the Gulf Coastal region of the U.S.A. Conservation Biology. 12(3): 533-548 Noss, R.F., E.T. LaRoe, and J.M. Scott. 1995. Endangered ecosystems of the United States: A preliminary assessment of loss and degradation. Biological Report 28. National Biological Service, Washington, D.C. Omernik, J.M. 1987. Ecoregions of the Conterminous United States. Annals of the Association of American Geographers. 77(1): 118-125. Outcalt, K.W., and R. M. Sheffield. 1996. The longleaf pine forest: trends and current conditions. Resource Bulletin SRS-9 Asheville, NC; U.S. Department of Agriculture, Forest Service, Southern Research Station. 23 pp. 106 Peet, R. K., and D. J. Allard. 1993. Longleaf Pine Vegetation of the Southern Atlantic and Eastern Gulf Coast Regions: A Preliminary Classification. In: The longleaf pine ecosystem: ecology, restoration and management. Proc. Tall Timbers Fire Ecol. Conf. No. 18. Perry, J.N., A.M. Liebhold, M.S. Rosenberg, J. Dungan, M. Miriti, A. Jakomulska, and S. Citron-Pousty. 2002. Illustrations and guidelines for selecting statistical methods for quantifying spatial pattern in ecological data. Ecography. 25: 578-600. Reese, H.M., T.M. Lillesand, D.E. Nagel, J.S. Stewart, R.A. Goldmann, T.E. Simmons, J.W. Chipman, and P.A. Tessar. 2002. Statewide land cover derived from multiseasonal Landsat TM data a retrospective of the WISCLAND project. Remote Sensing of Environment. 82: 224-237. SAS/STAT 2005 SAS OnlineDoc version eight URL: http://v8doc.sas.com/sashtml/ (last date accessed: 20 June 2005). Schmidt, K.S. and A.K. Skidmore. 2003. Spectral discrimination of vegetation types in a coastal wetland. Remote Sensing of Environment. 85: 92-108. Song, C., C.E. Woodcock, K.C. Seto, M.P. Lenney, and S.A. Macomber. 2001. Classification and change detection using landsat tm data: when and how to correct atmospheric effects? Remote Sensing of Environment. 75: 230-244. Tuldge, C. 1999. Plant (protecting endangered species). Index on Censorship. 28(6): 164- 166. U.S. Fish and Wildlife Service (USFWS). 2000. U.S. Fish and Wildlife Service region 4 ecosystem. URL: http://www.fws.gov/stand/standards/allshape.html (last date accessed: 25 July 2003). 107 Vogelmann, J.E., D. Helder, R. Morfitt, M. Choate, J.W. Merchant, and H. Bulley. 2001 Effects of landsat 5 thematic mapper and landsat 7 enhanced thematic mapper plus radiometric and geometric calibrations and corrections on landscape characterization. Remote Sensing of Environment. 78: 55-70. Weakley, A.S., K.D. Patterson, S. Landaal, M. Pyne, R.E. Evans, J.A. Teague, M.J. Russo, and others (compilers). 2000. International classification of ecological communities: terrestrial vegetation of the southeastern united states. Alabama subset for GAP Analysis Program report from biological conservation data systems and working draft of April 2000. Association for biodiversity information/The Nature Conservancy, southern resource office, community ecology group. Durham, North Carolina. Yang, L., S.V. Stehman, J.H. Smith, and J.D. Wickham. 2001. Thematic Accuracy of MRLC Land Cover for the Eastern United States. Remote Sensing of Environment, 76:418-422. Yu, B., I.M. Ostland, P. Gong, and R. Pu. 1999. Penalized discriminant analysis of in situ hyperspectral Data for conifer species recognition. Ieee Transactions on Geoscience and Remote Sensing. 37(5): 2569-2577. Yuan, D., and C.D. Elvidge. 1996. Comparison of relative radiometric normalization techniques. ISPRS Journal of Photogrammetry and Remote Sensing. 51: 117-126 108 CUMULATIVE BIBLIOGRAPHY Agresti, A. 1990. Categorical data analysis. Wiley-Interscience. New York. 558 pp. Agresti, A. 2002. Categorical data analysis. Wiley-Interscience. Hoboken, New Jersey. 710 pp. Akaike H. 1973. Information theory and an extension of the maximum likelihood principle In: B. Petrov and F. Cazakil editors, Proceedings of the 2 nd international Symposium on information Theory, Aakademiai Kidao Budapest. Bailey, N., T. Clemets, J. T. Lee, and S. Thompson, 2003. Modelling soil series data to facilitate targeted habitat restoration: a polytomous logistic regression approach. Journal of Environmental Management. 67: 395-407. Benz, U.C., P. Hofmann, G. Willhauck, I. Lingenfelder, and M. Heynen. 2004. Multi- resolution, object-oriented fuzzy analysis of remote sensing data for GIS-ready information. Journal of Photogrammetry and Remote Sensing. 58: 239-258. Brown de Colstoun, E. C., M. H. Story, C. Thompson, K. Commisso, T. G. Smith, J. R. Irons. 2003. National park vegetation mapping using multitemporal landsat 7 data and a decision tree classifier. Remote Sensing of Environment. 85: 316-327. Bull, S. B. and A. Donner. 1987. The efficiency of multinomial logistic regression compared with multiple group discriminant analysis. Journal of the American Statistical Association, 82(400): 1118-1122. 109 Canty, M.J., A.A. Nielsen, and M. Schmidt. 2004. Automatic radiometric normalization of multitemporal satellite imagery. Remote Sensing of Environment. 91: 441-451. Carter, Rober E., Mark D. MacKenzie, and Dean H. Gjerstad. 1999. Ecological land classification in the Southern Loam hills of south Alabama. Forest Ecology and Management 114: 395-404. Chavez, P.S. Jr. 1988. An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data. Remote Sensing of Environment. 24: 459-479. Chavez, P.S. Jr. 1989. Radiometric calibration of Landsat Thematic mapper multispectral images. Photogrammetric Engineering and Remote Sensing. 55(9):1285-1294. Chavez, P.S. Jr. 1996. Image-based atmospheric corrections revisited and improved. Photogrammetric Engineering and Remote Sensing. 66(9): 1025-1036. Chen, J.M. 1996. Canopy architecture and remote sensing of the fraction of photosynthetically active radiation absorbed by boreal conifer forests. Ieee transactions of geoscience and remote sensing. 34(6): 1353-1368. Du, Y., P.M. Teillet, and J. Cihlar. 2002. Radiometric normalization of multitemporal high-resolution satellite images with quality control for land cover change detection. Remote sensing of Environment. 82: 123-134. Efron, B. 1975. The efficiency of logistic regression compared to normal discriminant analysis. Journal of the American Statistical Association. 70(352): 892-898. 110 Elvidge C.D., D. Yuan, R.D. Weerackoon, and R.S. Lunneta. 1995. Relative radiometric normalization of Landsat Multispectral Scanner (MSS) data using an automatic scattergram-controlled regression. Photogrammetric Engineering and Remote Sensing. 61(10):1255-1260. ERDAS IMAGINE. 1997. ERDAS field guide; fourth edition. ERDAS Inc. Atlanta, Georgia. 656 pp. ESRI 2005. ArcGIS Desktop Help. URL: http://webhelp.esri.com/arcgisdesktop/9.1/index.cfm?TopicName=welcome (last date accessed: 7 June 2005). Foody, G.M., and R.A. Hill. 1996. Classification of tropical forest classes from Landsat TM data. International journal of remote sensing. 17(12): 2353-2367. Foody, G. M. 1996. Fuzzy modeling of vegetation from remotely sensed imagery. Ecological Modeling. 85: 3-12. Foody, G. M. 2002. Status of land cover classification accuracy assessment. Remote Sensing of Environment. 80: 185-201. Frost, C.C. 1993. Four centuries of changing landscape patterns in the longleaf pine ecosystem. In: The longleaf pine ecosystem: ecology, restoration and management. Proc. Tall Timbers Fire Ecol. Conf. No.18. 17-43. Galv?o, L. S., A. R. Formaggio, and D. A. Tistot. 2005. Discrimination of sugarcane varieties in Southeastern Brazil with EO-1 Hyperion data. Remote Sensing of Environment. 94: 523-534. Gillespie, A.J.R. 1999. Rationale for a national annual forest inventory program. Journal of Forestry. 97(12):16-20. 111 Gopal, S. and C. Woodcock. 1996. Remote sensing of forest change using artificial neural networks. IEEE Transactions on Geoscience and Remote Sensing. 34(2): 398-404. Goebel, P.C., B.J. Palik, K. Kirkman, and M.B. Drew. 2001. Forest ecosystems of a lower gulf coastal plain landscape: multifactor classification and analysis. Journal of the Torrey Botanical Society. 128(1): 47-75. Hall, F.G., D.E. Strebel, J.E. Nickeson, and S. J. Goetz. 1991. Radiometric rectification: toward a common radiometric response among multidate, multisensor images. Remote Sensing of Environment. 35: 11-27. Hall, F.G., Y.E. Shimabukuro, and K.F. Huemmrich. 1995. Remote sensing of forest biophysical structure using mixture decomposition and geometric reflectance models. Ecological applications. 5(4): 993-1013. Halpern, M., W. C. Blackwelder, and J. I. Verter. 1971. Estimation of the mujltivarite logistic risk function: a comparison of the discriminant function and maximum likelihood approaches. Journal of Chronic Disease. 24: 125-158. Hasegawa, O. and T. Kurita. 2002. Face and non-face classification by multinomial logit models and kernel feature compound vectors. Proceedings of the 9 th International Conference on Neural Information Processing. 2: 996-1000. Hedman, C.W., S.L. Grace, and S.E. King. 1999. Vegetation composition and structure of southern coastal plain pine forests: an ecological comparison. Forest Ecology and Management. 134: 233-247. Hogland, J.S., N. Billor, and M.D. MacKenzie. in progress. Comparing polytomous logistic regression and discriminant analysis: a remote sensing perspective. 112 Hogland, J.S., and M.D. MacKenzie. in progress. Bringing images to a common radiometric scale using Aggregate No Change Regression: a comparison between radiometric normalization techniques. Holm, R.G., R.D. Jackson, B. Yuan, M.S. Moran, P.N. Slater, and S.F. Bigger. 1989. Surface reflectance factor retrieval from Thematic Mapper data. Remote Sensing of Environment. 27: 47-57. Homer, C.,G., R.D. Ramsey, T.C. Edwards, Jr., and A. Falconer, 1997. Land cover-type modeling using a multi-scene Thematic Mapper mosaic, Photogrammetric Engineering and Remote Sensing. 63: 59?67. Homer, C., C. Huang, L. Yang, B. Wylie and M. Coan. 2004. Development of a 2001 National Landcover Database for the United States. Photogrammetric Engineering and Remote Sensing. 70(7) 829-840. Hosmer, T., D. W. Hosmer, and L. L. Fisher, 1983. A comparison of the maximum likelihood and discriminant function estimators of the coefficients of the logistic regression model for mixed continuous and discrete variables. Communications in Statistics. B12: 577-593. Hosmer, D. W., S. Lemeshow.1989. Applied logistic regression. Whiley-Interscience New York. 307 pp. Hosmer, D. W., S. Lemeshow. 2000. Applied logistic regression second addition. Whiley-Interscience New York. 375 pp. Hossain, M., S.Wright, and L. A. Petersen. 2002. Comparing performance of multinomial logistic regression and discriminant analysis for monitoring access to care for acute myocardial infarction. Journal of Clinical Epidemiology. 55: 400-406. 113 Huang, C., Z. Zhang, L. Yang, B. Wylie, and C. Homer, 2002. MRLC 2001 Image Preprocessing Procedure, USGS White Paper. Iqbal, M. 1983. An introduction to solar radiation. Academic Press; Toronto. 390 pp. Jakubauskas, M.E. 1996. Thematic mapper characterization of lodgepole pine seral stages in Yellowstone National Park, USA. Remote Sensing of Environment. 56: 118- 132. Jensen, J.R. (Editor) 1983. Urban/suburban land use analysis. In: Manual of Remote Sensing, 2 nd ed. American Society of Photogrammetry. Vol. 2, pp. 1571-1666. Jensen, J.R. 1986. Introductory digital image processing : a remote sensing perspective. : Prentice-Hall; Englewood Cliffs, N.J. 378 pp. Jensen, J.R. 2000. Remote sensing of the environment: an earth resource perspective. Prentice Hall; 2nd edition, Upper Saddle River, N.J. 544 pp. Johnson, R.A. and D.W. Wichern. 2002. Applied multivariate statistical data analysis. Prentice Hall; 5 th edition, Upper Saddle River, N.J. 767 pp. Kelly, J.F., and W.A. Bechtold. 1990. The longleaf pine resource. In: Proceedings of symposium on the management of longleaf pine, 1989. April 4-6; Long Beach, MS. Department of Agriculture, Forest Service, Southern Forest Experiment Station, 11-22. Keuchel, J., S. Naumann, M. Heiler, and A. Siegmund. 2003. Automatic land cover analysis for Tenerife by supervised classification using remotely sensed data. Remote Sensing of Environment. 86: 530-541. 114 Lear, D. H., W.D. Carroll, P.R. Kapeluck, and R. Johnson. 2005. History and restoration of the longleaf pine-grassland ecosystem: implications for species at risk. Forest Ecology and Management. 211: 150-165. Lillesand, T. M. and R. W. Kiefer. 1994. Remote sensing and image interpretation. Wiley and Sons; 3 rd edition, New York. 750 pp. Little, E.L., Jr. 1971. Atlas of United States trees. Volume 1, Conifers and Important Hardwoods: U.S. Department of Agriculture Miscellaneous Publication 1146, 9 pp., 200 maps. Mahopatra, K. and S. Kant. 2005. Tropical deforestation: a multinomial logistic model and some country-specific policy prescriptions. Forest Policy and Economics. 7: 1-24. Makela, H., and A. Pekkarinen. 2001. Estimation of timber volume at the sample plot level by means of image segmentation and Landsat TM imagery. Remote Sensing of Environment. 77: 66-75. Markham, B.L., and J.L. Barker, 1986. Landsat MSS and TM post-calibration dynamic ranges, exoatmospheric reflectances and at-satellite temperatures. EOSAT Landsat Technical Notes, 1: 3-8. McDonald, A.J., F.M. Gemmell, and P.E. Lewis. 1998. Investigation of the utility of spectral vegetation indicies for determining information on coniferous forests. Remote Sensing of Environment. 66: 250-272. Metternicht, G.I. 2003. Categorical fuzziness: a comparison between crisp and fuzzy class boundary modeling for mapping salt-affected soils using Landsat TM data and a classification based on anion ratios. Ecological Modelling. 168: 371-389. 115 Moran, M.S., R.D. Jackson, P.N. Slater, and P.M. Teillet. 1992. Evaluation of simplified procedures for retrieval of land surface reflectance factors from satellite sensor output. Remote Sensing of Environment. 41: 169-184. [NASA] National Aeronautics and Space Administration. Landsat 7 Science data Users Handbook, URL: http://ltpwww.gsfc.nasa.gov/IAS/handbook/handbook_toc.html (last date accessed: 5 January 2005). Neter, J., M. Kutner, C.J. Nachtsheim, and W. Wasserman. 1996. Applied linear statistical models. McGraw-Hill. Boston, Massachusetts. 1408 pp. Noel, J.M., W.J. Platt, and E.B. Moser. 1998. Structural characteristics of old-and second-growth stands of longleaf pine (Pinus palustris) in the Gulf Coastal region of the U.S.A. Conservation Biology. 12(3): 533-548. Noss, R.F., E.T. LaRoe, and J.M. Scott. 1995. Endangered ecosystems of the United States: A preliminary assessment of loss and degradation. Biological Report 28. National Biological Service, Washington, D.C. Olsson, H. 1993. Regression functions for multitemporal relative calibration of Thematic Mapper data over boreal forest. Remote Sensing of Environment. 46:89-102. Olthof, I., C. Butson, and R. Fraser. 2005a. Signature extension through space for northern landcover classification: a comparison of radiometric correction methods. Remote Sensing of Environment. 95: 290-302. Olthof, I., D. Pouliot, R. Fernandes, and R. Latifovic. 2005b. Landsat-7 ETM+ radiometric normalization comparison for northern mapping applications. Remote Sensing of Environment. 95: 388-398. 116 Omernik, J.M. 1987. Ecoregions of the Conterminous United States. Annals of the Association of American Geographers. 77(1): 118-125. Outcalt, K.W., and R. M. Sheffield. 1996. The longleaf pine forest: trends and current conditions. Resource Bulletin SRS-9 Asheville, NC; U.S. Department of Agriculture, Forest Service, Southern Research Station. 23 pp. Pal, M. and P. M. Mather. 2003. An assessment of the effectiveness of decision tree methods for land cover classification. Remote Sensing of Environment. 86: 554- 556. Pax-Lenney, M., C.E. Woodcock, S.A. Macomber, S. Gopal, and C. Song. 2001. Forest mapping with a generalized classifier and Landsat TM data. Remote Sensing of Environment. 77: 241-250. Peet, R. K., and D. J. Allard. 1993. Longleaf Pine Vegetation of the Southern Atlantic and Eastern Gulf Coast Regions: A Preliminary Classification. In: The longleaf pine ecosystem: ecology, restoration and management. Proc. Tall Timbers Fire Ecol. Conf. No. 18. Peng, C. J. and R. N. Nichols. 2003. Using multinomial logistic models to predict adolescent behavioral risk. Journal of Modern Applied Statistical Methods. 2(1): 1-13 Perry, J.N., A.M. Liebhold, M.S. Rosenberg, J. Dungan, M. Miriti, A. Jakomulska, and S. Citron-Pousty. 2002. Illustrations and guidelines for selecting statistical methods for quantifying spatial pattern in ecological data. Ecography. 25: 578-600. Press, S. J. and S. Wilson. 1978. Choosing between logistic regression and discriminant analysis. Journal of the American Statistical Association. 73(364): 699-705. 117 Ram?rez-Garc?a, P., J. L?pez-Blanco, and D. Oca?a. 1998. Mangrove vegetation assessment in the Santiago River Mouth, Mexico, by means of supervised classification using Landsat TM imagery. Forest Ecology and Management. 105: 217-229. Reese, H.M., T.M. Lillesand, D.E. Nagel, J.S. Stewart, R.A. Goldmann, T.E. Simmons, J.W. Chipman, and P.A. Tessar. 2002. Statewide land cover derived from multiseasonal Landsat TM data a retrospective of the WISCLAND project. Remote Sensing of Environment. 82: 224-237. SAS/STAT Software: Changes and Enhancements, Release 8.2, Cary, NC: SAS Institute Inc., 2001. SAS/STAT 2005 SAS OnlineDoc version eight URL: http://v8doc.sas.com/sashtml/ (last date accessed: 20 June 2005). Schmidt, K.S., A.K. Skidmore. 2003. Spectral discrimination of vegetation types in a coastal wetland. Remote Sensing of Environment. 85: 92-108. Schott, J.R., C. Salvaggio, and W.J. Volchok. 1988. Radiometric Scene Normalization using Pseudoinvarinat Features. Remote Sensing of Environment. 26: 1-16. Skole, D.L., C.O. Justice, J.R.G. Townshend, and A.C. Janetos. 1997. A land cover change monitoring program: strategy for an international effort. Mitigation Adaptation Strategies Global Change. 2:157-175. Song, C., C.E. Woodcock, K.C. Seto, M.P. Lenney, and S.A. Macomber. 2001. Classification and change detection using landsat tm data: when and how to correct atmospheric effects? Remote Sensing of Environment. 75: 230-244. 118 Teillet, P.M., and G. Fedosejevs. 1995. On the dark target approach to atmospheric correction of remotely sensed data. Canadian Journal of Remote Sensing. 21: 373-387. Tuldge, C. 1999. Plant (protecting endangered species). Index on Censorship. 28(6): 164- 166. U.S. Fish and Wildlife Service (USFWS). 2000. U.S. Fish and Wildlife Service region 4 ecosystem. URL: http://www.fws.gov/stand/standards/allshape.html (last date accessed: 25 July 2003). Vogelmann, J.E., D. Helder, R. Morfitt, M. Choate, J.W. Merchant, and H. Bulley. 2001 Effects of landsat 5 thematic mapper and landsat 7 enhanced thematic mapper plus radiometric and geometric calibrations and corrections on landscape characterization. Remote Sensing of Environment. 78: 55-70. Velloso, M.L., and F.J De Souza. Non-Parametric smoothing for relative radiometric correction on remotely sensed data. Proceedings of the XV Brazilian Symposium on Computer graphics and Image Processing. IEEE International 2002. Vogelmann, J.E., S.H. Howard, L. Yang, C.R. Larson, B.K. Wylie, and N. Van Driel. 2001. Completion of the 1990s national land cover data set for the conterminous United States from Landsat Thematic Mapper data and ancillary data sources. Photogrammetric Engineering and Remote Sensing. 67, 650-662. 119 Weakley, A.S., K.D. Patterson, S. Landaal, M. Pyne, R.E. Evans, J.A. Teague, M.J. Russo, and others (compilers). 2000. International classification of ecological communities: terrestrial vegetation of the southeastern united states. Alabama subset for GAP Analysis Program report from biological conservation data systems and working draft of April 2000. Association for biodiversity information/The Nature Conservancy, southern resource office, community ecology group. Durham, North Carolina. Wrigley, N., 1985.Categorical data analysis for geographers and environmental scientists. Longman New York. 385 pp. Yang, L., S.V. Stehman, J.H. Smith, and J.D. Wickham. 2001. Thematic Accuracy of MRLC Land Cover for the Eastern United States. Remote Sensing of Environment, 76:418-422. Yoshida, T. and S. Omatu. 1994. Neural network approach to land cover mapping. IEEE Transactions on Geoscience and Remote Sensing. 32(5): 1103-1109. Yu, B., I.M. Ostland, P. Gong, and R. Pu. 1999. Penalized discriminant analysis of in situ hyperspectral Data for conifer species recognition. Ieee transactions on Geoscience and remote sensing. 37(5): 2569-2577. Yuan, D., and C.D. Elvidge. 1996. Comparison of relative radiometric normalization techniques. ISPRS Journal of Photogrammetry and Remote Sensing. 51: 117-126 APPENDIX 1 Standardized slope and intercept estimates for each stage of level 1 in our hierarchical classification scheme (HCS). HCS stage Class W Band 1 W Band 2 W Band 3 W Band 4 W Band 5 W Band 7 S/F band 1 S/F band 2 S/F band 3 S/F band 4 S/F band 5 S/F band 7 LUC W, S/F -3.1161 -0.2767 -1.7509 3.4166 -2.0383 2.7979 5.3104 -3.334 -1.1294 3.9077 2.6607 2.4798 LUC W -1.2517 -1.5296 -4.0133 5.1994 -7.0932 6.8795 6.5879 -1.4274 -3.5158 4.9499 2.2323 1.2464 LUC S/F -0.4954 -2.198 1.3685 3.5964 -0.2254 -0.123 3.1765 -4.7399 -1.8856 4.5854 -3.4976 5.9296 Burnt Areas W 5.6347 -1.1848 -5.4059 6.7049 -6.3506 4.7961 -4.3463 -4.6395 -6.6306 5.0609 1.7377 5.9921 Burnt Areas S/F -1.1742 0.0216 -4.1659 5.339 -1.6305 0.7172 6.0852 -2.0864 -3.7443 2.132 -0.5385 4.6576 Smoke W -2.5441 -0.4788 -4.1895 4.1024 0.4597 -1.4046 7.4427 -5.2787 -1.8632 3.7648 2.4116 4.1279 Smoke S/F -2.7135 1.3773 -2.5504 4.7307 -2.172 1.1576 9.767 -0.9951 -3.9054 6.2497 -1.2994 4.2329 Clouds S/F -5.9394 4.1031 -6.8712 10.254 2.397 1.6716 11.7622 -17.2281 -1.8363 -2.0867 -7.8979 -4.6273 Shadow S/F -4.4426 0.1945 -2.1934 5.8292 -2.5155 2.9003 14.7463 -1.379 -9.9206 4.7731 -0.1159 4.4327 Fields -2.7271 1.7116 -3.1128 6.2186 -5.0633 4.2753 3.8808 -4.3941 -1.2421 6.9918 2.2538 3.2201 City -1.2263 1.0165 -1.8775 5.4121 -5.0909 3.434 6.5285 -2.6768 -0.1623 5.8117 -3.5973 5.4114 Deciduous -1.0227 -1.2741 -2.1514 5.7059 -2.3372 1.3384 4.643 -3.8568 -9.3526 5.2498 4.8278 0.4449 Evergreem -2.2348 0.9975 -1.6391 6.4748 -3.2161 -0.9606 1.0293 -2.7864 -3.9582 4.895 0.2099 2.5405 1 Wet Vegetated Areas 0.0306 -0.1799 -0.8438 2.908 -1.5621 -0.4837 4.2482 3.2286 -5.0042 3.1406 3.7382 -1.8564 LUC S/F - - - - - - 2.1235 -5.6794 1.6586 4.0858 4.577 2.2241 Fields - - - - - - 1.4017 -4.9152 1.676 7.341 2.5559 3.7537 City - - - - - - 2.5168 -2.2715 1.5472 5.5785 -0.8337 3.6504 Deciduous - - - - - - 1.7535 -4.4902 -2.614 4.7194 7.8462 -3.8358 Evergreem - - - - - - 1.369 -2.6778 -0.7064 4.7138 1.3701 -1.5604 2 Wet Vegetated Areas - - - - - - 2.1805 0.8503 -1.8881 2.0731 6.7956 -4.3077 LUC W 0.4676 -4.7482 1.0839 4.123 2.018 1.1106 - - - - - - Fields -0.5458 -0.8013 -1.8264 7.6573 -0.7847 3.4555 - - - - - - City 0.5427 -0.0548 -0.897 6.7515 -4.1107 3.9544 - - - - - - Deciduous 1.3189 -5.8125 -1.9078 6.4299 3.0132 -0.019 - - - - - - Evergreem -0.513 -4.6406 -0.3922 6.9721 2.0607 -3.3962 - - - - - - 3 Wet Vegetated Areas 2.5668 -2.9023 -0.7863 3.6629 3.0243 -2.5741 - - - - - - 120 121 Standardized slope estimates for each stage of level 2 in our hierarchical classification scheme. stage Ecosystem W band 1 W band 2 W band 3 W band 4 W band 5 W band 7 S/F band 1 S/F band 2 S/F band 3 S/F band 4 S/F band 5 S/F band 7 Elevation Slash -0.214 -0.117 -0.123 -0.251 0.177 - -0.421 0.711 -0.369 0.777 -0.144 - -0.742 Hardwood 0.799 0.735 -1.119 -1.500 2.206 - 0.873 -1.409 -1.399 1.065 0.588 - -0.192 Mixed Hardwood/Pine 0.946 0.473 -1.228 -0.864 2.079 - 0.301 -0.283 -0.691 0.297 -0.248 - 1.026 Mountain Longleaf 0.526 0.191 0.402 -0.124 -0.241 - -0.920 -0.197 -0.017 0.341 0.032 - 2.197 1 Loblolly 1.205 0.418 -1.174 -0.443 0.609 - 0.801 -0.330 -1.100 -0.840 0.591 - 2.106 Slash - - - - - -0.576 0.638 -0.369 0.573 -0.200 0.098 -1.033 Hardwood - - - - - 1.745 -1.740 -2.004 1.025 2.530 -0.519 1.045 Mixed Hardwood/Pine - - - - - 1.141 -0.370 -1.240 0.215 0.901 -0.027 1.641 Mountain Longleaf - - - - - -0.456 0.027 -0.033 0.301 0.507 -0.448 2.960 2 Loblolly - - - - - 1.577 -0.132 -1.354 -1.083 1.212 -0.544 2.616 Slash -0.027 0.111 -0.156 0.342 -0.810 0.811 - - - - - -0.906 Hardwood 1.648 0.480 -1.855 -1.584 3.099 -0.831 - - - - - 1.056 Mixed Hardwood/Pine 1.296 0.462 -1.789 -0.795 2.105 -0.556 - - - - - 1.311 Mountain Longleaf 0.367 0.010 0.343 0.155 0.146 -0.684 - - - - - 2.713 3 Loblolly 1.517 0.496 -1.438 -0.938 1.597 -1.559 - - - - - 2.080 121 122 APPENDIX 2 A) Class probabilities for each stage of our hierarchical classification can be calculated using tables 3 and 4 as follows: given that ( ) () ().-1 then ,11 ,1 1 1 hJ xx x ? ? ? = = == J h j j ,...,J-j ?? ? (1) Using the linear form of the logit ( ) () [] ij J j xx?x x x ?=??= ? ? ? ? ? ? ? ? ,1ln ? ? (2) and substituting () Jh for -1 ?? x ? ? = 1 1 J h (3) ? j and ? J are solved as follows; () ( ) () ? ? = ?+ ? = 1 1 1 J h j j j ?x ?x x exp exp ? (4) () () exp ? ? = ?+ = 1 1 1 1 J h j J ?x x? (5) Where ? j (x) = the mean probability of group j, with class means and variances = n? j ; n? j (1-? j ) respectively (Agresti 2002) given the values (x) and slope estimates (?) of explanatory variables in matrix notation. 123 B) Prior probabilities can be incorporated into stage 2 of our hierarchical classification by taking the natural log transformation of the odds of a given class, subtracting the natural log transformation of the ratio of that given class? original prior abundance estimates divided by the baseline class? prior abundance, and adding the natural log transformation of the ratio of that given class? prior probability by the baseline class?s prior probability as follows: Adjusted log odds j = ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? J prior j prior J j J j ln abundanceprior abundanceprior lnln ? ? (6) Manipulating equations 4 and 5 adjusted class probabilities can be estimated as follows; Adjusted ? ? = ? ? ? ? ? ? + ? ? ? ? ? ? = 1 1 1 J j j j odds log Adjustedexp j odds log Adjustedexp ? (7) Adjusted () odds log Adjustedexp j? ? = + = 1 1 1 1 J j J ? (8) APPENDIX 3 Standard error estimates for the explanatory variables of each stage of level 1 in our hierarchical classification scheme. W and S/F identify the period in which ETM+ bands were acquired and LUC and TF were identified for a winter leaf-off and summer/fall leaf-on season, respectively. HCS Stage Logits intercept W Band 1 W Band 2 W Band 3 W Band 4 W Band 5 W Band 7 S/F band 1 S/F band 2 S/F band 3 S/F band 4 S/F band 5 S/F band 7 LUC W, S/F 4.8229 0.1434 0.1463 0.086 0.0878 0.0984 0.1129 0.1185 0.1324 0.0846 0.0465 0.0818 0.1186 LUC W 4.5928 0.139 0.1438 0.0895 0.0848 0.0986 0.1129 0.116 0.129 0.0889 0.0452 0.0836 0.1231 LUC S/F 4.6715 0.142 0.1462 0.0866 0.0863 0.0988 0.1149 0.1205 0.1315 0.0865 0.045 0.0832 0.1224 Burnt W 9.244 0.2891 0.372 0.1904 0.1031 0.1208 0.1193 0.2842 0.4548 0.234 0.0937 0.152 0.2186 Burnt S/F 4.3171 0.1318 0.13 0.0841 0.0841 0.0971 0.1131 0.1078 0.1196 0.0829 0.0416 0.0799 0.1166 Smoke W 4.9291 0.1465 0.1608 0.0943 0.0864 0.1047 0.1245 0.1255 0.1473 0.0887 0.048 0.0823 0.1179 Smoke S/F 4.5308 0.1368 0.1403 0.0833 0.0847 0.0982 0.1135 0.1122 0.1237 0.0843 0.0439 0.0815 0.1177 Clouds S/F 72.5307 1.4017 1.8499 0.8682 0.4988 1.0174 1.6292 1.791 2.2977 2.1181 0.8093 1.728 1.8594 Shadow S/F 4.5393 0.1376 0.1414 0.0835 0.0849 0.0981 0.1133 0.1154 0.1265 0.0883 0.0441 0.0814 0.1169 Fields 4.6617 0.139 0.1423 0.0837 0.0844 0.0974 0.1116 0.1145 0.123 0.0825 0.0434 0.0819 0.1171 Urban / Trans. / Bare 4.605 0.1373 0.144 0.0832 0.0849 0.0976 0.1109 0.1143 0.1241 0.083 0.0442 0.0814 0.1168 Deciduous 4.4273 0.1338 0.1341 0.082 0.0842 0.0973 0.1129 0.111 0.119 0.0836 0.0433 0.0823 0.1208 Evergreen 4.4878 0.1355 0.1388 0.0874 0.0845 0.0995 0.1178 0.1149 0.1245 0.0888 0.0445 0.0839 0.1235 1 Wet Vegetated 4.3485 0.1286 0.1214 0.0753 0.0828 0.0957 0.1095 0.1047 0.1046 0.0726 0.0417 0.0786 0.115 LUC S/F 3.9467 - - - - - - 0.1069 0.1303 0.0745 0.0473 0.0797 0.1161 Fields 3.8998 - - - - - - 0.1056 0.1263 0.0737 0.046 0.0797 0.1162 Urban / Trans. / Bare 3.9505 - - - - - - 0.1056 0.124 0.0725 0.0469 0.0793 0.1158 Deciduous 3.6805 - - - - - - 0.1014 0.122 0.076 0.0451 0.0797 0.1191 Evergreen 3.5706 - - - - - - 0.1012 0.1206 0.0749 0.045 0.0794 0.1195 2 Wet Vegetated 3.5777 - - - - - - 0.0957 0.1091 0.0658 0.0421 0.076 0.1108 LUC W 2.4528 0.0699 0.0942 0.0526 0.0453 0.0534 0.0676 - - - - - - Fields 2.5061 0.0711 0.0957 0.0525 0.0441 0.0533 0.0674 - - - - - - Urban / Trans. / Bare 2.5085 0.0703 0.0965 0.0527 0.0447 0.0533 0.0671 - - - - - - Deciduous 2.3743 0.0676 0.0928 0.0524 0.0438 0.0539 0.0688 - - - - - - Evergreen 2.6171 0.0742 0.103 0.0648 0.0443 0.0596 0.0809 - - - - - - 3 Wet Vegetated 2.1213 0.056 0.0685 0.0435 0.0411 0.0491 0.0616 - - - - - - 124 Standard error estimates for the explanatory variables of each stage of level 2 in our HICS. HCS stage Logits intercept W band 1 W band 2 W band 3 W band 4 W band 5 W band 7 S/F band 1 S/F band 2 S/F band 3 S/F band 4 S/F band 5 S/F band 7 Elevation Slash 2.399 0.055 0.078 0.056 0.020 0.027 - 0.070 0.083 0.052 0.029 0.030 - 0.003 Hardwood 3.439 0.085 0.121 0.082 0.029 0.035 - 0.109 0.122 0.084 0.037 0.041 - 0.004 Mixed Hardwood/Pine 2.592 0.062 0.090 0.065 0.022 0.030 - 0.082 0.093 0.062 0.031 0.033 - 0.002 Mountain Longleaf 4.532 0.087 0.127 0.098 0.029 0.040 - 0.119 0.133 0.089 0.046 0.048 - 0.003 1 Loblolly 2.705 0.066 0.095 0.070 0.024 0.032 - 0.086 0.101 0.068 0.035 0.036 - 0.003 Slash 2.050 - - - - - - 0.053 0.075 0.051 0.022 0.030 0.053 0.003 Hardwood 2.677 - - - - - - 0.070 0.100 0.072 0.028 0.047 0.087 0.003 Mixed Hardwood/Pine 2.026 - - - - - - 0.052 0.076 0.055 0.022 0.033 0.059 0.002 Mountain Longleaf 3.812 - - - - - - 0.095 0.118 0.084 0.036 0.054 0.101 0.003 2 Loblolly 2.223 - - - - - - 0.058 0.087 0.061 0.026 0.037 0.065 0.002 Slash 1.619 0.046 0.070 0.049 0.017 0.033 0.049 - - - - - - 0.003 Hardwood 1.943 0.060 0.092 0.066 0.024 0.041 0.057 - - - - - - 0.002 Mixed Hardwood/Pine 1.657 0.050 0.076 0.056 0.019 0.036 0.051 - - - - - - 0.002 Mountain Longleaf 2.922 0.071 0.110 0.087 0.025 0.051 0.077 - - - - - - 0.002 3 Loblolly 1.745 0.052 0.080 0.059 0.020 0.038 0.056 - - - - - - 0.002 125