#layout_m <- c(2,2) # works if set to this?? ok set the resolution...

res_pix<-480

col_mfrow<-2

#row_mfrow<-2

row_mfrow<-2

png_file_name<- paste("Figure_14_paper_variogram_statistics_",out_prefix,".png", sep="")

png(filename=file.path(out_dir,png_file_name),

    width=col_mfrow*res_pix,height=row_mfrow*res_pix)

par(mfrow=c(row_mfrow,col_mfrow))

#png(paste("Figure14_paper_","_variogram_statistics","_",out_prefix,".png", sep=""),

#    height=960*layout_m[1],width=960*layout_m[2])

#mfrow(2,2)

p_hist<-histogram(tb_variogram$model,col="grey",

          main=list(label="Variogram model type",cex=1.8),

          xlab=list(label="Variogram model",cex=1.5))      

print(p_hist)

#hist(tt1$range)

#hist(tt1$psill)

#hist(tt1$Nug)

boxplot(tb_variogram$range~tb_variogram$month,outline=F)

boxplot(tb_variogram$Nug~tb_variogram$month,outline=F)

boxplot(tb_variogram$psill~tb_variogram$month,outline=F)

#boxplot(tb_variogram$model~tb_variogram$month,outline=F)

list_p_mae <- vector("list", length(names_mod))

for (k in 1:length(names_mod)){

  model_name <- names_mod[k]

  tt <- na.omit(data_v_mae)

  tt <- data_v_mae[data_v_mae[res_model_name],]

  p2 <- bubble(data_v_mae,res_model_name, main=paste("MAE per station for ",model_name,sep=""),

               col=matlab.like(25),na.rm=TRUE)

  print(p3)

png(paste("Figure13_paper_","_residuals_",date_selected,"_",out_prefix,".png", sep=""),

grid.arrange(list_p_mae[[1]],list_p_mae[[2]],list_p_mae[[5]],ncol=3)

infile_reg_outline <- "/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/OR83M_state_outline.shp"  #input region outline defined by polygon: Oregon

reg_outline <- readOGR(dsn=dirname(infile_reg_outline),layer=sub(".shp","",basename(infile_reg_outline)))

p + layer(sp.polygons(reg_outline,lwd=0.9,col="darkgray"))

p4 <- bubble(data_v_mae,"res_mod4",maxsize=4,col=c("red"),fill=FALSE)

p4 + layer(sp.polygons(reg_outline,lwd=0.9,col="darkgray"))

col_t <- colorRampPalette(c('blue', 'white', 'red'))

p_elev <-levelplot(subset(s_raster,"elev_s"),margin=FALSE)

p4 <-bubble(data_v_mae[data_v_mae$res_mod4>2.134,],"res_mod4",maxsize=4,col=c("blue"),fill=FALSE)

p_elev + p4 + layer(sp.polygons(reg_outline,lwd=0.9,col="green"))

title("mod4")

p_elev <-levelplot(subset(s_raster,"elev_s"))

p1 <-bubble(data_v_mae[data_v_mae$res_mod1>2.109,],"res_mod1",maxsize=4,col=c("blue"),fill=FALSE)

p_elev + p1 + layer(sp.polygons(reg_outline,lwd=0.9,col="green"))

#bubble(data_v_mae,"res_mod1")

#p<-spplot(data_v_mae,"res_mod1",maxsize=4,col=c("red"))

#p

#stations that are outliers in one but not the other...

id_setdiff<-setdiff(data_v_mae[data_v_mae$res_mod1>2.109,]$id,data_v_mae[data_v_mae$res_mod4>2.134,]$id)

data_id_setdiff <- data_v_mae[data_v_mae$id %in% id_setdiff,]

p_elev +layer(sp.polygons(reg_outline,lwd=0.9,col="green")) + layer(sp.points(data_id_setdiff,pch=4,cex=2,col="pink"))

list_data_v <- lapply(1:365,function(i){raster_prediction_obj_2$validation_mod_obj[[i]]$data_v} #testing with residuals

brks<-c(0,500,1000,1500,2000,2500,4000)

lab_brks<-1:6

boxplot(data_v_ag$res_mod2~elev_rcstat,ylim=c(-15,15))           

boxplot(data_v_ag$res_mod1~elev_rcstat,ylim=c(-15,15))                      

boxplot(data_v_ag$res_mod5~elev_rcstat,ylim=c(-15,15))                      

brks<-c(-20,-10,0,10,20,30,40)

lab_brks<-1:6

rcstat<-cut(data_v_ag$LST,breaks=brks,labels=lab_brks,right=F)

boxplot(data_v_ag$res_mod4 ~ rcstat,ylim=c(-15,15))           

y_range<-range(c(diff_fc))

x_range<-range(c(elev_rcstat))

plot(elev_rcstat,res_mod, ylab="diff_cf", xlab="ELEV_SRTM (m) ",

       ylim=y_range, xlim=x_range)

text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3)

grid(lwd=0.5,col="black")

title(paste("Testing stations residuals fusion vs Elevation",date_selected,sep=" "))

lat<- aggregate(lat~id,data=data_v_ag,FUN=mean)

lon<- aggregate(lon~id,data=data_v_ag,FUN=mean)                     

data_v_mae <- data.frame(lon=lon$lon,lat=lat$lat)

coordinates(data_v_mae)<- cbind(lon$lon,lat$lat)

mae_fun<-function(x){mean(abs(x))}                      

test<- aggregate(res_mod1~id,data=data_v_ag,FUN=mae_fun)

data_v_mae$mod1 <- test$res_mod1  

test<- aggregate(res_mod2~id,data=data_v_ag,FUN=mae_fun)

data_v_mae$mod2 <- test$res_mod2   

p<- bubble(data_v_mae,"mod1")

print(p)                    

p1 <- levelplot(elev,scales = list(draw = FALSE), colorkey = FALSE,par.settings = GrTheme)

hist(data_v_mae$mod1)

p<- bubble(data_v_mae,"mod1")

spplot()

p_elev <- levelplot(elev)

pc <-p+p1 + p

print(pc)

t<-melt(data_v_ag,

        measure=names_var, 

        id=c("res_mod1","res_mod2,"id"),

        na.rm=T)

n_tb<-cast(t,res_mod1_rc1+id~variable,mae_fun)

n_tb_tot <-cast(t,id~variable,length) #number of times the stations was used for validation

###################### END OF SCRIPT #######################