#MODIFIED ON: 12/18/2013            

function_analyses_paper2 <-"multi_timescales_paper_interpolation_functions_12182013.R"

#figure 1: Study area OR

#Figure 2: LST climatology production: daily mean compared to monthly mean

#figure 3: Prediction procedures: direct, CAI and FSS (figure created outside R)

#Figure 4. Accuracy and monthly hold out for FSS and CAI procedures and GWR, Kriging and GAM methods.

#Figure 5. Overtraining tendency, difference between training and testing

#Figure 7: Transect location (transect map)

#Figure 8: Transect profiles for one day 

#Figure 9: Image differencing and land cover: spatial patterns   

#Figure 10: LST and Tmax at stations, elevation and land cover.

#Figure 11: Spatial lag profiles and stations data  

#Figure 12: Daily deviation

#Figure 13: Spatial correlogram at stations, LST and elevation every 5 lags

ref_mod <- 3 #mod3

alt_mod <- 6 #mod6

LC_s <- subset(s_raster,LC_subset)

names(LC_s) <- LC_names

plot(LC_s)

#r_subset_name <- "elev_s"

#r_names <- c("elev_s")

#stat_list_elev <- extract_diff_by_landcover(r_stack_diff,s_raster,LC_subset,LC_names,avl)

#Figure 9: Tmax and LST averagees

names_tmp<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12")

LST_s<-subset(s_raster,names_tmp)

names_tmp<-c("nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08",

             "nobs_09","nobs_10","nobs_11","nobs_12")

LST_nobs<-subset(s_raster,names_tmp)

#LST_nobs<-mask(LST_nobs,LC_mask,filename="test2.tif")

#LST_s<-mask(LST_s,LC_mask,filename="test3.tif")

#c("Jan","Feb")

plot(LST_s)

plot(LST_nobs)

#Map 5: LST and TMax

#note differnces in patternin agricultural areas and 

min_values<-cellStats(LST_s,"min")

max_values<-cellStats(LST_s,"max")

mean_values<-cellStats(LST_s,"mean")

sd_values<-cellStats(LST_s,"sd")

#median_values<-cellStats(molst,"median") Does not extist

statistics_LST_s<-cbind(min_values,max_values,mean_values,sd_values) #This shows that some values are extremes...especially in October

LST_stat_data<-as.data.frame(statistics_LST_s)

rownames(LST_stat_data) <- month.abb

names(LST_stat_data)<-c("min","max","mean","sd")

# Statistics for number of valid observation stack

min_values<-cellStats(LST_nobs,"min")

max_values<-cellStats(LST_nobs,"max")

mean_values<-cellStats(LST_nobs,"mean")

sd_values<-cellStats(LST_nobs,"sd")

LST_nobs_stat_data<-cbind(min_values,max_values,mean_values,sd_values) #This shows that some values are extremes...especially in October

LST_nobs_stat_data <- as.data.frame(LST_nobs_stat_data)

rownames(LST_nobs_stat_data) <- month.abb

#X11(width=12,height=12)

#Plot statiscs (mean,min,max) for monthly LST images

plot(1:12,LST_stat_data$mean,type="b",ylim=c(-15,60),col="black",xlab="month",ylab="tmax (degree C)")

lines(1:12,LST_stat_data$min,type="b",col="blue")

lines(1:12,LST_stat_data$max,type="b",col="red")

text(1:12,LST_stat_data$mean,rownames(LST_stat_data),cex=1,pos=2)

legend("topleft",legend=c("min","mean","max"), cex=1.5, col=c("blue","black","red"),

       lty=1)

title(paste("LST statistics for Oregon", "2010",sep=" "))

#savePlot("lst_statistics_OR.png",type="png")

#Plot number of valid observations for LST

plot(1:12,LST_nobs_stat_data$mean,type="b",ylim=c(0,365),col="black",xlab="month",ylab="tmax (degree C)")

lines(1:12,LST_nobs_stat_data$min,type="b",col="blue")

lines(1:12,LST_nobs_stat_data$max,type="b",col="red")

text(1:12,LST_nobs_stat_data$mean,rownames(LST_stat_data),cex=1,pos=2)

legend("topleft",legend=c("min","mean","max"), cex=1.5, col=c("blue","black","red"),

       lty=1)

title(paste("LST number of valid observations for Oregon", "2010",sep=" "))

#savePlot("lst_nobs_OR.png",type="png")

#met_stations_obj <- load_obj(file.path(in_dir1,met_obj_file_1))

#met_stations_obj$monthly_query_ghcn_data

#

raster_obj <- load_obj(list_raster_obj_files$gam_CAI)

data_month <- (raster_obj$method_mod_obj[[1]]$data_month_s)

dates<-unique(raster_obj$tb_diagnostic_v$date)

#########################

#Figure

#LST for area with FOrest 50> and grass >50

########################

#Figure 

#get the  list of all prediction for a specific method

#list_lf <-extract_list_from_list_obj(load_obj(list_raster_obj_files$gam_CAI)$method_mod_obj,"dailyTmax") #getting objet

list_lf_gam_CAI <-extract_list_from_list_obj(load_obj(list_raster_obj_files[["gam_CAI"]])$method_mod_obj,"dailyTmax") #getting objet

list_lf_gam_fss <-extract_list_from_list_obj(load_obj(list_raster_obj_files[["gam_fss"]])$method_mod_obj,"dailyTmax") #getting objet

#nb_lag <- 10

list_filters<-lapply(1:nb_lag,FUN=autocor_filter_fun,f_type="queen") #generate lag 10 filters

list_param_stat_moran <- list(filter=list_filters[[10]],lf_list=list_lf_gam_CAI)

tt <- stat_moran_std_raster_fun(1,list_param=list_param_stat_moran)

tt <- mclapply(1:11, list_param=list_param_stat_moran, FUN=stat_moran_std_raster_fun,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement

list_param_stat_moran <- list(filter=list_filters[[10]],lf_list=list_lf_gam_fss)

tt_fss <- mclapply(1:365, list_param=list_param_stat_moran, FUN=stat_moran_std_raster_fun,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement

save(tt_fss,file=file.path(out_dir,"moran_std_tt_fss.RData"))

tx<- do.call(rbind,tt_fss)

mo<-as.integer(strftime(dates, "%m"))          # current month of the date being processed

tt_fss2 <- load_obj("moran_std_tt_fss2.RData")

tx<- do.call(rbind,tt_fss2)

dates_proc<-strptime(dates, "%Y%m%d")   # interpolation date being processed

dates_proc <- as.data.frame(dates_proc)

dates_proc$index <- 1:365

tx2 <- merge(tx,dates_proc,by="index")

tx2$month <- as.integer(strftime(tx2$dates_proc, "%m"))          # current month of the date being processed

names(tx2) <- c("index","moranI","std","pred_mod","date","month")

t<-melt(tx2,

          #measure=mod_var, 

          id=c("pred_mod","month"),

          na.rm=F)

#t$value<-as.numeric(t$value) #problem with char!!!

tb_mod_m_avg <-cast(t,pred_mod+month~variable,mean) #monthly mean for every model

#mo<-as.integer(strftime(date_proc, "%m"))          # current month of the date being processed

#day<-as.integer(strftime(date_proc, "%d"))

#year<-as.integer(strftime(date_proc, "%Y"))

# end of pasted

x<-subset(tb_mod_m_avg,pred_mod=="mod7")

plot(x$month,x$moranI,type="b",col="blue")

par(new=TRUE)              # key: ask for new plot without erasing old

plot(x$month,x$std,type="b",col="red",axes=F)

  #axis(4,xlab="",ylab="elevation(m)")  

axis(4,cex=1.2)

#now summarize by month...

#plot(data_month$TMax,add=TRUE)

#list_data_s <-extract_list_from_list_obj(load_obj(list_raster_obj_files$gam_CAI)$validation_mod_obj,"data_s")

#list_data_v <-extract_list_from_list_obj(load_obj(list_raster_obj_files$gam_CAI)$validation_mod_obj,"data_v")

#data_s <- list_data_s[[1]]

#res_m <- data_s$dailyTmax - data_s$mod3

#SSRES_m <- sum((res_m)^2,na.rm=T)

#res_dev <- data_s$dailyTmax - data_s$mod3_del

#SSRES_dev <- sum((res_dev)^2,na.rm=T)

#SST <- sum((data_s$dailyTmax - mean(data_s$dailyTmax))^2,na.rm=T)

#1- (SSRES_m/SST)

#1- (SSRES_dev/SST)

##############

lf_delta_gam_fss <-extract_list_from_list_obj(load_obj(list_raster_obj_files[["gam_fss"]])$method_mod_obj,"delta") #getting objet

test01 <-stack(lf_delta_gam_fss[[287]]) #Ot.14

test0 <-stack(lf_delta_gam_fss[[288]]) #Ot.15

test1 <-stack(lf_delta_gam_fss[[289]]) #Ot.16

test2<-stack(lf_delta_gam_fss[[290]]) #Ot.17

plot(test01)

plot(test0)

plot(test1)

plot(test2)

################################################

#Figure 10: Spatial lag profiles and stations data  

index<-244 #index corresponding to Sept 1 #For now create Moran's I for only one date...

index<-1

date_selected <- "20100901"

date_selected <- "20100101"

#lf

calculate_moranI_profile <- function(nb_lag,lf){

  list_filters<-lapply(1:nb_lag,FUN=autocor_filter_fun,f_type="queen") #generate lag 10 filters

  #moran_list <- lapply(list_filters,FUN=Moran,x=r)

  list_moran_df <- vector("list",length=length(lf))

  for (j in 1:length(lf)){

    r_stack <- stack(lf[[j]])

    list_param_moran <- list(list_filters=list_filters,r_stack=r_stack) #prepare parameters list for function

    #moran_r <-moran_multiple_fun(1,list_param=list_param_moran)

    nlayers(r_stack) 

    moran_I_df <-mclapply(1:nlayers(r_stack), list_param=list_param_moran, FUN=moran_multiple_fun,mc.preschedule=FALSE,mc.cores = 10) #This is the end bracket from mclapply(...) statement

    moran_df <- do.call(cbind,moran_I_df) #bind Moran's I value 10*nlayers data.frame

    moran_df$lag <-1:nrow(moran_df)

    list_moran_df[[j]] <- moran_df

  }

layout_m<-c(2,4) #one row two columns

par(mfrow=layout_m)

LST_s <- subset(s_raster,c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12"))

LST1 <- subset(LST_s,1)

test<-stack(LST_s,elev)

names(test) <- c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12","elev")

x_cor <-layerStats(test,stat="pearson",na.rm=T)

corr(values(LST1),values(elev))

plot(1:12,x_cor[[1]][1:12,13],type="b",ylab="corelation",xlab="month",main="correlation between elevation and LST")

#DATE MODIFIED: 12/16/2013            

#Modfiy...to allow any stat: min, max, mean,sd etc.

stat_moran_std_raster_fun<-function(i,list_param){

  f <-list_param$filter

  lf_list <- list_param$lf_list[[i]]

  list_var_stat<-vector("list",length(lf_list))

    raster_pred<-raster(unlist(lf_list[k])) 

    #tmp_rast<-mask(raster_pred,mask_rast)

    #raster_pred2<-tmp_rast

    t1<-cellStats(raster_pred,na.rm=TRUE,stat=sd)    #Calculating the standard deviation for the     

    m1 <- Moran(x=raster_pred,w=f) #Calculating Moran's I with window of 3 an default Queen's case

    #m1 <- Moran(raster_pred,w=3) #Calculating Moran's I with window of 3 an default Queen's case

  dat_var_stat$index<- i

## Utility function to quickly write out a stack or brick of rasterlayer to disk using names of layers