6

Using nlsLM from Package ‘minpack.lm’ is straightforward for a simple example like this:

MODEL:y=(exp(a*x+b*z+c)+d)^f
x=c(0.5,0.3,0.2,0.4)
z=c(0.1,0.6,1,0.9)
y=c(0.2,0.3,0.9,0.9)
fit=nlsLM(y~(exp(a*x+b*z+c)+d)^f,start = list(a = -.03, b = 0.5, c = 1,d=0.02,f=0.003))

And no problem with that. My real data are multiple rasters for one year.So I have one year of data for x,y,z but as rasters that contain pixels,I want to do a non-linear regression pixel by pixel(each regression for one pixel is not related to the next pixel). example data:

   r <- raster(nrows=10, ncols=10); r <- setValues(r, 1:ncell(r))
   r1 <- raster(nrows=10, ncols=10);r1 <- setValues(r1, 1:ncell(r))
   r2 <- raster(nrows=10, ncols=10);r2 <- setValues(r1, 1:ncell(r))
st1=stack(r,r1,r2)

   rl <- raster(nrows=10, ncols=10); r <- setValues(r, 1:ncell(r))
   rS <- raster(nrows=10, ncols=10);r1 <- setValues(r1, 1:ncell(r))
   rT <- raster(nrows=10, ncols=10);r2 <- setValues(r1, 1:ncell(r))
st2=stack(rl,rS,rT)

  re <- raster(nrows=10, ncols=10); r <- setValues(r, 1:ncell(r))
  ru <- raster(nrows=10, ncols=10);r1 <- setValues(r1, 1:ncell(r))
  rg <- raster(nrows=10, ncols=10);r2 <- setValues(r1, 1:ncell(r))
st3=stack(re,ru,rg)

so first x would be the first pixel of st1(for the whole temporal period)

so first y would be the first pixel of st2(for the whole temporal period)

so first z would be the first pixel of st3(for the whole temporal period)

do regression for this pixel then do the same for all other pixels

As a result I get a map (raster) 10*10 for each best parameter(a,b,c,d,f)

I hope it is clear and you have any help!

  • 1
    There may be some underlying problems to resolve before doing any coding. In particular, you are trying to fit a five-parameter model to just four data points (or only three, in the case of the rasters!). Another issue--which won't go away even with more data--is that your model is not identifiable, because exp(a*x+b*z+c)^d is algebraically the same as exp(d*a*x+d*b*z+d*c), which means only four of the five parameters can ever be determined. – whuber Mar 10 '15 at 20:45
  • Would the regression example at the bottom of this page help? – MickyT Mar 10 '15 at 22:27
  • @whuber Thanks for your comment.The model is not important for the moment but the principle of doing non-linear regression for multiple rasters is the most important. So you can please use whatever model to show we can do this? – usersam Mar 11 '15 at 8:21
2

This should do what you want.

The model I have used is not the one you posted and it probably doesn't make any sense at all, however it does demonstrate the principle.

# Set up the rasters
r1 <- r2 <- r3 <- r4 <- r5 <- r6 <- raster(nrows=10, ncols=10);
# Populate them with some values
r1 <- setValues(r1,runif(100,min=1,max=100));
r2 <- setValues(r2,runif(100,min=1,max=100));
r3 <- setValues(r3,runif(100,min=1,max=100));
r4 <- setValues(r4,runif(100,min=1,max=100));
r5 <- setValues(r5,runif(100,min=1,max=100));
r6 <- setValues(r6,runif(100,min=1,max=100));
# Stack them
st1 <- stack(r1,r2,r3);
st2 <- stack(r4,r5,r6);
# Set up the function
test <- function(r) {
    x <- r[1:3];
    y <- r[4:6];
    result <- c(NA,NA);
    try(result <- c(coef(nlsLM(y ~ a + b * x, start = list(a = 0.12345, b = 0.54321), na.action = na.omit))));
    result;
}
# Stack the stacks
s <- stack(st1,st2);
# Calculate a new Raster
rNLR <- calc(s, test);
  • Thanks for your answer. As I said my data are temporal for one year so st1 has 365 rasters and st2 has 365 rasters, etc...I wonder how we change this part x <- r[1:3]; y <- r[4:6] ?? – usersam Mar 11 '15 at 8:24
  • The model is not important but the way we do it is important.Thanks once again. – usersam Mar 11 '15 at 8:43
  • If you think of r as the pixel. At that point you would have 365 layers twice in the stack. So to put the pixel value into the vector x you specify it as x <- r[1~365]. Y is then set to the next 365 layer values in the stack. – MickyT Mar 11 '15 at 8:51
  • so Y would be y <- r[366~730]?? what If I have another variable z, what would be Z <- r[?~?] – usersam Mar 11 '15 at 8:55
  • 1
    To be honest, I wasn't trying to help with the actual model, rather the way to apply it to the raster. – MickyT Mar 12 '15 at 19:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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