Eddie,

Your stack consist of 24 layers but you used xcost<-cos(2*pi*time/12)
xsine<-sin(2*pi*time/12). Could this be the case?

I tried your code with my data and it took ages. ;(


On Wed, Nov 6, 2013 at 11:53 AM, Eddie Smith <eddieatr@gmail.com> wrote:

> Hi list,
>
> Thanks Jan. Yup, I think NA values are really the case here.
>
> I found a reply from Robert here
> http://r-sig-geo.2731867.n2.nabble.com/Problems-with-NA-values-and-lm-in-calc-in-raster-package-td7455990.html but
> I think I have mistakenly wrote my code.
>
> Could somebody kindly check this for me?
>
> #Create 24 dummy raster layers
> r <- raster(nrow=10, ncol=10)
> s <- stack( sapply(1:24, function(i) setValues(r, rnorm(ncell(r), i, 3) ))
> )
> time <- 1:nlayers(s)
>
> #Throw in NA values to the rasterstack
> s[1] <- NA
>
> #Create sine function that deals with NA values base on Robert
> fun1 <- function(x) {
> if (is.na(x[1])) {
> NA
> } else {
> xcost<-cos(2*pi*time/12)
> xsine<-sin(2*pi*time/12)
> m = lm(s~xcost+xsine)
> m$coefficients[2]
> }
> }
>
> fun2 <- function(x) {
> m <- NA
> try( m <- lm(s~xcost+xsine)$coefficients[2] ,silent=T)
> m
> }
>
> #run
> e1 <- calc(s, fun1)
> Error in `colnames<-`(`*tmp*`, value = "layer") :
>   length of 'dimnames' [2] not equal to array extent
>
> e2 <- calc(s, fun2)
>
> Thanks in advance
>
>
> On Mon, Nov 4, 2013 at 3:24 PM, Verbesselt, Jan <jan.verbesselt@wur.nl>wrote:
>
>> Dear Sofea,
>>
>> You could try something like this:
>>
>> fun3 <- function(x){
>>         time <- 1:length(x)
>>         xc <-cos(2*pi*t/365)
>>         xs <-sin(2*pi*t/365)
>>         fit <- lm(x~xc+xs)
>>         return(fit$coefficients[1])
>>         }
>>
>> s3 <- calc(s,  fun = fun3)
>>
>>
>> But check for NA's within your time series as these might trigger
>> errors/warnings. Maybe also check na.action options like na.approx(). Can
>> you provide a reproducible example in that case?
>>
>> Best,
>> Jan Verbesselt
>> Wageningen University
>> The Netherlands
>>
>>
>>
>> On 31/10/13 21:39, "Sofea Wright" <sofeawright@gmail.com> wrote:
>>
>> >Dear list,
>> >
>> >I am trying to fit a sine curve to my data. When dealing with non-raster
>> >data, it's quite straight forward:
>> >
>> ># Fit sine curve to describe the data
>> >temperature<-maxtemp$temperature
>> >time<-1:n
>> >xc<-cos(2*pi*time/365)
>> >xs<-sin(2*pi*time/365)
>> >fit<-lm(temperature~xc+xs)
>> >summary(fit)
>> >
>> >I am not sure how to do that with rasterstack as I keep on getting error
>> >message.
>> > > s
>> >class       : RasterStack
>> >dimensions  : 1000, 1621, 1621000, 365  (nrow, ncol, ncell, nlayers)
>> >resolution  : 0.04166667, 0.04166667  (x, y)
>> >extent      : 96.5, 165.5, -16.5, 26.5  (xmin, xmax, ymin, ymax)
>> >coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
>> >+towgs84=0,0,0
>> >
>> >s <- stack(tif)
>> >time <- 1:nlayers(s)
>> >fun1 <- function(x){
>> >  xc<-cos(2*pi*time/365)
>> >  xs<-sin(2*pi*time/365)
>> >          }
>> >
>> >fun2 <- function (x) {
>> >        lm(s~xc+xs)
>> >        }
>> >Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
>> >  cannot use this function
>> >
>> >s1 <- calc(s,  fun = fun1)
>> >s2 <- calc(s,  fun = fun2)
>> >
>> >Anybody?
>> >
>> >Thank you
>> >Sofea
>> >
>> >       [[alternative HTML version deleted]]
>> >
>> >_______________________________________________
>> >R-sig-Geo mailing list
>> >R-sig-Geo@r-project.org
>> >https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>

	[[alternative HTML version deleted]]

