Project

General

Profile

« Previous | Next » 

Revision b765ab2b

Added by Benoit Parmentier almost 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