library(ncf)

#

#

#ghcn_day_dat <- readOGR(dsn=dirname(met_stations_obj$daily_covar_ghcn_data),

#                    sub(".shp","",basename(met_stations_obj$daily_covar_ghcn_data)))

#ghcn_dayq_dat <- readOGR(dsn=dirname(met_stations_obj$daily_query_ghcn_data),

#                        sub(".shp","",basename(met_stations_obj$daily_query_ghcn_data)))

#y_lab_text <- "RMSE (C)"

#y_lab_text <- expression(paste("RMSE", ("^o","C")"), sep="")) 

y_lab_text <- "RMSE (°C)"

y_lab_text <- "MAE (°C)"

      ylab="MAE (°C)")

      ylab="RMSE (°C)")

      ylab="RMSE (°C)")

      ylab="RMSE (°C)")

      ylab="MAE (°C)")

#########Figure 6: Monthly RMSE averages for the three interpolation methods: GAM, GWR and Kriging.

ylab_text <- "MAE (°C)"

ylab_text<-"MAE (°C)"

index<-244 #index corresponding to Sept 1

#png(paste("Figure7__spatial_pattern_tmax_prediction_levelplot_",date_selected,out_prefix,".png", sep=""),

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

p<-levelplot(pred_temp_s,main="Interpolated Surfaces Method Comparison", ylab=NULL,xlab=NULL,

print(p)

#dev.off()

## FIGURE COMPARISON OF  MODELS COVARRIATES: Figure 7...

                 "mod5 = s(lat,long)+s(LST)","mod6 = s(lat,long)+s(DISTOC)","mod7 = s(lat,long)+s(FOREST)",

                 "mod8 = s(lat,long)+s(CANHGHT)","mod9 = s(lat,long)+s(LST,FOREST)","mod10 = s(lat,long)+s(LST,CANHGHT)")

p<- levelplot(pred_temp_s,main="Interpolated Surfaces Model Comparison baseline 1", ylab=NULL,xlab=NULL,

print(p)

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

#### Examining spatial correlograms for each model...use spedep but can be problematic when many datapoints

#The follwing method using Moran's I raster command works for now but does not provide CI/std dev...

r_stack <-pred_temp_s

r<- subset(r_stack,"mod1")

Moran(r) #with lag 1 and default rooks lag correlation

#generate filters for 10 lags: quick solution

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

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

list_param_moran <- list(list_filters=list_filters,r_stack=r_stack)

#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)

#prepare to automate the plotting of   all columns

mydata<-moran_df

dd <- do.call(make.groups, mydata[,-ncol(mydata)]) 

dd$lag <- mydata$lag 

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

png(paste("Figure_7b_spatial_correlogram_tmax_prediction_models_baseline1_gam_levelplot_",date_selected,out_prefix,".png", sep=""),

    height=480*layout_m[1],width=480*layout_m[2])

p<-xyplot(data ~ lag | which, dd,type="b",

          par.settings = list(axis.text = list(font = 2, cex = 1.3),layout=layout_m,

                              par.main.text=list(font=2,cex=2),strip.background=list(col="white")),par.strip.text=list(font=2,cex=1.5),

          strip=strip.custom(factor.levels=names_layers),

          xlab=list(label="Spatial lag neighbor", cex=2,font=2),

          ylab=list(label="Moran's I", cex=2, font=2))

print(p)

dev.off()

### median for first neighbour distances in data

ghcn_day_dat <- readOGR(dsn=dirname(met_stations_obj$daily_covar_ghcn_data),

                        sub(".shp","",basename(met_stations_obj$daily_covar_ghcn_data)))

loc_spdf <- readOGR(dsn=dirname(met_stations_obj$loc_stations),

                        sub(".shp","",basename(met_stations_obj$loc_stations)))

loc_2010_spdf<-subset(loc_spdf,loc_spdf$STAT_ID%in%ghcn_day_dat$station) #stations in 2010

loc_2010_spdf <- spTransform(loc_2010_spdf, CRS=CRS(proj4string(ghcn_day_dat)))

loc_knn1 <- knearneigh(coordinates(loc_2010_spdf), k=1) #lag1

class(loc_knn1) #knn object

loc_nb1 <- knn2nb(loc_knn1) #nb object

dsts_nb1 <-unlist(nbdists(loc_nb1,coordinates(loc_2010_spdf))) #vector with first nearest neighbour distance for each station

median(dsts_nb1) #median first nearest neighbour distance: 20023.03 m

mean(dsts_nb1) #mean first nearest neighbour distance : 22480.91 meters

sd(dsts_nb1) #

hist(dsts_nb1)

################ #FIGURE 8: difference image

p<-levelplot(pred_temp_s,main="Interpolated Surfaces Model Comparison", ylab=NULL,xlab=NULL,

print(p)

dev.off

### generate filter for Moran's I function in raster package

autocor_filter_fun <-function(no_lag=1,f_type="queen"){

  if(f_type=="queen"){

    no_rows <- 2*no_lag +1

    border_row <-rep(1,no_rows)

    other_row <- c(1,rep(0,no_rows-2),1)

    other_rows <- rep(other_row,no_rows-2)

    mat_data<- c(border_row,other_rows,border_row)

    autocor_filter<-matrix(mat_data,nrow=no_rows)

  }

  #if(f_type=="rook){} #add later

  return(autocor_filter)

}

#Now run Moran's I for raster image given a list of  filters for different lags and raster stack

moran_multiple_fun<-function(i,list_param){

  #Parameters:

  #list_filters: list of filters with different lags in the image

  #r_stack: stack of raster image, only the selected layer is used...

  list_filters <-list_param$list_filters

  r <- subset(list_param$r_stack,i)

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

  moran_v <-as.data.frame(unlist(moran_list))

  names(moran_v)<-names(r)

  return(moran_v)

}