Revision c05a4755
Added by Benoit Parmentier over 9 years ago
| climate/research/oregon/interpolation/global_run_scalingup_assessment_part1.R | ||
|---|---|---|
| 5 | 5 |
#Part 1 create summary tables and inputs files for figure in part 2 and part 3. |
| 6 | 6 |
#AUTHOR: Benoit Parmentier |
| 7 | 7 |
#CREATED ON: 03/23/2014 |
| 8 |
#MODIFIED ON: 04/27/2015
|
|
| 8 |
#MODIFIED ON: 05/12/2015
|
|
| 9 | 9 |
#Version: 4 |
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 | 11 |
#TO DO: |
| ... | ... | |
| 70 | 70 |
#reg1 (North Am), reg2(Europe),reg3(Asia), reg4 (South Am), reg5 (Africa), reg6 (Australia-Asia) |
| 71 | 71 |
#master directory containing the definition of tile size and tiles predicted |
| 72 | 72 |
#in_dir1 <- "/nobackupp6/aguzman4/climateLayers/output1000x3000_km/" |
| 73 |
in_dir1 <- "/nobackupp6/aguzman4/climateLayers/output1500x4500_km" #PARAM1
|
|
| 74 |
in_dir1b <- "/nobackupp6/aguzman4/climateLayers/output1500x4500_km/singles" #PARAM1, add for now in_dir1 can be a list... |
|
| 73 |
in_dir1 <- "/nobackupp6/aguzman4/climateLayers/out_15x45" #PARAM1
|
|
| 74 |
#in_dir1b <- "/nobackupp6/aguzman4/climateLayers/output1500x4500_km/singles" #PARAM1, add for now in_dir1 can be a list...
|
|
| 75 | 75 |
|
| 76 |
region_names <- c("reg5","reg5b") #selected region names, #PARAM2
|
|
| 76 | 77 |
|
| 77 |
region_names <- c("reg1","reg2","reg3","reg4","reg5","reg6") #selected region names, #PARAM2
|
|
| 78 |
region_namesb <- c("reg_1b","reg_1c","reg_2b","reg_3b","reg_6b") #selected region names, #PARAM2
|
|
| 78 |
#region_names <- c("reg1","reg2","reg3","reg4","reg5","reg6") #selected region names, #PARAM2
|
|
| 79 |
#region_namesb <- c("reg_1b","reg_1c","reg_2b","reg_3b","reg_6b") #selected region names, #PARAM2
|
|
| 79 | 80 |
|
| 80 | 81 |
y_var_name <- "dailyTmax" #PARAM3 |
| 81 | 82 |
interpolation_method <- c("gam_CAI") #PARAM4
|
| 82 |
out_prefix<-"run10_1500x4500_global_analyses_04172015" #PARAM5 |
|
| 83 |
out_prefix<-"run10_1500x4500_global_analyses_pred_2003_05122015" #PARAM5 |
|
| 84 |
#output_run10_1500x4500_global_analyses_pred_2003_04102015/ |
|
| 83 | 85 |
|
| 84 | 86 |
#out_dir<-"/data/project/layers/commons/NEX_data/" #On NCEAS Atlas |
| 85 | 87 |
#out_dir <- "/nobackup/bparmen1/" #on NEX |
| ... | ... | |
| 90 | 92 |
CRS_locs_WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84, #PARAM8
|
| 91 | 93 |
|
| 92 | 94 |
#day_to_mosaic <- c("20100101","20100901") #PARAM9
|
| 93 |
day_to_mosaic <- c("20100101","20100102","20100103","20100104","20100105",
|
|
| 94 |
"20100301","20100302","20100303","20100304","20100305",
|
|
| 95 |
"20100501","20100502","20100503","20100504","20100505",
|
|
| 96 |
"20100701","20100702","20100703","20100704","20100705",
|
|
| 97 |
"20100901","20100902","20100903","20100904","20100905",
|
|
| 98 |
"20101101","20101102","20101103","20101104","20101105")
|
|
| 95 |
day_to_mosaic <- c("20030101","20030102","20030103","20030104","20030105",
|
|
| 96 |
"20030301","20030302","20030303","20030304","20030305",
|
|
| 97 |
"20030501","20030502","20030503","20030504","20030505",
|
|
| 98 |
"20030701","20030702","20030703","20030704","20030705",
|
|
| 99 |
"20030901","20030902","20030903","20030904","20030905",
|
|
| 100 |
"20031101","20031102","20031103","20031104","20031105")
|
|
| 99 | 101 |
#day_to_mosaic <- NULL #if day to mosaic is null then mosaic all dates? |
| 100 | 102 |
|
| 101 | 103 |
file_format <- ".tif" #format for mosaiced files #PARAM10 |
| ... | ... | |
| 139 | 141 |
#list of shapefiles used to define tiles |
| 140 | 142 |
in_dir_shp_list <- list.files(in_dir_shp,".shp",full.names=T) |
| 141 | 143 |
|
| 142 |
## load problematic tiles |
|
| 144 |
## load problematic tiles or additional runs
|
|
| 143 | 145 |
|
| 144 |
in_dir_listb <- list.dirs(path=in_dir1b,recursive=FALSE) #get the list regions processed for this run |
|
| 146 |
#in_dir_listb <- list.dirs(path=in_dir1b,recursive=FALSE) #get the list regions processed for this run
|
|
| 145 | 147 |
#basename(in_dir_list) |
| 146 |
in_dir_listb<- lapply(region_namesb,FUN=function(x,y){y[grep(x,basename(y),invert=FALSE)]},y=in_dir_listb)
|
|
| 148 |
#in_dir_listb<- lapply(region_namesb,FUN=function(x,y){y[grep(x,basename(y),invert=FALSE)]},y=in_dir_listb)
|
|
| 147 | 149 |
|
| 148 |
in_dir_list_allb <- lapply(in_dir_listb,function(x){list.dirs(path=x,recursive=F)})
|
|
| 149 |
in_dir_listb <- unlist(in_dir_list_allb) |
|
| 150 |
#in_dir_list_allb <- lapply(in_dir_listb,function(x){list.dirs(path=x,recursive=F)})
|
|
| 151 |
#in_dir_listb <- unlist(in_dir_list_allb)
|
|
| 150 | 152 |
#in_dir_list <- in_dir_list[grep("bak",basename(basename(in_dir_list)),invert=TRUE)] #the first one is the in_dir1
|
| 151 |
in_dir_subsetb <- in_dir_listb[grep("subset",basename(in_dir_listb),invert=FALSE)] #select directory with shapefiles...
|
|
| 152 |
in_dir_shpb <- file.path(in_dir_subsetb,"shapefiles") |
|
| 153 |
#in_dir_subsetb <- in_dir_listb[grep("subset",basename(in_dir_listb),invert=FALSE)] #select directory with shapefiles...
|
|
| 154 |
#in_dir_shpb <- file.path(in_dir_subsetb,"shapefiles")
|
|
| 153 | 155 |
|
| 154 | 156 |
#select only directories used for predictions |
| 155 |
in_dir_regb <- in_dir_listb[grep(".*._.*.",basename(in_dir_listb),invert=FALSE)] #select directory with shapefiles...
|
|
| 157 |
#in_dir_regb <- in_dir_listb[grep(".*._.*.",basename(in_dir_listb),invert=FALSE)] #select directory with shapefiles...
|
|
| 156 | 158 |
#in_dir_reg <- in_dir_list[grep("july_tiffs",basename(in_dir_reg),invert=TRUE)] #select directory with shapefiles...
|
| 157 |
in_dir_listb <- in_dir_regb |
|
| 159 |
#in_dir_listb <- in_dir_regb
|
|
| 158 | 160 |
|
| 159 |
in_dir_listb <- in_dir_listb[grep("bak",basename(basename(in_dir_listb)),invert=TRUE)] #the first one is the in_dir1
|
|
| 161 |
#in_dir_listb <- in_dir_listb[grep("bak",basename(basename(in_dir_listb)),invert=TRUE)] #the first one is the in_dir1
|
|
| 160 | 162 |
#list of shapefiles used to define tiles |
| 161 |
in_dir_shp_listb <- list.files(in_dir_shpb,".shp",full.names=T) |
|
| 163 |
#in_dir_shp_listb <- list.files(in_dir_shpb,".shp",full.names=T)
|
|
| 162 | 164 |
|
| 163 | 165 |
|
| 164 | 166 |
#### Combine now... |
| 165 | 167 |
|
| 166 |
in_dir_list <- c(in_dir_list,in_dir_listb) |
|
| 167 |
in_dir_reg <- c(in_dir_reg,in_dir_regb) |
|
| 168 |
in_dir_shp <- c(in_dir_shp,in_dir_shpb) |
|
| 169 |
in_dir_shp_list <- c(in_dir_shp_list,in_dir_shp_listb) |
|
| 168 |
#in_dir_list <- c(in_dir_list,in_dir_listb)
|
|
| 169 |
#in_dir_reg <- c(in_dir_reg,in_dir_regb)
|
|
| 170 |
#in_dir_shp <- c(in_dir_shp,in_dir_shpb)
|
|
| 171 |
#in_dir_shp_list <- c(in_dir_shp_list,in_dir_shp_listb)
|
|
| 170 | 172 |
#in_dir_list <- c(in_dir_list,in_dir_listb) |
| 171 | 173 |
|
| 172 | 174 |
#system("ls /nobackup/bparmen1")
|
| ... | ... | |
| 288 | 290 |
file=file.path(out_dir,paste("tb_diagnostic_v_NA","_",out_prefix,".txt",sep="")),sep=",")
|
| 289 | 291 |
|
| 290 | 292 |
##Take where shutdown took place after pathcing |
| 291 |
summary_metrics_v_NA <- read.table(file=file.path(out_dir,paste("summary_metrics_v2_NA_",out_prefix,".txt",sep="")),sep=",")
|
|
| 293 |
#summary_metrics_v_NA <- read.table(file=file.path(out_dir,paste("summary_metrics_v2_NA_",out_prefix,".txt",sep="")),sep=",")
|
|
| 292 | 294 |
#fname <- file.path(out_dir,paste("summary_metrics_v2_NA_",out_suffix,".txt",sep=""))
|
| 293 |
tb_diagnostic_v_NA <- read.table(file=file.path(out_dir,paste("tb_diagnostic_v_NA","_",out_prefix,".txt",sep="")),sep=",")
|
|
| 295 |
#tb_diagnostic_v_NA <- read.table(file=file.path(out_dir,paste("tb_diagnostic_v_NA","_",out_prefix,".txt",sep="")),sep=",")
|
|
| 294 | 296 |
#tb_diagnostic_s_NA_run10_global_analyses_11302014.txt |
| 295 | 297 |
#tb_s <- read.table(file=file.path(out_dir,paste("tb_diagnostic_s_NA","_",out_suffix,".txt",sep="")),sep=",")
|
| 296 | 298 |
|
| ... | ... | |
| 375 | 377 |
# data_month <- extract_from_list_obj(robj1$clim_method_mod_obj,"data_month") |
| 376 | 378 |
# names(data_month) #this contains LST means (mm_1, mm_2 etc.) as well as TMax and other info |
| 377 | 379 |
# |
| 378 |
data_day_s_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(x$validation_mod_obj[["data_s"]])},mc.preschedule=FALSE,mc.cores = num_cores)
|
|
| 379 |
data_day_v_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(x$validation_mod_obj[["data_v"]])},mc.preschedule=FALSE,mc.cores = num_cores)
|
|
| 380 |
#data_day_s_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(x$validation_mod_obj[["data_s"]])},mc.preschedule=FALSE,mc.cores = num_cores)
|
|
| 381 |
#data_day_v_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(x$validation_mod_obj[["data_v"]])},mc.preschedule=FALSE,mc.cores = num_cores)
|
|
| 380 | 382 |
|
| 381 |
data_day_s_list <- mclapply(list_raster_obj_files[1:6],FUN=function(x){try(x<-load_obj(x));try(x$validation_mod_obj[["data_s"]])},mc.preschedule=FALSE,mc.cores = num_cores)
|
|
| 383 |
#data_day_s_list <- mclapply(list_raster_obj_files[1:6],FUN=function(x){try(x<-load_obj(x));try(x$validation_mod_obj[["data_s"]])},mc.preschedule=FALSE,mc.cores = num_cores)
|
|
| 382 | 384 |
|
| 383 |
data_day_v_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(extract_list_from_list_obj(x$validation_mod_obj,"data_v"))},mc.preschedule=FALSE,mc.cores = num_cores)
|
|
| 384 |
data_day_s_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(extract_list_from_list_obj(x$validation_mod_obj,"data_s"))},mc.preschedule=FALSE,mc.cores = num_cores)
|
|
| 385 |
#data_day_v_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(extract_list_from_list_obj(x$validation_mod_obj,"data_v"))},mc.preschedule=FALSE,mc.cores = num_cores)
|
|
| 386 |
#data_day_s_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(extract_list_from_list_obj(x$validation_mod_obj,"data_s"))},mc.preschedule=FALSE,mc.cores = num_cores)
|
|
| 385 | 387 |
|
| 386 |
list_data_day_v <- try(extract_list_from_list_obj(raster_obj$validation_mod_obj,"data_v")) |
|
| 387 |
list_data_day_s <- try(extract_list_from_list_obj(raster_obj$validation_mod_obj,"data_s")) |
|
| 388 |
sampling_dat_day <- extract_list_from_list_obj(raster_obj$method_mod_obj,"daily_dev_sampling_dat") |
|
| 388 |
#list_data_day_v <- try(extract_list_from_list_obj(raster_obj$validation_mod_obj,"data_v"))
|
|
| 389 |
#list_data_day_s <- try(extract_list_from_list_obj(raster_obj$validation_mod_obj,"data_s"))
|
|
| 390 |
#sampling_dat_day <- extract_list_from_list_obj(raster_obj$method_mod_obj,"daily_dev_sampling_dat")
|
|
| 389 | 391 |
#debug(pred_data_info_fun) |
| 390 | 392 |
#list_pred_data_day_s_info <- pred_data_info_fun(1,list_data=list_data_day_s,pred_mod=pred_mod,sampling_dat_info=sampling_dat_day) |
| 391 | 393 |
list_pred_data_day_s_info <- lapply(1:length(sampling_dat_day),FUN=pred_data_info_fun, |
| ... | ... | |
| 554 | 556 |
#MODULEPATH=$MODULEPATH:/nex/modules/files |
| 555 | 557 |
#module load pythonkits/gdal_1.10.0_python_2.7.3_nex |
| 556 | 558 |
|
| 557 |
for (j in 1:length(region_names)){
|
|
| 558 |
in_dir_mosaics <- file.path(in_dir1,region_names[j]) |
|
| 559 |
#recombine region first: |
|
| 560 |
region_names_mosaic <- list(region_names) |
|
| 561 |
names(region_names_mosaic) <- "reg5" |
|
| 562 |
in_dir_mosaics <- lapply(region_names_mosaic,FUN=function(x){file.path(in_dir1,x)})
|
|
| 563 |
|
|
| 564 |
for (j in 1:length(region_names_mosaic)){
|
|
| 565 |
#for(i in 1:length(region_names)) |
|
| 566 |
in_dir_mosaics <- lapply(region_names_mosaic,FUN=function(x){file.path(in_dir1,x)})
|
|
| 567 |
in_dir_mosaics <- paste(region_names_mosaic,collapse=" ") |
|
| 559 | 568 |
#out_dir_mosaics <- "/nobackupp6/aguzman4/climateLayers/output1000x3000_km/reg5/mosaicsMean" |
| 560 | 569 |
#Can be changed to have mosaics in different dir.. |
| 561 | 570 |
out_dir_mosaics <- out_dir |
| ... | ... | |
| 563 | 572 |
#tile_size <- basename(dirname(in_dir[[i]])) |
| 564 | 573 |
tile_size <- basename(in_dir1) |
| 565 | 574 |
|
| 566 |
prefix_str <- paste(region_names[j],"_",tile_size,sep="")
|
|
| 575 |
prefix_str <- paste(names(region_names_mosaic)[j],"_",tile_size,sep="")
|
|
| 567 | 576 |
|
| 568 | 577 |
mod_str <- "mod1" #use mod2 which corresponds to model with LST and elev |
| 569 | 578 |
|
| ... | ... | |
| 583 | 592 |
# Transform this into a function that takes in a list of files!!! We can skip the region stage to reduce the number of files.. |
| 584 | 593 |
|
| 585 | 594 |
#sh /nobackupp6/aguzman4/climateLayers/sharedCode/shMergeFromFile.sh list_mosaics_20100901.txt world_mosaics_1000x3000_20100901.tif |
| 586 |
|
|
| 595 |
in_dir_str <- in_dir_mosaics |
|
| 596 |
region <- "reg5" #can be the world... |
|
| 587 | 597 |
for (i in 1:length(day_to_mosaic)){
|
| 588 |
pattern_str <- paste("*.",day_to_mosaic[i],".*.tif",sep="")
|
|
| 589 |
lf_day_to_mosaic <- list.files(path=out_dir,pattern=pattern_str,full.names=T) |
|
| 598 |
pattern_str <- paste("*.","predicted_mod1",".*.",day_to_mosaic[i],".*.tif",sep="")
|
|
| 599 |
lf_day_to_mosaic <- lapply(1:length(unlist(in_dir_mosaics)),FUN=function(k){list.files(path=unlist(in_dir_mosaics)[k],pattern=pattern_str,full.names=T,recursive=T)})
|
|
| 600 |
lf_day_to_mosaic <- unlist(lf_day_to_mosaic) |
|
| 590 | 601 |
#write.table(lf_day_to_mosaic,file=file.path(out_dir,paste("list_to_mosaics_",day_to_mosaic[i],".txt",sep="")))
|
| 591 | 602 |
writeLines(lf_day_to_mosaic,con=file.path(out_dir,paste("list_to_mosaics_",day_to_mosaic[i],".txt",sep="")))
|
| 592 | 603 |
in_file_to_mosaics <- file.path(out_dir,paste("list_to_mosaics_",day_to_mosaic[i],".txt",sep=""))
|
| ... | ... | |
| 600 | 611 |
|
| 601 | 612 |
#prefix_str <- paste(region_names[i],"_",tile_size,sep="") |
| 602 | 613 |
mod_str <- "mod1" #use mod2 which corresponds to model with LST and elev |
| 603 |
out_mosaic_name <- paste("world_mosaics_",mod_str,"_",tile_size,"_",day_to_mosaic[i],"_",out_prefix,".tif",sep="")
|
|
| 614 |
out_mosaic_name <- paste(region,"_mosaics_",mod_str,"_",tile_size,"_",day_to_mosaic[i],"_",out_prefix,".tif",sep="")
|
|
| 604 | 615 |
module_path <- "/nobackupp6/aguzman4/climateLayers/sharedCode" #this should be a parameter for the function... |
| 605 | 616 |
cmd_str <- paste("sh", file.path(module_path,"shMergeFromFile.sh"),
|
| 606 | 617 |
in_file_to_mosaics, |
| ... | ... | |
| 663 | 674 |
|
| 664 | 675 |
Atlas_dir <- file.path("/data/project/layers/commons/NEX_data/",basename(out_dir),"mosaics")
|
| 665 | 676 |
Atlas_hostname <- "parmentier@atlas.nceas.ucsb.edu" |
| 666 |
lf_cp_f <- list.files(out_dir,full.names=T,pattern="*world.*.tif")#copy all files can filter later
|
|
| 677 |
lf_cp_f <- list.files(out_dir,full.names=T,pattern="*reg.*.tif")#copy all files can filter later
|
|
| 667 | 678 |
filenames_NEX <- paste(lf_cp_f,collapse=" ") #copy raster prediction object |
| 668 | 679 |
cmd_str <- paste("scp -p",filenames_NEX,paste(Atlas_hostname,Atlas_dir,sep=":"), sep=" ")
|
| 669 | 680 |
system(cmd_str) |
Also available in: Unified diff
global scaling up part 1, using additional tiles for reg5 Africa for year 2003