#system("wget -r --retr-symlinks ftp://ladsweb.nascom.nasa.gov/orders/500676499/ -P /home/adamw/acrobates/projects/interp/data/modis/MOD06_L2_hdf")

#tfs=fs2$path[which(fs2$month==vs$month[i]&fs2$type==vs$type[i])]

#do.call(c,mclapply(2:length(tfs),function(f) #identical(extent(raster(tfs[f-1])),extent(raster(tfs[f])))))

#raster(tfs[23])

library(raster); library(rasterVis)

out_prefix<-"20120723PRCP"

### load ROI

roi <- readOGR("../Test_sites/","Oregon")

roi=spTransform(roi,CRS("+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs"))

ghcn<-readOGR(".", filename)                  #reading shapefile

ghcn@data[,c("x","y")]=coordinates(ghcn)

mods=data.frame(

  formula=c(    

    "value ~ s(x_OR83M,y_OR83M)", 

#    "value ~ s(x_OR83M,y_OR83M) + s(distoc) + elev",

    "value ~ s(x_OR83M,y_OR83M) + s(distoc) + elev + ns + ew",

    "value ~ s(distoc) + s(CER_P20um) + elev + ns + ew",

    "value ~ s(distoc) + s(COT_mean) + elev + ns + ew",

    "value ~ s(distoc) + s(CLD_mean) + elev + ns + ew",

    "value ~ s(CLD_mean) + elev + ns + ew",

    "value ~ s(COT_mean) + elev + ns + ew",

    "value ~ s(CER_P20um) + elev + ns + ew",

    "value ~ s(CER_mean) + elev + ns + ew"

                                        #    "value ~ s(x_OR83M,y_OR83M) + s(distoc) + elev + ns + ew + s(CER_P20um)",

#    "value ~ s(x_OR83M,y_OR83M,CER_P20um) +s(x_OR83M,y_OR83M,CLD_mean) + elev + ns + ew",

#    "value ~ s(x_OR83M,y_OR83M) + s(CER_P20um,CLD_mean) + elev + ns + ew",

#    "value ~ x_OR83M + y_OR83M +s(CER_P20um) + elev + ns + ew",

#    "value ~ s(x_OR83M,y_OR83M,CER_P20um) + elev + ns + ew",

#    "value ~ s(x_OR83M,y_OR83M) + s(distoc,CER_P20um)+elev + ns + ew"

#    "value ~ s(x_OR83M,y_OR83M) + s(distoc) + elev + ns + ew + s(CLD_mean)",

#    "value ~ s(x_OR83M,y_OR83M) + s(distoc) + elev + ns + ew + s(COT_mean)"

    ),stringsAsFactors=F)

if(any(terms%in%layerNames(covar)))

i=1

m=3

savemodel=T;saveFullPrediction=T;scale=F;verbose=T

  results=do.call(rbind.data.frame,                   # Collect the results in a single data.frame

   lapply(1:length(dates),function(i,savemodel=T,saveFullPrediction=T,scale=F,verbose=T) {            # loop over dates

     if(verbose)      print(paste("Starting Date:",dates[i]))

     ## drop stations with no x coordinates (these are stations in the buffer region)

     tdata=tdata[!is.na(tdata$x_OR83M),]  ## REMOVE THIS WHEN FULL REGION IS INCLUDED

     #######

     if(saveFullPrediction){

       if(verbose) print(paste("Extracting covar information for this date (",dates[i],")"))

       ## subset full raster brick to include correct month of satellite data

       covarm=subset(covar,subset=which(covar@z$month=="00"|covar@z$month==month))

       covarm=setZ(covarm,getZ(covar)[which(covar@z$month=="00"|covar@z$month==month)])

       ## update layer names to match those in ghcn table

       layerNames(covarm)[getZ(covarm)!="00"]=sub(paste("_",month,sep=""),"",layerNames(covarm)[getZ(covarm)!="00"])

       ## convert to matrix for prediction later

       pred=extract(covarm,coordinates(covarm))

     }

     if(scale){

       ## get list of variables that need to be scaled to assess relative importance

       scalevars=!colnames(tdata)%in%c("station","month","count",

                                       "value","latitude","longitude","date","weights")

       tdata[,colnames(tdata)[scalevars]]=scale(tdata[,scalevars])

     }

     ## Create subsets for fitting and validation

     n = nrow(tdata)

     ns = round(n*prop)  #Create a sample from the data frame with 70% of the rows

     nv = n-ns             #create a sample for validation with prop of the rows

     ind.training = sample(nrow(tdata), size=ns, replace=FALSE) #This selects the index position for 70% of the rows taken randomly

     tdata$fit=F;tdata$fit[ind.training]=T   #add column to tdata with fit/training status

     tdata$class=as.factor(ifelse(tdata$fit,"Fitting","Validation"))

     ## lapply loops through models for the ith day, fits models, calculates the validation metrics, and saves them as R objects (if desired)

     if(verbose) print(paste("Starting the model fitting for this date (",dates[i],")"))

       lapply(1:nrow(mods),function(m,data=tdata[tdata$fit],...) {

         if(verbose) print(paste("Fitting Model:",mods$model[m]))

         err=try(gam(formula(mods$formula[m]),data=tdata[tdata$fit,]),silent=T)

         ## define the spatial covariance structure

#         library(geoR)

#         plot(variog(coords=as.matrix(data_s[,c("x","y")]),data=tdata$value[tdata$fit]))

#         plot(variofit(variog(coords=as.matrix(tdata[tdata$fit,c("x","y")]),data=tdata$value[tdata$fit]),cov.model="exponential"))

#         vg=likfit(geodata=list(coords=tdata[tdata$fit,c("x","y")],data=tdata$value[tdata$fit]),cov.model="exponential",ini.cov.pars=expand.grid(seq(0.1,2,len=5),seq(1000,500000,len=5)))

#         lines.variomodel(vg)

         corr=corExp(c(100000,0.1),form=~x+y,nugget=T,fixed=F)

#         corr=Initialize(corr,tdata[tdata$fit,])

         ## run the model

         mod=gamm(formula(mods$formula[m]),data=tdata[tdata$fit,],correlation=corr,weights=weights,method="ML")

#         mod=gam(formula(mods$formula[m]),data=tdata[tdata$fit,],weights=tdata$weights[tdata$fit])

         ### generate knots

#         knots=makegrid(bbox(ghcn),cellsize=100000)

#         coordinates(knots)=c("x1","x2")

#         proj4string(knots)=proj4string(ghcn)

#        niter=1500

#         t1=Sys.time()

#         f1=formula(value~CER_P20um+elevation+ns+ew)

#         m1=spLM(f1,data=tdata[tdata$fit,],coords=as.matrix(tdata[tdata$fit,c("x","y")]), 

#           starting=list("phi"=2.5,"sigma.sq"=4, "tau.sq"=1,"nu"=0.5),

#           sp.tuning=list("phi"=0.8, "sigma.sq"=0.1, "tau.sq"=0.2,"nu"=0.5),

#           priors=list("phi.Unif"=c(0.1, 10), "sigma.sq.IG"=c(2, 1),

#             "tau.sq.IG"=c(2, 1),"nu.unif"=c(0.01,3)),

#           cov.model="matern",

#           knots=coordinates(knots),

#           n.samples=niter, verbose=TRUE, n.report=100,sub.samples=c(500,niter,1))

#         t2=Sys.time()

#         y_mod<- spPredict(m1, tdata[,c("x","y")],pred.covars=model.matrix(lm(f1,data=tdata))) #Using the coeff to predict new values.

#         y_mod<- spPredict(m1, pred[,c("x_OR83M","y_OR83M")],pred.covars=model.matrix(lm(f1,data=data.frame(value=0,pred)))) #Using the coeff to predict new values.

#         y_mod$fit=apply(y_mod$y.pred,1,mean)

#         tdata$pred=y_mod$fit

         tpred=predict(mod$gam, newdata=tdata, se.fit = TRUE) #Using the coeff to predict new values.

         tdata$pred=as.vector(tpred$fit)

         tdata$pred.se=as.vector(tpred$se)

#         y_mod$fit[y_mod$fit<0]=NA

         tdata$resid<- itrans(tdata$value) - itrans(tdata$pred) #Residuals 

#           AIC=AIC(mod$gam),                # AIC

#           GCV=mod$gam$gcv.ubre,             # GCV

#           DEV=mod$gam$deviance,             # Deviance

           ME=mean(tdata$resid[!tdata$fit],na.rm=T),  # Mean error

           MAE=mean(abs(tdata$resid[!tdata$fit]),na.rm=T),  # Mean absolute error

           pME=mean((itrans(tdata$pred[!tdata$fit])/itrans(tdata$value[!tdata$fit]))[itrans(tdata$pred[!tdata$fit])/itrans(tdata$value[!tdata$fit])<Inf],na.rm=T),  # Mean % of total

           RMSE=sqrt(sum(tdata$resid[!tdata$fit]^2,na.rm=T)/nv),  # RMSE

           R2=summary(lm(itrans(tdata$pred[!tdata$fit])~itrans(tdata$value[!tdata$fit])))$r.squared,  # R^2

           R2w=summary(lm(itrans(tdata$pred[!tdata$fit])~itrans(tdata$value[!tdata$fit]),weights=tdata$weights[!tdata$fit]))$r.squared,

           R2lw=summary(lm(tdata$pred[!tdata$fit]~tdata$value[!tdata$fit],weights=tdata$weights[!tdata$fit]))$r.squared)

              if(verbose) print(tresults)

         write.csv(tresults,file=paste(path,"/",out_prefix,"_",gsub("-","",date),"_",mods$model[m],".csv",sep=""))

         if(savemodel)  save(mod,file= paste(path,"/",out_prefix,"_",gsub("-","",date),"_",mods$model[m],".Rdata",sep=""))

           if(verbose) print(paste("Doing predictions for model ",mods$model[m]," for this date (",dates[i],")"))

           ## do the prediction

           p1=predict(mod$gam,as.data.frame(pred),type="response",se.fit=T,block.size=10000)

           ## convert to raster stack (mean and se)

           p2=stack(SpatialPixelsDataFrame(coordinates(covarm),

             data=data.frame(

               pred=itrans(as.numeric(p1$fit)),

               se=itrans(as.numeric(p1$se.fit)))))

            writeRaster(p2,filename=paste(path,"/",out_prefix,"_",gsub("-","",date),"_",

                            mods$model[m],"_prediction.tif",sep=""),

                     format="GTiff",overwrite=T)

           ## draw map of full prediction

            ncols=100

           c.at=unique(quantile(ghcn$value,seq(0,1,len=ncols)))  #get quantile bins of raw data for all months

           c.ramp=colorRampPalette(c("brown","grey","lightgreen","darkgreen"))  #assign color ramp

           c.cols=c.ramp(length(c.at)-1)  #get ncol-1 colors to corresponding to c.at

           ## define se colors

           s.at=unique(quantile(subset(p2,"se"),seq(0,1,len=ncols),na.rm=T))  #get quantile bins of raw data for all months

           s.ramp=colorRampPalette(c("blue","red"))  #assign color ramp

           s.cols=s.ramp(length(s.at)-1)  #get ncol-1 colors to corresponding to c.at

           ## define residual color maps

           r.at=c(-Inf,-20,-5,-1,1,5,20,Inf)

           r.ramp=colorRampPalette(c("blue","darkblue","grey","darkred","red"))

           r.cols=r.ramp(length(r.at)-1)

           r.cut=cut(tdata$resid,breaks=r.at)

           ## plot the estimated values

           pdf(paste(path,"/",gsub("-","",date),"_",mods$model[m],"_prediction.pdf",sep=""),width=11,height=8.5)

           print(levelplot(subset(p2,"pred"),col.regions=c.ramp(length(c.at)),at=c.at,main=paste("Mean Predictions  (Month",m,")"),

                           sub=paste("Model:",mods$formula[m],"\n  Circles indicate training stations, triangles indicate validation stations"),margin=F)+

                 layer(panel.points(tdata$x_OR83M,tdata$y_OR83M,

                                    fill=as.character(cut(itrans(tdata$value),breaks=c.at,labels=c.cols)),col="black",

                                    pch=21,cex=tdata$weights),data=list(tdata=tdata[tdata$fit,],c.at=c.at,c.cols=c.cols))+

                 layer(panel.points(tdata$x_OR83M,tdata$y_OR83M,

                                    fill=as.character(cut(itrans(tdata$value),breaks=c.at,labels=c.cols)),col="black",pch=24,cex=tdata$weights),

                       data=list(tdata=tdata[!tdata$fit,],c.cols=c.cols,c.at=c.at))+

                 layer(sp.lines(as(roi,"SpatialLinesDataFrame"),col="black")))

           ## plot of prediction errors

           print(levelplot(subset(p2,"se"),col.regions=s.ramp(length(s.at)),at=s.at,main=paste("Prediction Standard Errors (Month",m,")"),

                     sub=paste("Model:",mods$formula[m],"\n  Circles indicate training stations, triangles indicate validation stations"),margin=F)+

                       layer(panel.points(tdata[,c("x_OR83M","y_OR83M")],fill=as.character(cut(itrans(tdata$pred.se),breaks=s.at,labels=s.cols)),col="black",pch=21),

                             data=list(tdata=tdata[tdata$fit,],s.at=s.at,s.cols=s.cols))+

               layer(panel.points(tdata[,c("x_OR83M","y_OR83M")],fill=as.character(cut(itrans(tdata$pred.se),breaks=s.at,labels=s.cols)),col="black",pch=24),

                     data=list(tdata=tdata[!tdata$fit,],s.cols=s.cols,s.at=s.at))+

                 layer(sp.lines(as(roi,"SpatialLinesDataFrame"),col="black")))

           ## plot of residuals

           print(levelplot(subset(p2,"pred"),col.regions=c.ramp(length(c.at)),at=c.at,main=paste("Mean Predictions with station residuals (month",m,")"),

                     sub=paste("Model:",mods$formula[m],"\n  Circles indicate training stations, triangles indicate validation stations"),margin=F,

                     key=list(text=list(rev(levels(r.cut)),col=rev(r.cols)),points=list(pch=21,fill=rev(r.cols)),space="left",title="Station \n Residuals"))+

             layer(panel.points(tdata[,c("x_OR83M","y_OR83M")],

                                fill=as.character(cut(itrans(tdata$resid),breaks=r.at,labels=r.cols)),col="black",pch=21,cex=tdata$weights),

                                        data=list(tdata=tdata[tdata$fit,],r.cols=r.cols,r.at=r.at))+

               layer(panel.points(tdata[,c("x_OR83M","y_OR83M")],

                                  fill=as.character(cut(itrans(tdata$resid),breaks=r.at,labels=r.cols)),col="black",pch=24,cex=tdata$weights),

                     data=list(tdata=tdata[!tdata$fit,],r.cols=r.cols,r.at=r.at))+

                 layer(sp.lines(as(roi,"SpatialLinesDataFrame"),col="black")))

          #         pdata=data.frame(Predicted=apply(itrans(as.matrix(y_mod$fit)),1,function(x) max(x,0)),Observed=itrans(data_v$value))

           print(xyplot(itrans(pred)~itrans(value),groups=class,ylab="Predicted",xlab="Observed",main=mods$formula[m],data=tdata,asp=.5,auto.key=T)+layer(panel.abline(0,1,col="red")))

 #          xyplot(pred~value,groups=class,ylab="Predicted",xlab="Observed",main=mods$formula[m],data=tdata,asp=.5,auto.key=T)+layer(panel.abline(0,1,col="red"))

           plot(mod$gam,pages=1,residuals=T,scale=-1,pers=T,shade=T,seWithMean=T); mtext(mods$formula[m],3,padj=-2)

           dev.off()

       }))  #end of loop for this time step

     ## identify model with minimum AIC

  write.csv(results, file= paste(path,"/",out_prefix,"_results.csv",sep=""))

## some plots

myTheme <- standard.theme(col = FALSE)

myTheme$superpose.symbol$pch <-1:20

myTheme$superpose.symbol$col=1:12

pdf(paste(path,"/",out_prefix,"_ModelSummary.pdf",sep=""),width=11,height=8.5)

for(i in c("CER_mean","CER_P20um","CLD","COT")) {

  maps=subset(covar,grep(i,layerNames(covar)))

  ncols=100

  c.at=unique(quantile(as.matrix(maps),seq(0,1,len=ncols)))  #get quantile bins of raw data for all months

  c.ramp=colorRampPalette(c("brown","grey","lightgreen","darkgreen"))  #assign color ramp

  print(levelplot(maps,col.regions=c.ramp(length(c.at)-1),at=c.at,main=i)+layer(sp.lines(as(roi,"SpatialLinesDataFrame"),col="black")))

}

## histogram of observed values by month

histogram(~value|format(as.Date(date),"%m"),data=ghcn@data,points=F,auto.key=list(space="right"),

          breaks=seq(0,50,len=50),fill="grey",col="grey",border="transparent",scales=list(x=list(log=F),y=list(relation="free",crt=90)),

          main="Histograms of observed mean daily precipitation, by month",xlab="Mean Daily Precipitation (mm)")

## plot scatterplots by month

for(i in sprintf("%02d",1:12))   print(splom(subset(covar,grep(paste("[CER|CLD|COT].*",i,sep=""),layerNames(covar),value=T)),main=paste("Month",i),

                                       sub="MOD06 cloud product climatologies for a single month"))

xyplot(value~date|variable,groups=form,data=resl,scales=list(y=list(relation="free")),auto.key=T,par.settings = myTheme)

resl2=cast(resl,date+month+variable~model)

myTheme$superpose.symbol$col=c("blue",rep("green",3),rep("red",3),rep("purple",3),rep("blue",2))

xyplot(mod6~mod8|variable,groups=month,data=resl2,scales=list(relation="free"),auto.key=list(space="right"),par.settings = myTheme)+layer(panel.abline(0,1))

dev.off()

print(xtable(mods),include.rownames=F,file=paste(path,"/",out_prefix,"_Models.tex",sep=""))

xyplot(form~value|variable,groups=month,data=resl,scales=list(x=list(relation="free"),y=list(draw=T, alternating=1)),auto.key=list(space="right"),par.settings = myTheme)