Using Landscape Spatial Relationships to Improve Estimates of Land-Cover Area from Coarse Resolution Remote Sensing

Using Landscape Spatial Relationships to Improve Estimates of Land-Cover Area from Coarse Resolution Remote Sensing Aaron Moody*

A

two-stage modeling strategy significantly improves land-cover area estimates from low spatial resolution remote sensing by correcting measurements of class proportions within large blocks of pixels. Vegetation classtype information is developed through supervised classification of Thematic Mapper spectral data at both fine (30 m) and coarse (1020 m) resolutions. Stage 1 models use measurements of landscape spatial properties to estimate the slopes and intercepts of proportion transition relationships between fine- and coarse-resolution classes within randomly located pixel blocks. Following this step, a Stage II model uses a linear estimator to predict true class proportions based on measured coarse-scale proportions and the slope and intercept estimates from the Stage I models. Model development and testing on a training site is followed by testing and inversion for a validation site. Model inversion involves using spatial variables measured only at the coarse resolution as input to the Stage I models. For both the training and the validation data, the procedure results in a statistically significant reduction in error when estimating land-cover area by class type within the sampling blocks. Elsevier Science Inc., 1998

INTRODUCTION It has been observed that, for a given land-cover taxonomy, spatial resolution will influence land-cover area estimates determined from remote sensing through image classification (Moody and Woodcock, 1994; Mayaux and

* Department of Geography, University of North Carolina at Chapel Hill Address correspondence to Aaron Moody, Dept. of Geography, Univ. of North Carolina at Chapel Hill, Chapel Hill, NC 27599-3220. E-mail: [email protected] Received 14 March 1997; revised 19 December 1997. REMOTE SENS. ENVIRON. 64:202–220 (1998) Elsevier Science Inc., 1998 655 Avenue of the Americas, New York, NY 10010

Lambin, 1995; Pax-Lenney and Woodcock, 1996). These resolution effects are attributed in part to a) a general tendency for dominant classes to increasingly dominate and secondary classes to increasingly dissipate at coarser resolutions (Turner et al., 1989); b) modulations of this basic effect that are related to the spatial organization of the classes in the landscape (Moody and Woodcock, 1996); and c) effects caused by the mixing of class spectra in the surface–atmosphere–sensor data stream (Moody and Woodcock, 1995; Starmer and Frohn, 1996). Models of these resolution effects can help improve estimates of land-cover area as well as changes in landcover area. There are at least two approaches to account for scale-dependent areal bias. Mixture models, incorporated into the image classification process, can provide estimates of subpixel composition assuming that linear or other mixing assumptions are adequate, that the classes are spectrally separable, and that pure class spectra are known or can be estimated (Oleson et al., 1995; Roberts et al., 1993; Adams et al., 1986). Alternatively, calibrated models, applied in a postclassification mode, can improve area estimates from coarse resolution data if the relationship between “true” and measured proportions can be modeled (Mayaux and Lambin, 1995; Kalkhan et al., 1995). This article introduces a two-stage modeling procedure that provides postclassification correction of landcover area estimates determined from coarse-resolution remote sensing. The basic model (Stage II model) is a linear model whose parameters characterize the relationship between true class proportions and coarse-scale proportion estimates. The parameters for the Stage II model are themselves statistically determined using variables that characterize the spatial organization of the landcover classes. These latter models are referred to below as Stage I models. 0034-4257/98/$19.00 PII S0034-4257(98)00003-0

Improved Area Estimates from Remote Sensing

Development and initial testing of the two-stage modeling approach proceeds using data for the Plumas National Forest in the northern Sierra Nevada Mountains in California. The trained models are then validated on the Stanislaus National Forest, located roughly 28 south of the Plumas. The approach presented in this paper is intended to support the following three objectives. • To incorporate into the correction procedure information about the landscape spatial organization which is considered to influence the degree of bias in coarse-resolution estimates of landcover area. • To develop an invertible area-correction procedure that can provide improved area estimates when only coarse-resolution data are available (that is, where there is no high-resolution training data). • To provide a model that is not cover-type specific and can, therefore, more easily generalize to a range of environmental settings. The goal of this article is to introduce a methodology that has the potential to fulfill the above criteria, and to present a first level of validation of the approach by applying it to a site that is independent from the training site, but of the same general type. The purpose of the methodology is to provide improved estimates of landcover area both on a grid-cell basis and on a regional basis. BACKGROUND Areal accuracy in land-cover data sets is important for two primary reasons. First, many models of surface processes, such as net primary productivity and surface– atmosphere fluxes of energy, carbon, and water, require land-cover data to provide state variables or parameters (Sellers et al., 1986; Foley, 1994; Dickinson, 1995). In these cases, land cover is a surrogate for more physiologically or structurally meaningful parameters like leaf area index or fraction of absorbed photosynthetically active radiation (Townshend et al., 1991). For global models that operate on a coarse grid, modeling strategies may employ area-integrated estimates of surface variables that can be inferred from low spatial resolution remote sensing. Recent research has focused on accommodating subgrid-cell effects through distributed parametrization of such models (Avissar, 1992; Koster and Suarez, 1992; Seth et al., 1994). However, for land cover data, even if subgrid-scale heterogeneity is explicitly handled, the original coarse resolution data from which the subgrid components are determined may contain significant bias in class proportions (Moody and Woodcock, 1994). At present, these errors may be within the error range of most coarse-scale models. However, as models of Earth system processes improve, and as finer grids are implemented, proportional bias may become increasingly important (Strahler et al., 1995).

203

Second, global land cover data are critical for monitoring the location, extent, and changes of major vegetation assemblages and ecological units. Monitoring surface cover supports a range of research related to modeling and understanding the response of vegetation to natural and anthropogenic forcings (Vitousek, 1994). Such forcings include climate change, pollution, fire, deforestation, cultivation, and watershed disruption. Although costly and logistically difficult, it is possible to perform global scale land cover assessments using high resolution data such as Landsat Multispectral Scanner (MSS) or Thematic Mapper (TM). For example, the Food and Agriculture Organization (FAO) 1990 Forest Resources Assessment is based in part on a stratified random sample of MSS scenes for the tropics (FAO, 1996). This effort is one of the most thorough and statistically rigorous attempts at broad-scale mapping of vegetation and would provide an excellent opportunity to test postclassification area correction methods such as that presented here. On an operational basis, however, it is impractical to perform as intensive an assessment of global vegetation as the FAO work if, for example, annual assessments are required. For the purposes of this work, it is assumed that global scale vegetation monitoring can be accomplished most realistically using high-temporalfrequency, remotely sensed data at low spatial resolution. There are limitations to accurately quantifying the area of existing vegetation types at the global scale using coarse-resolution remote sensing. In part these limitations relate to the quality and characteristics of data from the Advanced Very High Resolution Radiometer (AVHRR). The AVHRR is the primary sensor currently in use for global-scale remote sensing studies and suffers from a variety of quality constraints related to geodetic control, spectral resolution, radiometric calibration, and pixel distortion due to off-nadir viewing (Moody and Strahler, 1994). Many of these issues will be ameliorated with the availability of data from advanced instruments such as the Moderate Resolution Imaging Spectroradiometer (MODIS) (Running et al., 1994; Strahler et al., 1995). However, area estimates will still contain errors resulting from aggregation effects at the coarse scales of observation (Henderson-Sellers et al., 1985; Moody and Woodcock, 1994). The magnitude and scale-dependence of these aggregation effects depends on the true proportions of the component cover types in the landscape, the spatial organization of the classes, and the statistical and spatial properties of the sampled surface-leaving radiance (Moody and Woodcock, 1995). As noted above, there are two basic approaches for handling aggregation effects in coarse resolution landcover data. One approach is to extract and preserve subpixel information by integrating spectral mixture models into the image classification procedure. Models based on linear mixing assumptions have used an additive linear form (Quarmby et al., 1992), spectral factor analysis

204

Moody

(Huete, 1986), or multiple linear regression (Puyou-Lascassies et al., 1994). Spectral mixture models have also been developed based on physical principles. For example, Borel and Gerstl (1994) presented a model based on radiosity to account for nonlinear mixing due to scene geometric-optical effects. Other information-preserving approaches such as fuzzy classifiers also are used (Woodcock and Gopal, 1992; Woodcock et al., 1994). Another approach, and the one followed in this article, is to use post-classification procedures to correct biased area estimates. Such models are typically trained on sites for which both high and low spatial resolution data are available. Training is initially based on the relationship between class proportions measured at the two scales of representation. Several studies of this nature have employed a model that is derived from classical regression and referred to as the inverse estimator (Brown, 1982; Czaplewski and Catts, 1992; Walsh and Burk, 1993; Kalkhan et al., 1995). This resembles a Markovtype linear model in which the relationship between actual and observed class proportions is realized by a class transition matrix. The elements of the matrix represent class transition probabilities relating the reference (high resolution) proportions to the coarse resolution proportions for the calibration data. In its simplest form, the inverse estimator appears as: tc5Pc*r,

(1)

where tc is a k31 vector of reference (true) proportions for the k classes in the calibration data, r is a k31 vector of proportions from the coarse resolution data, and Pc is the k3k matrix of class transition probabilities calculated from the calibration data (Kalkhan et al., 1995). Pc and r are then used to estimate t when only r is known. Another approach relies on a linear model relating coarse-resolution area estimates to true, or reference, proportions (Mayaux and Lambin, 1995; Roesch et al., 1995; Moody and Woodcock, 1996). This is a form of the “classical” model in which known but incorrect values are estimated based on unknown correct values (Brown, 1982). It is sometimes known as the direct estimator (e.g., Yuan, 1996) and will be referred to as such below. Most simply, if Pi0 are the reference proportions for each class i (where i5l, . . . ,k), and Pir are the measured proportions at resolution r, then Pir5b01bl·Pi01e. (2) ˆ This is rearranged to estimate Pi0 in cases where only Pir is known. Separate linear models for each class i or nonlinear models can also be used to characterize these relationships. Additionally, the model coefficients themselves can be modeled and parametrized based on measurable scene spatial and statistical characteristics (Mayaux and Lambin, 1995; Moody and Woodcock, 1995). This is the approach used here, where Eq. (2) is effectively the Stage II model referred to below.

Czaplewski and Catts (1992) and Walsh and Burk (1993) found the inverse estimator [Eq. (1)] to outperform the direct estimator [Eq. (2)]. Moody and Woodcock (1996) also compared the two models for several independent test sites. Their results suggested that the inverse estimator performed better than the direct estimator at coarser scales, but that the opposite was true at finer scales. A weakness with the inverse model is that it may break down when used in an extrapolative sense, and produce spurious results if the training site either has very different interclass associations or contains different component classes from other areas in the region where the model is applied. This is because it is based on a transition matrix that maps confusion patterns between specific class pairs in the training site to the region for which improved area estimates are desired. Poor estimates are therefore likely to occur if the region, or areas within the region, do not have the same class composition, or the same confusion patterns as the training site. Other authors have commented on the instability of the inverse estimator (Yuan, 1996; Jupp, 1989). Several characteristics make the direct estimator a logical choice for the work presented here. First, based on earlier work by Mayaux and Lambin (1995) and Moody and Woodcock (1996) it appears that the performance of the direct estimator might be improved if variables that characterize the spatial characteristics of the landscape are incorporated into the model. Second, as formulated in this work, the direct estimator is independent of class type, the number of classes present, and the specific interclass transitions. This implies that the direct estimator might extrapolate better than the inverse estimator, if it can be made to accommodate spatial characteristics of landscape units regardless of their class. Finally, a central goal is to work toward the development of an invertible area correction model that can provide improved area estimates using only coarse resolution data and not requiring any 30 m training data. This is not feasible using a model such as the inverse estimator, where prior probabilities of interclass transitions across scales are required, necessitating the availability of both fine and coarse resolution data. As suggested above, two basic spatial effects contribute to biased area estimates. A first-order effect is the tendency of large classes to increasingly dominate the landscape when measured at coarser scales. Accordingly, small classes tend to diminish in size. Figure 1 illustrates this effect for the data used in this study. Second-order effects refer to modulations of this basic pattern due to characteristics of the landscape spatial organization (Moody and Woodcock, 1995). These effects result in the scatter about the smoothed fit notable in Figure 1. In either case, scale-dependent changes in the apparent area of classes result from class membership transitions across scales. This effect can be characterized in terms of proportion transition lines relating true and coarse-

Improved Area Estimates from Remote Sensing

205

to the Stage II model to predict 30 m proportions based only on information measured at 1020 m for the validation site. Inversion will require that the landscape pattern measures are relatively scale-invariant, or stable across scales. METHODS

Figure 1. Relationship between 30 m and 1020 m proportions. Each point represents an individual class type within a 51 km2 sampling cell.

resolution proportions. The slope of such a line summarizes the rate of transition. If transition rates depend on landscape spatial organization, it is sensible to try to model the coefficients of the proportion transition lines using measures of spatial pattern. Mayaux and Lambin (1995) used a similar approach to calibrate forest and nonforest area in the tropical belt. In this binary case, they improved results by adjusting the model coefficients based on scene spatial information (Mayaux and Lambin, 1995). In this research, landscape spatial properties (identified below) are used to predict the slopes and intercepts of the proportion transition lines for a large set of sampling units in independent training and validation sites. These are referred to as the Stage I models. A five-part strategy is employed as follows: a) identify a set of scale invariant spatial measures that characterize patterns in vegetation; b) use a subset of these measures (calculated at 30 m) to develop Stage I models for predicting the slope (bˆ 1) and intercept (bˆ 0) of the proportion transition lines for the training site; c) predict 30 m proportions (Pˆi30) for the training site by applying the Stage II model using the measured 1020 m proportions (Pi1020), and the slopes and intercepts estimated from the Stage I models; d) apply the calibrated Stage I models on the validation site, repeat step c, and evaluate the success of the procedure; e) invert and evaluate the procedure by running the Stage I models using the spatial variables as measured at 1020 m, and supply the results

Data The Plumas and Stanislaus National Forests are used as model training and validation sites, respectively. Both forests straddle the Sierra Nevada Mountains in northern California. These sites feature extreme elevational gradients from the foothills of the Central Valley up through treeline in the Sierra Nevada high country. There is also a strong moisture gradient between the more mesic front range and the xeric east side of the divide. The two sites are spatially separated by roughly 28 of latitude, but are similar in general character, and therefore provide a reasonable scenario for model validation. The major vegetation growth forms for the two sites are conifer, brush, hardwood, barren/grass, and water. These are the classes used in this research. Both sites have been studied as part of a project to develop vegetation mapping and timber inventory procedures for the U.S. Forest Service (Woodcock et al., 1994). Land-cover maps were produced using Landsat Thematic Mapper (TM) imagery and unsupervised image classification supported by air-photo and field validation (Woodcock et al., 1994). The accuracy of these maps has been assessed using methods based on fuzzy sets (Woodcock et al., 1994). Using a stringent MAX operator, the general land-cover classes are 84% and 79% accurate for the Plumas and Stanislaus sites, respectively. These values increase to 93% and 88% using a RIGHT operator (Woodcock and Gopal, 1992; Woodcock et al., 1994). Data Rescaling and Classification For each site, the TM spectral data are spatially resampled to 1020 m resolution, and reclassified at both 30 m and 1020 m. Resampling the TM data to 1020 m resolution is accomplished simply by averaging the spectral values within 34334 windows of 30 m pixels in each band. While it is considered adequate for the present purposes, this approach is somewhat unrealistic in that it assumes a square point spread function whose boundaries are the edges of the coarse pixel. Using training regions identified from the original class maps, statistical distributions of the spectral data for Bands 2–5 are estimated for each class using the 30 m TM data. These training data are used both to reclassify the 30 m data and to classify the 1020 m data using a maximum likelihood decision rule. This procedure is followed independently for the two sites. The resulting 30 m and 1020 m class maps are used for the remainder of the analysis.

206

Moody

Sampling and Development of Spatial Measures A set of noncontiguous, 51 km2 blocks of pixels are randomly distributed within each site. Each block contains 56,644 30-m pixels, and 49 1020-m pixels. The Plumas contains 55 blocks for model development and initial

evaluation. The Stanislaus contains 35 blocks for model validation and model inversion (see Fig. 2). The number of blocks represents approximately 30% of all possible such blocks from each site. Within each block, the following measurements are collected: 30 m and 1020 m

Figure 2. Stanislaus National Forest. Validation data set at 1020 m with sampling cells. From darker to lighter tones, the classes are hardwood, conifer, brush, and barren.

Improved Area Estimates from Remote Sensing

proportions for each class, a set of spatial measures at 30 m, and the same set of spatial measures at 1020 m. Additionally, within each sampled block, total proportion error and proportion error by class are measured both before and after the calibration model is applied in order to evaluate model performance. The spatial measures used to estimate the slope and intercept coefficients in the Stage I models are chosen based on three criteria. First, they must provide information about the spatial organization of the landscape that, based on existing research, is considered important in determining how land-cover proportions change with scale. Second, they must be reasonably self-correlated across scales to allow model inversion. Third, they must prove to be statistically significant for estimating the slope and intercept coefficients that will be used in the Stage II model. The spatial measures are calculated within each sampling block using the r.le software developed by Baker and Cai (1992). The first measure is the difference in size between the largest and the smallest classes in a block, or the range (rng). The second measure is entropy (ent) which is maximized when all pixels of a given class are as far away from one another as possible. Entropy is defined as k

ent52 o i5

k

o Pij·ln(Pij), 1 j 1

where k is the number of classes present in the block and Pij are elements of a k3k cooccurence matrix representing adjacency probabilities between classes i and j. The third measure is the variance (var) of the 1020 m NDVI data within each block, where NDVI is NDVI5

NIR2Red NIR1Red

(4)

and the variance is calculated as n

var5s210205

(NDVI2NDVI)2 o i 1 5

n21

.

k

i5

k

[(i2j)2·Pij]. o 1 j 1 5

(6)

The fifth measure is the Shannon index (shn) which responds both to the number of classes within a block (richness) and to the evenness with which their proportions are distributed. It is defined as k

H952 o (Pi·ln(Pi)), i 51

shn1020

ln(c)1020

ent1020

rng1020

0.74 — — —

— 0.46 — —

— — 0.84 —

— — — 0.67

shn30 ln(c)30 ent30 rng30

spatial measures used as the independent variables in the Stage I models. The Shannon Index (shn), entropy (ent) and range (rng) are all strongly self-correlated across scales, while contrast [ln(c)] is only moderately correlated. Since within-block variance in NDVI (var) is only determined for the 1020 m data, there is no cross-scale correlation for this variable. The correlations among variables are shown in Table 2. Note that there are some fairly strong correlations, especially between ent and shn, ent and rng, and shn and rng. Modeling Stage II Model The Stage II model is based on a linear relationship that approximates the transitions in land cover proportions between scales. For the two resolutions considered in this analysis, this appears simply as Pi10205b01b1·Pi30,

(7)

where Pi is the proportion of class i in the sampling cell. Table 1 provides the cross-scale correlations for the

(8)

where Pi30 is the proportion of class i at 30 m resolution and Pi1020 is the proportion of class i at 1020 m. Normally, this takes the form of a regression equation (with an error term), and the b coefficients are determined using least squares. However, the premise of this work is that the b coefficients are related to the spatial organization of the classes in the landscape. As such, the Stage I models (described below) are used to estimate b0 and b1 using the spatial variables described above. Once b0 and b1 are estimated, Eq. (8) can be rewritten as P 2b Pi305 i1020 0 b1

(5)

The fourth measure is the natural log of contrast [ln(c)], where c is an indication of the amount of local class-type variation present in a block defined as c5 o

Table 1. Cross-Scale Correlations of Independent Variables

(3)

5

207

(9)

to estimate class proportions at 30 m resolution. The b coefficients can also be estimated using spatial measures determined at 1020 m resolution. In this case the model is said to be inverted, where Pi30 is estimated using only information determined at 1020 m resolution. Model inversion will be tested below. Table 2. Correlations among Independent Variables

shn ln(c) ent rng var

shn

ln(c)

ent

rng

var

1.0 — — — —

0.33 1.0 — — —

0.91 0.34 1.0 — —

20.79 20.28 20.83 1.0 —

0.43 0.27 0.45 20.36 1.0

208

Moody

Stage I Models The Stage I models include a slope model and an intercept model. The slope model estimates the slope coefficient for the Stage II model using the previously described spatial measures as bˆ 15aˆ 01aˆ 1·rng1aˆ 2·ent1aˆ 3·var1aˆ 4·ln(c)1aˆ 5·shn. (10) Similarly, the intercept model estimates the intercept coefficient for the Stage II model as bˆ 05aˆ 01aˆ 1·rng1aˆ 2·ent1aˆ 3·var1aˆ 4·ln(c)1aˆ 5·shn, (11) where the estimated a-coefficients are not intended to be the same between the two models (10 and 11) but rather would differ depending on whether it is bˆ 0 or bˆ 1 that is being predicted. These models are developed on the training data for which least squares estimates of bˆ 0 and bˆ 1 are available based on measured values of Pi30 and Pi1020 for each sampled block. That is, for each of the 55 sampled blocks in the training site, a value for b0 and b1 is determined by regressing Pi30 on Pi1020. The estimated b coefficients are then used as the dependent variables in Equations (10) and (11). Parameter estimates for Models 10 and 11 are thus each based on an n of 55. As indicated above, the Stage I models are also run using the independent variables as determined at 1020 m resolution to test model inversion. Model Development, Validation, and Inversion The overall modeling procedure is as follows: First, the Stage I models are used to estimate the slope and intercept parameters needed for the Stage II model. The Stage I models are developed as described above using the spatial data and proportion data measured within the 55 blocks sampled on the training site. Second, using the slope and intercept parameters predicted from the Stage I models, the Stage II model is used to predict class proportions at 30 m based on measured class proportions at 1020 m. This second step is performed once for the training site as an initial evaluation, and again for the validation site as an independent evaluation of the area correction procedure. Third, model inversion for the validation site proceeds in the same fashion as the forward modeling described above, except that the values for the spatial measures are determined from the 1020 m data rather than the 30 m data. In the inverted mode the procedure requires only coarse resolution data. If successful, inversion of the procedure will obviate the need for high resolution training or calibration data, thus making the method more practical and operational for global scale land-cover mapping. Model Evaluation Model performance is evaluated in several ways. These methods focus on the error within sampling blocks, rather than for the overall sites because it is the ability of the procedure to improve area estimates within blocks that is of greatest interest here. Also, it is presumed that

if area estimates are significantly improved within a majority of individual blocks, then region-wide area estimates will also be improved. First, 1020 m proportions for each class in each of the sampled blocks are plotted against the 30 m proportions both pre- and postcorrection. If the modeling procedure works, then the points on these plots will fall closer to the zero–one line after the Stage II model is applied (e.g., see Figs. 4a and b). Second, total error within each block g is determined as k

TEg5 o |Pi10202Pi30|. i 51

(12)

For each block, TEg can be determined before and after correction. These values can be plotted against each other to visualize how the correction procedure is affecting the error within the blocks. This approach also allows a t-test to determine whether any reductions in mean total proportion error are statistically significant (e.g., Figs. 4c and d). Finally, individual sample blocks corresponding to cases where the modeling procedure performed either extremely well or extremely poorly are examined. This provides insight regarding the sensitivity of the area correction procedure to specific surface characteristics. RESULTS Figure 1 illustrates the general relationship between measured class proportions at the two resolutions of interest for the entire set of sample cells. Notice that the smoothed fit through these data suggests a nearly linear relationship between class proportions at the two resolutions. Linearity is violated near the origin where class proportions are low. It is also evident by comparing the smoothed fit line with the zero–one line that proportions for large classes are overestimated at 1020 m while proportions for small classes are underestimated. The crossover point is near 0.35. If no error resulted from estimating the class proportions at the two different scales, all of the points would fall along the zero–one line. As such, a successful correction procedure should result in a new graph where the points all fall closer to the zero–one line. Because of the basic nature of the proportion scaling behavior, the slope of the relationship between fine- and coarse-scale proportions will always be greater than 1. The closer the slope value is to 1, the less error will result when estimating class proportions at 1020 m. Similarly, the intercept value will always be less than 0. Note that the Stage II model involves the subtraction of the intercept from the coarse-resolution area estimate [Eq. (9)]. Hence, corrected area estimates cannot take on negative values. Regression summaries for the Stage I models are presented in Table 3. The five spatial variables are significant at the 0.05 level for both the slope model and

Improved Area Estimates from Remote Sensing

Table 3. Regression Summary for Stage I Modelsa

Slope model b0 shn ln(c) ent rng var

Coefficient

Standard Error

t-Value

P.|t|

20.135 20.411 20.076 0.819 0.897 20.001

0.212 0.145 0.025 0.127 0.173 0.0003

20.638 22.843 23.038 6.426 5.201 23.357

0.527 0.007 0.004 0.000 0.000 0.002

0.051 0.014 0.003 20.027 20.031 0.000

0.007 0.005 0.001 0.004 0.006 0.000

7.601 3.026 3.554 26.696 25.589 3.844

0.000 0.004 0.001 0.000 0.000 0.0003

Intercept model b0 shn ln(c) ent rng var a

Slope model R2adj50.57. Intercept model R2adj50.61.

the intercept model. The adjusted-R2 values are 0.57 and 0.61, respectively. Interpretation of slope coefficients in a multiple regression can be problematic due to variable interactions and multicollinearity. In this case, the signs of the model coefficients correspond to the signs of simple correlations between the spatial variables and the dependent variable (slope of the transition line). It is thus reasonable to interpret the signs of the coefficients from the Stage I models as indicative of the basic relationships between the independent variables and the b1 values. To interpret the coefficients for the variables, it is helpful to keep in mind that as the slope value decreases toward 1, the class proportion estimates become more stable, or change less, across scales. As an example, the negative slope coefficient for shn suggests that as the number of classes within a sample block increases, and/or as the evenness with which those classes are distributed within the block increases, the slope of the scale transition line reduces toward 1.0. As such there will be less error within blocks whose Shannon index is high. Figure 3 is a close look at two Plumas sample blocks, one for which shn is low (0.34) (block 54, left) and one for which shn is high (1.22) (block 28, right). For each case, 30 m NDVI data, 30 m class maps, and 1020 m class maps are shown. In block 54 the total error is relatively large (36%). In block 28, the total error is relatively low (26%). The summary of the Stage I slope model (Table 3) suggests that the stability of class proportions across scales (or their tendency to remain unchanged) increases with richness and evenness (shn), within-patch cohesiveness [ln(c)], and spectral variance (var). Stability decreases with between-patch disaggregation (ent) and range of class sizes (rng). The results for the overall area correction procedure for the training data (Plumas) are presented in Figure 4. Figure 4a shows the relationship between 30 m and 1020

209

m proportions before correction. Figure 4b shows the same relationship after the area-correction model is applied to the 1020 m proportion estimates. As desired, the area estimates for the smaller classes (below about 0.4 on the x-axis) are brought up toward the line and the estimates for the larger classes (above 0.4) are brought down toward the line. In Figures 4a and 4b, each point represents one class within a sample block. In Figure 4c, each point represents the total error [TE, Eq. (12)] within each block. Blocks whose post-correction TE is greater than their precorrection TE are indicated by points falling above the zero–one line. Blocks whose postcorrection TE value is less than their precorrection TE value are indicated by points falling below the zero–one line. Notice that a) most points fall below the line, indicating that an improvement in area estimates has been achieved within the sample blocks, and b) the model has a tendency to perform poorly for blocks within which there is little total error initially. This effect is illustrated in Figure 5 which shows NDVI, 30 m class data, and 1020 m class data for Plumas sample blocks #13 (left) and #3 (right). For block 13, the total error is quite high (52%) and is considerably reduced (to 32%) using the correction procedure. For block 3, the total error is low (8%) and is actually made worse (18%) by the correction procedure. The points corresponding to these blocks are indicated in Figure 4c. Boxplots are used (Fig. 4d) to compare the distribution of the total error values before and after correction. The upper and lower boundary of each box represents the interquartile distance and the midline is the median value. Based on these error distributions, the mean value of the postcorrection errors (0.18) is significantly lower than the mean error before correction (0.24) (t52.68; P.t50.009). Comparison of model results is presented in Table 4. The results for the area correction model applied to the validation data (Stanislaus) are presented in Figure 6. As seen in Figure 6a, the relationship between 30 m proportions and uncorrected 1020 m proportions follows the expected pattern and exhibits the same basic behavior as the training site. Figure 6b shows the same relationship after the 1020 m proportions have been corrected and, as desired, the values above about 0.4 are pulled down toward the zero–one line, and the values below 0.4 are pulled up toward the line. The withinblock total error (TE) values before and after correction are compared in Figure 6c. Again, most of the points fall below the zero–one line, indicating that the total error within most blocks has been reduced. As with the training data, the model works better for sample blocks that have a large amount of error than it does for cells that have little error initially. Figure 7 shows the NDVI data, 30 m class data, and 1020 m class data for Stanislaus sample blocks #4 (left) for which the model performed well, and #18 (right) for which the model performed

210

Moody

Figure 3. Comparison of two sample cells from the Plumas with different shn values. The panels on the left (top to bottom) are NDVI data, 30 m proportions, and 1020 m proportions for cell #54 which has a low shn value (0.34). The panels on the right are the same data fields for cell #28 which has a high shn value (1.22). The total error in cell #54 is 36%. The total error in cell #28 is 26%. The locations of these points are identified in Figure 4c. From lightest to darkest, the classes are: barren, brush, hardwood, conifer, and water.

Improved Area Estimates from Remote Sensing

211

Figure 4. Model results for the training data. Relationship between 30 m proportions and 1020 m proportions: a) precorrection; b) postcorrection; c) relationship between total-error values before and after correction; d) comparison of total error distributions pre- and postcorrection. The points identified in 4c correspond to the sample cells shown in Figures 3 and 5.

poorly. For block 4 the original total error is 43% and is reduced to 15%. For block 18 the original total error is 18% and is increased to 21%. The points for these blocks are indicated on Figure 6c. Figure 6d compares the error distributions before and after correction. Again, the mean total error value for the postcorrection data (0.22) is significantly lower than the mean value for the precorrection data (0.31) (t54.20; P.t50.000). Finally, the procedure is run in the inverted mode for the validation data. These results are presented in Figure 8. The same general patterns are observed as for the two cases described above. The inverted correction procedure does not work as well for the validation data as the for-

ward method does; however, it still produces a significant reduction in mean total error (from 0.31 to 0.24) as illustrated in Figures 8c and d (t53.13; P.t50.003). A strong intercept effect is evident in all three model applications. This effect can cause classes that are not actually present in a cell to be introduced in proportions equivalent to the magnitude of the intercept used to correct them. This effect is seen in Figures 4b, 6b, and 8b, where the intercept produces a stacking of points above the 0-value on the x-axis. This is problematic. However, when the intercepts are constrained to zero, the model performs considerably worse in all cases. Although it is not tested here, it is possible to include a

212

Moody

Figure 5. Comparison of two sample cells from the Plumas with different total-error (TE) values. The panels on the left are NDVI data, 30 m proportions, and 1020 m proportions for cell #13 for which TE is quite high (52%) and is reduced by the correction procedure (to 32%). The panels on the right are the same data fields for cell #3 for which TE is low (8%) and is made worse by the correction procedure (to 18%). The locations of these points are identified in Figure 4c. The greyscale scheme is the same as for Figure 3.

Improved Area Estimates from Remote Sensing

213

Table 4. Comparison of Mean Pre- and Postcorrection Total-Error Values for the Three Modelsa

Forward: PL Forward: ST Inverted: ST Inverse: ST Inverse: ST.NO

Precorrected

Postcorrected

t-Value

P.|t|

Cells Improved (%)

0.24 0.31 0.31 0.31 0.31

0.18 0.22 0.24 0.26 0.21

2.68 4.20 3.13 1.12 3.37

0.009 0.000 0.003 0.270 0.001

80 91 89 77 83

a The results from the inverse procedure determined both with and without (NO) outliers are included for comparison.

Figure 6. Forward model results for the validation data. Relationship between 30 m proportions and 1020 m proportions: a) precorrection; b) postcorrection; c) relationship between total-error values before and after correction; d) comparison of total error distributions pre- and postcorrection. The points identified in 6c correspond to the sample cells shown in Figure 7.

214

Moody

Figure 7. Comparison of two sample cells from the Stanislaus for which the model performed well (left, #4) and poorly (right, #18). The points for these cells are identified in Figure 6c. The greyscale scheme is the same as for Figure 3.

Improved Area Estimates from Remote Sensing

215

Figure 8. Inverted model results for the validation data. Relationship between 30 m proportions and 1020 m proportions: a) precorrection; b) postcorrection; c) relationship between total-error values before and after correction; d) comparison of total error distributions pre- and postcorrection.

decision rule that will not allow a class to be introduced within a sample block if its initial estimated proportion at 1020 m is zero. This approach, of course, would fail to correct for classes that are present in small quantities but disappear at coarse scales of observation. Since several authors have shown that the inverse method [Eq. (1)] outperforms the direct method of area correction used here (Czaplewski and Catts, 1992; Walsh and Burk, 1993; Moody and Woodcock, 1996), a comparison is carried out to determine if this is still the case when spatial information is incorporated into the direct method. Figure 9 illustrates the results when the inverse method is applied to the 35 sample blocks from the

Stanislaus site using a transition matrix that is calibrated on the 55 sample blocks from the Plumas site [see Eq. (1)]. These results should be compared to the results from the forward application of the model to the Stanislaus data presented in Figure 6. Figure 9b exhibits the stacking of points above the zero value on the x-axis which is analogous to the intercept effect in the direct method. While the inverse method works well overall, numerous points exhibit large over or underestimations. When total error is determined for each sample block there are three blocks whose total error values are quite large (near 1.0). The inverse method performs considerably worse than the direct method and produces no sig-

216

Moody

Figure 9. Inverse method results for the validation data. Relationship between 30 m proportions and 1020 m proportions; a) precorrection; b) postcorrection; c) relationship between total-error values before and after correction; d) comparison of total error distributions pre- and postcorrection.

nificant overall improvement in area estimates (t51.12; P.t50.27) when the outliers are included in the test. When the two largest outliers are removed, a significant error reduction is achieved (t53.37; P.t50.002). See Table 4 for a comparison of the performance for all the models. Figure 10 shows class-specific proportions for six cases representing selected sample blocks. Figures 10a, b, and c represent cases where the model performs poorly for the training data, validation, and validation in the inverted mode, respectively. Notice that, in Figures 10a, b, and c, there is either very little difference between the postcorrection and precorrection values or the

postcorrection values actually fall further away from the 30 m lines than the precorrection data. In these cases, the errors for the individual classes were small initially. Figures 10d, e, and f represent cases where the model performs well. Notice that in these figures, the postcorrection values typically fall closer to the 30 m values than do the precorrection data. In these cases, the total error for each block was large initially. DISCUSSION The two-stage modeling procedure provides significantly improved estimates of land-cover area in both the for-

Figure 10. Class specific proportions for six sample cells representing select cases. Figures 10a, b, and c represent cases where the model performs poorly for the training data, validation, and validation in the inverted mode, respectively. Figures 10d, e, and f represent cases where the model performs well.

Improved Area Estimates from Remote Sensing

217

218

Moody

ward and inverted modes for the training and validation sites considered here. However, several key issues remain to be discussed. These include the limitations, assumptions and weaknesses of the model and analysis. First, it is assumed that 30 m and 1020 m class area estimates are linearly related. This appears to be true within most of the range of class proportions considered in this analysis. However, since class proportions necessarily range between 0 and 1, linearity will tend to be violated in cases where proportions are either very small or very large. Consequently, the linear model will probably not work as well where very few classes are involved. For example, if the goal was to correct area estimates for forested and nonforested area (a two-class case), linear assumptions would be violated for sample blocks that were either highly forest-dominated or highly non-forestdominated. For the data presented here, examination of Figure 1 shows that the relationship between measured proportions at 30 m and 1020 m resolution is weak for classes that make up less than about 20% of a sample block (i.e., points falling to the left of 0.2 on the x-axis in Fig. 1). In fact, it is for these points that the modeling procedure performs worst in both training and validation cases. To account for these problems, nonlinear or piecewise functions might be more appropriate. However, depending on the specific approach, these models may be more complex and require estimates for a greater number of parameters, thereby giving up the simplicity of the approach presented in this paper. Modification of the basic approach presented here might also include structural features that constrain corrected proportions to sum to 1. In this analysis it is assumed that spatially averaged spectral values from Thematic Mapper data are, for the present purposes, reasonable approximations of coarse spatial resolution data that might come from a sensor with similar spectral characteristics. Therefore, the spectral mixing of the land-cover classes considered is not simulated realistically. It is unclear how this would affect a) the apparent stability of area estimates across resolutions and b) the ability of the model to correct error in these area estimates. The convolution of spectral and spatial effects needs to be considered in terms of its impact on area estimates across resolutions. Although the training and validation landscapes used are spatially separated, they do represent two similar environments with the same classes present and the same types of environmental controls on the spatial distribution of those classes. While this provides a valuable firstcase test of the model, it remains to be seen how stable the modeling procedure is in landscapes with a different number of component classes and different spatial (as well as spectral) class distributions. Further, it is not known whether the set of spatial measures used here will provide the same power for model training in other landscapes. The sensitivities of the modeling procedure to a)

Table 5. Regression Summary for Stage I Slope Model Applied to the Validation Dataa Slope Model

Coefficient

Standard Error

t-Value

P.|t|

b0 shn ln(c) ent rng var

20.200 20.646 20.122 0.998 0.874 20.0004

0.645 0.211 0.039 0.283 0.503 0.001

20.311 23.066 23.124 3.526 1.7383 20.569

0.758 0.006 0.004 0.001 0.093 0.573

a

Slope model R2adj50.42.

the spatial extent of the area over which area estimates are to be corrected, and b) the size of the blocks used to train the model also need to be evaluated. The Stage I models are based on 55 blocks of data from the Plumas National Forest. Several of the predictor variables have fairly high correlations and there may be spatial autocorrelation among the sampled blocks. Furthermore, since Eqs. (9) and (10) each require the estimation of six coefficients, it is reasonable to question whether the Stage I models are overfit to the training data. Since the 35 cells from the Stanislaus Forest are reserved for validation of the overall model, they are obviously not available for training. However, the Stage I model parameters were reestimated using data from the Stanislaus National Forest to see whether the coefficients appear to have been overfit to the training site. The results of the slope model are presented in Table 5. Although the larger standard errors for rng and var prohibit these two variables from being significant in the new model, the signs and relative magnitudes of the coefficients are the same as for the training data (compare to Table 3). The similarities between the two sets of parameter estimates suggest that the Stage I model is probably not overfit to the training data. There are two concerns regarding the spatial measures used. First, there is only empirical justification for the use of this particular set of measures. Although they make sense, and behave as expected in the model, it would be valuable to have either a theoretical justification, or a better developed empirical justification, for these particular measures. In this regard, simulation studies that are designed to evaluate the independent, as well as interactive effects of different spatial characteristics on scale-dependent error, would be valuable. Further, in order to be useful for model inversion, it is necessary to find a set of measures that are strongly and consistently (i.e., in various landscapes) self-correlated across scales. It is not known whether the set of measures used here best fits this criterion. Finally, there is considerable uncertainty regarding the actual sensitivities of ecosystem and climate models to errors in land-cover area estimates. In cases where land-cover classes are used to parameterize such models,

Improved Area Estimates from Remote Sensing

a look-up table approach is typically employed to translate cover classes into either mean biophysical parameter values or to call up parameter distribution functions for each class type. An assumption made in this article is that bias in area estimates at the subgrid scale of ecosystem models can lead to biases in model parameterization and consequently in model outputs. Further, it is assumed that any biases introduced by errors in land-cover area will lead to model output errors that are significant with respect to the overall noise and error in the model as a whole. Simulation studies which run ecosystem models under different levels of bias in land-cover area would help quantify and better understand these potential effects. CONCLUSION Earlier work demonstrated that large proportion errors arise as land-cover areas are estimated at progressively coarser scales (Moody and Woodcock, 1994). In subsequent analyses, these errors were shown to be driven largely by a few specific spatial characteristics of the landscape. Specifically, these are: a) the size of the largest and smallest classes in an area; b) the number of classes in the area; c) the evenness with which those classes are distributed; and d) the cohesiveness with which those classes are distributed within the area (Moody and Woodcock, 1995; 1996). The model presented here shows that, in the landscapes considered, spatial measures that characterize these basic properties are useful in predicting the rate at which errors in area estimates will arise in transitioning between 30 m and 1020 m class maps derived from multispectral satellite data. Further, the parameters that indicate this rate of change in class size can be used to improve coarse-resolution area estimates in a postclassification mode. Although there are several unresolved issues, the model results are encouraging for the prospect of reducing bias in area estimates that are derived from coarseresolution remotely sensed data. While the methods presented here have the advantage of simplicity, they rely on a set of assumptions that need to be more thoroughly explored. In particular, the methods and assumptions should be tested under a wider variety of landscape scenarios in order to evaluate their operational utility. Similarly well focused sensitivity studies based on simulated data can go a long distance in improving understanding of the basic effects observed in this analysis. These activities are part of ongoing research. This work was supported by the National Aeronautics and Space Administration under Contract NAGW-5158. I extend my thanks to Ray Czaplewski and two other anonymous reviewers for their helpful critique of an earlier draft of this manuscript.

219

REFERENCES Adams, J., Smith, M., and Johnson, P. (1986), Spectral mixture modeling: A new analysis of rock and soil types at the Viking Lander 1 site. J. Geophys. Res. 91(B8):8098–8112. Avissar, R. (1992), Conceptual aspects of a statistical-dynamical approach to represent landscape subgrid-scale heterogeneities in atmospheric models. J. Geophys. Res. 97(D3):2729–2742. Baker, W. L., and Cai, Y. (1992), The r.le programs for multiscale analysis of landscape structure using the GRASS geographical information system. Landscape Ecol. 7(4): 291–302. Borel, C., and Gerstl, S. (1994), Nonlinear spectral mixing models for vegetative and soil surfaces. Remote Sens. Environ. 47:403–416. Brown, P. (1982), Multivariate calibration. J. Roy. Stat. Soc. B 44:287–321. Czaplewski, R., and Catts, G. (1992), Calibration of remotely sensed proportion or area estimates for misclassification error. Remote Sens. Environ. 39:29–43. Dickinson, R. (1995), Land processes in climate models. Remote Sens. Environ. 51:27–38. Foley, J. (1994), Net primary productivity in the terrestrial biosphere: the application of a global model. J. Geophys. Res. 99(D10):20,773–20,783. Food and Agriculture Organization (FAO) (1996), Forest Resources Assessment 1990—Survey of Tropical Forest Cover and Study of Change Processes. Forestry Paper No. 130, Rome, 154 pp. Henderson-Sellers, A., Wilson, M., and Thomas, G. (1985), The effect of spatial resolution on archives of land cover type. Clim. Change 7:391–402. Huete, A. (1986), Separation of soil-plant mixtures by factor analysis. Remote Sens. Environ. 19:237–251. Jupp, D. L. B. (1989), The stability of global estimates from confusion matrices. Int. J. Remote Sens. 10:1563–1569. Kalkhan, M., Reich, R., and Czaplewski, R. (1995), Evaluation of statistical properties of the inverse estimator for remotely sense areal estimates using simple random sampling. In Proc. Am. Soc. Photogramm. and Remote Sensing Conf., 27 Feb.–2 Mar., Charlotte, NC, Vol. 3, pp. 258–270. Koster, R., and Suarez, M. (1992), Modeling the land surface boundary in climate models as a composite of independent vegetation stands. J. Geophys. Res. 97(D3):2697–2715. Mayaux, P., and Lambin, E. (1995), Estimation of tropical forest area from coarse spatial resolution data: a two step correction function for proportional errors. Remote Sens. Environ. 53:1–16. Moody, A., and Strahler, A. H. (1994), Characteristics of composited AVHRR data and problems in their classification. Int. J. Remote Sens. 15(17):3473–3491. Moody, A., and Woodcock, C. E. (1994), Scale-dependent errors in the estimation of land-cover proportions—implications for global land-cover datasets. Photogramm. Eng. Remote Sens. 60(5):585–594. Moody, A., and Woodcock, C. E. (1995), The influence of scale and the spatial characteristics of landscapes on land-cover mapping using remote sensing. Landscape Ecol. 10(6): 363–379. Moody, A., and Woodcock, C. E. (1996), Calibration-based mod-

220

Moody

els for correction of area estimates derived from coarse resolution land-cover data. Remote Sens. Environ. 58:225–241. Oleson, K., Sarlin, S., Garrison, J., Smith, S., Privette, J., and Emery, W. (1995), Unmixing multiple land-cover type reflectances from coarse spatial resolution satellite data. Remote Sens. Environ. 54:98–112. Pax-Lenney, M., and Woodcock, C. E. (1996), The effect of spatial resolution on the ability to monitor the status of agricultural lands. Remote Sens. Environ. 61:210–220. Puyou-Lascassies, P., Flouzat, G., Gay, M., and Vignolles, C. (1994), Validation of the use of multiple linear regression as a tool for unmixing coarse spatial resolution images. Remote Sens. Environ. 49:155–166. Quarmby, N., Townshend, J., Settle, J., and White, K. (1992), Linear mixture modeling applied to AVHRR data for crop area estimation. Int. J. Remote Sens. 13:415–425. Roberts, D., Smith, M., and Adams, J. (1993), Green vegetation, nonphotosynthetic vegetation, and soils in AVIRIS data. Remote Sens. Environ. 44:222–269. Roesch, F., Van Deusen, P., and Zhu, Z. (1995), A comparison of various estimators for updating forest area coverage using AVHRR and forest inventory data. Photogramm. Eng. Remote Sens. 61(3):307–311. Running, S., Justice, C., Hall, D., et al. (1994), Terrestrial remote sensing science and algorithms planned for EOS/MODIS. Int. J. Remote Sens. 15:3587–3620. Sellers, P., Mintz, Y., Sud, Y., and Dalcher, A. (1986), A simple biosphere model (SiB) for use within general circulation models. J. Atmos. Sci. 43(6):505–531. Seth, A., Giorgio, F., and Dickinson, R. (1994), Simulating fluxes from heterogeneous land surfaces: explicit subgrid

method employing the biosphere–atmosphere transfer scheme (BATS). J. Geophys. Res. 99(D9):18,651–18,667. Starmer, W. J., and Frohn, R. C. (1996), The effects of changing spatial resolution and spatial pattern on IGBP landcover classes in the Washington, D.C. Area. In Proc. Am. Soc. Photogramm. and Remote Sens. Conf. 21–25 April, Baltimore, MD, Vol. 1, pp. 169–178. Strahler, A., Moody, A., Lambin, E., et al. (1995), MODIS land cover product: algorithm technical basis document, NASA Document for MODIS Product No. 11, Parameter Nos. 2669 and 2671, NASA, Greenbelt, MD. Townshend, J., Justice, C., Li, W., Gurney, C., and McManus, J. (1991), Global land cover classification by remote sensing: present capabilities and future possibilities. Remote Sens. Environ. 35:243–255. Turner, M., O’Neill, R., Gardner, R., and Milne, B. (1989), Effects of changing spatial scale on the analysis of landscape pattern. Landscape Ecol. 3(3/4):153–162. Vitousek, P. (1994), Beyond global warming: ecology and global change. Ecology 75(7):1861–1876. Walsh, T., and Burk, T. (1993), Calibration of satellite classifications of land area. Remote Sens. Environ. 46:281–290. Woodcock, C. E., and Gopal, S. (1992), Accuracy assessment of the Stanislaus Forest vegetation map using fuzzy sets. In Proc. Fourth Biennial Remote Sens. Appl. Conf., Orlando, FL, April, pp. 378–394. Woodcock, C. E., Gopal, S., Macomber, S. A., and Jakabhazy, V. D. (1994), Accuracy assessment of the vegetation map of the Plumas National Forest, Technical Report, Center for Remote Sensing, Boston University, Boston, MA, 19 pp. Yuan, D. (1996), Natural constraints for inverse area estimate corrections. Photogramm. Eng. Remote Sens. 62(4):413–417.