Landsat-7
ETM+
images. We selected imagery from a three-year
window surrounding the nominal year. Within this three-year
window, we identified three images covering different pheno-
logical stages (Sexton
et al
., 2013): a leaf-off image (November-
February), an early growing season image (March–June), and a
late growing season image (July–October). We only used imag-
ery with no more than 10% cloud cover over land areas and no
more than 10% snow cover, as identified by the Fmask-based
pixel
QA
layer (Zhu
et al
., 2015; Zhu and Woodcock, 2012). For
classifications from the 1986-1996 period, there were ten cases
where we could not find images of a given nominal classifica-
tion year that met these criteria in one of the phenological
periods (typically during the leaf-off stage), so only two images
were used for these classifications. All images were reprojected
from their native
UTM
projections to the Albers Equal Area pro-
jection used by the
NLCD
(Homer et al., 2015) and resampled to
30 m resolution using a bilinear method. We converted surface
reflectance of each image to the tasseled cap brightness, green-
ness, and wetness indices using coefficients provided by Crist
and Cicone (1984). While the spectral response of the
ETM+
sensor differs from the
TM
sensor, using the same tasseled cap
coefficients for both sensors is not problematic since
AASG
em-
pirically adapts the spectral signatures of each land cover class
to the properties of each image. The multi-season tasseled cap
indices were used as predictors for the classification along
with the brightness temperature from each image, with stable
sites identified based on image difference histograms derived
from the tasseled cap brightness index during the growing sea-
son. We also included a single binary water layer (based on the
pixel
QA
layer) as an additional predictor in the classifications.
We masked pixels that were identified as cloud-covered
(with medium or high confidence) in the pixel
QA
layer. To
avoid missing pixels in the final classified images, we filled
these masked pixels using two approaches designed to capital-
ize on close correspondence between reflectance of a given pix-
el and its spatial and temporal neighbors (Figure S1). First, we
fit ordinary least squares linear regression models between each
tasseled cap index of a given image and the same layer of the
other image(s) of the same nominal year, and we used this fitted
regression relationship to fill cloud-contaminated pixels (Figure
S1). If two images were available for fitting in a given year, we
selected the image that provided the best fit (highest R
2
). Any
remaining cloud-contaminated pixels (i.e., due to overlapping
clouds between two image dates) were filled using the lower-
performing regression model. Second, in very rare cases where
cloud-masked pixels were still present in the image, we used
two low-pass moving window filters to estimate missing values,
first a 3-by-3 pixel moving window then a 7-by-7 pixel moving
window for any remaining unfilled pixels (Figure S1).
Topographic Data
In addition to multi-season Landsat imagery, we also included
four topographic metrics derived from the Shuttle Radar Topog-
raphy Mission digital elevation model as predictors in the clas-
sifications. These included elevation, slope, upslope accumu-
lated area (
UAA
), and topographic wetness index (
TWI
).
UAA
and
TWI
represent the upslope surface area that could contribute
runoff to a given point and the tendency of water to accumulate
in a given area, respectively (Beven and Kirkby, 1979). They
are particularly useful for distinguishing wetlands from other
vegetated land cover classes (Dannenberg
et al
., 2016). Slope,
UAA
, and
TWI
were all calculated using the multiple flow direc-
tion D-Infinity method in
TauDEM
(Tarboton, 1997), after which
all topographic data were reprojected to an Albers Equal Area
projection with 30 meter spatial resolution. Missing values for
each topographic metric were filled using the low-pass moving-
window filters described in the previous section, with remain-
ing missing values filled with the median value of the scene.
There were two cases where large areas in northern Florida
and along the border between North and South Carolina were
missing from the
SRTM
elevation dataset (Figure S2), which we
had to fill using the median value. However, these areas are
quite flat, so the regional median value is a reasonable approxi-
mation. Additionally, the nonparametric machine learning
classifier used here should treat the topographic data as being
non-informative predictors in these missing regions and would
therefore give higher weight to the remotely sensed imagery.
Post-Classification Processing and Assessment
Images of different dates from the same path and row do not
have precisely the same spatial extents, and differences in the
image extents of the multi-season imagery resulted in clear
edge effects on the margins of each classified scene where
the images do not perfectly overlap. We therefore visually
trimmed each classified image to remove these outer mar-
gins. For consistency with the
NLCD
(Homer
et al
., 2015), we
filtered each classified image using class-specific minimum
mapping units (
MMUs
) to reduce “salt-and-pepper” effects and
to improve accuracy (Figure S3). These included twelve pixel
MMUs
for the two agricultural classes and five pixel
MMUs
for all remaining classes except the three Developed classes,
which were not filtered. Pixels within clusters below these
MMU
thresholds were assigned the class label of the nearest
unfiltered pixel using the “nibble” function in ArcGIS
®
10.5.1.
We assessed both the consistency of our
AASG
classification
products with the
NLCD
as well as the accuracy of the classifi-
cation products using independent reference data. To assess
consistency with the
NLCD
, which is of primary concern for
multi-temporal assessment of land cover change, we conduct-
ed a cross-comparison of our 2001 and 2006
AASG
classifica-
tions to the 2001 and 2006
NLCD
classifications (which were
not used in the training process). We compared these classi-
fication products on both an areal and a pixel-by-pixel basis,
and assessed agreement between the
NLCD
and
AASG
classifica-
tions using overall agreement (
OA*
), “producer’s” agreement
(
PA*
), and “user’s” agreement (
UA*
), with
NLCD
as the refer-
ence. Here, we use the term “agreement” rather than “accu-
racy” to reflect the fact that the
NLCD
is not a true reference
dataset and we are therefore only assessing consistency, not
accuracy, between two land cover datasets that both contain
error. To distinguish the agreement metrics from the accuracy
metrics, we denote the agreement metrics with an asterisk (*).
In addition to assessing the consistency of
AASG
with
NLCD
, we also assessed the accuracy of the
AASG
classifica-
tions using best practice recommendations (Olofsson
et al
.,
2014). Since many of the classes used in the Anderson Level
2 classification scheme are not easily distinguished visually,
we assessed classification accuracy using a simplified seven
class system: water, developed (which includes all developed
intensities in a single class), barren, woody vegetation (which
combines the three forest classes, shrub/scrub, and woody
wetland), herbaceous vegetation (including both grassland/
herbaceous and hay/pasture), cultivated crops, and emergent
herbaceous wetland. We sampled a total of roughly 550 pix-
els, which yields an estimated standard error of approximate-
ly 1.5% on the final overall accuracy (Olofsson
et al
., 2014).
We used a stratified random sampling design to select pixels
for assessment, with a mixed sample allocation that achieves
a minimum allocation of 50 pixels per class and then allo-
cates the remaining pixels proportionally to the prevalence
of each class (Table S1). The sample allocation and selection
was defined based on the 2006
AASG
classification, with the
same points used for the remaining assessments.
We assigned reference labels for these sample points using
a combination of fine resolution aerial imagery and Landsat
imagery in Google Earth Pro
™
. Across the study region, there
was very little fine-resolution imagery available prior to 1993,
so we were only able to assess the 1996, 2001, and 2006
562
September 2018
PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING
