Project

General

Profile

« Previous | Next » 

Revision 0b599be0

Added by Benoit Parmentier over 8 years ago

moving function to script and fixing temporal profiles extracted

View differences:

climate/research/oregon/interpolation/NASA2016_conference_temperature_predictions.R
5 5
#Figures and data for the AAG conference are also produced.
6 6
#AUTHOR: Benoit Parmentier 
7 7
#CREATED ON: 05/01/2016  
8
#MODIFIED ON: 05/02/2016            
8
#MODIFIED ON: 05/03/2016            
9 9
#Version: 1
10 10
#PROJECT: Environmental Layers project     
11 11
#COMMENTS: Initial commit, script based on part 2 of assessment, will be modified further for overall assessment 
......
51 51
library(colorRamps)
52 52
library(zoo)
53 53
library(xts)
54
library(lubridate)
54 55

  
55 56
###### Function used in the script #######
56 57
  
57 58
#script_path <- "/nobackupp8/bparmen1/env_layers_scripts" #path to script
58 59
script_path <- "/home/parmentier/Data/IPLANT_project/env_layers_scripts" #path to script
59 60

  
60
#Mosaic related
61
## NASA poster and paper related
62
source(file.path(script_path,"NASA2016_conference_temperature_predictions_function_05032016b.R"))
63

  
64
#Mosaic related on NEX
61 65
#script_path <- "/home/parmentier/Data/IPLANT_project/env_layers_scripts"
62 66
function_mosaicing_functions <- "global_run_scalingup_mosaicing_function_04232016.R" #PARAM12
63 67
function_mosaicing <-"global_run_scalingup_mosaicing_05012016.R"
64 68
source(file.path(script_path,function_mosaicing)) #source all functions used in this script 
65 69
source(file.path(script_path,function_mosaicing_functions)) #source all functions used in this script 
66 70

  
67
#Assessment
71
#Assessment on NEX
68 72
function_assessment_part1_functions <- "global_run_scalingup_assessment_part1_functions_02112015.R" #PARAM12
69 73
function_assessment_part1a <-"global_run_scalingup_assessment_part1a_01042016.R"
70 74
function_assessment_part2 <- "global_run_scalingup_assessment_part2_02092016.R"
......
76 80
source(file.path(script_path,function_assessment_part2_functions)) #source all functions used in this script 
77 81
source(file.path(script_path,function_assessment_part3)) #source all functions used in this script 
78 82

  
79
pre_process_raster_mosaic_fun <- function(i,list_param){
80
  
81
  
82
  ## Extract parameters
83
  
84
  lf <- list_param$lf
85
  python_bin <- list_param$python_bin
86
  infile_mask <- list_param$infile_mask
87
  scaling <- list_param$scaling
88
  mask_pred <- list_param$mask_pred
89
  NA_flag_val <- list_param$NA_flag_val
90
  out_suffix <- list_param$out_suffix
91
  out_dir <- list_param$out_dir
92
  
93
  raster_name_in <- lf[i]
94
  
95
  #Step 1: match extent and resolution
96
  
97
  lf_files <- c(raster_name_in) #match to mask
98
  rast_ref <- infile_mask
99
  ##Maching resolution is probably only necessary for the r mosaic function
100
  #Modify later to take into account option R or python...
101
  list_param_raster_match <- list(lf_files,rast_ref,file_format,python_bin,out_suffix,out_dir)
102
  names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","python_bin","out_suffix","out_dir_str")
103
  r_pred_matched <- raster_match(1,list_param_raster_match)
104
  raster_name_in <- c(r_pred_matched)
105

  
106
  #Step 2: mask
107
  if(mask_pred==TRUE){
108
    r_mask <- raster(infile_mask)
109
    extension_str <- extension(raster_name_in)
110
    raster_name_tmp <- gsub(extension_str,"",basename(raster_name_in))
111
    raster_name <- file.path(out_dir,paste(raster_name_tmp,"_masked.tif",sep = ""))
112
    r_pred <- mask(raster(raster_name_in),r_mask,filename = raster_name,overwrite = TRUE)
113
  }
114
  
115
  NAvalue(r_pred) <- NA_flag_val
116
  #r_pred <- setMinMax(r_pred)
117

  
118
  #Step 3: remove scaling factor
119
  raster_name_in <- filename(r_pred)
120
  extension_str <- extension(raster_name_in)
121
  raster_name_tmp <- gsub(extension_str,"",basename(filename(r_pred)))
122
  raster_name_out <- file.path(out_dir,paste(raster_name_tmp,"_rescaled.tif",sep = ""))
123
  #r_pred <- overlay(r_pred, fun=function(x){x*1/scaling},filename=raster_name,overwrite=T)
124

  
125
  #raster_name_in <- "comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990102_reg4_1999_m_gam_CAI_dailyTmax_19990102_reg4_1999_m__meeting_NASA_reg4_04292016_masked.tif"
126
  python_cmd <- file.path(python_bin,"gdal_calc.py")
127
  cmd_str3 <- paste(python_cmd, 
128
                     paste("-A ", raster_name_in,sep=""),
129
                     paste("--outfile=",raster_name_out,sep=""),
130
                     #paste("--type=","Int32",sep=""),
131
                     "--co='COMPRESS=LZW'",
132
                     paste("--NoDataValue=",NA_flag_val,sep=""),
133
                     paste("--calc='(A)*",scaling,"'",sep=""),
134
                     "--overwrite",sep=" ") #division by zero is problematic...
135
  #system(cmd_str3)
136
  system(cmd_str3)
137
  #NAvalue(r_pred) <- NA_flag_val
138
  #r_pred <- 
139
  r_pred <- setMinMax(raster(raster_name_out))
140
  
141
  return(raster_name_out)
142
}
143

  
144
plot_raster_mosaic <- function(i,list_param){
145
  #Function to plot mosaic for poster
146
  #
147
  l_dates <- list_param$l_dates
148
  r_mosaiced_scaled <- list_param$r_mosaiced_scaled
149
  NA_flag_val <- list_param$NA_flag_val
150
  out_dir <- list_param$out_dir
151
  out_suffix <- list_param$out_suffix
152
  region_name <- list_param$region_name
153
  variable_name <- list_param$variable_name
154

  
155
#for (i in 1:length(nlayers(r_mosaic_scaled))){
156
  
157
  date_proc <- l_dates[i]
158
  r_pred <- subset(r_mosaic_scaled,i)
159
  NAvalue(r_pred)<- NA_flag_val 
160
 
161
  date_proc <- l_dates[i]
162
  date_val <- as.Date(strptime(date_proc,"%Y%m%d"))
163
  #month_name <- month.name(date_val)
164
  month_str <- format(date_val, "%b") ## Month, char, abbreviated
165
  year_str <- format(date_val, "%Y") ## Year with century
166
  day_str <- as.numeric(format(date_val, "%d")) ## numeric month
167
  date_str <- paste(month_str," ",day_str,", ",year_str,sep="")
168
  
169
  res_pix <- 1200
170
  #res_pix <- 480
171
  col_mfrow <- 1
172
  row_mfrow <- 1
173
  
174
  png_filename <-  file.path(out_dir,paste("Figure4_clim_mosaics_day_","_",date_proc,"_",region_name,"_",out_suffix,".png",sep =""))
175
  title_str <-  paste("Predicted ",variable_name, " on ",date_str , " ", sep = "")
176
  
177
  png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix)
178
  plot(r_pred,main =title_str,cex.main =1.5,col=matlab.like(255),zlim=c(-50,50),
179
       legend.shrink=0.8,legend.width=0.8)
180
       #axis.args = list(cex.axis = 1.6), #control size of legend z
181
       #legend.args=list(text='dNBR', side=4, line=2.5, cex=2.2))
182
       #legend.args=list(text='dNBR', side=4, line=2.49, cex=1.6))
183
  dev.off()
184

  
185
  return(png_filename)
186
}
187

  
188
extract_date <- function(i,x,item_no=NULL){
189
  y <- unlist(strsplit(x[[i]],"_"))
190
  if(is.null(item_no)){
191
    date_str <- y[length(y)-2] #count from end
192
  }else{
193
    date_str <- y[item_no]
194
  }
195
  return(date_str)
196
}
197

  
198 83
###############################
199 84
####### Parameters, constants and arguments ###
200 85

  
......
252 137
#python_bin <- "/nobackupp6/aguzman4/climateLayers/sharedModules2/bin" #PARAM 30
253 138
python_bin <- "/usr/bin" #PARAM 30
254 139

  
255
day_start <- "19840101" #PARAM 12 arg 12
256
day_end <- "20021231" #PARAM 13 arg 13
140
day_start <- "1986101" #PARAM 12 arg 12
141
day_end <- "19981231" #PARAM 13 arg 13
257 142

  
258 143
#infile_mask <- "/nobackupp8/bparmen1/NEX_data/regions_input_files/r_mask_LST_reg4.tif"
259 144
infile_mask <- "/data/project/layers/commons/NEX_data/regions_input_files/r_mask_LST_reg4.tif"
......
275 160
l_dates <- c("19920101","19920102","19920103","19920701","19920702","19990703")
276 161

  
277 162
df_points_extracted_fname <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/data_points_extracted.txt"
163
NA_flag_val_mosaic <- -3399999901438340239948148078125514752.000
278 164

  
279 165
##################### START SCRIPT #################
280 166

  
......
326 212
#paste(raster_name[1:7],collapse="_")
327 213
#add filename option later
328 214

  
329
NA_flag_val_mosaic <- -3399999901438340239948148078125514752.000
215
#NA_flag_val_mosaic <- -3399999901438340239948148078125514752.000
330 216

  
331 217
list_param_plot_raster_mosaic <- list(l_dates,r_mosaic_scaled,NA_flag_val_mosaic,out_dir,out_suffix,
332 218
                                      region_name,variable_name)
......
335 221

  
336 222
lf_mosaic_plot_fig <- mclapply(1:length(lf_mosaic_scaled),FUN=plot_raster_mosaic,list_param=list_param_plot_raster_mosaic,mc.preschedule=FALSE,mc.cores = num_cores)                         
337 223

  
338

  
339

  
340 224
############### PART2: temporal profile #############
341 225
#### Extract time series
342 226
###
......
379 263

  
380 264
station_id <- 8
381 265
var_name <-paste0("ID_",station_id)
382
#aggregate(sdf_tmp
266

  
267
##Screen for unique date values
383 268
if(max(unique_date_tb)>1){
384 269
#  formula_str <- paste(var_name," ~ ","TRIP_START_DATE_f",sep="")
385 270
   var_pix <- aggregate(ID_8 ~ date, data = df_points, mean) #aggregate by date
386 271
}
387 272

  
388
 
389
#start_date <-input$dates[1]
390
#end_date <-input$dates[2]
273
var_pix$ID_8 <- var_pix$ID_8*scaling
391 274

  
392
d_z_tmp <- zoo(df_points[,station_id],df_points$date)
275
d_z_tmp <- zoo(var_pix$ID_8,var_pix$date)
276
names(d_z_tmp)<-"ID_8"
277
min(d_z_tmp$ID_8)
278
max(d_z_tmp$ID_8)
393 279

  
280
plot(d_z_tmp)
281

  
282
day_start <- "1986-01-01" #PARAM 12 arg 12
283
day_end <- "1998-12-31" #PARAM 13 arg 13
284

  
285
start_date <- as.Date(day_start)
286
end_date <- as.Date(day_end)
287
start_year <- year(start_date)
288
end_year <- year(end_date)
289

  
290
d_z <- window(d_z_tmp,start=start_date,end=end_date)
291
#d_z2 <- window(d_z_tmp2,start=start_date,end=end_date)
292

  
293
title_str <- paste("Predicted daily ",variable_name," for the ", start_year,"-",end_year," time period",sep="")
294
plot(d_z,ylab="tmax in deg C",xlab="daily time steps",
295
     main=title_str,
296
     lty=3)
394 297

  
395
d_z <- window(d_z_tmp,start=start_date,end=end_date)   
396 298

  
397
  
398 299
#data_pixel <- data_df[id_selected,]
399 300
#data_pixel$rainfall <- as.numeric(data_pixel$rainfall)
400 301
#d_z_tmp <-zoo(data_pixel$rainfall,as.Date(data_pixel$date))
......
402 303
#data_pixel <- as.data.frame(data_pixel)
403 304
#d_z_tmp2 <- zoo(data_pixel[[var_name]],as.Date(data_pixel$date))
404 305
    
405
#start_date <-input$dates[1]
406
#end_date <-input$dates[2]
407
    
408
#d_z <- window(d_z_tmp,start=start_date,end=end_date)
409
#d_z2 <- window(d_z_tmp2,start=start_date,end=end_date)
410
#df2 <- as.data.frame(d_z2)
411
#names(df2)<- var_name
412

  
413 306
#df_tmp <- subset(data_var,data_var$ID_stat==id_name)
414 307
#if(da)
415 308

  

Also available in: Unified diff