Revision 0fe1b6fb
Added by Benoit Parmentier over 10 years ago
| climate/research/oregon/interpolation/global_run_scalingup_assessment_part2.R | ||
|---|---|---|
| 5 | 5 |
#Analyses, figures, tables and data are also produced in the script. |
| 6 | 6 |
#AUTHOR: Benoit Parmentier |
| 7 | 7 |
#CREATED ON: 03/23/2014 |
| 8 |
#MODIFIED ON: 08/16/2014
|
|
| 8 |
#MODIFIED ON: 08/28/2014
|
|
| 9 | 9 |
#Version: 3 |
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 |
#COMMENTS: analyses for run 3 global using 2 specific tiles
|
|
| 11 |
#COMMENTS: analyses for run 5 global using 6 specific tiles
|
|
| 12 | 12 |
################################################################################################# |
| 13 | 13 |
|
| 14 | 14 |
### Loading R library and packages |
| ... | ... | |
| 69 | 69 |
#on ATLAS |
| 70 | 70 |
#in_dir1 <- "/data/project/layers/commons/NEX_data/test_run1_03232014/output" #On Atlas |
| 71 | 71 |
#parent output dir : contains subset of the data produced on NEX |
| 72 |
in_dir1 <- "/data/project/layers/commons/NEX_data/output_run4_global_analyses_08142014/output20Deg"
|
|
| 72 |
in_dir1 <- "/data/project/layers/commons/NEX_data/output_run5_global_analyses_08252014/output20Deg"
|
|
| 73 | 73 |
# parent output dir for the curent script analyes |
| 74 | 74 |
#out_dir <-"/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/" #On NCEAS Atlas |
| 75 |
out_dir <-"/data/project/layers/commons/NEX_data/output_run4_global_analyses_08142014/"
|
|
| 75 |
out_dir <-"/data/project/layers/commons/NEX_data/output_run5_global_analyses_08252014/"
|
|
| 76 | 76 |
# input dir containing shapefiles defining tiles |
| 77 |
in_dir_shp <- "/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/output/subset/shapefiles"
|
|
| 77 |
#in_dir_shp <- "/data/project/layers/commons/NEX_data/output_run5_global_analyses_08252014/output/subset/shapefiles"
|
|
| 78 | 78 |
|
| 79 | 79 |
#On NEX |
| 80 | 80 |
#contains all data from the run by Alberto |
| ... | ... | |
| 85 | 85 |
|
| 86 | 86 |
y_var_nay_var_name <- "dailyTmax" |
| 87 | 87 |
interpolation_method <- c("gam_CAI")
|
| 88 |
out_prefix<-"run4_global_analyses_08142014"
|
|
| 88 |
out_prefix<-"run5_global_analyses_08252014"
|
|
| 89 | 89 |
|
| 90 | 90 |
|
| 91 | 91 |
#out_dir <-paste(out_dir,"_",out_prefix,sep="") |
| ... | ... | |
| 110 | 110 |
#lf_tables <- list.files(out_dir,) |
| 111 | 111 |
summary_metrics_v <- read.table(file=file.path(out_dir,paste("summary_metrics_v2_NA_",out_prefix,".txt",sep="")),sep=",")
|
| 112 | 112 |
tb <- read.table(file=file.path(out_dir,paste("tb_diagnostic_v_NA","_",out_prefix,".txt",sep="")),sep=",")
|
| 113 |
tb_s <- read.table(file=file.path(out_dir,paste("tb_diagnostic_s_NA","_",out_prefix,".txt",sep="")),sep=",")
|
|
| 114 |
|
|
| 115 |
tb_month_s_<- read.table(file=file.path(out_dir,paste("tb_month_diagnostic_s_NA","_",out_prefix,".txt",sep="")),sep=",")
|
|
| 116 |
tb_month_s <- read.table("tb_month_diagnostic_s_NA_run5_global_analyses_08252014.txt",sep=",")
|
|
| 113 | 117 |
pred_data_month_info <- read.table(file=file.path(out_dir,paste("pred_data_month_info_",out_prefix,".txt",sep="")),sep=",")
|
| 114 | 118 |
pred_data_day_info <- read.table(file=file.path(out_dir,paste("pred_data_day_info_",out_prefix,".txt",sep="")),sep=",")
|
| 115 | 119 |
df_tile_processed <- read.table(file=file.path(out_dir,paste("df_tile_processed_",out_prefix,".txt",sep="")),sep=",")
|
| ... | ... | |
| 129 | 133 |
|
| 130 | 134 |
### First get background map to display where study area is located |
| 131 | 135 |
#can make this more general later on.. |
| 132 |
if region_name=="USA"{
|
|
| 136 |
if (region_name=="USA"){
|
|
| 133 | 137 |
usa_map <- getData('GADM', country='USA', level=1) #Get US map
|
| 134 | 138 |
#usa_map <- getData('GADM', country=region_name,level=1) #Get US map, this is not working right now
|
| 135 | 139 |
usa_map <- usa_map[usa_map$NAME_1!="Alaska",] #remove Alaska |
| 136 | 140 |
reg_layer <- usa_map[usa_map$NAME_1!="Hawaii",] #remove Hawai |
| 137 | 141 |
} |
| 138 | 142 |
|
| 139 |
if region_name=="world"{
|
|
| 143 |
if (region_name=="world"){
|
|
| 140 | 144 |
#http://www.diva-gis.org/Data |
| 141 | 145 |
countries_shp <-"/data/project/layers/commons/NEX_data/countries.shp" |
| 142 | 146 |
path_to_shp <- dirname(countries_shp) |
| ... | ... | |
| 150 | 154 |
shps_tiles <- vector("list",length(list_shp_reg_files))
|
| 151 | 155 |
#collect info: read in all shapfiles |
| 152 | 156 |
#This is slow...make a function and use mclapply?? |
| 157 |
|
|
| 153 | 158 |
for(i in 1:length(list_shp_reg_files)){
|
| 154 | 159 |
path_to_shp <- dirname(list_shp_reg_files[[i]]) |
| 160 |
path_to_shp <- in_dir1 |
|
| 155 | 161 |
layer_name <- sub(".shp","",basename(list_shp_reg_files[[i]]))
|
| 156 | 162 |
shp1 <- readOGR(path_to_shp, layer_name) |
| 157 | 163 |
#shp1<-readOGR(dirname(list_shp_reg_files[[i]]),sub(".shp","",basename(list_shp_reg_files[[i]])))
|
| ... | ... | |
| 285 | 291 |
|
| 286 | 292 |
|
| 287 | 293 |
##### Diagnostic gam |
| 294 |
#### |
|
| 295 |
gam_diag_10x10 <- read.table("gam_diagnostic_10x10_df_run4_global_analyses_08142014.txt",sep=",")
|
|
| 296 |
gam_diag_10x10$month <- as.factor(gam_diag_10x10$month) |
|
| 288 | 297 |
|
| 289 |
gam_diagnostic_df |
|
| 290 |
|
|
| 291 |
boxplot(rmse~k,data=subset(gam_diagnostic_df,pred_mod=="mod1"),main="mod1 and term=s(lat,lon)",ylab="RMSE_f",xlab="k") |
|
| 292 |
boxplot(rmse~k,data=subset(gam_diagnostic_df,pred_mod=="mod2"),main="mod2 and term=s(lat,lon)",ylab="RMSE_f",xlab="k") |
|
| 293 |
|
|
| 294 |
boxplot(rmse~k,data=subset(gam_diagnostic_df,pred_mod=="mod1"),main="mod1 and term=s(elev_s)",ylab="RMSE_f",xlab="k") |
|
| 295 |
boxplot(rmse~k,data=subset(gam_diagnostic_df,pred_mod=="mod2"),main="mod2 and term=s(elev_s)",ylab="RMSE_f",xlab="k") |
|
| 296 |
|
|
| 297 |
boxplot(rmse~k,data=subset(gam_diagnostic_df,pred_mod=="mod2"),main="mod2 and term=s(LST)",ylab="RMSE_f",xlab="k") |
|
| 298 |
|
|
| 299 |
#boxplot(rmse~pred_mod,data=tb,ylim=c(0,5),outline=FALSE)#,names=tb$pred_mod) |
|
| 300 |
#title("RMSE per model over all tiles")
|
|
| 301 |
#bwplot(rmse~k | term + month,data=subset(gam_diagnostic_df,pred_mod=="mod2"),) |
|
| 302 |
xyplot(rmse~k | term + month,data=subset(gam_diagnostic_df,pred_mod=="mod2"),type="p") |
|
| 303 |
xyplot(rmse~k | term + month,data=subset(gam_diagnostic_df,pred_mod=="mod1"),type="p") |
|
| 304 |
xyplot(rmse~k | term + month,data=subset(gam_diagnostic_df,pred_mod=="mod1" & tile_id=="tile_7"),type="b") |
|
| 305 |
xyplot(rmse~k | term + month,data=subset(gam_diagnostic_df,pred_mod=="mod2" & tile_id=="tile_7"),type="b") |
|
| 306 |
xyplot(rmse~k | term + month,data=subset(gam_diagnostic_df,pred_mod=="mod1" & tile_id=="tile_8"),type="b") |
|
| 307 |
xyplot(rmse~k | term + month,data=subset(gam_diagnostic_df,pred_mod=="mod2" & tile_id=="tile_8"),type="b") |
|
| 308 |
|
|
| 309 |
xyplot(rmse~k | term + month,group=tile_id ,data=subset(gam_diagnostic_df,pred_mod=="mod1"),type="b", |
|
| 310 |
auto.key=list(space = "top", cex=1.0,columns=8)) |
|
| 311 |
dev.off() |
|
| 312 |
|
|
| 313 |
|
|
| 298 |
xyplot(rmse~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 299 |
pred_mod!="mod_kr"),type="b") |
|
| 314 | 300 |
|
| 315 |
plot(n~tile_id,data=gam_diagnostic_df,type="h") |
|
| 301 |
xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 302 |
pred_mod!="mod_kr"),type="h") |
|
| 316 | 303 |
|
| 317 |
## Experimenting |
|
| 318 |
list_mod <- load_obj(lf_diagnostic_obj$tile_1[2])$list_mod |
|
| 319 |
data_training_lf <- as.data.frame(list_mod[[1]]$model) |
|
| 320 |
dim(data_training_lf) #is 13 |
|
| 321 |
mod_t1 <-gam(y_var ~ s(lat,lon,k=6) + s(elev_s,k=4) + s(LST,k=4) , data= data_training_lf) |
|
| 322 |
mod_t2 <-gam(y_var ~ s(lat,lon,k=5) + s(elev_s,k=5) + s(LST,k=5) , data= data_training_lf) |
|
| 304 |
xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 305 |
pred_mod!="mod_kr"),type="h") |
|
| 323 | 306 |
|
| 324 |
list_mod <- load_obj(lf_diagnostic_obj$tile_7[4])$list_mod |
|
| 325 |
data_training_lf <- as.data.frame(list_mod[[1]]$model) |
|
| 326 |
dim(data_training_lf) #is 13 |
|
| 327 |
mod_t1 <-gam(y_var ~ s(lat,lon,k=6) + s(elev_s,k=4) + s(LST,k=4) , data= data_training_lf) |
|
| 328 |
mod_t2 <-gam(y_var ~ s(lat,lon,k=5) + s(elev_s,k=5) + s(LST,k=5) , data= data_training_lf) |
|
| 307 |
xyplot(n~month | tile_id + pred_mod,data=subset(as.data.frame(tb_month_s), |
|
| 308 |
pred_mod!="mod_kr"),type="h") |
|
| 329 | 309 |
|
| 330 | 310 |
# |
| 331 | 311 |
## Figure 3b |
Also available in: Unified diff
run5 assessment NEX part2 with 3 set of k combination