-1

I want to fit a polynomial function (max. 3rd order) on each raster cell over all my spectral bands (Landsat 1-7) creating a new raster(stack) representing the coefficients. I got my data (including NA values) in a stack with 6 Layer (Landsat Band 1-7[excluding 6]).

I guess somehow I should tell the polynomial function on which spectral wavelength the bands are located

  1. Landsat7 Wavelength (micrometers)
    • Band 1 0.45-0.52
    • Band 2 0.52-0.60
    • Band 3 0.63-0.69
    • Band 4 0.77-0.90
    • Band 5 1.55-1.75
    • Band 7 2.09-2.35

so that it can fit it properly. Has anyone an idea how to do that polynomial fitting of each cell and extracting the coefficients in R? Thanks for any help!

user2978751
  • 57
  • 1
  • 10
  • edit: my NA-values are on the same position in each band (rasterlayer) - so exact the same amount on the same position. – user2978751 May 25 '15 at 11:39

1 Answers1

1

Your question is not very clear, as you do not specify what you are fitting. I am guessing it is band number. You can do something like this. 

library(raster)
b <- brick(system.file("external/rlogo.grd", package="raster"))
b[[2]][125:225] <- NA
s <- stack(b, flip(b, 'y'))
names(s) <- paste0('b', 1:6)
bands <- 1:6
f <- function(x) {
    # in case of NAs; match the number of coefficients returned
    if (any(is.na(x))) return(c(NA, NA, NA))
    m <- lm(x ~ bands + I(bands^2))
    coefficients(m)
}
z <- calc(s, f)
z
plot(z)

If you need to speed this up you can follow the example here: https://gis.stackexchange.com/questions/144211/want-cell-linear-regression-values-for-a-netcdf-or-multi-band-raster/144408#144408

Community
  • 1
  • 1
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Indeed i want to fit band number - to **get the spectral curve or signature for each pixel value**. Thanks very much for the snippet - but **isn't it necessary to specify the wavelength for the bands as they differ in their "distance" to one another?** (wavelength of Band 1 leads directly into Band 2 - whereas between Band 4 and Band 5 are about 0.6 micrometers difference in wavelength (that is not covered by measurement)). I guess this should be taken into account in my fitting? (The 6 Bands on my X-Axis [Y representing intensity of reflectance] are not of equal distance to one another...) – user2978751 May 25 '15 at 20:55
  • You can certainly replace bands by a vector of wavelenghts. Whether that is sensible or not I cannot know, as I do not have any context or why you want to do this. – Robert Hijmans May 26 '15 at 01:56
  • I tried to apply the function on a rasterstack with NA values and unfortunately it did not work. (I also tried: `d <- na.omit(data.frame(x))# if (nrow(d) < 3) return(NA)# m <- lm(x ~ bands + I(bands^2)#`). Should it work with NAs? Concerning the context: I want to add additional info to my data. So for each pixel not only the observed values in the 6 spectral bands but also an approximation of a curve representing the spectral signature of the pixel will be basis for further mapping. Unfortunately so far I wasn't able to fit a polynomial function matching good the observed value. – user2978751 May 27 '15 at 15:10
  • I have updated my answer to show how to deal with NAs – Robert Hijmans May 27 '15 at 15:52