Running regression in each sample separately?

Question

I'd like to smooth an NDVI image by regression analysis; response of time series NDVI image to Julian date to represent the seasonal changes in NDVI as a function of Julian day. The time series NDVI image has 26 bands (8-day images during April to October) and 12618/4144 col/row. I have chosen 30 point samples in each land cover (30 x 5 LC=150) from all band of time series NDVI image using stratified random sampling method. I need to discuss my further steps with you all.

Do I need to run the regression for each sample point (150 samples) and calculate regression equation extracted from step, 1 in each point sample (one by one) again to reproduce predicted NDVI image?

Is there any method to simplify these processes to save time?

Answer 1

Using the raster package in R you could apply a pixel-wise regression estimate of NDVI ~ time. Here is an example for a linear model, locally-weighted polynomial regression and regression coefficients.

library(raster)

# Create some example data
r <- raster(nrow=100, ncol=100)
  r[] <- runif(ncell(r),-1,1)
    rt <- stack(r)      
      for(i in 2:26) {
        r <- rt[[1]] 
          r[] <- runif(ncell(r),-1,1)     
            rt <- addLayer(rt, r)  
      }

# Create a time vector to act as x
time <- sort(sample(1:365,nlayers(rt))) 
        
# linear (lm) regression estimate(s) of ndvi ~ time
t.lm.predict <- function(x) {if (is.na(x[1])) {NA} else {predict(lm(x ~ time))}} 
f.pred <- calc(rt, t.lm.predict)
  plot(f.pred)

# locally-weighted polynomial regression of ndvi ~ time
t.lowess <- function(x,...) { if (is.na(x[1])) { NA } else { lowess(x,y,...)$y } } 
f.pred <- calc(rt, t.lowess)
  plot(f.pred)
      
# slope and intercept of ndvi ~ time
t.lm.coef <- function(x) {
  if (is.na(x[1])) { NA } else { lm(x ~ time)$coefficients }
  }
f.coef <- calc(rt, t.lm.coef)