# NDVISE2.Rcm
# analyze vegetation/temperature from MODIS satellite

#place <- "Khao Leam"
#place <- "Mayo"
place <- "Bajo"
#place <- "Java"

setwd("D:/bogorWorkshop")

pl <- substr(place,1,1)
if (pl=="K") read.csv("TDkhaoleam.csv",header=FALSE,as.is=TRUE) -> T1
if (pl=="M") read.csv("TDmayo.csv",header=FALSE,as.is=TRUE) -> T1
if (pl=="B") read.csv("TDbajo.csv",header=FALSE,as.is=TRUE) -> T1
if (pl=="J") read.csv("TDbogban.csv",header=FALSE,as.is=TRUE) -> T1

dim(T1)

names(T1)[1:7] <- c("rID","fn1","yrDay","latLon","fn4","var","t0") 
T1$yrDay <- substr(T1$yrDay,2,8)

str(T1[,c(1:7)])

T1$yrDay <- as.integer(T1$yrDay)	# convert character to integer

nc <- ncol(T1)      			# number of columns

names(T1)[8:nc] <- paste("t",c(1:(nc-7)),sep="")
T1$rID <- as.integer(c(1:nrow(T1)))

lat <- substr(T1$latLon[1],4,10)
for (i in c(2:7)) {
 if (substr(lat,i,i)=="L") break
}
lat <- substr(lat,1,i-1)
lat <- as.numeric(lat)
lon <- substr(T1$latLon[1],5,20)
for (i in c(1:8) ) {
 if (substr(lon,i,i+2)=="Lon") break
}
lon <- substr(lon,i+3,i+10)
for (i in c(2:8)) {
 if (substr(lon,i,i)=="S") break
}
lon <- substr(lon,1,i-1)
lon <- as.numeric(lon)

td <- T1[,c(3,8:nc)]

year <- as.integer(substr(td$yrDay,1,4))
day <- as.integer(substr(td$yrDay,5,7))

tab <- addmargins(table(day,year))
tab

ncd <- ncol(td)      			# number of columns

if (pl %in% c("B","K","M"))
 days <- 2000000+c(1169,1177,13105,13209,13217,13225,13241)
if (pl=="J")
 days <- 2000000+c(1169,1177,13105,13209,13217,13225,13233,13241)

 for (d in days) {
 td.a <- td[td$yrDay==2001161,]
 td.a$yrDay <- d
 td.a[,c(2:ncd)] <- NA
 td <- rbind(td,td.a)
}

td <- td[order(td$yrDay),]

td <- subset(td,td$yrDay>2000169 & td$yrDay<2015177)
year <- as.integer(substr(td$yrDay,1,4))
day <- as.integer(substr(td$yrDay,5,7))

nYrs <- length(unique(year))-1

tab <- addmargins(table(day,year))
tab

td <- td[,c(2:ncd)]
td[td=="F"] <- 0
for (j in c(1:(ncd-1))) {
 td[,j] <- as.numeric(td[,j])
}
td <- as.matrix(td) 
td <- ifelse(td==0,NA,td)
td <- td-273.15

np <- ncol(td)				# number of pixels for LST_day
wid <- sqrt(np)				# width of area (assumed square) in km
nt <- nrow(td)				# number of days
td <- t(td)				# transpose
dim(td) <- c(np*nt,1)			# stack as single column

pixlT <- c(1:wid^2)			# number of pixels in each row for LST
pixlT <- as.vector(pixlT)
yt <- as.data.frame(rep(pixlT,nt))
names(yt) <- "pixlT"
yt$pixlT <- as.integer(yt$pixlT)

wv <- 4*(wid-1)+1			# number of pixels in each row for NDVI
pixlV <- c(1:wv^2)
dim(pixlV) <- c(wv,wv)
pixlV <- t(pixlV)
pixlV <- t(pixlV)
dim(pixlV) <- c(wv^2,1)
pixlV <- as.vector(pixlV)

yrDay <- 1000*year+day
yt$yrDay <- rep(yrDay-2000000,each=np)
yt$td <- as.numeric(td)

yt$year <- floor(yt$yrDay/1000)
yt$day <- yt$yrDay-1000*yt$year

yt$tOK <- ifelse(is.na(yt$td),0,1)
ytOK <- subset(yt,tOK==1)

addmargins(tab <- table(ytOK$day,ytOK$year+2000))
tab
round(100*addmargins(tab)[,17]/(nYrs*np),1)
100-mean(100*tab/np)
#---------------------------
# plot of day temperatures for central pixel

pixel0 <- (np+1)/2

windows(7,5)
par(las=1,mar=c(0,0.5,0.5,0.3),oma=c(3,2,3,0.8),tcl=0,mgp=c(1.1,0.1,0))

ylab <- expression(paste("LST.Day (",degree,"C)",sep=""))

ymax <- max(yt$td,na.rm=T)
ymin <- min(yt$td,na.rm=T)

pxl <- pixel0 

tdk <- subset(yt,pixlT==pxl)
plot(tdk$day,tdk$td,type="n",ylim=c(ymin,ymax),xaxt="n")

xlab <- c(1+16*c(0:22))
abline(v=xlab,col="grey")
axis(side=1,at=xlab,lab=xlab,padj=0)
axis(side=1,at=182.5,lab="Day of Year",padj=1.6,cex.axis=1.2,tcl=0)
mtext(side=3,line=0.2,adj=-0.07,ylab,cex=1.2)

points(tdk$day,tdk$td,pch=20)

n <- sum(ifelse(is.na(tdk$td),0,1))

text(185,ymax,adj=c(0.5,0.5),"Central Pixel",cex=1.2)
text(365,ymax,adj=c(1,0.5),paste("n: ",n,sep=""),cex=1.2)

title1 <- paste("Lat:",lat," Lon:",lon," ",place,sep=" ")
mtext(outer=TRUE,side=3,adj=0.5,line=0.5,title1,cex=1.4)

table(tdk$day,tdk$tOK)

#-----------------------------------
# ... with spline function fitted

days <- unique(yt$day)
days <- days[order(days)]
days <- as.data.frame(days)
names(days) <- "day"

windows(7,5)
par(las=1,mar=c(0,0.5,0.5,0.3),oma=c(3,2,3,0.8),tcl=0,mgp=c(1.1,0.1,0))

ylab <- expression(paste("LST.Day (",degree,"C)",sep=""))

ymax <- max(yt$td,na.rm=T)
ymin <- min(yt$td,na.rm=T)

pxl <- pixel0 

tdk <- subset(yt,pixlT==pxl)
plot(tdk$day,tdk$td,type="n",ylim=c(ymin,ymax),xaxt="n")

xlab <- c(1+16*c(0:22))
abline(v=xlab,col="grey")
axis(side=1,at=xlab,lab=xlab,padj=0)
axis(side=1,at=182.5,lab="Day of Year",padj=1.6,cex.axis=1.2,tcl=0)
mtext(side=3,line=0.2,adj=-0.07,ylab,cex=1.2)

points(tdk$day,tdk$td,pch=20)

n <- sum(ifelse(is.na(tdk$td),0,1))

 kn <- c(20,50,80,130,180,230,280,320)
 text(kn,ymin,adj=c(0.5,0),"+",col="blue",cex=1.2)
 p <- length(kn)		# number of spline knots
 set <- !is.na(tdk$td)
 s <- tdk[set,c(3:5)]
 s <- s[order(s$day),]
 x <- s$day 
 yy <- as.data.frame(s[,1:2])
 names(yy) <- c("y1","yr")
 yy$day <- s$day
 yy$x <- x
 d1 <- kn[p]-kn[p-1]; d2 <- kn[p]-kn[p-2]; d3 <- kn[p]-kn[p-3]
 d4 <- kn[p-1]-kn[p-2]; d5 <- kn[p-1]-kn[p-3]; d6 <- kn[p-2]-kn[p-3]
 for (j in c(1:(p-4))) {
  sj <- ifelse(x>kn[j],(x-kn[j])^3,0)
  sj <- sj-((kn[p]-kn[j])*(kn[p-1]-kn[j])*(kn[p-2]-kn[j])/(d3*d5*d6))*ifelse(x>kn[p-3],(x-kn[p-3])^3,0)
  sj <- sj+((kn[p]-kn[j])*(kn[p-1]-kn[j])*(kn[p-3]-kn[j])/(d2*d4*d6))*ifelse(x>kn[p-2],(x-kn[p-2])^3,0)
  sj <- sj-((kn[p]-kn[j])*(kn[p-2]-kn[j])*(kn[p-3]-kn[j])/(d1*d4*d5))*ifelse(x>kn[p-1],(x-kn[p-1])^3,0)
  sj <- sj+((kn[p-1]-kn[j])*(kn[p-2]-kn[j])*(kn[p-3]-kn[j])/(d1*d2*d3))*ifelse(x>kn[p],(x-kn[p])^3,0)
  yy[,(j+4)] <- sj
  names(yy)[j+4] <- paste("s",j,sep="")
 }
 x1 <- x-365*ifelse(x>kn[p-3],(x-kn[p-3])^3,0)/(d3*d5*d6)
 x1 <- x1+365*ifelse(x>kn[p-2],(x-kn[p-2])^3,0)/(d2*d4*d6)
 x1 <- x1-365*ifelse(x>kn[p-1],(x-kn[p-1])^3,0)/(d1*d4*d5)
 x1 <- x1+365*ifelse(x>kn[p],(x-kn[p])^3,0)/(d1*d2*d3)
 mod1 <- lm(data=yy,y1~x1+s1+s2+s3+s4)
 points(x,mod1$fit,type="l",col=2,lwd=2)
 res <- summary(mod1)
 rsq <- round(res$adj.r,3)

text(0,ymax,adj=c(0,1),paste("r-sq: ",rsq,sep=""),cex=1.2)

text(185,ymax,adj=c(0.5,1),"Central Pixel",cex=1.2)
text(365,ymax,adj=c(1,1),paste("n: ",n,sep=""),cex=1.2)

title1 <- paste("Lat:",lat," Lon:",lon," ",place,sep=" ")
mtext(outer=TRUE,side=3,adj=0.5,line=0.5,title1,cex=1.4)

#-----------------------------------
# Plots of LST_Day by day of year for seven locations

ymax1 <- -10
ymin1 <- 10

pixel0 <- (np+1)/2

days <- unique(yt$day)
days <- days[order(days)]
days <- as.data.frame(days)
names(days) <- "day"
seasonal <- NULL

locLab <- c("NW","NE","West","Centre","East","SW","SE")
if (pl=="N") locLab <- c("NW","North","West","Centre","East","South","SE")

windows(14,10)
par(mfrow=c(3,3),las=1,mar=c(0,0.5,0.5,0.3),oma=c(3,2,3,0.8),tcl=0,mgp=c(1.1,0.1,0))

ylab <- expression(paste("LST.Day (",degree,"C)",sep=""))

ymax <- max(yt$td,na.rm=T)
ymin <- min(yt$td,na.rm=T)

for (k in c(1:7)) {
 if (k==1) pxl <- pixel0-198
 if (k==2) pxl <- pixel0-189
 if (k==3) pxl <- pixel0-9
 if (k==4) pxl <- pixel0 
 if (k==5) pxl <- pixel0+9
 if (k==6) pxl <- pixel0+189 
 if (k==7) pxl <- pixel0+198

 tdk <- subset(yt,pixlT==pxl)
 if (k %in% c(3,6)) plot(1,  type="n",axes=F)
 plot(tdk$day,tdk$td,type="n",ylim=c(ymin,ymax),xaxt="n",yaxt="n")
 if (k %in% c(3,6,7)) {
  xlab <- c(1+16*c(0:22))
  axis(side=1,at=xlab,lab=xlab,padj=0)
  axis(side=1,at=182.5,lab="Day of Year",padj=1.6,cex.axis=1.4,tcl=0,font=2)
 }
 if (k %in% c(1,3,6)) axis(side=2,cex.axis=1.4)
 if (k==1) mtext(side=3,line=0.2,adj=-0.11,ylab,cex=1.2)

 points(tdk$day,tdk$td,pch=20,cex=1.2)

 n <- sum(ifelse(is.na(tdk$td),0,1))
 if (n>10) {
  kn <- c(20,50,80,130,180,230,280,320)
  text(kn,ymin,adj=c(0.5,0),"+",col="blue",cex=1.2)
  p <- length(kn)		# number of spline knots
  set <- !is.na(tdk$td)
  s <- tdk[set,c(3:5)]
  s <- s[order(s$day),]
  x <- s$day 
  yy <- as.data.frame(s[,1:2])
  names(yy) <- c("y1","yr")
  yy$day <- s$day
  yy$x <- x
  d1 <- kn[p]-kn[p-1]; d2 <- kn[p]-kn[p-2]; d3 <- kn[p]-kn[p-3]
  d4 <- kn[p-1]-kn[p-2]; d5 <- kn[p-1]-kn[p-3]; d6 <- kn[p-2]-kn[p-3]
  for (j in c(1:(p-4))) {
   sj <- ifelse(x>kn[j],(x-kn[j])^3,0)
   sj <- sj-((kn[p]-kn[j])*(kn[p-1]-kn[j])*(kn[p-2]-kn[j])/(d3*d5*d6))*ifelse(x>kn[p-3],(x-kn[p-3])^3,0)
   sj <- sj+((kn[p]-kn[j])*(kn[p-1]-kn[j])*(kn[p-3]-kn[j])/(d2*d4*d6))*ifelse(x>kn[p-2],(x-kn[p-2])^3,0)
   sj <- sj-((kn[p]-kn[j])*(kn[p-2]-kn[j])*(kn[p-3]-kn[j])/(d1*d4*d5))*ifelse(x>kn[p-1],(x-kn[p-1])^3,0)
   sj <- sj+((kn[p-1]-kn[j])*(kn[p-2]-kn[j])*(kn[p-3]-kn[j])/(d1*d2*d3))*ifelse(x>kn[p],(x-kn[p])^3,0)
   yy[,(j+4)] <- sj
   names(yy)[j+4] <- paste("s",j,sep="")
  }
  x1 <- x-365*ifelse(x>kn[p-3],(x-kn[p-3])^3,0)/(d3*d5*d6)
  x1 <- x1+365*ifelse(x>kn[p-2],(x-kn[p-2])^3,0)/(d2*d4*d6)
  x1 <- x1-365*ifelse(x>kn[p-1],(x-kn[p-1])^3,0)/(d1*d4*d5)
  x1 <- x1+365*ifelse(x>kn[p],(x-kn[p])^3,0)/(d1*d2*d3)
  mod1 <- lm(data=yy,y1~x1+s1+s2+s3+s4)
  points(x,mod1$fit,type="l",col=2,lwd=2)
  ss <- as.data.frame(unique(cbind(yy$day,yy$yr,yy$y1,round(mod1$fit,4))))
  names(ss) <- c("day","yr","y1","ss")
  ss <- merge(days,ss,by.x="day",by.y="day",all=T)
  ss$loc <- locLab[k]
  seasonal <- rbind(seasonal,ss)
  ymax1 <- max(ymax1,max(tdk$td,na.rm=T))
  ymin1 <- min(ymin1,min(tdk$td,na.rm=T))
 }
 text(175,ymax,adj=c(1,1),locLab[k],cex=1.5)
 text(345,ymax,adj=c(1,1),paste("n: ",n,sep=""),cex=1.5)
}
title1 <- paste("Lat:",lat," Lon:",lon," ",place,sep=" ")
mtext(outer=TRUE,side=3,adj=0.5,line=0.5,title1,cex=1.4)
#-----------------------------------
# plot of LST_Day by yrDay for central pixel

ymax <- ymax1+1
ymin <- ymin

windows(8,5)
par(las=1,mar=c(0,0.5,0.5,0.3),oma=c(3.2,2.5,3,0.8),tcl=0,mgp=c(1.1,0.1,0))

 plot(1,type="n",ylim=c(ymin,ymax),xlim=c(30,702),xaxt="n",yaxt="n")
 axis(side=1,at=c(22+46*c(0:15)),lab=c(0:15),padj=0.2,tcl=0,cex.axis=1.4)
 axis(side=1,at=22+46*c(11,13,15),lab=c(11,13,15),padj=0.2,tcl=0,cex.axis=1.4)
 axis(side=1,at=345,lab="Year after 2000",padj=1.8,tcl=0,cex.axis=1.5)
 axis(side=2,cex.axis=1.5,tcl=0.2)
 tds <- subset(seasonal,loc==locLab[4])
 abline(h=5*c(2:6),col=8)
 abline(v=c(0,46*c(0:15)),col=8)
 text(600,ymax,adj=c(1,0.5),locLab[4],cex=1.5)
 n <- nrow(tds)
 text(690,ymax,adj=c(1,0.5),paste("n: ",n,sep=""),cex=1.5)
 ylab <- expression(paste("LST.Day (",degree,"C)",sep=""))
 mtext(side=3,line=0.2,adj=-0.07,ylab,cex=1.4)

  tds <- subset(tds,!is.na(tds$y1))
  tds$obs <- 46*tds$yr+(tds$day-1)/8
  tds1 <- tds[order(tds$yr,tds$day),]
  tdsa <- tds1$y1-tds1$ss
  tdsa <- tdsa+mean(tds1$y1,na.rm=T)-mean(tdsa,na.rm=T)
  points(tds1$obs,tds1$y1,pch=20,cex=1.3)

  lm(data=tds1,tdsa~obs) -> mod
  summary(mod) -> rez
  range <- max(tds1$obs)-min(tds1$obs)
  incPerDec <- rez$coef[2,1:2]*range/1.5
  fit <- mod$fit+tds1$ss-mean(tds1$ss,na.rm=T)
  points(tds1$obs,fit,type="l",col=2,lwd=2)
  points(tds1$obs,mod$fit,type="l",col=4,lwd=2)

  inc <- round(incPerDec[1],3)
  pval <- rez$coef[2,4]
  pval <- ifelse(pval<0.0005,"<0.0005",round(rez$coef[2,4],3))
  text(10,ymax,adj=c(0,0.5),paste("Inc/dec: ",inc," p-val: ",pval,sep=""),cex=1.5)
  box()

mtext(outer=TRUE,side=3,adj=0.5,line=0.5,title1,cex=1.4)

#-----------------------------------
# Plots of LST_Day by yrDay

seasonal <- subset(seasonal,!is.na(y1))
seasonal$obs <- 0
seasonal$anomaly <- 0
seasonal$id <- 0

ymax <- ymax1+2

windows(16,10)
par(mfrow=c(3,3),las=1,mar=c(0,0.5,0.5,0.3),oma=c(3.2,2.5,3,0.8),tcl=0,mgp=c(1.1,0.1,0))
for (k in c(1:7)) {

 if (k %in% c(3,6)) plot(1,type="n",axes=F)
 plot(1,type="n",ylim=c(ymin,ymax),xlim=c(30,702),xaxt="n",yaxt="n")
  if (k %in% c(3,6,7)) {
   axis(side=1,at=c(22+46*c(0:15)),lab=c(0:15),padj=0.2,tcl=0,cex.axis=1.4)
   axis(side=1,at=22+46*c(10,12,14),lab=c(10,12,14),padj=0.2,tcl=0,cex.axis=1.4)
   axis(side=1,at=345,lab="Year after 2000",padj=1.8,tcl=0,cex.axis=1.5)
  }
 if (k %in% c(1,3,6)) axis(side=2,cex.axis=1.5)
 tds <- subset(seasonal,loc==locLab[k])
 abline(h=5*c(2:6),col=8)
 abline(v=c(0,46*c(0:15)),col=8)
 text(560,ymax,adj=c(1,0.5),locLab[k],cex=1.5)
 n <- nrow(tds)
 text(690,ymax,adj=c(1,0.5),paste("n: ",n,sep=""),cex=1.5)
 ylab <- expression(paste("LST.Day (",degree,"C)",sep=""))
 if (k==1) mtext(side=3,line=0.2,adj=-0.11,ylab,cex=1.2)

 if (n>10) {
  tds <- subset(tds,!is.na(tds$y1))
  tds1 <- tds[order(tds$yr,tds$day),]
  tds1$obs <- 46*tds1$yr+(tds1$day-1)/8
  tdsa <- tds1$y1-tds1$ss
  tdsa <- tdsa+mean(tds1$y1,na.rm=T)-mean(tdsa,na.rm=T)
  points(tds1$obs,tds1$y1,pch=20,cex=1.3)

  lm(data=tds1,tdsa~obs) -> mod
  summary(mod) -> rez
  range <- max(tds1$obs)-min(tds1$obs)
  incPerDec <- rez$coef[2,1:2]*range/1.5
  fit <- mod$fit+tds1$ss-mean(tds1$ss,na.rm=T)
  points(tds1$obs,fit,type="l",col=2,lwd=2)
  points(tds1$obs,mod$fit,type="l",col=4,lwd=2)
  seasonal$obs[seasonal$loc==locLab[k]] <- tds1$obs
  seasonal$id[seasonal$loc==locLab[k]] <- k
  seasonal$anomaly[seasonal$loc==locLab[k]] <- tdsa

  inc <- round(incPerDec[1],3)
  pval <- rez$coef[2,4]
  pval <- ifelse(pval<0.0005,"<0.0005",round(rez$coef[2,4],3))
  text(10,ymax,adj=c(0,0.5),paste("Inc/dec: ",inc," p-val: ",pval,sep=""),cex=1.5)
  box()
 }
}
mtext(outer=TRUE,side=3,adj=0.5,line=0.5,title1,cex=1.4)

library(geepack)

options(scipen=7)

CI <- NULL
ssl1 <- subset(seasonal,yr<5)
ssl1 <- ssl1[order(ssl1$id),]
geeglm(data=ssl1,anomaly~obs,corstr="exchangeable",id=id) -> mod1
summary(mod1) -> rez1
summary(mod1)
CI <- cbind(CI,mod1$coef[2]+1.96*c(-1,0,1)*rez1$coef[2,2])

ssl2 <- subset(seasonal,yr>4 & yr<10)
ssl2 <- ssl2[order(ssl2$id),]
geeglm(data=ssl2,anomaly~obs,corstr="exchangeable",id=id) -> mod2
summary(mod2) -> rez2
summary(mod2)
CI <- cbind(CI,mod2$coef[2]+1.96*c(-1,0,1)*rez2$coef[2,2])

ssl3 <- subset(seasonal,yr>9)
ssl3 <- ssl3[order(ssl3$id),]
geeglm(data=ssl3,anomaly~obs,corstr="exchangeable",id=id) -> mod3
summary(mod3) -> rez3
summary(mod3)
CI <- cbind(CI,mod3$coef[2]+1.96*c(-1,0,1)*rez3$coef[2,2])

CI <- 460*CI			# time scale 460 observations/decade	

ymin <- min(CI)
ymax <- max(CI)

windows(5,4)
par(las=1,mar=c(0,0,0,0),oma=c(2.2,3.5,3.2,0.8),tcl=0.2,mgp=c(1.1,0.1,0))
plot(1,type="n",ylim=c(ymin,ymax),xlim=c(0.5,3.5),xaxt="n",xlab="",ylab="")
abline(h=0,col="red")
abline(v=c(1.5,2.5),col="dimgrey",lty=2)
for (j in c(1:3)) {
 points(c(j,j),CI[c(1,3),j],type="l",lwd=2)
}
points(c(1:3),CI[2,c(1:3)],pch=20,cex=1.4)
axis(side=1,at=c(1:3),lab=c("2000-2004","2005-2009","2010-2015"))
ylab1 <- expression(paste("LST.Day Increase/decade (",degree,"C)",sep=""))
mtext(side=3,adj=-0.25,line=0.1,ylab1,cex=1.1)
mtext(side=3,adj=0.5,line=1.6,title1,cex=1.2)

#-----------------------------------------
# plot NDVI

if (pl=="K") read.csv("NVkhaoleam.csv",header=FALSE,as.is=TRUE) -> V1
if (pl=="M") read.csv("NVmayo.csv",header=FALSE,as.is=TRUE) -> V1
if (pl=="B") read.csv("NVbajo.csv",header=FALSE,as.is=TRUE) -> V1
if (pl=="J") read.csv("NVbogban.csv",header=FALSE,as.is=TRUE) -> V1

dim(V1)

names(V1)[1:7] <- c("rID","fn1","yrDay","latLon","fn4","var","v0") 
V1$yrDay <- substr(V1$yrDay,2,8)

str(V1[,c(1:7)])

V1$yrDay <- as.integer(V1$yrDay)	# convert character to integer

nc <- ncol(V1)      			# number of columns

names(V1)[8:nc] <- paste("v",c(1:(nc-7)),sep="")

vi <- V1[,c(3,8:nc)]
rm(V1)

year <- as.integer(substr(vi$yrDay,1,4))
day <- as.integer(substr(vi$yrDay,5,7))

tab <- addmargins(table(day,year))
tab

if (pl=="K")  days <- 2000000+c(13209,15049,15065)
if (pl=="M")  days <- 2000000+c(13209,13225)
if (pl=="B")  days <- 2000000+c(13209,13225)
if (pl=="J")  days <- 2000000+c(13209,15065)

for (d in days) {
 vi.a <- vi[vi$yrDay==2013193,]
 vi.a$yrDay <- d
 vi.a[,c(2:(nc-6))] <- NA
 vi <- rbind(vi,vi.a)
}

vi <- vi[order(vi$yrDay),]
vi <- subset(vi,vi$yrDay>2000169 & vi$yrDay<2015177)

year <- as.integer(substr(vi$yrDay,1,4))
day <- as.integer(substr(vi$yrDay,5,7))

tab <- addmargins(table(day,year))
tab

vi <- vi[,c(2:(nc-6))]
np <- ncol(vi)				# number of pixels

vi[vi=="F"] <- "-1"

for (j in c(1:np)) {
 vi[,j] <- as.numeric(vi[,j])
}

nt <- nrow(vi)				# number of days
vi <- t(vi)				# transpose
dim(vi) <- c(np*nt,1)			# stack as single column
pixl <- c(1:np)
y <- as.data.frame(rep(pixl,nt))
names(y) <- "pixel"
y$pixel <- as.integer(y$pixel)
yrDay <- 1000*year+day
y$yrDay <- rep(yrDay-2000000,each=np)
y$vi <- as.numeric(vi)
y$vi <- ifelse(y$vi==(-1),NA,y$vi)

y$year <- floor(y$yrDay/1000)
y$day <- y$yrDay-1000*y$year

#---------------------------
# plot of NDVI for central pixel

pixel0 <- (np+1)/2

windows(7,5)
par(las=1,mar=c(0,0.5,0.5,0.3),oma=c(3,2,3,0.8),tcl=0,mgp=c(1.1,0.1,0))

ylab <- "NDVI"

ymax <- 1.1
ymin <- 0

pxl <- pixel0 

vik <- subset(y,pixel==pxl)

plot(vik$day,vik$vi,type="n",ylim=c(ymin,ymax),xaxt="n")

xlab <- c(1+16*c(0:22))
abline(v=xlab,col="grey")
abline(h=0.2*c(0:5),col="grey")
axis(side=1,at=xlab,lab=xlab,padj=0)
axis(side=1,at=182.5,lab="Day of Year",padj=1.6,cex.axis=1.2,tcl=0)
mtext(side=3,line=0.2,adj=-0.07,ylab,cex=1.2)

points(vik$day,vik$vi,pch=20)

n <- sum(ifelse(is.na(vik$vi),0,1))

text(185,ymax,adj=c(0.5,1),"Central Pixel",cex=1.2)
text(365,ymax,adj=c(1,1),paste("n: ",n,sep=""),cex=1.2)

title <- paste("Lat:",lat," Lon:",lon," ",place,sep=" ")
mtext(outer=TRUE,side=3,adj=0.5,line=0.5,title,cex=1.4)

#------------------------------------------------------------
# merge pixels

if (pl=="K") key <- 1
if (pl=="M") key <- 11
if (pl=="B") key <- 6
if (pl=="J") key <- 8

sq.np <- sqrt(np)
px <- rep(1:wid,each=4)[1:sq.np]
px <- rep(px,4)

if (wid==7) {
 px <- c(px,7+px,14+px,21+px,28+px,35+px,42+px)[1:np]
 px <- rep(px,wid)[1:np]
 dim(px) <- c(sq.np,sq.np)
}
if (wid==21) {
 px <- c(px,21+px,42+px,63+px,84+px,105+px,126+px,147+px,168+px,189+px,
	210+px,231+px,252+px,273+px,294+px,315+px,336+px,357+px,378+px,
	399+px,420+px)[1:np]
 px <- rep(px,wid)[1:np]
 dim(px) <- c(sq.np,sq.np)
}
if (key==1) pixlsVT <- px
if (key==2) pixlsVT <- px[c(2:sq.np,sq.np),]
if (key==3) pixlsVT <- px[c(3:sq.np,sq.np,sq.np),]
if (key==4) pixlsVT <- px[c(4:sq.np,sq.np,sq.np,sq.np),]
if (key==5) pixlsVT <- px[,c(2:sq.np,sq.np)]
if (key==6) pixlsVT <- px[c(2:sq.np,sq.np),c(2:sq.np,sq.np)]
if (key==7) pixlsVT <- px[c(3:sq.np,sq.np,sq.np),c(2:sq.np,sq.np)]
if (key==8) pixlsVT <- px[c(4:sq.np,sq.np,sq.np,sq.np),c(2:sq.np,sq.np)]
if (key==9) pixlsVT <- px[,c(3:sq.np,sq.np,sq.np)]
if (key==10) pixlsVT <- px[c(2:sq.np,sq.np),c(3:sq.np,sq.np,sq.np)]
if (key==11) pixlsVT <- px[c(3:sq.np,sq.np,sq.np),c(3:sq.np,sq.np,sq.np)]
if (key==12) pixlsVT <- px[c(4:sq.np,sq.np,sq.np,sq.np),c(3:sq.np,sq.np,sq.np)]
if (key==13) pixlsVT <- px[,c(4:sq.np,sq.np,sq.np,sq.np)]
if (key==14) pixlsVT <- px[c(2:sq.np,sq.np),c(4:sq.np,sq.np,sq.np,sq.np)]
if (key==15) pixlsVT <- px[c(3:sq.np,sq.np,sq.np),c(4:sq.np,sq.np,sq.np,sq.np)]
if (key==16) pixlsVT <- px[c(4:sq.np,sq.np,sq.np,sq.np),c(4:sq.np,sq.np,sq.np,sq.np)]

dim(pixlsVT) <- c(np,1)

y$pixelT <- rep(pixlsVT,nt)

yyt <- merge(y,yt,by.x=c("pixelT","yrDay"),by.y=c("pixlT","yrDay"),all.x=T)

#---------------------------
# plot of NDVI for central pixel with doubtful values displayed

pixel0 <- (np+1)/2

windows(7,5)
par(las=1,mar=c(0,0.5,0.5,0.3),oma=c(3,2,3,0.8),tcl=0,mgp=c(1.1,0.1,0))

ylab <- "NDVI"

ymax <- 1.1
ymin <- 0

pxl <- pixel0 

vik <- subset(yyt,pixel==pxl)
vik <- vik[order(vik$yrDay),]

vik$clr <- ifelse(is.na(vik$td) & !is.na(vik$vi),"cyan3",1)
vik$siz <- ifelse(is.na(vik$td) & !is.na(vik$vi),2,1)

plot(vik$day.x,vik$vi,type="n",ylim=c(ymin,ymax),xaxt="n")

xlab <- c(1+16*c(0:22))
abline(v=xlab,col="grey")
abline(h=0.2*c(0:5),col="grey")
axis(side=1,at=xlab,lab=xlab,padj=0)
axis(side=1,at=182.5,lab="Day of Year",padj=1.6,cex.axis=1.2,tcl=0)
mtext(side=3,line=0.2,adj=-0.07,ylab,cex=1.2)

points(vik$day.x,vik$vi,pch=20,col=vik$clr,cex=vik$siz)

points(vik$day.x,vik$vi,pch=20,cex=1)

vOK <- ifelse(!is.na(vik$vi),1,0)
tOK <- ifelse(!is.na(vik$td),1,0)
n <- table(vOK,tOK)[2,2]

text(185,ymax,adj=c(0.5,1),"Central Pixel",cex=1.2)
text(365,ymax,adj=c(1,1),paste("n: ",n,sep=""),cex=1.2)

title <- paste("Lat:",lat," Lon:",lon," ",place,sep=" ")
mtext(outer=TRUE,side=3,adj=0.5,line=0.5,title,cex=1.4)

#------------------------------------------------------------
# ... with spline function fitted

kn <- c(20,90,150,210,280,350)	# spline function knots
text(kn,0,"+",col="blue",cex=1.2)
   
p <- length(kn)			# number of spline knots
set <- vik$tOK==1 & !is.na(vik$vi)
s <- vik[set,c(4:6)]
s <- s[order(s$day.x),]
x <- s$day.x
yy <- as.data.frame(s[,1:2])
names(yy) <- c("y1","yr")
yy$day <- s$day.x
yy$x <- x
d1 <- kn[p]-kn[p-1]; d2 <- kn[p]-kn[p-2]; d3 <- kn[p]-kn[p-3]
d4 <- kn[p-1]-kn[p-2]; d5 <- kn[p-1]-kn[p-3]; d6 <- kn[p-2]-kn[p-3]
for (j in c(1:(p-4))) {
 sj <- ifelse(x>kn[j],(x-kn[j])^3,0)
 sj <- sj-((kn[p]-kn[j])*(kn[p-1]-kn[j])*(kn[p-2]-kn[j])/(d3*d5*d6))*ifelse(x>kn[p-3],(x-kn[p-3])^3,0)
 sj <- sj+((kn[p]-kn[j])*(kn[p-1]-kn[j])*(kn[p-3]-kn[j])/(d2*d4*d6))*ifelse(x>kn[p-2],(x-kn[p-2])^3,0)
 sj <- sj-((kn[p]-kn[j])*(kn[p-2]-kn[j])*(kn[p-3]-kn[j])/(d1*d4*d5))*ifelse(x>kn[p-1],(x-kn[p-1])^3,0)
 sj <- sj+((kn[p-1]-kn[j])*(kn[p-2]-kn[j])*(kn[p-3]-kn[j])/(d1*d2*d3))*ifelse(x>kn[p],(x-kn[p])^3,0)
 yy[,(j+4)] <- sj
 names(yy)[j+4] <- paste("s",j,sep="")
}
mod1 <- lm(data=yy,y1~s1+s2)		# fit cubic spline model

points(x,mod1$fit,type="l",col=2,lwd=2)
#---------------------------------------------------------------
# ... and using weighted regression

c1 <- -0.02			# zap data with residuals less than c1 in rainy season
c2 <- 0.2			# zap data with residuals greater in magnitude than c2 outside rainy season
e <- 0				# minimum weight
t1 <- 40			# end of rainly season
t2 <- 275			# beginning of rainy season

options(scipen=8)		# show numbers to 8 decimals (not scientific notation)

windows(7,5)
par(las=1,mar=c(0,0.5,0.5,0.3),oma=c(3,2,3,0.8),tcl=0,mgp=c(1.1,0.1,0))

ylab <- "NDVI"

ymax <- 1.1
ymin <- 0

pxl <- pixel0 

vik <- subset(yyt,pixel==pxl)
vik$clr <- ifelse(is.na(vik$td),"cyan3",1)

plot(vik$day.x,vik$vi,type="n",ylim=c(ymin,ymax),xaxt="n")

xlab <- c(1+16*c(0:22))
abline(v=xlab,col="grey")
abline(h=0.2*c(0:5),col="grey")
axis(side=1,at=xlab,lab=xlab,padj=0)
axis(side=1,at=182.5,lab="Day of Year",padj=1.6,cex.axis=1.2,tcl=0)
mtext(side=3,line=0.2,adj=-0.05,ylab,cex=1.2)

points(vik$day.x[!vik$clr==1],vik$vi[!vik$clr==1],pch=20,col="cyan3",cex=2)
points(vik$day.x[vik$clr==1],vik$vi[vik$clr==1],pch=20,cex=1)

kn <- c(20,90,150,210,280,350)		# spline function knots
text(kn,0,"+",col="blue",cex=1.2)
p <- length(kn)				# number of spline knots
set <- vik$tOK==1 & !is.na(vik$vi)
s <- vik[set,c(4:6)]
s <- s[order(s$day.x),]
x <- s$day.x 
yy <- as.data.frame(s[,1:2])
names(yy) <- c("y1","yr")
yy$day <- s$day.x
yy$x <- x
d1 <- kn[p]-kn[p-1]; d2 <- kn[p]-kn[p-2]; d3 <- kn[p]-kn[p-3]
d4 <- kn[p-1]-kn[p-2]; d5 <- kn[p-1]-kn[p-3]; d6 <- kn[p-2]-kn[p-3]
for (j in c(1:(p-4))) {
 sj <- ifelse(x>kn[j],(x-kn[j])^3,0)
 sj <- sj-((kn[p]-kn[j])*(kn[p-1]-kn[j])*(kn[p-2]-kn[j])/(d3*d5*d6))*ifelse(x>kn[p-3],(x-kn[p-3])^3,0)
 sj <- sj+((kn[p]-kn[j])*(kn[p-1]-kn[j])*(kn[p-3]-kn[j])/(d2*d4*d6))*ifelse(x>kn[p-2],(x-kn[p-2])^3,0)
 sj <- sj-((kn[p]-kn[j])*(kn[p-2]-kn[j])*(kn[p-3]-kn[j])/(d1*d4*d5))*ifelse(x>kn[p-1],(x-kn[p-1])^3,0)
 sj <- sj+((kn[p-1]-kn[j])*(kn[p-2]-kn[j])*(kn[p-3]-kn[j])/(d1*d2*d3))*ifelse(x>kn[p],(x-kn[p])^3,0)
 yy[,(j+4)] <- sj
 names(yy)[j+4] <- paste("s",j,sep="")
}
mod1 <- lm(data=yy,y1~s1+s2)		# fit cubic spline model

points(x,mod1$fit,type="l",col=2,lwd=1)	# unweighted model

 wts <- ifelse(mod1$res<c1 & (x<t1 | x>t2),e,ifelse(abs(mod1$res)>c2 & x>t1 & x<t2,e,1))
 mod1 <- lm(data=yy,y1~s1+s2,weights=wts)
 wts <- ifelse(mod1$res<c1 & (x<t1 | x>t2),e,ifelse(abs(mod1$res)>c2 & x>t1 & x<t2,e,1))
 mod1 <- lm(data=yy,y1~s1+s2,weights=wts)
 wts <- ifelse(mod1$res<c1 & (x<t1 | x>t2),e,ifelse(abs(mod1$res)>c2 & x>t1 & x<t2,e,1))
 mod1 <- lm(data=yy,y1~s1+s2,weights=wts)
 wts <- ifelse(mod1$res<c1 & (x<t1 | x>t2),e,ifelse(abs(mod1$res)>c2 & x>t1 & x<t2,e,1))
 mod1 <- lm(data=yy,y1~s1+s2,weights=wts)
 wts <- ifelse(mod1$res<c1 & (x<t1 | x>t2),e,ifelse(abs(mod1$res)>c2 & x>t1 & x<t2,e,1))
 mod1 <- lm(data=yy,y1~s1+s2,weights=wts)
 wts <- ifelse(mod1$res<c1 & (x<t1 | x>t2),e,ifelse(abs(mod1$res)>c2 & x>t1 & x<t2,e,1))
 yy$wts <- wts
 mod1 <- lm(data=yy,y1~s1+s2,weights=wts)

points(x,mod1$fit,type="l",col=2,lwd=2)

ch <- ifelse(wts<1,4,20)
siz <- ifelse(wts<1,1.6,1)
points(yy$day,yy$y1,pch=ch,cex=siz)

n <- sum(ifelse(wts>0.99,1,0))

text(185,ymax,adj=c(0.5,1),"Central Pixel",cex=1.2)
text(365,ymax,adj=c(1,1),paste("n: ",n,sep=""),cex=1.2)

title <- paste("Lat:",lat," Lon:",lon," ",place,sep=" ")
mtext(outer=TRUE,side=3,adj=0.5,line=0.5,title,cex=1.4)

#---------------------------
# Plots of NDVI by day of year

c1 <- -0.02			# zap data with residuals less than c1 in rainy season
c2 <- 0.2			# zap data with residuals greater in magnitude than c2 outside rainy season
e <- 0				# minimum weight

pixels0 <- c(3034,3040,3522,3528)

ymax <- 1.12
ymin <- 0

days <- unique(yyt$day.x)
days <- days[order(days)]
days <- as.data.frame(days)
names(days) <- "day"
seasonal <- NULL

locLab <- c("NW","NE","West","Centre","East","SW","SE")

# need to store spline fits and weights

for (loc in c(1:7)) {
 if (loc==1) pixels <- pixels0-2952
 if (loc==2) pixels <- pixels0-2916
 if (loc==3) pixels <- pixels0-36 
 if (loc==4) pixels <- pixels0
 if (loc==5) pixels <- pixels0+36 
 if (loc==6) pixels <- pixels0+2916 
 if (loc==7) pixels <- pixels0+2952

 v0 <- subset(yyt,pixel %in% pixels)
 v0$vi <- ifelse(is.na(v0$td),NA,v0$vi)

 windows(8,6)
 par(mfrow=c(2,2),las=1,mar=c(0,0.5,0.5,0.3),oma=c(3,2,2,0.6),tcl=0,mgp=c(1.1,0.1,0))
 for (k in c(1:4)) {
  vk <- subset(v0,pixel==pixels[k])
  ymin <- min(0,min(vk$vi,na.rm=TRUE))
  plot(vk$day.x,vk$vi,type="n",ylim=c(ymin,ymax),xaxt="n",yaxt="n")
  if (k>2) {
   xlab <- 1+16*c(0:22)
   axis(side=1,at=xlab,lab=xlab,padj=0)
   axis(side=1,at=182.5,lab="Day of Year",padj=1.6,cex.axis=1.2,tcl=0,font=2)
  }
  if (k %in% c(1,3)) axis(side=2,cex.axis=1.2)
  abline(h=c(ymin,1),col=8)
  if (k==1) mtext(side=3,line=0.2,adj=-0.09,"NDVI",cex=1.2)
  text(1,1.1,adj=c(0,0.5),paste("Pixel ",pixels[k],sep=""),cex=1.2)

  vks <- subset(vk,!is.na(vi))		# omit records where NDVI is unknown

  n <- nrow(vks)
  if (n>100) {				# fit model only if 100+ known NDVI values
   kn <- c(20,90,150,210,280,350)
   text(kn,ymin,"+",col="blue",cex=1.2)
   p <- length(kn)			# number of spline knots
   set <- !is.na(vk$vi)
   s <- vk[set,c(4:6)]
   s <- s[order(s$day.x),]
   x <- s$day.x 
   yy <- as.data.frame(s[,1:2])
   names(yy) <- c("y1","yr")
   yy$day <- s$day.x
   yy$x <- x
   d1 <- kn[p]-kn[p-1]; d2 <- kn[p]-kn[p-2]; d3 <- kn[p]-kn[p-3]
   d4 <- kn[p-1]-kn[p-2]; d5 <- kn[p-1]-kn[p-3]; d6 <- kn[p-2]-kn[p-3]
   for (j in c(1:(p-4))) {
    sj <- ifelse(x>kn[j],(x-kn[j])^3,0)
    sj <- sj-((kn[p]-kn[j])*(kn[p-1]-kn[j])*(kn[p-2]-kn[j])/(d3*d5*d6))*ifelse(x>kn[p-3],(x-kn[p-3])^3,0)
    sj <- sj+((kn[p]-kn[j])*(kn[p-1]-kn[j])*(kn[p-3]-kn[j])/(d2*d4*d6))*ifelse(x>kn[p-2],(x-kn[p-2])^3,0)
    sj <- sj-((kn[p]-kn[j])*(kn[p-2]-kn[j])*(kn[p-3]-kn[j])/(d1*d4*d5))*ifelse(x>kn[p-1],(x-kn[p-1])^3,0)
    sj <- sj+((kn[p-1]-kn[j])*(kn[p-2]-kn[j])*(kn[p-3]-kn[j])/(d1*d2*d3))*ifelse(x>kn[p],(x-kn[p])^3,0)
    yy[,(j+4)] <- sj
    names(yy)[j+4] <- paste("s",j,sep="")
   }
   mod1 <- lm(data=yy,y1~s1+s2)		# fit cubic spline model

   points(x,mod1$fit,type="l",col=2,lwd=1)

    wts <- ifelse(mod1$res<c1 & (x<t1 | x>t2),e,ifelse(abs(mod1$res)>c2 & x>t1 & x<t2,e,1))
    mod1 <- lm(data=yy,y1~s1+s2,weights=wts)
    wts <- ifelse(mod1$res<c1 & (x<t1 | x>t2),e,ifelse(abs(mod1$res)>c2 & x>t1 & x<t2,e,1))
    mod1 <- lm(data=yy,y1~s1+s2,weights=wts)
    wts <- ifelse(mod1$res<c1 & (x<t1 | x>t2),e,ifelse(abs(mod1$res)>c2 & x>t1 & x<t2,e,1))
    mod1 <- lm(data=yy,y1~s1+s2,weights=wts)
    wts <- ifelse(mod1$res<c1 & (x<t1 | x>t2),e,ifelse(abs(mod1$res)>c2 & x>t1 & x<t2,e,1))
    mod1 <- lm(data=yy,y1~s1+s2,weights=wts)
    wts <- ifelse(mod1$res<c1 & (x<t1 | x>t2),e,ifelse(abs(mod1$res)>c2 & x>t1 & x<t2,e,1))
    mod1 <- lm(data=yy,y1~s1+s2,weights=wts)
    wts <- ifelse(mod1$res<c1 & (x<t1 | x>t2),e,ifelse(abs(mod1$res)>c2 & x>t1 & x<t2,e,1))
    mod1 <- lm(data=yy,y1~s1+s2,weights=wts)
    points(x,mod1$fit,type="l",col=2,lwd=2)
   ch <- ifelse(wts<1,4,20)
   siz <- ifelse(wts<1,1.4,1)
   points(yy$day,yy$y1,pch=20)
   points(yy$day,yy$y1,pch=ch,cex=siz)

   n <- sum(ifelse(wts>0.99,1,0))
   text(350,1.1,adj=c(1,0.5),paste("n: ",n,sep=""),cex=1.2)

   ss <- as.data.frame(unique(cbind(x,yy$yr,yy$y1,round(mod1$fit,4),wts)))
   names(ss) <- c("day","yr","y1","ss","wts")
   ss <- merge(days,ss,by.x="day",by.y="day",all=T)
   ss$loc <- loc
   ss$pxl <- pixels[k]
   seasonal <- rbind(seasonal,ss)
  }
 }
 title <- paste("Lat: ",lat," Lon: ",lon," Location: ",locLab[loc],sep="")
 mtext(outer=TRUE,side=3,adj=0.5,line=0,title,cex=1.4)
}
#---------------------------
# Plots of NDVI by yrDay

locLab <- c("NW","NE","West","Centre","East","SW","SE")

seasonal <- subset(seasonal,!is.na(y1))
seasonal$obs <- 0
seasonal$anomaly <- 0
seasonal$id <- 0

for (lc in c(1:7)) {
 if (lc==1) pixels <- pixels0-2952
 if (lc==2) pixels <- pixels0-2916
 if (lc==3) pixels <- pixels0-36 
 if (lc==4) pixels <- pixels0
 if (lc==5) pixels <- pixels0+36 
 if (lc==6) pixels <- pixels0+2916 
 if (lc==7) pixels <- pixels0+2952

 vis <- subset(seasonal,loc==lc)

 if (nrow(vis)>0) {
  windows(10,7)
  par(mfrow=c(2,2),las=1,mar=c(0,0.5,0.5,0.3),oma=c(3,2.5,2,0.6),tcl=0,mgp=c(1.1,0.1,0))
  for (k in c(1:4)) {
   vk <- subset(vis,pxl==pixels[k])
   vk <- subset(vk,!is.na(y1))
   vk$day.x <- vk$day
   vk$obs <- 23*vk$yr+(vk$day.x-1)/16
   vk1 <- vk[order(vk$yr,vk$day.x),]
   vk1$visa <- vk1$y1-vk1$ss
   vk1$visa <- vk1$visa+mean(vk1$y1,na.rm=T)-mean(vk1$visa,na.rm=T)
   plot(vk1$obs,vk1$y1,type="n",ylim=c(ymin,ymax),xaxt="n",yaxt="n")
   abline(v=c(0,23*c(0:15)),col=8)
   if (k %in% c(1,3)) axis(side=2,cex.axis=1.3)
   if (k>2) {
    axis(side=1,at=14+23*c(0:15),lab=c(0:15),padj=0,tcl=0,cex.axis=1.2)
    axis(side=1,at=14+23*c(11,13,15),lab=c(11,13,15),padj=0,tcl=0,cex.axis=1.2)
    axis(side=1,at=175,lab="Year after 2000",padj=1.6,tcl=0,cex.axis=1.4)
   }
   abline(h=c(0,1),col="dimgrey")

   ch <- ifelse(vk1$wts<1,4,20)
   siz <- ifelse(vk1$wts<1,1.4,1)
   points(vk1$obs,vk1$y1,pch=20)
   points(vk1$obs,vk1$y1,pch=ch,cex=siz)

   seasonal$obs[seasonal$pxl==pixels[k]] <- vk1$obs
   seasonal$id[seasonal$pxl==pixels[k]] <- pixels[k]
   seasonal$anomaly[seasonal$pxl==pixels[k]] <- vk1$visa

   lm(data=vk1,visa~obs,weights=wts) -> mod
   summary(mod) -> rez
   nObs <- length(mod$fit)
   incPerDec <- (mod$fit[nObs]-mod$fit[1])/1.5
   incPerObs <- rez$coef[2,1:2]
   se.incPerDec <- incPerObs[2]*incPerDec/incPerObs[1]
   fit <- mod$fit+vk1$ss-mean(vk1$ss,na.rm=T)
   points(vk1$obs,mod$fit,type="l",col=3,lwd=2)
   points(vk1$obs,fit,type="l",col=2,lwd=2)
   inc <- round(incPerDec,3)
   pval <- rez$coef[2,4]
   pval <- ifelse(pval<0.0005,"<0.0005",round(rez$coef[2,4],3))
   nOK <- sum(vk1$wts)
   lab <- paste("Pixel ",pixels[k],"   Inc/dec: ",inc,"  p-val: ",pval,"  n: ",nOK,sep="")
   legend("top",inset=0.005,leg=lab,bg="ivory",cex=1.3,y.intersp=0.2,x.intersp=0)
   if (k==1) mtext(side=3,line=0.2,adj=-0.09,"NDVI",cex=1.2)
   box()
  }
  title <- paste("Lat: ",lat," Lon: ",lon," Location: ",locLab[lc],sep="")
  mtext(outer=TRUE,side=3,adj=0.5,line=0,title,cex=1.4)
 }
}

CI <- NULL
ssl1 <- subset(seasonal,yr<5)
ssl1 <- ssl1[order(ssl1$id),]
geeglm(data=ssl1,anomaly~obs,corstr="exchangeable",id=id) -> mod1
summary(mod1) -> rez1
summary(mod1)
CI <- cbind(CI,mod1$coef[2]+1.96*c(-1,0,1)*rez1$coef[2,2])

ssl2 <- subset(seasonal,yr>4 & yr<10)
ssl2 <- ssl2[order(ssl2$id),]
geeglm(data=ssl2,anomaly~obs,corstr="exchangeable",id=id) -> mod2
summary(mod2) -> rez2
summary(mod2)
CI <- cbind(CI,mod2$coef[2]+1.96*c(-1,0,1)*rez2$coef[2,2])

ssl3 <- subset(seasonal,yr>9)
ssl3 <- ssl3[order(ssl3$id),]
geeglm(data=ssl3,anomaly~obs,corstr="exchangeable",id=id) -> mod3
summary(mod3) -> rez3
summary(mod3)
CI <- cbind(CI,mod3$coef[2]+1.96*c(-1,0,1)*rez3$coef[2,2])

CI <- 230*CI			# time scale 230 observations/decade	

ymin <- min(CI)
ymax <- max(CI)

windows(5,4)
par(las=1,mar=c(0,0,0,0),oma=c(2.2,3.5,3.2,0.8),tcl=0.2,mgp=c(1.1,0.1,0))
plot(1,type="n",ylim=c(ymin,ymax),xlim=c(0.5,3.5),xaxt="n",xlab="",ylab="")
abline(h=0,col="red")
abline(v=c(1.5,2.5),col="dimgrey",lty=2)
for (j in c(1:3)) {
 points(c(j,j),CI[c(1,3),j],type="l",lwd=2)
}
points(c(1:3),CI[2,c(1:3)],pch=20,cex=1.4)
axis(side=1,at=c(1:3),lab=c("2000-2004","2005-2009","2010-2015"))
ylab1 <- "NDVI Increase/decade"
mtext(side=3,adj=-0.25,line=0.1,ylab1,cex=1.1)
mtext(side=3,adj=0.5,line=1.6,title1,cex=1.2)

#-----------------------------------------end of program
