[R-sig-Geo] Fit a sine curve to raster stack

[Sofea Wright]

Thanks a lot Arnaud Mosnier!

The code is running however, the final output is strange. The minimum
values and maximum values is NA. I would expect it to give me the
slope (coeefficients[2] is the slope). I just could not figure out
where I went wrong.
> e2
class       : RasterLayer
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : layer
values      : NA, NA  (min, max)

The full code(based on Eddie) is now:

#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(x~xcost+xsine)
m$coefficients[2]
}
}

fun2 <- function(x) {
m <- NA
try( m <- lm(x~xcost+xsine)$coefficients[2] ,silent=T)
m
}

#run
e1 <- calc(s, fun1)
e2 <- calc(s, fun2)



On Thu, Nov 7, 2013 at 2:42 PM, Arnaud Mosnier <a.mosnier at gmail.com> wrote:
> Just change the s for a x in
>
> 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]
> }
> }
>
> this should work !
>
> fun1 <- function(x) {
> if (is.na(x[1])) {
> NA
> } else {
> xcost<-cos(2*pi*time/12)
> xsine<-sin(2*pi*time/12)
> m = lm(x~xcost+xsine)
> m$coefficients[2]
> }
> }
>
>
>
> Message: 2
> Date: Wed, 6 Nov 2013 11:53:26 +0000
> From: Eddie Smith <eddieatr at gmail.com>
> To: "Verbesselt, Jan" <jan.verbesselt at wur.nl>
> Cc: "r-sig-geo at r-project.org" <r-sig-geo at r-project.org>
> Subject: Re: [R-sig-Geo] Fit a sine curve to raster stack
> Message-ID:
>         <
> CABaJH78P4TPbG4pR4a1kTGABE7jZv-o0C=TuQvnW7LVSgfHDew at mail.gmail.com>
> Content-Type: text/plain
>
> 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
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo