Project

General

Profile

« Previous | Next » 

Revision 42c025ec

Added by Benoit Parmentier over 9 years ago

global assessement part 1, changes to accomadate additional tiles for 1500x4500km tiles

View differences:

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/15/2015
8 
#MODIFIED ON: 04/24/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/out_15x45/" #PARAM1
74 
i#n_dir1b <- "/nobackupp6/aguzman4/climateLayers/output1500x4500_km/singles" #PARAM1, add for now in_dir1 can be a list...
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...
75 75 

  		




  76
  
  
    
region_names <- c("reg5") #selected region names, #PARAM2


  		




  77
  
  
    
#region_namesb <- c("reg_1b","reg_2b","reg_6b") #selected region names, #PARAM2


  		




  
  76
  
    

  		




  
  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
  79
  
    

  		




  79
  80
  
    
y_var_name <- "dailyTmax" #PARAM3


  		




  80
  81
  
    
interpolation_method <- c("gam_CAI") #PARAM4


  		




  81
  
  
    
out_prefix<-"run10_1500x4500_global_analyses_pred_2003_04102015" #PARAM5


  		




  
  82
  
    
out_prefix<-"run10_1500x4500_global_analyses_04172015" #PARAM5


  		




  82
  83
  
    

  		




  83
  84
  
    
#out_dir<-"/data/project/layers/commons/NEX_data/" #On NCEAS Atlas


  		




  84
  85
  
    
#out_dir <- "/nobackup/bparmen1/" #on NEX


  		




  ......




  89
  90
  
    
CRS_locs_WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84, #PARAM8


  		




  90
  91
  
    

  		




  91
  92
  
    
#day_to_mosaic <- c("20100101","20100901") #PARAM9


  		




  92
  
  
    

  		




  93
  
  
    
day_to_mosaic <- c("20030101","20030102","20030103","20030104","20030105",


  		




  94
  
  
    
                   "20030301","20030302","20030303","20030304","20030305",


  		




  95
  
  
    
                   "20030501","20030502","20030503","20030504","20030505",


  		




  96
  
  
    
                   "20030701","20030702","20030703","20030704","20030705",


  		




  97
  
  
    
                   "20030901","20030902","20030903","20030904","20030905",


  		




  98
  
  
    
                   "20031101","20031102","20031103","20031104","20031105") #PARAM7


  		




  99
  
  
    

  		




  100
  
  
    
#day_to_mosaic <- NULL #if day to mosaic is null then mosaic all dates?


  		




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


  		




  
  99
  
    
day_to_mosaic <- NULL #if day to mosaic is null then mosaic all dates?


  		




  101
  100
  
    

  		




  102
  101
  
    
file_format <- ".tif" #format for mosaiced files #PARAM10


  		




  103
  102
  
    
NA_flag_val <- -9999  #No data value, #PARAM11


  		




  ......




  141
  140
  
    

  		




  142
  141
  
    
## load problematic tiles


  		




  143
  142
  
    

  		




  144
  
  
    
# in_dir_listb <- list.dirs(path=in_dir1b,recursive=FALSE) #get the list regions processed for this run


  		




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


  		




  147
  
  
    
# 


  		




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


  		




  154
  
  
    
# #select only directories used for predictions


  		




  155
  
  
    
# in_dir_regb <- in_dir_listb[grep(".*._.*.",basename(in_dir_listb),invert=FALSE)] #select directory with shapefiles...


  		




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


  		




  158
  
  
    
#     


  		




  159
  
  
    
# in_dir_listb <- in_dir_listb[grep("bak",basename(basename(in_dir_listb)),invert=TRUE)] #the first one is the in_dir1


  		




  160
  
  
    
# #list of shapefiles used to define tiles


  		




  161
  
  
    
# in_dir_shp_listb <- list.files(in_dir_shpb,".shp",full.names=T)


  		




  
  143
  
    
in_dir_listb <- list.dirs(path=in_dir1b,recursive=FALSE) #get the list regions processed for this run


  		




  
  144
  
    
#basename(in_dir_list)


  		




  
  145
  
    
in_dir_listb<- lapply(region_namesb,FUN=function(x,y){y[grep(x,basename(y),invert=FALSE)]},y=in_dir_listb) 


  		




  
  146
  
    

  		




  
  147
  
    
in_dir_list_allb  <- lapply(in_dir_listb,function(x){list.dirs(path=x,recursive=F)})


  		




  
  148
  
    
in_dir_listb <- unlist(in_dir_list_allb)


  		




  
  149
  
    
#in_dir_list <- in_dir_list[grep("bak",basename(basename(in_dir_list)),invert=TRUE)] #the first one is the in_dir1


  		




  
  150
  
    
in_dir_subsetb <- in_dir_listb[grep("subset",basename(in_dir_listb),invert=FALSE)] #select directory with shapefiles...


  		




  
  151
  
    
in_dir_shpb <- file.path(in_dir_subsetb,"shapefiles")


  		




  
  152
  
    

  		




  
  153
  
    
#select only directories used for predictions


  		




  
  154
  
    
in_dir_regb <- in_dir_listb[grep(".*._.*.",basename(in_dir_listb),invert=FALSE)] #select directory with shapefiles...


  		




  
  155
  
    
#in_dir_reg <- in_dir_list[grep("july_tiffs",basename(in_dir_reg),invert=TRUE)] #select directory with shapefiles...


  		




  
  156
  
    
in_dir_listb <- in_dir_regb


  		




  
  157
  
    
    


  		




  
  158
  
    
in_dir_listb <- in_dir_listb[grep("bak",basename(basename(in_dir_listb)),invert=TRUE)] #the first one is the in_dir1


  		




  
  159
  
    
#list of shapefiles used to define tiles


  		




  
  160
  
    
in_dir_shp_listb <- list.files(in_dir_shpb,".shp",full.names=T)


  		




  162
  161
  
    

  		




  163
  162
  
    

  		




  164
  163
  
    
#### Combine now...


  		




  165
  164
  
    

  		




  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)


  		




  
  165
  
    
in_dir_list <- c(in_dir_list,in_dir_listb)


  		




  
  166
  
    
in_dir_reg <- c(in_dir_reg,in_dir_regb)


  		




  
  167
  
    
in_dir_shp <- c(in_dir_shp,in_dir_shpb)


  		




  
  168
  
    
in_dir_shp_list <- c(in_dir_shp_list,in_dir_shp_listb)


  		




  170
  169
  
    
#in_dir_list <- c(in_dir_list,in_dir_listb)


  		




  171
  170
  
    

  		




  172
  171
  
    
#system("ls /nobackup/bparmen1")


  		




  ......




  225
  224
  
    
#Get the number of models predicted


  		




  226
  225
  
    
nb_mod <- length(unique(robj1$tb_diagnostic_v$pred_mod))


  		




  227
  226
  
    
list_formulas <- (robj1$clim_method_mod_obj[[1]]$formulas)


  		




  
  227
  
    
dates_predicted <- (unique(robj1$tb_diagnostic_v$date))


  		




  
  228
  
    

  		




  228
  229
  
    

  		




  229
  230
  
    
#list_tb_diagnostic_v <- mclapply(lf_validation_obj,FUN=function(x){try( x<- load_obj(x)); try(extract_from_list_obj(x,"metrics_v"))},mc.preschedule=FALSE,mc.cores = 6)                           


  		




  230
  231
  
    
#names(list_tb_diagnostic_v) <- list_names_tile_id


  		




  ......




  466
  467
  
    

  		




  467
  468
  
    
#dates_l <- unique(robj1$tb_diagnostic_s$date) #list of dates to query tif


  		




  468
  469
  
    
#create date!!!


  		




  469
  
  
    
idx <- seq(as.Date('2010-01-01'), as.Date('2010-12-31'), 'day')


  		




  470
  
  
    
#idx <- seq(as.Date('20100101'), as.Date('20101231'), 'day')


  		




  471
  
  
    
#date_l <- strptime(idx[1], "%Y%m%d") # interpolation date being processed


  		




  472
  
  
    
dates_l <- format(idx, "%Y%m%d") # interpolation date being processed


  		




  
  470
  
    
if(is.null(day_to_mosaic)){


  		




  
  471
  
    
  


  		




  
  472
  
    
  #idx <- seq(as.Date('2010-01-01'), as.Date('2010-12-31'), 'day')


  		




  
  473
  
    
  #idx <- seq(as.Date('20100101'), as.Date('20101231'), 'day')


  		




  
  474
  
    
  #date_l <- strptime(idx[1], "%Y%m%d") # interpolation date being processed


  		




  
  475
  
    
  #dates_l <- format(idx, "%Y%m%d") # interpolation date being processed


  		




  
  476
  
    
  day_to_mosaic <- dates_predicted #should be 365 days...


  		




  
  477
  
    
  #l_dates <- day_to_mosaic


  		




  
  478
  
    
}


  		




  
  479
  
    
#else{


  		




  
  480
  
    
#  l_dates <- paste(day_to_mosaic,collapse=",")


  		




  
  481
  
    
#}


  		




  473
  482
  
    

  		




  474
  483
  
    
## make this a function? report on number of tiles used for mosaic...


  		




  475
  484
  
    

  		




  ......




  482
  491
  
    
#system("MODULEPATH=$MODULEPATH:/nex/modules/files")


  		




  483
  492
  
    
#system("module load /nex/modules/files/pythonkits/gdal_1.10.0_python_2.7.3_nex")


  		




  484
  493
  
    

  		




  485
  
  
    
module_path <- ""


  		




  
  494
  
    
#module_path <- ""


  		




  486
  495
  
    
module_path <- "/nobackupp6/aguzman4/climateLayers/sharedCode/"


  		




  487
  496
  
    
#/nobackupp6/aguzman4/climateLayers/sharedCode/mosaicUsingGdalMerge.py


  		




  488
  
  
    
#l_dates <- paste(day_to_mosaic,collapse=",",sep=" ")


  		




  489
  
  
    
l_dates <- paste(day_to_mosaic,collapse=",")


  		




  
  497
  
    
l_dates <- paste(day_to_mosaic,collapse=",",sep=" ")


  		




  490
  498
  
    
## use region 2 first


  		




  491
  499
  
    

  		




  492
  500
  
    
### FIRST mosaics by processing region


  		




  ......




  503
  511
  
    
#MODULEPATH=$MODULEPATH:/nex/modules/files


  		




  504
  512
  
    
#module load pythonkits/gdal_1.10.0_python_2.7.3_nex


  		




  505
  513
  
    

  		




  506
  
  
    
for (i in 1:length(region_names)){


  		




  507
  
  
    
  in_dir_mosaics <- file.path(in_dir1,region_names[i])


  		




  
  514
  
    
for (j in 1:length(region_names)){


  		




  
  515
  
    
  in_dir_mosaics <- file.path(in_dir1,region_names[j])


  		




  508
  516
  
    
  #out_dir_mosaics <- "/nobackupp6/aguzman4/climateLayers/output1000x3000_km/reg5/mosaicsMean"


  		




  509
  517
  
    
  #Can be changed to have mosaics in different dir..


  		




  510
  518
  
    
  out_dir_mosaics <- out_dir


  		




  ......




  512
  520
  
    
  #tile_size <- basename(dirname(in_dir[[i]]))


  		




  513
  521
  
    
  tile_size <- basename(in_dir1)


  		




  514
  522
  
    

  		




  515
  
  
    
  prefix_str <- paste(region_names[i],"_",tile_size,sep="")


  		




  
  523
  
    
  prefix_str <- paste(region_names[j],"_",tile_size,sep="")


  		




  516
  524
  
    

  		




  517
  525
  
    
  mod_str <- "mod1" #use mod2 which corresponds to model with LST and elev


  		




  518
  526
  
    

  		












Also available in:  Unified diff