I have Landsat images for the same area over 20 years (one photo per year in same season). I would like to detect changes occured in forested area like clear-cut and bark beetle outbreak over time. I would like to use Maximum likelihood supervised classification.

My questions are:

  1. As the signature of change (clear-cut, bark beetle outbreak) is the same over time, can I use the same signature file from one year to classify all 20 years?

    • in this case I suppose that firstly I have to normalize all scenes over time using Dark Object subtraction or other technique - which one should be the best to detect changes in forested area?


  1. Do I have to create signature file for every year and then run classification?

    • in that case, there is no problem with changing number of possibly detected classes of change - before outbreak you can find only one class = forest, after outbreak you can find "forest", "clear-cut", "regrowth"...

I can use ERDAS and ArcGIS as software.

  • Thank you for your answer @WhiteboxDev. Images are already normalized but ok, I will not use the same signature file. What do you think about the changing number of classes over years? thank you very much again.
    – maycca
    Oct 8, 2014 at 14:44
  • Thanks again. I was really thinking that "normalisation" is normalising both - atmospheric conditions and phenological phase of forest, varying moisture conditions, etc... I didnt look on that in this way. Good luck to you :)
    – maycca
    Oct 8, 2014 at 15:00
  • Have you looked at "How to represent trend over time?": gis.stackexchange.com/questions/52502/…. Also, have you looked at Erdas's change detection program, DeltaCue?
    – Aaron
    Oct 8, 2014 at 16:06
  • Hi @Aaron, thank you for response, I have read the "how-to-represent-trend-over-time". However I dont think that it is possible to express my trends using simple linear regression. P.ex. if I am plotting raster cell value (Y) to time (X) it is not growing or declining for whole period, it is mostly like a curve. Pixel values Y should firstly be similar(representing undisturbed forest), after few years could decline (as bark beetle outbreak), then there could be clearcut = steep growtth or regrowth = slow decline). For now I dont have a licence for DeltaCue. Can you see the other possibility?
    – maycca
    Oct 9, 2014 at 10:08
  • as a remark, the full Landsat archive have already been processed for forest disturbance . Maybe you will make a better job by adjusting locally, but probably this will be enough for your needs earthenginepartners.appspot.com/science-2013-global-forest
    – radouxju
    Nov 24, 2017 at 8:33

1 Answer 1


Responding directly to your question: applying the same signature file / spectral response to each image is somewhat risky, due to the potential for variation in how the clear-cut area looks, and how the health forest looks, due to phenological variation. As such, I'd classify each image on its own, and then go on to do change analysis on the output. Note that using a sub-pixel analysis will most likely provide better results than a hard classification. This approach is very operator dependent, as each image / classification will require a set of endmembers.

If you feel more adventurous, then I'd suggest that you look into LandTrendr. It is a tool developed to work with temporal stacks of Landsat images, and is designed to evaluate disturbances. Here is an example view taken from the landtrendr website: LandTrendr example view

See more about LandTrendr at http://landtrendr.forestry.oregonstate.edu/

And in the article:

Kennedy, Yang, Cohen; Detecting trends in forest disturbance and recovery using yearly Landsat time series 1. LandTrendr — Temporal segmentation algorithms; Remote Sensing of Environment 114 (2010) 2897–2910

Another option is BFAST in R: http://cran.r-project.org/web/packages/bfast/index.html.

BFAST splits the timeseries into two components, a linear trend and a harmonic seasonal. Below is an illustration of a BFAST analysis.

Slightly modified BFAST result graph

On this illustration, I have pointed out the two test parameters that the script outlined below will extract. The 'Magnitude of Change' is implemented in BFAST, while the 'st_diff' is something that I have added. Another possible addition that could be implemented is looking at the change in slope between the two components.

Below is an example of how BFAST could be implemented for a regular timeseries of a single index (NDVI, EVI, NDSI, TC-Greeness etc.), which is good for MODIS. To use it with Landsat it may be necessary to modify the timeseries-structure.

#bfast works on 1 parameter only. In this case I use NDVI, but EVI and Tasseled Cap Greeness are also good for vegetation studies.
#the code can is useful for any type of timeseries of indices.
vegndvi<-brick("ndvi_subset_brick.grd") #loading a raster brick of data. Each timestep is a layer in the brick.

#Establishing a helper-function - the function applies BFAST to 1 pixel at a time
xbfast <- function(data) {  
  ndvi <- ts(data, frequency=46, start=1) #Here a regular timeseries is established. This works for MODIS, due to the high revisit time. Landsat may need an irregularly spaced time series.
  result <- bfast(ndvi, season="harmonic", max.iter=2, breaks=1) #bfast is applied to the data series
  niter <- length(result$output) #extracting results
  out <- result$output[[niter]] #more extraction
  bp <- out$Wt.bp #breakpoint of the seasonality component
  st <- out$St #the seasonality component
  st_a <- st[1:bp] #seasonality untill the breakpoint
  st_b <- st[bp:NROW(ndvi)] #seasonality after the breakpoint
  st_amin <- min(st_a) #determines the smallest value in the seasonality-component before the break
  st_amax <- max(st_a) #determines the largest value in the seasonality-component before the break
  st_bmin <- min(st_b) #determines the smallest value in the seasonality-component after the break
  st_bmax <- max(st_b) #determines the largest value in the seasonality-component after the break
  st_adif <- st_amax - st_amin #determines the range of the seasonality-component before the break
  st_bdif <- st_bmax - st_bmin #determines the range of the seasonality-component after the break
  st_dif <- st_bdif - st_adif #determines the change in the seasonality across the break
  Magni<-result$Magnitude #magnitude of the biggest change detected in the trend component
  Timing<-result$Time #timing of the biggest change detected in the trend component

#Establishing a function to apply BFAST to a raster-stack
bfastfun <- function(y) {
  percNA <- apply(y, 1, FUN=function(x) (sum(is.na(x))/length(x)) ) #determines the % of NA-values in the given pixel
  i <- (percNA<0.2) #tests the percentage of NA against a stop/go barrier
  res <- matrix(NA, length(i), 4) #creating an empty result matrix
  if (sum(i) > 0) { #if go, then apply.
    res <- t(apply(y[i,], 1, xbfast)) #applying the BFAST helper function to the data
  return(res) #providing the result matrix.

processed <- calc(dec_subset, fun=bfastfun)

diff=subset(processed, 1)
time=subset(processed, 2)
magn=subset(processed, 3)
magn_time=subset(processed, 4)

If a more detailed analysis of BFAST processing and results is desired, I wrote my masters thesis on the topic and can provide that for additional reading.

  • HI Mikkel, thank you for your answer. Yes, I have found this approach, but unfortunately I am not very familiar with ENVI and IDL. Do you think that I can run the script in ERDAS or in R? or it is possible only in ENVI? thank you, maria
    – maycca
    Jan 12, 2015 at 9:29
  • The LandTrendr code unfortunately only functions with IDL. In R, you have the option of looking at BFAST: cran.r-project.org/web/packages/bfast/index.html - I have only used this on MODIS data, but it should be possible to use on Landsat as well. Jan 12, 2015 at 10:13
  • If you want, I can provide a code example from one of my MODIS projects that you can attempt to adjust into something that you can use for Landsat. Jan 12, 2015 at 10:49
  • Do you think you can provide it to me? If it is not a probleme, I think it will be very helpful... I have my e-mail adress on my stack.exchange profile. Thank you Mikkel.
    – maycca
    Jan 12, 2015 at 11:13
  • I have added a description of BFAST and a code-example to my original post to allow others to use the ideas. Jan 12, 2015 at 12:16

Not the answer you're looking for? Browse other questions tagged or ask your own question.