Revision d9ebfb1a
Added by Benoit Parmentier about 9 years ago
| climate/research/oregon/interpolation/global_run_scalingup_assessment_part2_functions.R | ||
|---|---|---|
| 1 |
############################## INTERPOLATION OF TEMPERATURES ####################################### |
|
| 2 |
####################### Script for assessment of scaling up on NEX : part2 ############################## |
|
| 3 |
#This script uses the worklfow code applied to the globe. Results currently reside on NEX/PLEIADES NASA. |
|
| 4 |
#Accuracy methods are added in the the function scripts to evaluate results. |
|
| 5 |
#Analyses, figures, tables and data are also produced in the script. |
|
| 6 |
#AUTHOR: Benoit Parmentier |
|
| 7 |
#CREATED ON: 03/23/2014 |
|
| 8 |
#MODIFIED ON: 09/23/2015 |
|
| 9 |
#Version: 4 |
|
| 10 |
#PROJECT: Environmental Layers project |
|
| 11 |
#COMMENTS: analyses for run 10 global analyses,all regions 1500x4500km with additional tiles to increase overlap |
|
| 12 |
#TODO: |
|
| 13 |
#1) Split functions and master script |
|
| 14 |
#2) Make this is a script/function callable from the shell/bash |
|
| 15 |
#3) Check image format for tif |
|
| 16 |
|
|
| 17 |
################################################################################################# |
|
| 18 |
|
|
| 19 |
### Loading R library and packages |
|
| 20 |
#library used in the workflow production: |
|
| 21 |
library(gtools) # loading some useful tools |
|
| 22 |
library(mgcv) # GAM package by Simon Wood |
|
| 23 |
library(sp) # Spatial pacakge with class definition by Bivand et al. |
|
| 24 |
library(spdep) # Spatial pacakge with methods and spatial stat. by Bivand et al. |
|
| 25 |
library(rgdal) # GDAL wrapper for R, spatial utilities |
|
| 26 |
library(gstat) # Kriging and co-kriging by Pebesma et al. |
|
| 27 |
library(fields) # NCAR Spatial Interpolation methods such as kriging, splines |
|
| 28 |
library(raster) # Hijmans et al. package for raster processing |
|
| 29 |
library(gdata) # various tools with xls reading, cbindX |
|
| 30 |
library(rasterVis) # Raster plotting functions |
|
| 31 |
library(parallel) # Parallelization of processes with multiple cores |
|
| 32 |
library(maptools) # Tools and functions for sp and other spatial objects e.g. spCbind |
|
| 33 |
library(maps) # Tools and data for spatial/geographic objects |
|
| 34 |
library(reshape) # Change shape of object, summarize results |
|
| 35 |
library(plotrix) # Additional plotting functions |
|
| 36 |
library(plyr) # Various tools including rbind.fill |
|
| 37 |
library(spgwr) # GWR method |
|
| 38 |
library(automap) # Kriging automatic fitting of variogram using gstat |
|
| 39 |
library(rgeos) # Geometric, topologic library of functions |
|
| 40 |
#RPostgreSQL # Interface R and Postgres, not used in this script |
|
| 41 |
library(gridExtra) |
|
| 42 |
#Additional libraries not used in workflow |
|
| 43 |
library(pgirmess) # Krusall Wallis test with mulitple options, Kruskalmc {pgirmess}
|
|
| 44 |
library(colorRamps) |
|
| 45 |
library(zoo) |
|
| 46 |
library(xts) |
|
| 47 |
|
|
| 48 |
#### FUNCTION USED IN SCRIPT |
|
| 49 |
|
|
| 50 |
function_analyses_paper1 <-"contribution_of_covariates_paper_interpolation_functions_07182014.R" #first interp paper |
|
| 51 |
function_analyses_paper2 <-"multi_timescales_paper_interpolation_functions_08132014.R" |
|
| 52 |
|
|
| 53 |
load_obj <- function(f) |
|
| 54 |
{
|
|
| 55 |
env <- new.env() |
|
| 56 |
nm <- load(f, env)[1] |
|
| 57 |
env[[nm]] |
|
| 58 |
} |
|
| 59 |
|
|
| 60 |
create_dir_fun <- function(out_dir,out_suffix){
|
|
| 61 |
if(!is.null(out_suffix)){
|
|
| 62 |
out_name <- paste("output_",out_suffix,sep="")
|
|
| 63 |
out_dir <- file.path(out_dir,out_name) |
|
| 64 |
} |
|
| 65 |
#create if does not exists |
|
| 66 |
if(!file.exists(out_dir)){
|
|
| 67 |
dir.create(out_dir) |
|
| 68 |
} |
|
| 69 |
return(out_dir) |
|
| 70 |
} |
|
| 71 |
|
|
| 72 |
#Remove models that were not fitted from the list |
|
| 73 |
#All modesl that are "try-error" are removed |
|
| 74 |
remove_errors_list<-function(list_items){
|
|
| 75 |
|
|
| 76 |
#This function removes "error" items in a list |
|
| 77 |
list_tmp<-list_items |
|
| 78 |
if(is.null(names(list_tmp))){
|
|
| 79 |
names(list_tmp) <- paste("l",1:length(list_tmp),sep="_")
|
|
| 80 |
names(list_items) <- paste("l",1:length(list_tmp),sep="_")
|
|
| 81 |
} |
|
| 82 |
|
|
| 83 |
for(i in 1:length(list_items)){
|
|
| 84 |
if(inherits(list_items[[i]],"try-error")){
|
|
| 85 |
list_tmp[[i]]<-0 |
|
| 86 |
}else{
|
|
| 87 |
list_tmp[[i]]<-1 |
|
| 88 |
} |
|
| 89 |
} |
|
| 90 |
cnames<-names(list_tmp[list_tmp>0]) |
|
| 91 |
x <- list_items[match(cnames,names(list_items))] |
|
| 92 |
return(x) |
|
| 93 |
} |
|
| 94 |
|
|
| 95 |
#turn term from list into data.frame |
|
| 96 |
#name_col<-function(i,x){
|
|
| 97 |
#x_mat<-x[[i]] |
|
| 98 |
#x_df<-as.data.frame(x_mat) |
|
| 99 |
#x_df$mod_name<-rep(names(x)[i],nrow(x_df)) |
|
| 100 |
#x_df$term_name <-row.names(x_df) |
|
| 101 |
#return(x_df) |
|
| 102 |
#} |
|
| 103 |
#Function to rasterize a table with coordinates and variables...,maybe add option for ref image?? |
|
| 104 |
rasterize_df_fun <- function(data_tb,coord_names,proj_str,out_suffix,out_dir=".",file_format=".rst",NA_flag_val=-9999,tolerance_val= 0.000120005){
|
|
| 105 |
data_spdf <- data_tb |
|
| 106 |
coordinates(data_spdf) <- cbind(data_spdf[[coord_names[1]]],data_spdf[[coord_names[2]]]) |
|
| 107 |
proj4string(data_spdf) <- proj_str |
|
| 108 |
|
|
| 109 |
data_pix <- try(as(data_spdf,"SpatialPixelsDataFrame")) |
|
| 110 |
#tolerance_val <- 0.000120005 |
|
| 111 |
#tolerance_val <- 0.000856898 |
|
| 112 |
if(inherits(data_pix,"try-error")){
|
|
| 113 |
data_pix <- SpatialPixelsDataFrame(data_spdf, data=data_spdf@data, tolerance=tolerance_val) |
|
| 114 |
} |
|
| 115 |
|
|
| 116 |
#test <- as(data_spdf,"SpatialPixelsDataFrame") |
|
| 117 |
|
|
| 118 |
# set up an 'empty' raster, here via an extent object derived from your data |
|
| 119 |
#e <- extent(s100[,1:2]) |
|
| 120 |
#e <- e + 1000 # add this as all y's are the same |
|
| 121 |
|
|
| 122 |
#r <- raster(e, ncol=10, nrow=2) |
|
| 123 |
# or r <- raster(xmn=, xmx=, ... |
|
| 124 |
|
|
| 125 |
data_grid <- as(data_pix,"SpatialGridDataFrame") #making it a regural grid |
|
| 126 |
r_ref <- raster(data_grid) #this is the ref image |
|
| 127 |
rast_list <- vector("list",length=ncol(data_tb))
|
|
| 128 |
|
|
| 129 |
for(i in 1:(ncol(data_tb))){
|
|
| 130 |
field_name <- names(data_tb)[i] |
|
| 131 |
var <- as.numeric(data_spdf[[field_name]]) |
|
| 132 |
data_spdf$var <- var |
|
| 133 |
#r <-rasterize(data_spdf,r_ref,field_name) |
|
| 134 |
r <-rasterize(data_spdf,r_ref,"var",NAflag=NA_flag_val,fun=mean) #prolem with NA in NDVI!! |
|
| 135 |
|
|
| 136 |
data_name<-paste("r_",field_name,sep="") #can add more later...
|
|
| 137 |
#raster_name<-paste(data_name,out_names[j],".tif", sep="") |
|
| 138 |
raster_name<-paste(data_name,out_suffix,file_format, sep="") |
|
| 139 |
|
|
| 140 |
writeRaster(r, NAflag=NA_flag_val,filename=file.path(out_dir,raster_name),overwrite=TRUE) |
|
| 141 |
#Writing the data in a raster file format... |
|
| 142 |
rast_list[i] <-file.path(out_dir,raster_name) |
|
| 143 |
} |
|
| 144 |
return(unlist(rast_list)) |
|
| 145 |
} |
|
| 146 |
|
|
| 147 |
plot_raster_tb_diagnostic <- function(reg_layer,tb_dat,df_tile_processed,date_selected,mod_selected,var_selected,out_suffix){
|
|
| 148 |
|
|
| 149 |
test <- subset(tb_dat,pred_mod==mod_selected & date==date_selected,select=c("tile_id",var_selected))
|
|
| 150 |
|
|
| 151 |
test_data_tb <- merge(df_tile_processed,test,by="tile_id",all=T) #keep all |
|
| 152 |
test_r <- subset(test_data_tb,select=c("lat","lon","tile_id",var_selected))
|
|
| 153 |
out_suffix_str <- paste(var_selected,mod_selected,date_selected,out_suffix,sep="_") |
|
| 154 |
coord_names <- c("lon","lat")
|
|
| 155 |
l_rast <- rasterize_df_fun(test_r,coord_names,proj_str,out_suffix_str,out_dir=".",file_format=".tif",NA_flag_val=-9999,tolerance_val=0.000120005) |
|
| 156 |
#mod_kr_stack <- stack(mod_kr_l_rast) |
|
| 157 |
d_tb_rast <- stack(l_rast) |
|
| 158 |
names(d_tb_rast) <- names(test_r) |
|
| 159 |
#plot(d_tb_rast) |
|
| 160 |
r <- subset(d_tb_rast,"rmse") |
|
| 161 |
names(r) <- paste(mod_selected,var_selected,date_selected,sep="_") |
|
| 162 |
#plot info: with labels |
|
| 163 |
|
|
| 164 |
res_pix <- 1200 |
|
| 165 |
col_mfrow <- 1 |
|
| 166 |
row_mfrow <- 1 |
|
| 167 |
|
|
| 168 |
png(filename=paste("Figure9_",names(r),"_map_processed_region_",region_name,"_",out_suffix,".png",sep=""),
|
|
| 169 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 170 |
|
|
| 171 |
#plot(reg_layer) |
|
| 172 |
#p1 <- spplot(reg_layer,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color? |
|
| 173 |
title_str <- paste(names(r),"for ", region_name,sep="") |
|
| 174 |
|
|
| 175 |
p0 <- levelplot(r,col.regions=matlab.like(25),margin=F,main=title_str) |
|
| 176 |
p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black')) |
|
| 177 |
|
|
| 178 |
p <- p0 + p_shp |
|
| 179 |
print(p) |
|
| 180 |
|
|
| 181 |
dev.off() |
|
| 182 |
|
|
| 183 |
} |
|
| 184 |
|
|
| 185 |
create_raster_from_tb_diagnostic <- function(i,list_param){
|
|
| 186 |
#create a raster image using tile centroids and given fields from tb diagnostic data |
|
| 187 |
tb_dat <- list_param$tb_dat |
|
| 188 |
df_tile_processed <- list_param$df_tile_processed |
|
| 189 |
date_selected <- list_param$date_selected[i] |
|
| 190 |
mod_selected <- list_param$mod_selected |
|
| 191 |
var_selected <- list_param$var_selected |
|
| 192 |
out_suffix <- list_param$out_suffix |
|
| 193 |
|
|
| 194 |
test <- subset(tb_dat,pred_mod==mod_selected & date==date_selected,select=c("tile_id",var_selected))
|
|
| 195 |
|
|
| 196 |
test_data_tb <- merge(df_tile_processed,test,by="tile_id",all=T) #keep all |
|
| 197 |
test_r <- subset(test_data_tb,select=c("lat","lon","tile_id",var_selected))
|
|
| 198 |
out_suffix_str <- paste(var_selected,mod_selected,date_selected,out_suffix,sep="_") |
|
| 199 |
coord_names <- c("lon","lat")
|
|
| 200 |
l_rast <- rasterize_df_fun(test_r,coord_names,proj_str,out_suffix_str,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005) |
|
| 201 |
#mod_kr_stack <- stack(mod_kr_l_rast) |
|
| 202 |
#d_tb_rast <- stack(l_rast) |
|
| 203 |
#r <- subset(d_tb_rast,var_selected) |
|
| 204 |
#names(d_tb_rast) <- names(test_r) |
|
| 205 |
return(l_rast[4]) |
|
| 206 |
} |
|
| 207 |
|
|
| 208 |
assign_FID_spatial_polygons_df <-function(list_spdf,ID_str=NULL){
|
|
| 209 |
list_spdf_tmp <- vector("list",length(list_spdf))
|
|
| 210 |
if(is.null(ID_str)){
|
|
| 211 |
nf <- 0 #number of features |
|
| 212 |
#for(i in 1:length(spdf)){
|
|
| 213 |
# shp1 <- list_spdf[[i]] |
|
| 214 |
# f <- nrow(shp1) |
|
| 215 |
# nf <- nf + f |
|
| 216 |
#} |
|
| 217 |
#This assumes that the list has one feature per item list |
|
| 218 |
nf <- length(list_spdf) |
|
| 219 |
ID_str <- as.character(1:nf) |
|
| 220 |
} |
|
| 221 |
for(i in 1:length(list_spdf)){
|
|
| 222 |
#test=spRbind(shps_tiles[[1]],shps_tiles[[2]]) |
|
| 223 |
shp1 <- list_spdf[[i]] |
|
| 224 |
shp1$FID <- ID_str |
|
| 225 |
shp1<- spChFIDs(shp1, as.character(shp1$FID)) #assign ID |
|
| 226 |
list_spdf_tmp[[i]] <-shp1 |
|
| 227 |
} |
|
| 228 |
return(list_spdf_tmp) |
|
| 229 |
} |
|
| 230 |
|
|
| 231 |
combine_spatial_polygons_df_fun <- function(list_spdf_tmp,ID_str=NULL){
|
|
| 232 |
if(is.null(ID_str)){
|
|
| 233 |
#call function |
|
| 234 |
list_spdf_tmp <- assign_FID_spatial_polygons_df |
|
| 235 |
} |
|
| 236 |
combined_spdf <- list_spdf_tmp[[1]] |
|
| 237 |
for(i in 2:length(list_spdf_tmp)){
|
|
| 238 |
combined_spdf <- rbind(combined_spdf,list_spdf_tmp[[i]]) |
|
| 239 |
#sapply(slot(shps_tiles[[2]], "polygons"), function(x) slot(x, "ID")) |
|
| 240 |
#rownames(as(alaska.tract, "data.frame")) |
|
| 241 |
} |
|
| 242 |
return(combined_spdf) |
|
| 243 |
} |
|
| 244 |
|
|
| 245 |
plot_daily_mosaics <- function(i,list_param_plot_daily_mosaics){
|
|
| 246 |
#Purpose: |
|
| 247 |
#This functions mask mosaics files for a default range (-100,100 deg). |
|
| 248 |
#It produces a masked tif in a given dataType format (FLT4S) |
|
| 249 |
#It creates a figure of mosaiced region being interpolated. |
|
| 250 |
#Author: Benoit Parmentier |
|
| 251 |
#Parameters: |
|
| 252 |
#lf_m: list of files |
|
| 253 |
#reg_name:region name with tile size included |
|
| 254 |
#To do: |
|
| 255 |
#Add filenames |
|
| 256 |
#Add range |
|
| 257 |
#Add output dir |
|
| 258 |
#Add dataType_val option |
|
| 259 |
|
|
| 260 |
##### BEGIN ######## |
|
| 261 |
|
|
| 262 |
#Parse the list of parameters |
|
| 263 |
lf_m <- list_param_plot_daily_mosaics$lf_m |
|
| 264 |
reg_name <- list_param_plot_daily_mosaics$reg_name |
|
| 265 |
out_dir_str <- list_param_plot_daily_mosaics$out_dir_str |
|
| 266 |
out_suffix <- list_param_plot_daily_mosaics$out_suffix |
|
| 267 |
l_dates <- list_param_plot_daily_mosaics$l_dates |
|
| 268 |
|
|
| 269 |
|
|
| 270 |
#list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str) |
|
| 271 |
|
|
| 272 |
|
|
| 273 |
#rast_list <- vector("list",length=length(lf_m))
|
|
| 274 |
r_test<- raster(lf_m[i]) |
|
| 275 |
|
|
| 276 |
m <- c(-Inf, -100, NA, |
|
| 277 |
-100, 100, 1, |
|
| 278 |
100, Inf, NA) #can change the thresholds later |
|
| 279 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
|
| 280 |
rc <- reclassify(r_test, rclmat,filename=paste("rc_tmp_",i,".tif",sep=""),dataType="FLT4S",overwrite=T)
|
|
| 281 |
file_name <- unlist(strsplit(basename(lf_m[i]),"_")) |
|
| 282 |
|
|
| 283 |
#date_proc <- file_name[7] #specific tot he current naming of files |
|
| 284 |
date_proc <- l_dates[i] |
|
| 285 |
#paste(raster_name[1:7],collapse="_") |
|
| 286 |
#add filename option later |
|
| 287 |
extension_str <- extension(filename(r_test)) |
|
| 288 |
raster_name_tmp <- gsub(extension_str,"",basename(filename(r_test))) |
|
| 289 |
raster_name <- file.path(out_dir_str,paste(raster_name_tmp,"_masked.tif",sep="")) |
|
| 290 |
r_pred <- mask(r_test,rc,filename=raster_name,overwrite=TRUE) |
|
| 291 |
|
|
| 292 |
res_pix <- 1200 |
|
| 293 |
#res_pix <- 480 |
|
| 294 |
|
|
| 295 |
col_mfrow <- 1 |
|
| 296 |
row_mfrow <- 1 |
|
| 297 |
|
|
| 298 |
png(filename=paste("Figure9_clim_mosaics_day_test","_",date_proc,"_",reg_name,"_",out_suffix,".png",sep=""),
|
|
| 299 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 300 |
|
|
| 301 |
plot(r_pred,main=paste("Predicted on ",date_proc , " ", reg_name,sep=""),cex.main=1.5)
|
|
| 302 |
dev.off() |
|
| 303 |
|
|
| 304 |
return(raster_name) |
|
| 305 |
|
|
| 306 |
} |
|
| 307 |
|
|
| 308 |
plot_screen_raster_val<-function(i,list_param){
|
|
| 309 |
##USAGE### |
|
| 310 |
#Screen raster list and produced plot as png. |
|
| 311 |
fname <-list_param$lf_raster_fname[i] |
|
| 312 |
screenRast <- list_param$screenRast |
|
| 313 |
l_dates<- list_param$l_dates |
|
| 314 |
out_dir_str <- list_param$out_dir_str |
|
| 315 |
prefix_str <-list_param$prefix_str |
|
| 316 |
out_suffix_str <- list_param$out_suffix_str |
|
| 317 |
|
|
| 318 |
### START SCRIPT #### |
|
| 319 |
date_proc <- l_dates[i] |
|
| 320 |
|
|
| 321 |
if(screenRast==TRUE){
|
|
| 322 |
r_test <- raster(fname) |
|
| 323 |
|
|
| 324 |
m <- c(-Inf, -100, NA, |
|
| 325 |
-100, 100, 1, |
|
| 326 |
100, Inf, NA) #can change the thresholds later |
|
| 327 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
|
| 328 |
rc <- reclassify(r_test, rclmat,filename=paste("rc_tmp_",i,".tif",sep=""),dataType="FLT4S",overwrite=T)
|
|
| 329 |
#file_name <- unlist(strsplit(basename(lf_m[i]),"_")) |
|
| 330 |
extension_str <- extension(filename(r_test)) |
|
| 331 |
raster_name_tmp <- gsub(extension_str,"",basename(filename(r_test))) |
|
| 332 |
raster_name <- file.path(out_dir_str,paste(raster_name_tmp,"_masked.tif",sep="")) |
|
| 333 |
|
|
| 334 |
r_pred <- mask(r_test,rc,filename=raster_name,overwrite=TRUE) |
|
| 335 |
}else{
|
|
| 336 |
r_pred <- raster(fname) |
|
| 337 |
} |
|
| 338 |
|
|
| 339 |
#date_proc <- file_name[7] #specific tot he current naming of files |
|
| 340 |
#date_proc<- "2010010101" |
|
| 341 |
|
|
| 342 |
#paste(raster_name[1:7],collapse="_") |
|
| 343 |
#add filename option later |
|
| 344 |
|
|
| 345 |
res_pix <- 960 |
|
| 346 |
#res_pix <- 480 |
|
| 347 |
|
|
| 348 |
col_mfrow <- 1 |
|
| 349 |
row_mfrow <- 1 |
|
| 350 |
|
|
| 351 |
# png(filename=paste("Figure10_clim_world_mosaics_day_","_",date_proc,"_",tile_size,"_",out_suffix,".png",sep=""),
|
|
| 352 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 353 |
png(filename=paste(prefix_str,"_",date_proc,"_",tile_size,"_",out_suffix_str,".png",sep=""), |
|
| 354 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 355 |
|
|
| 356 |
plot(r_pred,main=paste("Predicted on ",date_proc , " ", tile_size,sep=""),cex.main=1.5)
|
|
| 357 |
dev.off() |
|
| 358 |
|
|
| 359 |
} |
|
| 360 |
|
|
| 361 |
############################################ |
|
| 362 |
#### Parameters and constants |
|
| 363 |
|
|
| 364 |
#on ATLAS |
|
| 365 |
#in_dir1 <- "/data/project/layers/commons/NEX_data/test_run1_03232014/output" #On Atlas |
|
| 366 |
#parent output dir : contains subset of the data produced on NEX |
|
| 367 |
#in_dir1 <- "/data/project/layers/commons/NEX_data/output_run6_global_analyses_09162014/output20Deg2" |
|
| 368 |
# parent output dir for the curent script analyes |
|
| 369 |
#out_dir <-"/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/" #On NCEAS Atlas |
|
| 370 |
# input dir containing shapefiles defining tiles |
|
| 371 |
#in_dir_shp <- "/data/project/layers/commons/NEX_data/output_run5_global_analyses_08252014/output/subset/shapefiles" |
|
| 372 |
|
|
| 373 |
#On NEX |
|
| 374 |
#contains all data from the run by Alberto |
|
| 375 |
#in_dir1 <- "/nobackupp4/aguzman4/climateLayers/output4" #On NEX |
|
| 376 |
#parent output dir for the current script analyes |
|
| 377 |
#out_dir <- "/nobackup/bparmen1/" #on NEX |
|
| 378 |
#in_dir_shp <- "/nobackupp4/aguzman4/climateLayers/output4/subset/shapefiles/" |
|
| 379 |
|
|
| 380 |
y_var_name <- "dailyTmax" #PARAM1 |
|
| 381 |
interpolation_method <- c("gam_CAI") #PARAM2
|
|
| 382 |
#out_suffix<-"run10_global_analyses_01282015" |
|
| 383 |
#out_suffix <- "output_run10_1000x3000_global_analyses_02102015" |
|
| 384 |
out_suffix <- "run10_1500x4500_global_analyses_pred_1982_09152015" #PARAM3 |
|
| 385 |
out_dir <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1982_09152015" #PARAM4 |
|
| 386 |
create_out_dir_param <- FALSE #PARAM 5 |
|
| 387 |
|
|
| 388 |
mosaic_plot <- FALSE #PARAM6 |
|
| 389 |
|
|
| 390 |
#if daily mosaics NULL then mosaicas all days of the year |
|
| 391 |
|
|
| 392 |
day_to_mosaic <- c("19820101","19820102","19820103","19820104","19820105",
|
|
| 393 |
"19820106","19820107","19820108","19820109","19820110", |
|
| 394 |
"1982011") |
|
| 395 |
|
|
| 396 |
CRS_WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84 #CONSTANT1
|
|
| 397 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
|
| 398 |
|
|
| 399 |
proj_str<- CRS_WGS84 #PARAM 8 #check this parameter |
|
| 400 |
file_format <- ".rst" #PARAM 9 |
|
| 401 |
NA_value <- -9999 #PARAM10 |
|
| 402 |
NA_flag_val <- NA_value |
|
| 403 |
|
|
| 404 |
tile_size <- "1500x4500" #PARAM 11 |
|
| 405 |
multiple_region <- TRUE #PARAM 12 |
|
| 406 |
|
|
| 407 |
region_name <- "world" #PARAM 13 |
|
| 408 |
plot_region <- TRUE |
|
| 409 |
num_cores <- 6 #PARAM 14 |
|
| 410 |
reg_modified <- TRUE |
|
| 411 |
region <- c("reg4") #reference region to merge if necessary #PARAM 16
|
|
| 412 |
|
|
| 413 |
########################## START SCRIPT ############################## |
|
| 414 |
|
|
| 415 |
|
|
| 416 |
####### PART 1: Read in data ######## |
|
| 417 |
|
|
| 418 |
if(create_out_dir_param==TRUE){
|
|
| 419 |
out_dir <- create_dir_fun(out_dir,out_suffix) |
|
| 420 |
setwd(out_dir) |
|
| 421 |
}else{
|
|
| 422 |
setwd(out_dir) #use previoulsy defined directory |
|
| 423 |
} |
|
| 424 |
|
|
| 425 |
setwd(out_dir) |
|
| 426 |
|
|
| 427 |
###Table 1: Average accuracy metrics |
|
| 428 |
###Table 2: daily accuracy metrics for all tiles |
|
| 429 |
|
|
| 430 |
summary_metrics_v <- read.table(file=file.path(out_dir,paste("summary_metrics_v2_NA_",out_suffix,".txt",sep="")),sep=",")
|
|
| 431 |
#fname <- file.path(out_dir,paste("summary_metrics_v2_NA_",out_suffix,".txt",sep=""))
|
|
| 432 |
tb <- read.table(file=file.path(out_dir,paste("tb_diagnostic_v_NA","_",out_suffix,".txt",sep="")),sep=",")
|
|
| 433 |
#tb_diagnostic_s_NA_run10_global_analyses_11302014.txt |
|
| 434 |
tb_s <- read.table(file=file.path(out_dir,paste("tb_diagnostic_s_NA","_",out_suffix,".txt",sep="")),sep=",")
|
|
| 435 |
|
|
| 436 |
tb_month_s <- read.table(file=file.path(out_dir,paste("tb_month_diagnostic_s_NA","_",out_suffix,".txt",sep="")),sep=",")
|
|
| 437 |
pred_data_month_info <- read.table(file=file.path(out_dir,paste("pred_data_month_info_",out_suffix,".txt",sep="")),sep=",")
|
|
| 438 |
pred_data_day_info <- read.table(file=file.path(out_dir,paste("pred_data_day_info_",out_suffix,".txt",sep="")),sep=",")
|
|
| 439 |
df_tile_processed <- read.table(file=file.path(out_dir,paste("df_tile_processed_",out_suffix,".txt",sep="")),sep=",")
|
|
| 440 |
|
|
| 441 |
#add column for tile size later on!!! |
|
| 442 |
|
|
| 443 |
tb$pred_mod <- as.character(tb$pred_mod) |
|
| 444 |
summary_metrics_v$pred_mod <- as.character(summary_metrics_v$pred_mod) |
|
| 445 |
summary_metrics_v$tile_id <- as.character(summary_metrics_v$tile_id) |
|
| 446 |
df_tile_processed$tile_id <- as.character(df_tile_processed$tile_id) |
|
| 447 |
|
|
| 448 |
tb_month_s$pred_mod <- as.character(tb_month_s$pred_mod) |
|
| 449 |
tb_month_s$tile_id<- as.character(tb_month_s$tile_id) |
|
| 450 |
tb_s$pred_mod <- as.character(tb_s$pred_mod) |
|
| 451 |
tb_s$tile_id <- as.character(tb_s$tile_id) |
|
| 452 |
|
|
| 453 |
|
|
| 454 |
#multiple regions? |
|
| 455 |
if(multiple_region==TRUE){
|
|
| 456 |
df_tile_processed$reg <- basename(dirname(as.character(df_tile_processed$path_NEX))) |
|
| 457 |
|
|
| 458 |
tb <- merge(tb,df_tile_processed,by="tile_id") |
|
| 459 |
tb_s <- merge(tb_s,df_tile_processed,by="tile_id") |
|
| 460 |
tb_month_s<- merge(tb_month_s,df_tile_processed,by="tile_id") |
|
| 461 |
summary_metrics_v <- merge(summary_metrics_v,df_tile_processed,by="tile_id") |
|
| 462 |
|
|
| 463 |
} |
|
| 464 |
|
|
| 465 |
tb_all <- tb |
|
| 466 |
|
|
| 467 |
summary_metrics_v_all <- summary_metrics_v |
|
| 468 |
#deal with additional tiles... |
|
| 469 |
# if(reg_modified==T){
|
|
| 470 |
# |
|
| 471 |
# summary_metrics_v_tmp <- summary_metrics_v |
|
| 472 |
# #summary_metrics_v_tmp$reg[summary_metrics_v_tmp$reg=="reg_1b"] <- "reg1" |
|
| 473 |
# #summary_metrics_v_tmp$reg[summary_metrics_v_tmp$reg=="reg_1c"] <- "reg1" |
|
| 474 |
# #summary_metrics_v_tmp$reg[summary_metrics_v_tmp$reg=="reg_3b"] <- "reg3" |
|
| 475 |
# summary_metrics_v_tmp$reg[summary_metrics_v_tmp$reg=="reg5b"] <- "reg5" |
|
| 476 |
# |
|
| 477 |
# summary_metrics_v_tmp$reg_all <- summary_metrics_v$reg |
|
| 478 |
# ### |
|
| 479 |
# summary_metrics_v<- summary_metrics_v_tmp |
|
| 480 |
# |
|
| 481 |
# ### |
|
| 482 |
# tb_tmp <- tb |
|
| 483 |
# #tb_tmp$reg[tb_tmp$reg=="reg_1b"] <- "reg1" |
|
| 484 |
# #tb_tmp$reg[tb_tmp$reg=="reg_1c"] <- "reg1" |
|
| 485 |
# #tb_tmp$reg[tb_tmp$reg=="reg_3b"] <- "reg3" |
|
| 486 |
# tb_tmp$reg[tb_tmp$reg=="reg5b"] <- "reg5" |
|
| 487 |
# |
|
| 488 |
# ### |
|
| 489 |
# tb <- tb_tmp |
|
| 490 |
# } |
|
| 491 |
|
|
| 492 |
table(summary_metrics_v_all$reg) |
|
| 493 |
table(summary_metrics_v$reg) |
|
| 494 |
table(tb_all$reg) |
|
| 495 |
table(tb$reg) |
|
| 496 |
|
|
| 497 |
############ PART 2: PRODUCE FIGURES ################ |
|
| 498 |
|
|
| 499 |
########################### |
|
| 500 |
### Figure 1: plot location of the study area with tiles processed |
|
| 501 |
|
|
| 502 |
#df_tiled_processed <- na.omit(df_tile_processed) #remove other list of folders irrelevant |
|
| 503 |
#list_shp_reg_files <- df_tiled_processed$shp_files |
|
| 504 |
list_shp_reg_files<- as.character(df_tile_processed$shp_files) |
|
| 505 |
#list_shp_reg_files <- file.path("/data/project/layers/commons/NEX_data/",out_dir,
|
|
| 506 |
# as.character(df_tile_processed$tile_coord),"shapefiles",basename(list_shp_reg_files)) |
|
| 507 |
list_shp_reg_files <- file.path("/data/project/layers/commons/NEX_data/",out_dir,
|
|
| 508 |
"shapefiles",basename(list_shp_reg_files)) |
|
| 509 |
|
|
| 510 |
#table(summary_metrics_v$reg) |
|
| 511 |
#table(summary_metrics_v$reg) |
|
| 512 |
|
|
| 513 |
### Potential function starts here: |
|
| 514 |
#function(in_dir,out_dir,list_shp_reg_files,title_str,region_name,num_cores,out_suffix,out_suffix) |
|
| 515 |
|
|
| 516 |
### First get background map to display where study area is located |
|
| 517 |
#can make this more general later on..should have this already in a local directory on Atlas or NEX!!!! |
|
| 518 |
|
|
| 519 |
if (region_name=="USA"){
|
|
| 520 |
usa_map <- getData('GADM', country='USA', level=1) #Get US map
|
|
| 521 |
#usa_map <- getData('GADM', country=region_name,level=1) #Get US map, this is not working right now
|
|
| 522 |
usa_map <- usa_map[usa_map$NAME_1!="Alaska",] #remove Alaska |
|
| 523 |
reg_layer <- usa_map[usa_map$NAME_1!="Hawaii",] #remove Hawai |
|
| 524 |
} |
|
| 525 |
|
|
| 526 |
if (region_name=="world"){
|
|
| 527 |
#http://www.diva-gis.org/Data |
|
| 528 |
countries_shp <-"/data/project/layers/commons/NEX_data/countries.shp" |
|
| 529 |
path_to_shp <- dirname(countries_shp) |
|
| 530 |
layer_name <- sub(".shp","",basename(countries_shp))
|
|
| 531 |
reg_layer <- readOGR(path_to_shp, layer_name) |
|
| 532 |
#proj4string(reg_layer) <- CRS_locs_WGS84 |
|
| 533 |
#reg_shp<-readOGR(dirname(list_shp_reg_files[[i]]),sub(".shp","",basename(list_shp_reg_files[[i]])))
|
|
| 534 |
} |
|
| 535 |
|
|
| 536 |
centroids_pts <- vector("list",length(list_shp_reg_files))
|
|
| 537 |
shps_tiles <- vector("list",length(list_shp_reg_files))
|
|
| 538 |
#collect info: read in all shapfiles |
|
| 539 |
#This is slow...make a function and use mclapply?? |
|
| 540 |
#/data/project/layers/commons/NEX_data/output_run6_global_analyses_09162014/shapefiles |
|
| 541 |
|
|
| 542 |
for(i in 1:length(list_shp_reg_files)){
|
|
| 543 |
#path_to_shp <- dirname(list_shp_reg_files[[i]]) |
|
| 544 |
path_to_shp <- file.path(out_dir,"/shapefiles") |
|
| 545 |
layer_name <- sub(".shp","",basename(list_shp_reg_files[[i]]))
|
|
| 546 |
shp1 <- try(readOGR(path_to_shp, layer_name)) #use try to resolve error below |
|
| 547 |
#shp_61.0_-160.0 |
|
| 548 |
#Geographical CRS given to non-conformant data: -186.331747678 |
|
| 549 |
|
|
| 550 |
#shp1<-readOGR(dirname(list_shp_reg_files[[i]]),sub(".shp","",basename(list_shp_reg_files[[i]])))
|
|
| 551 |
if (!inherits(shp1,"try-error")) {
|
|
| 552 |
pt <- gCentroid(shp1) |
|
| 553 |
centroids_pts[[i]] <- pt |
|
| 554 |
}else{
|
|
| 555 |
pt <- shp1 |
|
| 556 |
centroids_pts[[i]] <- pt |
|
| 557 |
} |
|
| 558 |
shps_tiles[[i]] <- shp1 |
|
| 559 |
#centroids_pts[[i]] <- centroids |
|
| 560 |
} |
|
| 561 |
|
|
| 562 |
#fun <- function(i,list_shp_files) |
|
| 563 |
#coord_names <- c("lon","lat")
|
|
| 564 |
#l_ras#t <- rasterize_df_fun(test,coord_names,proj_str,out_suffix=out_suffix,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005) |
|
| 565 |
|
|
| 566 |
#remove try-error polygons...we loose three tiles because they extend beyond -180 deg |
|
| 567 |
tmp <- shps_tiles |
|
| 568 |
shps_tiles <- remove_errors_list(shps_tiles) #[[!inherits(shps_tiles,"try-error")]] |
|
| 569 |
#shps_tiles <- tmp |
|
| 570 |
length(tmp)-length(shps_tiles) #number of tiles with error message |
|
| 571 |
|
|
| 572 |
tmp_pts <- centroids_pts |
|
| 573 |
centroids_pts <- remove_errors_list(centroids_pts) #[[!inherits(shps_tiles,"try-error")]] |
|
| 574 |
#centroids_pts <- tmp_pts |
|
| 575 |
|
|
| 576 |
#plot info: with labels |
|
| 577 |
res_pix <-1200 |
|
| 578 |
col_mfrow <- 1 |
|
| 579 |
row_mfrow <- 1 |
|
| 580 |
|
|
| 581 |
png(filename=paste("Figure1_tile_processed_region_",region_name,"_",out_suffix,".png",sep=""),
|
|
| 582 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 583 |
|
|
| 584 |
plot(reg_layer) |
|
| 585 |
#Add polygon tiles... |
|
| 586 |
for(i in 1:length(shps_tiles)){
|
|
| 587 |
shp1 <- shps_tiles[[i]] |
|
| 588 |
pt <- centroids_pts[[i]] |
|
| 589 |
if(!inherits(shp1,"try-error")){
|
|
| 590 |
plot(shp1,add=T,border="blue") |
|
| 591 |
#plot(pt,add=T,cex=2,pch=5) |
|
| 592 |
label_id <- df_tile_processed$tile_id[i] |
|
| 593 |
text(coordinates(pt)[1],coordinates(pt)[2],labels=i,cex=1.3,font=2,col=c("red"))
|
|
| 594 |
} |
|
| 595 |
} |
|
| 596 |
#title(paste("Tiles ", tile_size,region_name,sep=""))
|
|
| 597 |
|
|
| 598 |
dev.off() |
|
| 599 |
|
|
| 600 |
#unique(summaty_metrics$tile_id) |
|
| 601 |
#text(lat-shp,) |
|
| 602 |
#union(list_shp_reg_files[[1]],list_shp_reg_files[[2]]) |
|
| 603 |
|
|
| 604 |
############### |
|
| 605 |
### Figure 2: boxplot of average accuracy by model and by tiles |
|
| 606 |
|
|
| 607 |
|
|
| 608 |
tb$tile_id <- factor(tb$tile_id, levels=unique(tb$tile_id)) |
|
| 609 |
|
|
| 610 |
model_name <- as.character(unique(tb$pred_mod)) |
|
| 611 |
|
|
| 612 |
|
|
| 613 |
## Figure 2a |
|
| 614 |
|
|
| 615 |
for(i in 1:length(model_name)){
|
|
| 616 |
|
|
| 617 |
res_pix <- 480 |
|
| 618 |
col_mfrow <- 1 |
|
| 619 |
row_mfrow <- 1 |
|
| 620 |
|
|
| 621 |
png(filename=paste("Figure2a_boxplot_with_oultiers_by_tiles_",model_name[i],"_",out_suffix,".png",sep=""),
|
|
| 622 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 623 |
|
|
| 624 |
boxplot(rmse~tile_id,data=subset(tb,tb$pred_mod==model_name[i])) |
|
| 625 |
title(paste("RMSE per ",model_name[i]))
|
|
| 626 |
|
|
| 627 |
dev.off() |
|
| 628 |
} |
|
| 629 |
|
|
| 630 |
## Figure 2b |
|
| 631 |
#wtih ylim and removing trailing... |
|
| 632 |
for(i in 1:length(model_name)){
|
|
| 633 |
|
|
| 634 |
res_pix <- 480 |
|
| 635 |
col_mfrow <- 1 |
|
| 636 |
row_mfrow <- 1 |
|
| 637 |
png(filename=paste("Figure2b_boxplot_scaling_by_tiles","_",model_name[i],"_",out_suffix,".png",sep=""),
|
|
| 638 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 639 |
|
|
| 640 |
model_name <- unique(tb$pred_mod) |
|
| 641 |
boxplot(rmse~tile_id,data=subset(tb,tb$pred_mod==model_name[i]) |
|
| 642 |
,ylim=c(0,4),outline=FALSE) |
|
| 643 |
title(paste("RMSE per ",model_name[i]))
|
|
| 644 |
dev.off() |
|
| 645 |
} |
|
| 646 |
#bwplot(rmse~tile_id, data=subset(tb,tb$pred_mod=="mod1")) |
|
| 647 |
|
|
| 648 |
############### |
|
| 649 |
### Figure 3: boxplot of average RMSE by model acrosss all tiles |
|
| 650 |
|
|
| 651 |
## Figure 3a |
|
| 652 |
res_pix <- 480 |
|
| 653 |
col_mfrow <- 1 |
|
| 654 |
row_mfrow <- 1 |
|
| 655 |
|
|
| 656 |
png(filename=paste("Figure3a_boxplot_overall_region_with_oultiers_",model_name[i],"_",out_suffix,".png",sep=""),
|
|
| 657 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 658 |
|
|
| 659 |
boxplot(rmse~pred_mod,data=tb)#,names=tb$pred_mod) |
|
| 660 |
title("RMSE per model over all tiles")
|
|
| 661 |
|
|
| 662 |
dev.off() |
|
| 663 |
|
|
| 664 |
## Figure 3b |
|
| 665 |
png(filename=paste("Figure3b_boxplot_overall_region_scaling_",model_name[i],"_",out_suffix,".png",sep=""),
|
|
| 666 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 667 |
|
|
| 668 |
boxplot(rmse~pred_mod,data=tb,ylim=c(0,5),outline=FALSE)#,names=tb$pred_mod) |
|
| 669 |
title("RMSE per model over all tiles")
|
|
| 670 |
|
|
| 671 |
dev.off() |
|
| 672 |
|
|
| 673 |
|
|
| 674 |
################ |
|
| 675 |
### Figure 4: plot predicted tiff for specific date per model |
|
| 676 |
|
|
| 677 |
#y_var_name <-"dailyTmax" |
|
| 678 |
#index <-244 #index corresponding to Sept 1 |
|
| 679 |
|
|
| 680 |
# if (mosaic_plot==TRUE){
|
|
| 681 |
# index <- 1 #index corresponding to Jan 1 |
|
| 682 |
# date_selected <- "20100901" |
|
| 683 |
# name_method_var <- paste(interpolation_method,"_",y_var_name,"_",sep="") |
|
| 684 |
# |
|
| 685 |
# pattern_str <- paste("mosaiced","_",name_method_var,"predicted",".*.",date_selected,".*.tif",sep="")
|
|
| 686 |
# lf_pred_list <- list.files(pattern=pattern_str) |
|
| 687 |
# |
|
| 688 |
# for(i in 1:length(lf_pred_list)){
|
|
| 689 |
# |
|
| 690 |
# |
|
| 691 |
# r_pred <- raster(lf_pred_list[i]) |
|
| 692 |
# |
|
| 693 |
# res_pix <- 480 |
|
| 694 |
# col_mfrow <- 1 |
|
| 695 |
# row_mfrow <- 1 |
|
| 696 |
# |
|
| 697 |
# png(filename=paste("Figure4_models_predicted_surfaces_",model_name[i],"_",name_method_var,"_",data_selected,"_",out_suffix,".png",sep=""),
|
|
| 698 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 699 |
# |
|
| 700 |
# plot(r_pred) |
|
| 701 |
# title(paste("Mosaiced",model_name[i],name_method_var,date_selected,sep=" "))
|
|
| 702 |
# dev.off() |
|
| 703 |
# } |
|
| 704 |
# #Plot Delta and clim... |
|
| 705 |
# |
|
| 706 |
# ## plotting of delta and clim for later scripts... |
|
| 707 |
# |
|
| 708 |
# } |
|
| 709 |
|
|
| 710 |
|
|
| 711 |
###################### |
|
| 712 |
### Figure 5: plot accuracy ranked |
|
| 713 |
|
|
| 714 |
#Turn summary table to a point shp |
|
| 715 |
|
|
| 716 |
list_df_ac_mod <- vector("list",length=length(model_name))
|
|
| 717 |
for (i in 1:length(model_name)){
|
|
| 718 |
|
|
| 719 |
ac_mod <- summary_metrics_v[summary_metrics_v$pred_mod==model_name[i],] |
|
| 720 |
### Ranking by tile... |
|
| 721 |
df_ac_mod <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("pred_mod","rmse","mae","tile_id")]
|
|
| 722 |
list_df_ac_mod[[i]] <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("rmse","mae","tile_id")]
|
|
| 723 |
|
|
| 724 |
res_pix <- 480 |
|
| 725 |
col_mfrow <- 1 |
|
| 726 |
row_mfrow <- 1 |
|
| 727 |
|
|
| 728 |
png(filename=paste("Figure5_ac_metrics_ranked_",model_name[i],"_",out_suffix,".png",sep=""),
|
|
| 729 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 730 |
x<- as.character(df_ac_mod$tile_id) |
|
| 731 |
barplot(df_ac_mod$rmse, names.arg=x) |
|
| 732 |
#plot(ac_mod1,cex=sqrt(ac_mod1$rmse),pch=1,add=T) |
|
| 733 |
#plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T) |
|
| 734 |
title(paste("RMSE ranked by tile for ",model_name[i],sep=""))
|
|
| 735 |
|
|
| 736 |
dev.off() |
|
| 737 |
|
|
| 738 |
} |
|
| 739 |
|
|
| 740 |
###################### |
|
| 741 |
### Figure 6: plot map of average RMSE per tile at centroids |
|
| 742 |
|
|
| 743 |
#Turn summary table to a point shp |
|
| 744 |
|
|
| 745 |
# coordinates(summary_metrics_v) <- cbind(summary_metrics_v$lon,summary_metrics_v$lat) |
|
| 746 |
# proj4string(summary_metrics_v) <- CRS_WGS84 |
|
| 747 |
# #lf_list <- lf_pred_list |
|
| 748 |
# list_df_ac_mod <- vector("list",length=length(lf_pred_list))
|
|
| 749 |
# for (i in 1:length(lf_list)){
|
|
| 750 |
# |
|
| 751 |
# ac_mod <- summary_metrics_v[summary_metrics_v$pred_mod==model_name[i],] |
|
| 752 |
# r_pred <- raster(lf_list[i]) |
|
| 753 |
# |
|
| 754 |
# res_pix <- 480 |
|
| 755 |
# col_mfrow <- 1 |
|
| 756 |
# row_mfrow <- 1 |
|
| 757 |
# |
|
| 758 |
# png(filename=paste("Figure6_ac_metrics_map_centroids_tile_",model_name[i],"_",out_suffix,".png",sep=""),
|
|
| 759 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 760 |
# |
|
| 761 |
# plot(r_pred) |
|
| 762 |
# |
|
| 763 |
# #plot(ac_mod1,cex=sqrt(ac_mod1$rmse),pch=1,add=T) |
|
| 764 |
# plot(ac_mod,cex=(ac_mod$rmse^2)/10,pch=1,add=T) |
|
| 765 |
# #plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T) |
|
| 766 |
# title(paste("Averrage RMSE per tile and by ",model_name[i]))
|
|
| 767 |
# |
|
| 768 |
# dev.off() |
|
| 769 |
# |
|
| 770 |
# ### Ranking by tile... |
|
| 771 |
# #df_ac_mod <- |
|
| 772 |
# list_df_ac_mod[[i]] <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("rmse","mae","tile_id")]
|
|
| 773 |
# } |
|
| 774 |
|
|
| 775 |
#quick kriging... |
|
| 776 |
#autokrige(rmse~1,r2,) |
|
| 777 |
|
|
| 778 |
|
|
| 779 |
### Without |
|
| 780 |
|
|
| 781 |
#list_df_ac_mod <- vector("list",length=length(lf_pred_list))
|
|
| 782 |
list_df_ac_mod <- vector("list",length=2)
|
|
| 783 |
|
|
| 784 |
for (i in 1:length(model_name)){
|
|
| 785 |
|
|
| 786 |
ac_mod <- summary_metrics_v[summary_metrics_v$pred_mod==model_name[i],] |
|
| 787 |
#r_pred <- raster(lf_list[i]) |
|
| 788 |
|
|
| 789 |
res_pix <- 1200 |
|
| 790 |
#res_pix <- 480 |
|
| 791 |
|
|
| 792 |
col_mfrow <- 1 |
|
| 793 |
row_mfrow <- 1 |
|
| 794 |
|
|
| 795 |
png(filename=paste("Figure6_ac_metrics_map_centroids_tile_",model_name[i],"_",out_suffix,".png",sep=""),
|
|
| 796 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 797 |
|
|
| 798 |
#plot(r_pred) |
|
| 799 |
#plot(reg_layer) |
|
| 800 |
#plot(ac_mod1,cex=sqrt(ac_mod1$rmse),pch=1,add=T) |
|
| 801 |
#plot(ac_mod,cex=(ac_mod$rmse^2)/10,pch=1,col="red",add=T) |
|
| 802 |
|
|
| 803 |
#coordinates(ac_mod) <- ac_mod[,c("lon","lat")]
|
|
| 804 |
coordinates(ac_mod) <- ac_mod[,c("lon.x","lat.x")] #solve this later
|
|
| 805 |
p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black')) |
|
| 806 |
#title("(a) Mean for 1 January")
|
|
| 807 |
p <- bubble(ac_mod,"rmse",main=paste("Averrage RMSE per tile and by ",model_name[i]))
|
|
| 808 |
p1 <- p+p_shp |
|
| 809 |
print(p1) |
|
| 810 |
#plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T) |
|
| 811 |
#title(paste("Averrage RMSE per tile and by ",model_name[i]))
|
|
| 812 |
|
|
| 813 |
dev.off() |
|
| 814 |
|
|
| 815 |
### Ranking by tile... |
|
| 816 |
#df_ac_mod <- |
|
| 817 |
list_df_ac_mod[[i]] <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("rmse","mae","tile_id")]
|
|
| 818 |
} |
|
| 819 |
|
|
| 820 |
## Number of tiles with information: |
|
| 821 |
sum(df_tile_processed$metrics_v) #26,number of tiles with raster object |
|
| 822 |
length(df_tile_processed$metrics_v) #26,number of tiles in the region |
|
| 823 |
sum(df_tile_processed$metrics_v)/length(df_tile_processed$metrics_v) #80 of tiles with info |
|
| 824 |
|
|
| 825 |
#coordinates |
|
| 826 |
#coordinates(summary_metrics_v) <- c("lon","lat")
|
|
| 827 |
coordinates(summary_metrics_v) <- c("lon.y","lat.y")
|
|
| 828 |
|
|
| 829 |
threshold_missing_day <- c(367,365,300,200) |
|
| 830 |
|
|
| 831 |
nb<-nrow(subset(summary_metrics_v,model_name=="mod1")) |
|
| 832 |
sum(subset(summary_metrics_v,model_name=="mod1")$n_missing)/nb #33/35 |
|
| 833 |
|
|
| 834 |
## Make this a figure... |
|
| 835 |
|
|
| 836 |
#plot(summary_metrics_v) |
|
| 837 |
#Make this a function later so that we can explore many thresholds... |
|
| 838 |
|
|
| 839 |
j<-1 #for model name 1 |
|
| 840 |
for(i in 1:length(threshold_missing_day)){
|
|
| 841 |
|
|
| 842 |
#summary_metrics_v$n_missing <- summary_metrics_v$n == 365 |
|
| 843 |
#summary_metrics_v$n_missing <- summary_metrics_v$n < 365 |
|
| 844 |
summary_metrics_v$n_missing <- summary_metrics_v$n < threshold_missing_day[i] |
|
| 845 |
summary_metrics_v_subset <- subset(summary_metrics_v,model_name=="mod1") |
|
| 846 |
|
|
| 847 |
#res_pix <- 1200 |
|
| 848 |
res_pix <- 960 |
|
| 849 |
|
|
| 850 |
col_mfrow <- 1 |
|
| 851 |
row_mfrow <- 1 |
|
| 852 |
|
|
| 853 |
png(filename=paste("Figure7a_ac_metrics_map_centroids_tile_",model_name[j],"_","missing_day_",threshold_missing_day[i],
|
|
| 854 |
"_",out_suffix,".png",sep=""), |
|
| 855 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 856 |
|
|
| 857 |
model_name[j] |
|
| 858 |
|
|
| 859 |
p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black')) |
|
| 860 |
#title("(a) Mean for 1 January")
|
|
| 861 |
p <- bubble(summary_metrics_v_subset,"n_missing",main=paste("Missing per tile and by ",model_name[j]," for ",
|
|
| 862 |
threshold_missing_day[i])) |
|
| 863 |
p1 <- p+p_shp |
|
| 864 |
print(p1) |
|
| 865 |
#plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T) |
|
| 866 |
#title(paste("Averrage RMSE per tile and by ",model_name[i]))
|
|
| 867 |
|
|
| 868 |
dev.off() |
|
| 869 |
} |
|
| 870 |
|
|
| 871 |
###################### |
|
| 872 |
### Figure 7: Number of predictions: daily and monthly |
|
| 873 |
|
|
| 874 |
#xyplot(rmse~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 875 |
# pred_mod!="mod_kr"),type="b") |
|
| 876 |
|
|
| 877 |
#xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 878 |
# pred_mod!="mod_kr"),type="h") |
|
| 879 |
|
|
| 880 |
#cor |
|
| 881 |
|
|
| 882 |
# |
|
| 883 |
## Figure 7a |
|
| 884 |
png(filename=paste("Figure7a_number_daily_predictions_per_models","_",out_suffix,".png",sep=""),
|
|
| 885 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 886 |
|
|
| 887 |
xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 888 |
pred_mod!="mod_kr"),type="h") |
|
| 889 |
dev.off() |
|
| 890 |
|
|
| 891 |
table(tb$pred_mod) |
|
| 892 |
table(tb$index_d) |
|
| 893 |
table(subset(tb,pred_mod!="mod_kr")) |
|
| 894 |
table(subset(tb,pred_mod=="mod1")$index_d) |
|
| 895 |
#aggregate() |
|
| 896 |
tb$predicted <- 1 |
|
| 897 |
test <- aggregate(predicted~pred_mod+tile_id,data=tb,sum) |
|
| 898 |
xyplot(predicted~pred_mod | tile_id,data=subset(as.data.frame(test), |
|
| 899 |
pred_mod!="mod_kr"),type="h") |
|
| 900 |
|
|
| 901 |
test |
|
| 902 |
|
|
| 903 |
as.character(unique(test$tile_id)) #141 tiles |
|
| 904 |
|
|
| 905 |
dim(subset(test,test$predicted==365 & test$pred_mod=="mod1")) |
|
| 906 |
histogram(subset(test, test$pred_mod=="mod1")$predicted) |
|
| 907 |
unique(subset(test, test$pred_mod=="mod1")$predicted) |
|
| 908 |
table((subset(test, test$pred_mod=="mod1")$predicted)) |
|
| 909 |
|
|
| 910 |
#LST_avgm_min <- aggregate(LST~month,data=data_month_all,min) |
|
| 911 |
histogram(test$predicted~test$tile_id) |
|
| 912 |
#table(tb) |
|
| 913 |
## Figure 7b |
|
| 914 |
#png(filename=paste("Figure7b_number_daily_predictions_per_models","_",out_suffix,".png",sep=""),
|
|
| 915 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 916 |
|
|
| 917 |
#xyplot(n~month | tile_id + pred_mod,data=subset(as.data.frame(tb_month_s), |
|
| 918 |
# pred_mod!="mod_kr"),type="h") |
|
| 919 |
#xyplot(n~month | tile_id,data=subset(as.data.frame(tb_month_s), |
|
| 920 |
# pred_mod="mod_1"),type="h") |
|
| 921 |
#test=subset(as.data.frame(tb_month_s),pred_mod="mod_1") |
|
| 922 |
#table(tb_month_s$month) |
|
| 923 |
#dev.off() |
|
| 924 |
# |
|
| 925 |
|
|
| 926 |
|
|
| 927 |
########################################################## |
|
| 928 |
##### Figure 8: Breaking down accuracy by regions!! ##### |
|
| 929 |
|
|
| 930 |
#summary_metrics_v <- merge(summary_metrics_v,df_tile_processed,by="tile_id") |
|
| 931 |
|
|
| 932 |
## Figure 8a |
|
| 933 |
res_pix <- 480 |
|
| 934 |
col_mfrow <- 1 |
|
| 935 |
row_mfrow <- 1 |
|
| 936 |
|
|
| 937 |
png(filename=paste("Figure8a_boxplot_overall_separated_by_region_with_oultiers_",model_name[i],"_",out_suffix,".png",sep=""),
|
|
| 938 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 939 |
|
|
| 940 |
p<- bwplot(rmse~pred_mod | reg, data=tb, |
|
| 941 |
main="RMSE per model and region over all tiles") |
|
| 942 |
print(p) |
|
| 943 |
dev.off() |
|
| 944 |
|
|
| 945 |
## Figure 8b |
|
| 946 |
png(filename=paste("Figure8b_boxplot_overall_separated_by_region_scaling_",model_name[i],"_",out_suffix,".png",sep=""),
|
|
| 947 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 948 |
|
|
| 949 |
boxplot(rmse~pred_mod,data=tb,ylim=c(0,5),outline=FALSE)#,names=tb$pred_mod) |
|
| 950 |
title("RMSE per model over all tiles")
|
|
| 951 |
p<- bwplot(rmse~pred_mod | reg, data=tb,ylim=c(0,5), |
|
| 952 |
main="RMSE per model and region over all tiles") |
|
| 953 |
print(p) |
|
| 954 |
dev.off() |
|
| 955 |
|
|
| 956 |
##################################################### |
|
| 957 |
#### Figure 9: plotting mosaics of regions ########### |
|
| 958 |
## plot mosaics... |
|
| 959 |
|
|
| 960 |
#First collect file names |
|
| 961 |
|
|
| 962 |
|
|
| 963 |
#names(lf_mosaics_reg) <- l_reg_name |
|
| 964 |
|
|
| 965 |
#This part should be automated... |
|
| 966 |
#plot Australia |
|
| 967 |
#lf_m <- lf_mosaics_reg[[2]] |
|
| 968 |
#out_dir_str <- out_dir |
|
| 969 |
#reg_name <- "reg6_1000x3000" |
|
| 970 |
#lapply() |
|
| 971 |
#list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str,out_suffix=out_suffix) |
|
| 972 |
#list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str,out_suffix=out_suffix,l_dates=day_to_mosaic) |
|
| 973 |
|
|
| 974 |
#lf_m_mask_reg4_1500x4500 <- mclapply(1:2,FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 6) |
|
| 975 |
#debug(plot_daily_mosaics) |
|
| 976 |
#lf_m_mask_reg6_1000x3000 <- plot_daily_mosaics(1,list_param=list_param_plot_daily_mosaics) |
|
| 977 |
|
|
| 978 |
#lf_m_mask_reg6_1000x3000 <- mclapply(1:length(lf_m),FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 10) |
|
| 979 |
|
|
| 980 |
|
|
| 981 |
################## WORLD MOSAICS NEEDS MAJOR CLEAN UP OF CODE HERE |
|
| 982 |
##make functions!! |
|
| 983 |
###Combine mosaics with modified code from Alberto |
|
| 984 |
|
|
| 985 |
#use list from above!! |
|
| 986 |
|
|
| 987 |
# test_list <-list.files(path=file.path(out_dir,"mosaics"), |
|
| 988 |
# pattern=paste("^world_mosaics.*.tif$",sep=""),
|
|
| 989 |
# ) |
|
| 990 |
# #world_mosaics_mod1_output1500x4500_km_20101105_run10_1500x4500_global_analyses_03112015.tif |
|
| 991 |
# |
|
| 992 |
# #test_list<-lapply(1:30,FUN=function(i){lapply(1:x[[i]]},x=lf_mosaics_mask_reg)
|
|
| 993 |
# #test_list<-unlist(test_list) |
|
| 994 |
# #mosaic_list_mean <- vector("list",length=1)
|
|
| 995 |
# mosaic_list_mean <- test_list |
|
| 996 |
# out_rastnames <- "world_test_mosaic_20100101" |
|
| 997 |
# out_path <- out_dir |
|
| 998 |
# |
|
| 999 |
# list_param_mosaic<-list(mosaic_list_mean,out_path,out_rastnames,file_format,NA_flag_val,out_suffix) |
|
| 1000 |
# names(list_param_mosaic)<-c("mosaic_list","out_path","out_rastnames","file_format","NA_flag_val","out_suffix")
|
|
| 1001 |
# #mean_m_list <-mclapply(1:12, list_param=list_param_mosaic, mosaic_m_raster_list,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement |
|
| 1002 |
# |
|
| 1003 |
# lf <- mosaic_m_raster_list(1,list_param_mosaic) |
|
| 1004 |
# |
|
| 1005 |
# debug(mosaic_m_raster_list) |
|
| 1006 |
# mosaic_list<-list_param$mosaic_list |
|
| 1007 |
# out_path<-list_param$out_path |
|
| 1008 |
# out_names<-list_param$out_rastnames |
|
| 1009 |
# file_format <- list_param$file_format |
|
| 1010 |
# NA_flag_val <- list_param$NA_flag_val |
|
| 1011 |
# out_suffix <- list_param$out_suffix |
|
| 1012 |
|
|
| 1013 |
##Now mosaic for mean: should reorder files!! |
|
| 1014 |
#out_rastnames_mean<-paste("_",lst_pat,"_","mean",out_suffix,sep="")
|
|
| 1015 |
#list_date_names<-c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
|
|
| 1016 |
#mosaic_list<-split(tmp_str1,list_date_names) |
|
| 1017 |
#new_list<-vector("list",length(mosaic_list))
|
|
| 1018 |
# for (i in 1:length(list_date_names)){
|
|
| 1019 |
# j<-grep(list_date_names[i],mosaic_list,value=FALSE) |
|
| 1020 |
# names(mosaic_list)[j]<-list_date_names[i] |
|
| 1021 |
# new_list[i]<-mosaic_list[j] |
|
| 1022 |
# } |
|
| 1023 |
# mosaic_list_mean <-new_list #list ready for mosaicing |
|
| 1024 |
# out_rastnames_mean<-paste(list_date_names,out_rastnames_mean,sep="") |
|
| 1025 |
|
|
| 1026 |
### Now mosaic tiles...Note that function could be improved to use less memory |
|
| 1027 |
|
|
| 1028 |
|
|
| 1029 |
################## PLOTTING WORLD MOSAICS ################ |
|
| 1030 |
|
|
| 1031 |
#lf_world_pred <-list.files(path=file.path(out_dir,"mosaics"), |
|
| 1032 |
# pattern=paste("^world_mosaics.*.tif$",sep=""),full.names=T)
|
|
| 1033 |
|
|
| 1034 |
lf_world_pred <-list.files(path=file.path(out_dir,"mosaics"), |
|
| 1035 |
pattern=paste("^reg5.*.",".tif$",sep=""),full.names=T)
|
|
| 1036 |
l_reg_name <- unique(df_tile_processed$reg) |
|
| 1037 |
lf_world_pred <-list.files(path=file.path(out_dir,l_reg_name[[i]]), |
|
| 1038 |
pattern=paste(".tif$",sep=""),full.names=T)
|
|
| 1039 |
|
|
| 1040 |
#mosaic_list_mean <- test_list |
|
| 1041 |
#out_rastnames <- "world_test_mosaic_20100101" |
|
| 1042 |
#out_path <- out_dir |
|
| 1043 |
|
|
| 1044 |
#lf_world_pred <- list.files(pattern="world.*2010090.*.tif$") |
|
| 1045 |
#lf_raster_fname <- list.files(pattern="world.*2010*.*02162015.tif$",full.names=T) |
|
| 1046 |
lf_raster_fname <- lf_world_pred |
|
| 1047 |
prefix_str <- "Figure10_clim_world_mosaics_day_" |
|
| 1048 |
l_dates <-day_to_mosaic |
|
| 1049 |
screenRast=FALSE |
|
| 1050 |
list_param_plot_screen_raster <- list(lf_raster_fname,screenRast,l_dates,out_dir,prefix_str,out_suffix) |
|
| 1051 |
names(list_param_plot_screen_raster) <- c("lf_raster_fname","screenRast","l_dates","out_dir_str","prefix_str","out_suffix_str")
|
|
| 1052 |
|
|
| 1053 |
debug(plot_screen_raster_val) |
|
| 1054 |
|
|
| 1055 |
world_m_list1<- plot_screen_raster_val(1,list_param_plot_screen_raster) |
|
| 1056 |
#world_m_list <- mclapply(11:30, list_param=list_param_plot_screen_raster, plot_screen_raster_val,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement |
|
| 1057 |
world_m_list <- mclapply(1:length(l_dates), list_param=list_param_plot_screen_raster, plot_screen_raster_val,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement |
|
| 1058 |
|
|
| 1059 |
s_pred <- stack(lf_raster_fname) |
|
| 1060 |
|
|
| 1061 |
res_pix <- 1500 |
|
| 1062 |
col_mfrow <- 3 |
|
| 1063 |
row_mfrow <- 2 |
|
| 1064 |
|
|
| 1065 |
png(filename=paste("Figure10_levelplot_combined_",region_name,"_",out_suffix,".png",sep=""),
|
|
| 1066 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 1067 |
|
|
| 1068 |
levelplot(s_pred,layers=1:6,col.regions=rev(terrain.colors(255)),cex=4) |
|
| 1069 |
|
|
| 1070 |
dev.off() |
|
| 1071 |
|
|
| 1072 |
# blues<- designer.colors(6, c( "blue", "white") ) |
|
| 1073 |
# reds <- designer.colors(6, c( "white","red") ) |
|
| 1074 |
# colorTable<- c( blues[-6], reds) |
|
| 1075 |
# breaks with a gap of 10 to 17 assigned the white color |
|
| 1076 |
# brks<- c(seq( 1, 10,,6), seq( 17, 25,,6)) |
|
| 1077 |
# image.plot( x,y,z,breaks=brks, col=colorTable) |
|
| 1078 |
# |
|
| 1079 |
|
|
| 1080 |
#lf_world_mask_reg <- vector("list",length=length(lf_world_pred))
|
|
| 1081 |
#for(i in 1:length(lf_world_pred)){
|
|
| 1082 |
# |
|
| 1083 |
# lf_m <- lf_mosaics_reg[i] |
|
| 1084 |
# out_dir_str <- out_dir |
|
| 1085 |
#reg_name <- paste(l_reg_name[i],"_",tile_size,sep="") #make this automatic |
|
| 1086 |
#lapply() |
|
| 1087 |
# list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=tile_size,out_dir_str=out_dir_str,out_suffix=out_suffix,l_dates=day_to_mosaic) |
|
| 1088 |
#lf_m_mask_reg4_1500x4500 <- mclapply(1:2,FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 6) |
|
| 1089 |
|
|
| 1090 |
#lf_world_mask_reg[[i]] <- mclapply(1:length(lf_m),FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 10) |
|
| 1091 |
} |
|
| 1092 |
|
|
| 1093 |
############# NEW MASK AND DATA |
|
| 1094 |
## Plot areas and day predicted as first check |
|
| 1095 |
|
|
| 1096 |
l_reg_name <- unique(df_tile_processed$reg) |
|
| 1097 |
#(subset(df_tile_processed$reg == l_reg_name[i],date) |
|
| 1098 |
|
|
| 1099 |
for(i in 1:length(l_reg_name)){
|
|
| 1100 |
lf_world_pred<-list.files(path=file.path(out_dir,l_reg_name[[i]]), |
|
| 1101 |
pattern=paste(".tif$",sep=""),full.names=T)
|
|
| 1102 |
|
|
| 1103 |
#mosaic_list_mean <- test_list |
|
| 1104 |
#out_rastnames <- "world_test_mosaic_20100101" |
|
| 1105 |
#out_path <- out_dir |
|
| 1106 |
|
|
| 1107 |
#lf_world_pred <- list.files(pattern="world.*2010090.*.tif$") |
|
| 1108 |
#lf_raster_fname <- list.files(pattern="world.*2010*.*02162015.tif$",full.names=T) |
|
| 1109 |
lf_raster_fname <- lf_world_pred |
|
| 1110 |
prefix_str <- paste("Figure10_",l_reg_name[i],sep="")
|
|
| 1111 |
|
|
| 1112 |
l_dates <- basename(lf_raster_fname) |
|
| 1113 |
tmp_name <- gsub(".tif","",l_dates)
|
|
| 1114 |
tmp_name <- gsub("gam_CAI_dailyTmax_predicted_mod1_0_1_","",tmp_name)
|
|
| 1115 |
#l_dates <- tmp_name |
|
| 1116 |
l_dates <- paste(l_reg_name[i],"_",tmp_name,sep="") |
|
| 1117 |
|
|
| 1118 |
screenRast=TRUE |
|
| 1119 |
list_param_plot_screen_raster <- list(lf_raster_fname,screenRast,l_dates,out_dir,prefix_str,out_suffix) |
|
| 1120 |
names(list_param_plot_screen_raster) <- c("lf_raster_fname","screenRast","l_dates","out_dir_str","prefix_str","out_suffix_str")
|
|
| 1121 |
|
|
| 1122 |
#undebug(plot_screen_raster_val) |
|
| 1123 |
|
|
| 1124 |
#world_m_list1<- plot_screen_raster_val(1,list_param_plot_screen_raster) |
|
| 1125 |
#world_m_list <- mclapply(11:30, list_param=list_param_plot_screen_raster, plot_screen_raster_val,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement |
|
| 1126 |
world_m_list <- mclapply(1:length(l_dates), list_param=list_param_plot_screen_raster, plot_screen_raster_val,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement |
|
| 1127 |
|
|
| 1128 |
#s_pred <- stack(lf_raster_fname) |
|
| 1129 |
|
|
| 1130 |
#res_pix <- 1500 |
|
| 1131 |
#col_mfrow <- 3 |
|
| 1132 |
#row_mfrow <- 2 |
|
| 1133 |
|
|
| 1134 |
#png(filename=paste("Figure10_levelplot_combined_",region_name,"_",out_suffix,".png",sep=""),
|
|
| 1135 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 1136 |
|
|
| 1137 |
#levelplot(s_pred,layers=1:6,col.regions=rev(terrain.colors(255)),cex=4) |
|
| 1138 |
|
|
| 1139 |
#dev.off() |
|
| 1140 |
} |
|
| 1141 |
|
|
| 1142 |
|
|
| 1143 |
|
|
| 1144 |
##################### END OF SCRIPT ###################### |
|
Also available in: Unified diff
assessment part2 splitting functions into new script