Project

General

Profile

« Previous | Next » 

Revision b765ab2b

Added by Benoit Parmentier about 9 years ago

mosaicing function adding function to produce accuracy metric surface from tile RMSE

View differences:

climate/research/oregon/interpolation/global_run_scalingup_mosaicing_function.R
4 4 
#Different options to explore mosaicing are tested. This script only contains functions.
5 5 
#AUTHOR: Benoit Parmentier
6 6 
#CREATED ON: 04/14/2015
7 
#MODIFIED ON: 11/12/2015
7 
#MODIFIED ON: 12/02/2015
8 8 
#Version: 1
9 9 
#PROJECT: Environmental Layers project
10 10 
#COMMENTS: first commit of function script to test mosaicing using 1500x4500km and other tiles
......
99 99 
  #metric_name: accuracy metric selected to be mapped, RMSE, MAE etc.
100 100 
  #pred_mod_name: model selected, such as mod1, mod_kr etc.
101 101 
  #y_var_name: variable being modeled e.g."dailyTmax", dailyTmin, precip
102 
  #interpolation_method <- list_param$interpolation_method #c("gam_CAI") #PARAM3
102 
  #interpolation_method: names of the interpolation/modeling method
103 103 
  #date_processed: day being processed , e.g. 19920101
104 
  #num_cores : number of cores used in the parallelization
104 105 
  #NA_flag_val: value used as flag in the raster
105 106 
  #file_format: e.g. tif., .rst
106 107 
  #out_dir_str: output directory
107 
  #out_suffix_str: <- list_param$out_suffix
108 
  #out_suffix_str: output suffix
108 109 
  #Outputs:
109 110 
  #raster list of weights and product of wegihts and inuts
110 111 
  #TODO:
......
113 114 
  #
114 115 
  ############
115 116
117 
  #Functions
118 
  create_raster_df_centroids_fun <- function(j,list_param){
119 
    #This function generates raster images from metrics values from a data.frame.
120 
    #The raster layer is assigned a unique value from the pixel at every location.
121 
    #Input Parameters:
122 
    #  #lf: list of raster files
123 
    #df_centroids: data.frame table #fitting or validation table with all days
124 
    #metric_name: accuracy metric selected to be mapped, RMSE, MAE etc.
125 
    #date_processed: day being processed , e.g. 19920101
126 
    #num_cores : number of cores used in the parallelization
127 
    #NA_flag_val: value used as flag in the raster
128 
    #file_format: e.g. tif., .rst
129 
    #out_dir_str: output directory
130 
    #out_suffix_str: output suffix
131 
    #Outputs:
132
133 
    #### PARSE arguments
134
135 
    df_centroids <- list_param$df_centroids
136 
    metric_name <- list_param$metric_name
137 
    #interpolation_method <- list_param$interpolation_method #c("gam_CAI") #PARAM3
138 
    #date_processed <- list_param$metric_name
139 
    #num_cores <- list_param$num_cores
140 
    NA_flag_val <- list_param$NA_flag_val
141 
    file_format <- list_param$file_format
142 
    out_dir_str <- list_param$out_dir
143 
    out_suffix_str <- list_param$out_suffix
144 

  		




  
  145
  
    
    ####### START SCRIPT #####    


  		




  
  146
  
    
    


  		




  
  147
  
    
    inFilename <- df_centroids$files[j]


  		




  
  148
  
    
    r1 <- raster(inFilename)


  		




  
  149
  
    
    r1[] <- df_centroids[[metric_name]][j] #improve this


  		




  
  150
  
    
    #set1f <- function(x){rep(NA, x)}


  		




  
  151
  
    
    #r_init <- init(r_in, fun=set1f)


  		




  
  152
  
    
    lf_tmp <- gsub(file_format,"",lf)


  		




  
  153
  
    
  


  		




  
  154
  
    
    extension_str <- extension(inFilename)


  		




  
  155
  
    
    raster_name_tmp <- gsub(extension_str,"",basename(inFilename))


  		




  
  156
  
    
    outFilename <- file.path(out_dir,paste(raster_name_tmp,"_",metric_name,"_",out_suffix,file_format,sep="")) #for use in function later...


  		




  
  157
  
    
  


  		




  
  158
  
    
    writeRaster(r1, NAflag=NA_flag_val,filename=outFilename,overwrite=TRUE)  


  		




  
  159
  
    
    #list_raster_name[[j]] <- outFilename


  		




  
  160
  
    
    return(outFilename)


  		




  
  161
  
    
  }


  		




  
  162
  
    
  


  		




  
  163
  
    
  ####### PARSE ARGUMENTS


  		




  
  164
  
    
  


  		




  116
  165
  
    
  lf <- list_param$lf[[i]] #list of files to mosaic


  		




  117
  166
  
    
  tb <- list_param$tb #fitting or validation table with all days


  		




  118
  167
  
    
  metric_name <- list_param$metric_name #RMSE, MAE etc.


  		




  ......




  120
  169
  
    
  y_var_name <- list_param$y_var_name #"dailyTmax" #PARAM2


  		




  121
  170
  
    
  interpolation_method <- list_param$interpolation_method #c("gam_CAI") #PARAM3


  		




  122
  171
  
    
  date_processed <- list_param$days_to_process[i]


  		




  
  172
  
    
  num_cores <- list_param$num_cores #number of cores used


  		




  123
  173
  
    
  NA_flag_val <- list_param$NA_flag_val


  		




  124
  174
  
    
  #NAflag,file_format,out_suffix etc...


  		




  125
  175
  
    
  file_format <- list_param$file_format


  		




  ......




  128
  178
  
    
   


  		




  129
  179
  
    
  ####### START SCRIPT #####


  		




  130
  180
  
    
  


  		




  131
  
  
    
  #r_in <- raster(lf[i]) #input image


  		




  132
  
  
    

  		




  133
  
  
    
  #date_processed <- day_to_mosaic[i]


  		




  134
  
  
    
  #lf_to_mosaic <-list.files(path=file.path(in_dir_tiles),    


  		




  135
  
  
    
  #         pattern=paste(".*.",date_processed,".*.tif$",sep=""),full.names=T) #choosing date 2...20100901


  		




  136
  
  
    

  		




  137
  181
  
    
  lf_tmp <- gsub(file_format,"",lf)


  		




  138
  182
  
    
  tx<-strsplit(as.character(lf_tmp),"_")


  		




  139
  183
  
    
  #deal with the fact that we have number "1" attached to the out_suffix (centroids of tiles)


  		




  ......




  156
  200
  
    
  tb_date$tile_coord <- as.character(tb_date$tile_coord)


  		




  157
  201
  
    
  df_centroids <- merge(df_centroids,tb_date,by="tile_coord")


  		




  158
  202
  
    

  		




  159
  
  
    
  #r1 <- raster(lf[i])


  		




  160
  
  
    
  


  		




  161
  
  
    
  #function(j,df_centroids,metric_name,date_processed,interpolation,)


  		




  162
  
  
    
  #loop over files, make this a function later, works for now


  		




  163
  203
  
    
  #use mclapply  


  		




  164
  
  
    
  list_raster_name <- vector("list",length=length(lf))


  		




  165
  
  
    
  


  		




  166
  
  
    
  for(j in 1:length(lf)){ 


  		




  167
  
  
    
    


  		




  168
  
  
    
    inFilename <- df_centroids$files[j]


  		




  169
  
  
    
    r1 <- raster(inFilename)


  		




  170
  
  
    
    r1[] <- df_centroids[[metric_name]][j] #improve this


  		




  171
  
  
    
    #set1f <- function(x){rep(NA, x)}


  		




  172
  
  
    
    #r_init <- init(r_in, fun=set1f)


  		




  173
  
  
    
    lf_tmp <- gsub(file_format,"",lf)


  		




  174
  
  
    
  


  		




  175
  
  
    
    extension_str <- extension(inFilename)


  		




  176
  
  
    
    raster_name_tmp <- gsub(extension_str,"",basename(inFilename))


  		




  177
  
  
    
    outFilename <- file.path(out_dir,paste(raster_name_tmp,"_",metric_name,"_",out_suffix,file_format,sep="")) #for use in function later...


  		




  178
  
  
    
  


  		




  179
  
  
    
    #raster_name <- file.path(out_dir_str,paste(raster_name_tmp,"_",out_suffix,file_format,sep=""))#output file


  		




  
  204
  
    
  #list_raster_name <- vector("list",length=length(lf))


  		




  
  205
  
    
  list_param_raster_df_centroids <- list(df_centroids,metric_name,NA_flag_val,file_format,out_dir,out_suffix)


  		




  
  206
  
    
  names(list_param_raster_df_centroids) <- c("df_centroids","metric_name","NA_flag_val","file_format","out_dir","out_suffix")


  		




  180
  207
  
    

  		




  181
  
  
    
    #raster_name_tmp <- paste("r_",metric_name,"_",interpolation_method,"_",date_processed,"_",out_suffix_str,sep="")


  		




  182
  
  
    
    #raster_name <- file.path(out_dir_str,raster_name_tmp)


  		




  183
  
  
    
    writeRaster(r1, NAflag=NA_flag_val,filename=outFilename,overwrite=TRUE)  


  		




  184
  
  
    
    list_raster_name[[j]] <- outFilename


  		




  185
  
  
    
  }


  		




  
  208
  
    
  #undebug(create_raster_df_centroids_fun)


  		




  
  209
  
    
  #test_lf <- lapply(1,FUN=create_raster_df_centroids_fun,list_param=list_param_raster_df_centroids)                           


  		




  186
  210
  
    
  


  		




  
  211
  
    
  list_raster_name <- mclapply(1:length(lf),FUN=create_raster_df_centroids_fun,list_param=list_param_raster_df_centroids,mc.preschedule=FALSE,mc.cores = num_cores)                           


  		




  
  212
  
    

  		




  187
  213
  
    
  raster_created_obj <- list(list_raster_name,df_centroids)


  		




  188
  214
  
    
  names(raster_created_obj) <- c("list_raster_name","df_centroids")


  		




  189
  215
  
    
  return(raster_created_obj)


  		




  190
  
  
    

  		




  191
  216
  
    
}


  		




  192
  217
  
    

  		




  193
  218
  
    
#### end of function


  		




  ......




  381
  406
  
    
  rast_list<-file.path(out_path,raster_name)


  		




  382
  407
  
    
  


  		




  383
  408
  
    
  ## The Raster and rgdal packages write temporary files on the disk when memory is an issue. This can potential build up


  		




  384
  
  
    
  ## in long  loops and can fill up hard drives resulting in errors. The following  sections removes these files 


  		




  385
  
  
    
  ## as they are created in the loop. This code section  can be transformed into a "clean-up function later on


  		




  386
  
  
    
  ## Start remove


  		




  387
  
  
    
  #tempfiles<-list.files(tempdir(),full.names=T) #GDAL transient files are not removed


  		




  388
  
  
    
  #files_to_remove<-grep(out_suffix,tempfiles,value=T) #list files to remove


  		




  389
  
  
    
  #if(length(files_to_remove)>0){


  		




  390
  
  
    
  #  file.remove(files_to_remove)


  		




  391
  
  
    
  #}


  		




  392
  
  
    
  #now remove temp files from raster package located in rasterTmpDir


  		




  
  409
  
    

  		




  393
  410
  
    
  removeTmpFiles(h=0) #did not work if h is not set to 0


  		




  394
  411
  
    
  ## end of remove section


  		




  395
  412
  
    
  


  		




  ......




  771
  788
  
    
                     paste("-A ", r_prod_sum_raster_name,sep=""),


  		




  772
  789
  
    
                     paste("-B ", r_weights_sum_raster_name,sep=""),


  		




  773
  790
  
    
                     paste("--outfile=",r_m_weighted_mean_raster_name,sep=""),


  		




  
  791
  
    
                     paste("NoDataValue=",NA_flag_val,sep=""),


  		




  774
  792
  
    
                     "--calc='A/B'","--overwrite",sep=" ") #division by zero is problematic...


  		




  775
  793
  
    
    system(cmd_str)


  		




  776
  
  
    
    


  		




  
  794
  
    
    cmd_mosaic_logfile <- file.path(out_dir,paste("cmd_mosaic_",out_suffix,".txt",sep=""))


  		




  
  795
  
    
    writeLines(cmd_str1,con=cmd_mosaic_logfile) #weights files to mosaic 


  		




  
  796
  
    
    #writeLines(cmd_str2,con=file.path(out_dir,paste("cmd_mosaic_",out_suffix,".txt",sep=""))) #weights files to mosaic 


  		




  
  797
  
    
    cat(cmd_str2, file=cmd_mosaic_logfile, append=TRUE, sep = "\n")


  		




  
  798
  
    

  		




  777
  799
  
    
    #cmd_str <- "/nobackupp6/aguzman4/climateLayers/sharedModules/bin/gdal_calc.py -A r_prod_weights_sum_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920101_reg4_run10_1500x4500_global_analyses_pred_1992_10052015.tif -B r_weights_sum_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920101_reg4_run10_1500x4500_global_analyses_pred_1992_10052015.tif --outfile='test2.tif' --calc='A/B' --overwrite"


  		




  778
  800
  
    
    


  		




  779
  801
  
    
    #writeRaster(r_m_weighted_mean, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE)  


  		




  ......




  993
  1015
  
    
}


  		




  994
  1016
  
    

  		




  995
  1017
  
    
#Use this instead of daily mosaic plot function


  		




  
  1018
  
    
#Add NAvalue flag!!


  		




  996
  1019
  
    
plot_screen_raster_val<-function(i,list_param){


  		




  997
  1020
  
    
  ##USAGE###


  		




  998
  1021
  
    
  #Screen raster list and produced plot as png.


  		












Also available in:  Unified diff