Change Detection Methods
and the Boundary Waters Canoe Area
Thomas Juntunen
For FR5262 – Fall Semester, 2010
Objectives
As far back as 1975, it was predicted that "even low-resolution data from the Landsat MSS scanner, if combined and enhanced, will disclose 80 to 90% of the exchanges of land use between forest and non-forest categories. In addition, such data will show 25 to 90% of the less distinct features in the forest, depending on category." (Aldrich, 1975, after Coppin & Bauer, 1996). Coppin and Bauer (1994) eventually studied whether any of the prevailing change detection algorithms -- at least those as of the late '90s -- were capable of reliably detecting relatively subtle changes in land cover, particularly the sort natural resource managers would be interested in. To that end they compared two change detection methods, with a variety of generated change features, on Landsat Thematic Mapper data of an area in northern Minnesota, over three different time periods. They concluded their methodology would be sufficient as a first pass to let managers safely ignore stands classified as unchanged, but better methods and sensors would be required before anything more useful could be accomplished.
I wished to compare the basic change detection methods covered in class against more statistically oriented methods, such as those inventoried in Coppin and Bauer (1996), to see if spectrally subtle differences -- such as revegetation after wildfires -- could be reliably detected. I gathered Landsat 5 Thematic Mapper scenes from two periods, four years apart, of areas in the BWCA that underwent two major wildfires during the intervening time. These areas were selected to include large-scale, dramatic change events that should be highly noticeable in the data. The satellite data were gathered and prepared as much as possible according to the methodologies laid out in Coppin and Bauer (1994).
Materials, Tools and Concepts
Satellite imagery – I tried to acquire Landsat 7 ETM+ imagery, but the post-2003 Scan Line Module (SLM) problem made many scenes useless, and many of the scenes in my chosen time frame contained far too much cloud cover. I was able to acquire suitable scenes from Landsat 5 TM imagery, albeit with non-zero but under 10 percent cloud cover. Rows, paths and other metadata are listed in Table 1 in the appendix.
Digital aerial photography – My original plan included performing an accuracy assessment using NAIP data (and indeed, I acquired county-wide mosaics for Cook and Lake Counties in Minnesota, for 2003 and 2008). While I did not perform such an assessment, the NAIP images, which are created as 3.75-minute quadrangles, defined the area of interest, which consisted of 15 such quadrangles combined into two distinct polygons. The quadrangle metadata is listed in Table 2 in the appendix.
Boundary polygons of two wildfires were obtained from the USGS Monitored Trends in Burn Severity (MTBS) database for the Cavity Lake fire (FS-0909-043-20060714), which started on 2006-07-14 and affected more than 26,000 acres of the BWCA; and the Famine Lake fire (FS-0909-089-20060908), which started on 2006-09-08 and affected almost 4,000 acres. Both areas are entirely within the boundary of Cook County in Minnesota.
All image processing, including differencing, linear transformations and change detection was performed with ERDAS IMAGINE 2010. Cartographic output and some statistics were produced with ESRI ArcGIS Desktop 9.3.1. Line graphs were produced in Excel 2007.
Procedures
Data Acquisition, Radiometric and Geometric Rectification
Both classification of and change detection within remote sensing imagery require calibration and rectification steps prior to analysis, particularly when using multi-temporal data. If data sets in a time series are not in precide alignment and using the same spectral baseline, some of what's detected may simply be noise from error. It isn't possible to get perfect results, but a reasonably high confidence level is acheivable with care during pre-processing. I relied on the preprocessing EROS performs for Landsat data for radiometric calibration according to the sensor, for orthorectification of imagery to the terrain, as well as for georectification across scenes. As a check, the extents and world files were found to be identical:
West: -91.134419, East: -90.676035, North: 48.192138, South: 47.934734
Left: 639315.0, Right: 672735.0, Top: 5339355.0, Bottom: 5311635.0
30.00000000000, 0.0, 0.0, -30.00000000000, 639330.0000000000, 5339340.0000000000
Atmospheric Correction
Coppin and Bauer (1994) asserted the lack of sufficiently detailed atmospheric data for remote wilderness areas frequently left dark subtraction (of spectrally stable features from across the time series of images) as the most viable means of atmospheric correction. The images used in this study, located within a protected wilderness area, had few suitable features outside of several large bodies of deep water. Ten points were selected that were distributed within such lakes across both images and the averages computed. The images were then normalized by subtracting the mean water pixel values from all pixels. Table 3 in the appendix lists the values used.
Pixel-based methods versus Stands
It is worth noting here that Coppin and Bauer's 1994 study focused on detecting a range of forest canopy events that are of interest to natural resource managers. They noted that while most algorithms and techniques available were pixel based, managers operated in terms of “stands”. To accommodate the difference, classification was performed in advance of change detection. Between their reference data and the resolution of the available imagery, they were reduced to just four classes to meet change detection requirements of exhaustivity, mutual exclusion but relevant to resource management objectives: net canopy loss, net canopy gain, no change and “recent storm damage” (structural/textural canopy changes).
Anniversary Window and Data Intervals
Anniversary dates or windows are used to minimize differences in reflectance from seasonal and sun-angle changes. Coppin and Bauer (1994) believed mid-summer imagery at peak green was best for disturbance monitoring in northern Minnesota. They studied changes in two-, four- and six-year windows and found a two-year cycle was “optimal” to study storm damage (among other things) in the Great Lakes region (at least with Landsat TM imagery), and four- or six-year cycles were best for “human-induced and natural canopy disturbances such as thinning, cutting and dieback.” Revegetation after wildfires most resembles regrowth after clear-cutting, so I went with a minimum four-year interval for our study. There were no 2006 Landsat 5 scenes within our July-August anniversary window that were sufficiently cloud-free to be of use, so I eschewed a two-year interval as part of the study.
I chose a July-August window because Coppin and Bauer (1994) used the same window for the following reasons (quoted from their paper):
- Forest cover phenological stability
- Lowest monthly percent cloud cover, based on Landsat imagery available for [their] area of interest
- Lowest seasonal soil moisture content
- Minimal sun-angle effects (minimal variation in solar zenith angles between dates in window)
Data Enhancement for Interpretability
The thermal band (TM6) was eliminated at this stage for two primary reasons. First, Coppin and Bauer (1994) advised that “other investigators have shown that, for identification of surface types, thermal identification is not readily associated with that in the reflective part of the spectrum...”. Secondly, some tools in IMAGINE require that all bands have the same spatial resolution. The IMAGINE subset tool was used again to generate images with only six bands: TM1 through 5 and TM7.
Coppin and Bauer (1994) investigated seven vegetative indices“deemed promising by other investigators” because while the six reflective bands in TM data could be transformed into three (or four) dimensions in feature space, “no consensus can be found in the literature on which indexes or data transformations represent these features best in the context of green vegetation assessment.” For my purposes, the only vegetation index outputs were created for each of our datasets (henceforth, the 2004 and 2008 datasets): what Coppin and Bauer called “Crippen's NDVI”. Crippen's NDVI takes the form: TM4/(TM4 + TM3) which yields numbers from 0 to 1 instead of -1 to 1 for standard NDVI. The Raster > Scientific > Two Image functions were used to create the TM4 + TM3 output and then the final Crippen's output for both 2004 and 2008.
In the end, Coppin and Bauer (1994) generated 14 change features, half as standardized differences and half as second principal components, and applied them to their images across three time intervals in combinations of two up through six features at a time. From these results, they used Jeffries-Matusita (J-M) distances for best and mean separability distances among the change signatures. They found the top five change features (in terms of frequency) were the: standardized difference (STD) of brightness, second principal component (PC2) of greenness, PC2 of brightness, PC2 of the green ration (TM4/TM2) and STD of greenness. While six change features produced the largest average J-M distance, using only five or even four features decreased the Kappa values of the change detection process a very small amount. Accordingly, I compared the five features that produced the best J-M distance for a four-year interval. These are summarized in table 4 in the Appendix.
Within IMAGINE, the Tasseled Cap (TC) components were generated for 2004 and 2008 using the six-band subset images previously created, with the Raster > Spectral > Tasseled Cap command. This produced a six-component image that, according to ERDAS documentation, has: Brightness, Greenness, Wetness, Haze, Fifth and Sixth and the six respective output bands.
Principal Components images were generated by first making a two-layer stack of the appropriate components (TC Greenness, Red Bands, etc.) and then using Raster > Spectral > Principal Components, selecting two output components.
Standardized differences, as defined by Coppin and Bauer (1994), were created by using Raster > Scientific > Two Image functions to take the appropriate 2004 and 2008 components and either subtracting or adding them, and then once more with those two outputs, using division for the final result.
In all these operations, Ignore Zeros in Stats was checked and an AOI was defined, using the AOI file previously described.
Results
Coppin and Bauer (1994) generated 14 change features, half as standardized differences and half as second principal components, and applied them to their images across three time intervals in combinations of two up through six features at a time. From these results, they used Jeffries-Matusita (J-M) distances for best and mean separability distances among the change signatures.
Incorporating six change features produced the largest average J-M distances, using only five or even four features decreased accuracy (as measured by Kappa) only slightly while retaining good J-M distances.
I generated the ten change features as described in Table 4 in the appendix, and overlaid the burn boundary polygons on each output. In most cases it was not difficult to see differences between the burn area and the rest of the area of interest. Using one of the Hawth’s Tools (Thematic Raster Index) I gathered pixel counts of the classes generated by the IMAGINE Image Difference tool, at 10-, 20- and 40-percent thresholds. In most cases, somewhere between 20- and 40-percent differences, the class counts of pixels inside the burn areas taper off to the point of being almost indistinguishable from how all pixels in the area of interest are distributed.
References
Coppin, Pol R., and Marvin E. Bauer. 1994. Processing of Multitemporal Landsat TM Imagery to Optimize Extraction of Forest Cover Change Features, IEEE Transactions on Geoscience and Remote Sensing, 32(4):918-927.
Coppin, Pol R., and Marvin E. Bauer. 1996. Change Detection in Forest Ecosystems with Remote Sensing Digital Imagery, Remote Sensing Reviews, 13:207-234.
Aldrich, R.C. 1975. Detecting disturbances in a forest environment. Photogramm. Eng. Remote Sens., 41:39-48.
Appendix
Table 1: Landsat 5 Thematic Mapper Scenes - Level 1T (Terrain Corrected)
Level 1T Standard Terrain Correction provides systematic radiometric and geometric accuracy by incorporating ground control points while employing a DEM from the GLS2005 dataset for topographic accuracy. Scenes were projected into UTM Zone 15 North with the WGS84 datum.
Scenes acquired were from path 026 row 027, covering an area of northern Minnesota including the area of interest near the border with Canada in Cook County. Scenes were acquired on 2004-07-14 and 2008-08-10 and each contain some clouds, but less than 10% scene coverage. Grid cells are 30m in size.
Corner Latitude and Longitudes from scene metadata:
PRODUCT_UL_CORNER_LAT = 48.4252163
PRODUCT_UL_CORNER_LON = -92.4674259
PRODUCT_UR_CORNER_LAT = 48.3611961
PRODUCT_UR_CORNER_LON = -89.1264665
PRODUCT_LL_CORNER_LAT = 46.4709197
PRODUCT_LL_CORNER_LON = -92.4867970
PRODUCT_LR_CORNER_LAT = 46.4111048
PRODUCT_LR_CORNER_LON = -89.2669583
Table 2: USGS 3.75" Quadrangles Comprising the Area of Interest
TILETILE NAMEACQUISITION DATE (NAIP Images)
q1051seMunker Island - SE2003-06-16
q1051swMunker Island - SW2003-06-16
q1050seEster Lake - SE2004-08-05
q1050swEster Lake - SW2004-08-05
q1151neGillis Lake - NE2003-06-16
q1151nwGillis Lake - NW2003-06-16
q1150neOgishkemuncie Lake - NE2004-08-05
q1150nwOgishkemuncie Lake - NW2004-08-05
q1153swGunflint Lake - SW2004-08-05
q1152seLong Island Lake - SE2004-08-05
q1151seGillis Lake - SE2003-06-16
q1151swGillis Lake - SW2003-06-16
q1150seOgishkemuncie Lake - SE2004-08-05
q1253nwBrule Lake - NW2004-08-05
q1252neCherokee Lake - NE2004-08-05
Table 3: Coordinates and Spectral Values Used for Dark Subtraction Atmospheric Correction
For the image acquired on 2004-07-14:
X / Y / Band 1 / Band 2 / Band 3 / Band 4 / Band 5 / Band 6 / Band 7652184.7789 / 5331962.877 / 62 / 24 / 20 / 19 / 10 / 124 / 6
643050.6066 / 5327613.271 / 64 / 23 / 18 / 18 / 10 / 125 / 6
646916.9229 / 5324568.547 / 62 / 21 / 18 / 16 / 9 / 125 / 5
658032.5824 / 5323022.021 / 65 / 24 / 19 / 17 / 13 / 125 / 8
665571.8992 / 5320750.56 / 61 / 22 / 19 / 17 / 11 / 125 / 5
670791.4263 / 5312389.651 / 65 / 24 / 19 / 17 / 12 / 121 / 8
667843.3601 / 5328966.482 / 62 / 23 / 18 / 16 / 10 / 122 / 3
646771.9361 / 5342595.247 / 64 / 23 / 18 / 15 / 9 / 121 / 6
656776.0296 / 5343465.168 / 62 / 21 / 17 / 13 / 9 / 125 / 6
669873.1762 / 5344431.747 / 63 / 21 / 17 / 14 / 8 / 124 / 4
63 / 22.6 / 18.3 / 16.2 / 10.1 / 123.7 / 5.7
For the image acquired on 2008-08-10:
X / Y / Band 1 / Band 2 / Band 3 / Band 4 / Band 5 / Band 6 / Band 7652184.7789 / 5331962.877 / 47 / 16 / 11 / 8 / 3 / 127 / 3
643050.6066 / 5327613.271 / 47 / 15 / 11 / 9 / 5 / 128 / 5
646916.9229 / 5324568.547 / 49 / 16 / 10 / 8 / 6 / 127 / 4
658032.5824 / 5323022.021 / 49 / 17 / 11 / 8 / 6 / 127 / 5
665571.8992 / 5320750.56 / 47 / 16 / 11 / 10 / 6 / 126 / 3
670791.4263 / 5312389.651 / 44 / 15 / 11 / 8 / 7 / 125 / 4
667843.3601 / 5328966.482 / 46 / 15 / 11 / 8 / 5 / 126 / 3
646771.9361 / 5342595.247 / 48 / 16 / 11 / 9 / 5 / 125 / 4
656776.0296 / 5343465.168 / 48 / 16 / 11 / 7 / 4 / 126 / 4
669873.1762 / 5344431.747 / 47 / 15 / 11 / 8 / 5 / 128 / 4
47.2 / 15.7 / 10.9 / 8.3 / 5.2 / 126.5 / 3.9
Table 4: Change Features Compared
IndexFormula
STD Brightness(2004 TC Brightness - 2008 TC Brightness)/(2004 TC Brightness + 2008 TC Brightness)
STD Greenness(2004 TC Greenness - 2008 TC Greenness)/(2004 TC Greenness + 2008 TC Greenness)
STD Crippen's NDVI(2004 IPVI - 2008 IPVI)/(2004 IPVI + 2008 IPVI)
STD Red Bands(2004 TM3 - 2008 TM3)/(2004 TM3 + 2008 TM3)
PC2 BrightnessSecond Principal Component(2004 TC Brightness, 2008 TC Brightness)
PC2 GreennessSecond Principal Component(2004 TC Greenness, 2008 TC Greenness)
PC2 WetnessSecond Principal Component(2004 TC Wetness, 2008 TC Wetness)
PC2 Green RatiosSecond Principal Component(2004 Green Ratio, 2008 Green Ratio)
PC2 MIR RatiosSecond Principal Component(2004 MIR Ratio, 2008 MIR Ratio)
PC2 Red BandsSecond Principal Component(2004 Red Bands, 2-008 Red Bands)
STD = Standardized Difference: (Time 1 - Time 2)/(Time 1 + Time 2)
TC = Tasseled Cap component
IPVI = Crippen's NDVI: TM4/(TM4 + TM3)
PC2 = Second Principal Component
Green Ratio = TM4/TM2
MIR Ratio = TM4/TM5
Red Band = TM3