#MODIFIED ON: 10/05/2014            

library(zoo)

library(xts)

assign_FID_spatial_polygons_df <-function(list_spdf,ID_str=NULL){

  list_spdf_tmp <- vector("list",length(list_spdf))

  if(is.null(ID_str)){

    nf <- 0 #number of features

    #for(i in 1:length(spdf)){

    #    shp1 <- list_spdf[[i]]

    #    f <- nrow(shp1)

    #    nf <- nf + f

    #}

    #This assumes that the list has one feature per item list

    nf <- length(list_spdf)

    ID_str <- as.character(1:nf)

  }

  for(i in 1:length(list_spdf)){

    #test=spRbind(shps_tiles[[1]],shps_tiles[[2]])

    shp1 <- list_spdf[[i]]

    shp1$FID <- ID_str

    shp1<- spChFIDs(shp1, as.character(shp1$FID)) #assign ID

    list_spdf_tmp[[i]]  <-shp1

  }

  return(list_spdf_tmp)

}

combine_spatial_polygons_df_fun <- function(list_spdf_tmp,ID_str=NULL){

  if(is.null(ID_str)){

    #call function

    list_spdf_tmp <- assign_FID_spatial_polygons_df

  }

  combined_spdf <- list_spdf_tmp[[1]]

  for(i in 2:length(list_spdf_tmp)){

    combined_spdf <- rbind(combined_spdf,list_spdf_tmp[[i]])

    #sapply(slot(shps_tiles[[2]], "polygons"), function(x) slot(x, "ID"))

    #rownames(as(alaska.tract, "data.frame"))

  }

  return(combined_spdf)

}

out_dir <-"/data/project/layers/commons/NEX_data/output_run7_global_analyses_10042014/"

out_prefix<-"run7_global_analyses_10042014"

  path_to_shp <- file.path(out_dir,"/shapefiles")

  coordinates(ac_mod) <- ac_mod[,c("lon","lat")] 

#idx <- seq(as.Date('2010-01-01'),by="day", length=365)

r_mod2_rmse <- stack(l_rast)

r_mod2_n <- stack(l_rast_n)

#r_mod2_n <- stack(list.files(pattern="r_nn_mod2.*.rst"))

#r_mod1_n <- stack(list.files(pattern="r_nn_mod1.*.rst"))

l_rast_n[[1]]

#"./r_nn_mod2_20101231_run7_global_analyses_10042014.rst"

# Create a raster stack with time series tag

r_mod2_rmse <- setZ(r_mod2_rmse, idx)

names(r_mod2_rmse) <- l_date_selected

r_mod2_n <- setZ(r_mod2_n, idx)

names(r_mod2_n) <- l_date_selected

writeRaster(r_mod2_n,by=F)

levelplot(r_mod2_rmse,layers=1:12,panel=panel.levelplot.raster, col.regions=matlab.like(25))+p_shp

p_hist <- histogram(r_mod2_rmse,FUN=as.yearmon)

p_hist2 <- p_hist

print(p_hist)

p_hist$x.limits <- rep(list(c(-2,6)),12)

print(p_hist)

p_bw <- bwplot(r_mod2_rmse,FUN=as.yearmon)

print(p_bw)

p_bw2 <- p_bw

p_bw2$y.limits <- c(-2,6)

print(p_bw2)

r_mod2_rmse_m <- zApply(r_mod2_rmse,by=as.yearmon,fun="mean")

p_lev_rmse_m <- levelplot(r_mod2_rmse_m,panel=panel.levelplot.raster,col.regions=matlab.like(25)) + p_shp

print(p_lev_rmse_m)

## analysis of the daily nmber of validation stations per tile

r_mod2_n_m <- zApply(r_mod2_n,by=as.yearmon,fun="mean")

p_lev_n_m <- levelplot(r_mod2_n_m,panel=panel.levelplot.raster,col.regions=matlab.like(25))+p_shp

print(p_lev_n_m)

p_bw_n2 <- bwplot(r_mod2_n,FUN=as.yearmon,do.out=F)

print(p_bw_n)

p_bw_n$y.limits <- c(0,16)

print(p_bw_n)

p_ho <- hovmoller(r_mod2_rmse)

p_lev_n <- levelplot(r_mod2_n,layers=1:12,panel=panel.levelplot.raster, col.regions=matlab.like(25))+p_shp

print(p_lev_n)

### mod1

#as(r_mod2_rmse)

# dd <- merge(df_tile_processed,pred_data_month_info,by="tile_id")

# coordinates(dd) <- c(dd$x,dd$y)

# 

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

### Figure 9: Plot the number of stations in a processing tile

## Get ID from tile number...

#ID_str <- unlist(lapply(1:nrow(df_tile_processed),function(i){unlist(strsplit(as.character(df_tile_processed$tile_id[i]),"_"))[2]}))

#combined_shp$tile_id <- df_tile_processed$tile_id

#r <- raster(lf_pred_list[i])

#plot(combined_shp)

#p0 <- spplot(combined_shp, "Stations",col.regions=matlab.like(100))

#p1 <- spplot(usa_map,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color?

#p2 <- spplot(can_map,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color?

#p3 <- spplot(mex_map,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color?

#p0 +p1+p2+p3

#df_dat <- subset(pred_data_month_info, pred_mod == "mod1" & date =="20100115")

#shp_dat <-merge(x=combined_shp,y=df_dat,by="tile_id",all.x=T) #if tile is missing then add rows with NA

#pol <- SpatialPolygons(combined_shp@polygons,proj4string=CRS(CRS_WGS84))

#spp <- SpatialPolygonsDataFrame(pol,data=shp_dat)

#p0 <- spplot(spp, "n_mod",col.regions=matlab.like(100))

#p1 <- spplot(reg_layer,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color?

#p0 + p1 + p2 + p3