#####################################  METHODS COMPARISON ##########################################

#This script is not aimed at producing new interpolation surfaces. It utilizes the R ojbects created 

# during the interpolation phase.                       #

# At this stage the script produces figures of various accuracy metrics and compare methods:       #

#- spatial density of station network and accuracy metric 

#- visualization of maps of prediction and difference for comparison 

#DATE: 10/30/2012                                                                                 #

#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491 --                                    #

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

#### PART II EXAMINIG PREDICTIONS AND RESIDUALS TEMPORAL PROFILES ##########

ghcn_value<-cbind(ghcn_locs,ghcn_value[,1:365]) #This data frame contains values for 365 days

                                                #for 186 stations of year 2010...

#id<-c("USW00024284","USC00354126","USC00358536","USC00354835",

stat_list<-vector("list",3 )  #list containing the selected stations...

  y1<-as.numeric(ghcn_fus_pred[m1,6:ncol(ghcn_fus_pred)]) #vector containing fusion time series of predictecd tmax

  y2<-as.numeric(ghcn_cai_pred[m2,6:ncol(ghcn_cai_pred)]) #vector containing CAI time series of predictecd

  y3<-as.numeric(ghcn_value[m3,6:ncol(ghcn_value)])  #vector containing observed time series of predictecd

  res2<-y2-y3 #CAI time series residuals

  res1<-y1-y3 #fusion time series residuals

  plot(x,y1,type="l",col="red",ylab="tmax (C)",xlab="Day of year")

  lines(x,y3,col="green")

  legend("topright",legend=c("fus","CAI","OBS"), cex=1.2, col=c("red","blue","green"),

  ### RESIDUALS PLOT

  plot(x,res2,type="l",col="blue", ylab="tmax (C)",xlab="Day of year") #res2 contains residuals from cai

  lines(x,res1,col="red")        #res1 contains fus

  lines(x,zero,col="green")      

  legend("topright",legend=c("fus","CAI"), cex=1.2, col=c("red","blue"),

         lty=1)

  title(paste("temporal profile for station ", id[i],sep=""))

### RESIDUALS FOR EVERY STATION...############

########### TRANSECT-- DAY OF YEAR PLOT...#########

id<-c("USW00024284","USC00354126","USC00358536","USC00354835",

"USC00356252","USC00359316","USC00358246","USC00350694",

"USC00350699","USW00024230","USC00353542")

id_order<-1:11

m<-match(id,ghcn_locs$station)

dat_id<-ghcn_locs[m,]  #creating new subset

#dat_id<-subset(ghcn_locs[gj])

filename<-sub(".shp","",infile6)             #Removing the extension from file.

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

X11()

s.range <- c(min(minValue(mm_01)), max(maxValue(mm_01)))

col.breaks <- pretty(s.range, n=50)

lab.breaks <- pretty(s.range, n=5)

temp.colors <- colorRampPalette(c('blue', 'white', 'red'))

plot(mm_01, breaks=col.breaks, col=temp.colors(length(col.breaks)-1),

     axis=list(at=lab.breaks, labels=lab.breaks))

plot(reg_outline, add=TRUE)

plot(dat_id,cex=1.5,add=TRUE)

title("Selected stations for comparison",line=3)

title("(Background: mean January LST)", cex=0.5, line=2)

coords<-coordinates(dat_id)

text(x=coords[,1],y=coords[,2],labels=as.character(id_order),cex=1.5, adj=c(0,1),offset=2) #c(0,1) for lower right corner!

savePlot(paste("temporal_profile_station_locations_map",out_prefix,".png", sep=""), type="png")

dev.off()

date_selection<-c("01-01-2010","01-09-2010")

date_str<-c("20100101","20100901")

  doy<-243+5

  plot(1:11,stat_subset$value,type="b",col="green",ylab="tmax",xlab="station transtect number")

  lines(1:11,stat_subset$fus,type="b",col="red")

  lines(1:11,stat_subset$cai,type="b",col="blue")

  legend("bottomright",legend=c("obs","fus","cai"), cex=1.2, col=c("green","red","blue"),

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

########## USING TEMPORAL IMAGES...############

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

#For all stations

ghcn$station

# For validation and training stations...

unique(ghcn$station)

############ GRID BOX AVERAGING ####################

####### Create the averaged grid box...##############

ghcn_sub<-as.data.frame(subset(ghcn, select=c("station","x_OR83M","y_OR83M")))

ghcn_sub_melt<-melt(ghcn_sub,

                    measure=c("x_OR83M","y_OR83M"), 

                    id=c("station"),

                    na.rm=F)

ghcn_stations<-as.data.frame(cast(ghcn_sub_melt,station~variable,mean))

coords<- ghcn_stations[,c('x_OR83M','y_OR83M')]

coordinates(ghcn_stations)<-coords

proj4string(ghcn_stations)<-proj_str  #Need to assign coordinates...

oc_all<-vector("numeric",nrow(ghcn_stations))

oc_all<-oc_all+1

ghcn_stations$oc_all<-oc_all

rast_oc_all<-rasterize(ghcn_stations,rast_agg,"oc_all",na.rm=TRUE,fun=sum)

ac_agg50$oc_all<-values(rast_oc_all)

plot(rast_oc_all, main="Number of stations in coarsened 50km grid")

plot(reg_outline, add=TRUE)

fdens_all50<-as.data.frame(freq(rast_oc_all))

tot50<-sum(fdens_all50$count[1:(nrow(fdens_all50)-1)])

percent<-(fdens_all50$count/tot50)*100

percent[length(percent)]<-NA

fdens_all50$percent<-percent

  #ghcn_sub<-as.data.frame(subset(ghcn, select=c("station","x_OR83M","y_OR83M")))

  #ghcn_sub_melt<-melt(ghcn_sub,

  #                    measure=c("x_OR83M","y_OR83M"), 

  #                    id=c("station"),

  #                    na.rm=F)

  #ghcn_stations<-as.data.frame(cast(ghcn_sub_melt,station~variable,mean))

  #coords<- ghcn_stations[,c('x_OR83M','y_OR83M')]

  #coordinates(ghcn_stations)<-coords

  #proj4string(ghcn_stations)<-proj_str  #Need to assign coordinates...

  #oc_all<-vector("numeric",nrow(ghcn_stations))

  #oc_all<-oc_all+1

  #ghcn_stations$oc_all<-oc_all

  #rast_oc_all<-rasterize(ghcn_stations,rast_agg,"oc_all",na.rm=TRUE,fun=sum)

  #ac_agg50$oc_all<-values(rast_oc_all)

  plot(rast_oc_t, main="Number of training stations in coarsened 50km grid")

names_statistics<-c("mean","sd","min","max")

stat_val_m<-matrix(NA,nrow=length(date_list2),ncol=length(names_statistics))

colnames(stat_val_m)<-names_statistics

rownames(stat_val_m)<-date_list2

stat_val_m<-as.data.frame(stat_val_m)

  raster_pred1<-raster(l_f[[index]])  # Fusion image

  raster_pred2<-raster(l_f2[[index]]) # CAI image

  stat_val_m$mean[i]<-cellStats(raster_pred1,na.rm=TRUE,stat=mean)    #Calculating the standard deviation for the 

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

  t2<-cellStats(raster_pred2,na.rm=TRUE,stat=mean)    #Calculating the standard deviation for the 

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

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

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

  t2<-cellStats(raster_pred2,na.rm=TRUE,stat=min)    #Calculating the standard deviation for the 

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

  t2<-cellStats(raster_pred2,na.rm=TRUE,stat=max)    #Calculating the standard deviation for the 

  hist(predictions,freq=FALSE,maxpixels=ncells(predictions),xlabel="tmax (C)")

  X11(6,12)

  title(paste("Difference between CAI and fusion for ",date_list2[[k]],sep=""))

  savePlot(paste("comparison_diff_CAI_fusion_tmax_prediction_",date_list2[[k]],out_prefix,".png", sep=""), type="png")

  dev.off()

}  

symbols(c(2e5, 4e5), c(2e5, 3e5), circles=rep(2e4, 2), inches=FALSE, add=TRUE)

text(c(2e5, 4e5), c(2e5, 3e5), labels=1:2,cex=0.5)

points(coordinates(pts), type="c")

text(coordinates(pts), labels=9:11, cex=0.8)

points(coordinates(pts), type="b", pch=as.character(1:length(pts))

       points(coordinates(pts), type="b", pch=as.character(9:11)

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