Revision 6db687b5
Added by Benoit Parmentier over 8 years ago
| 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/01/2016
|
|
| 8 |
#MODIFIED ON: 05/02/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 |
| ... | ... | |
| 101 | 101 |
y_var_month <- "TMin" |
| 102 | 102 |
} |
| 103 | 103 |
|
| 104 |
|
|
| 105 |
##Add for precip later... |
|
| 106 |
if (var == "TMAX") {
|
|
| 107 |
variable_name <- "maximum temperature" |
|
| 108 |
} |
|
| 109 |
if (var == "TMIN") {
|
|
| 110 |
variable_name <- "minimum temperature" |
|
| 111 |
} |
|
| 112 |
|
|
| 104 | 113 |
#interpolation_method<-c("gam_fusion") #other otpions to be added later
|
| 105 | 114 |
interpolation_method<-c("gam_CAI") #param 2
|
| 106 | 115 |
CRS_interp <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0" #param 3 |
| ... | ... | |
| 133 | 142 |
plotting_figures <- TRUE #running part2 of assessment to generate figures... # param 14 |
| 134 | 143 |
#num_cores <- args[8] #number of cores used # param 13, arg 8 |
| 135 | 144 |
num_cores <- 11 #number of cores used # param 13, arg 8 |
| 145 |
#python_bin <- "/nobackupp6/aguzman4/climateLayers/sharedModules2/bin" #PARAM 30 |
|
| 146 |
python_bin <- "/usr/bin" #PARAM 30 |
|
| 136 | 147 |
|
| 137 | 148 |
day_start <- "19990101" #PARAM 12 arg 12 |
| 138 | 149 |
day_end <- "19990103" #PARAM 13 arg 13 |
| ... | ... | |
| 160 | 171 |
|
| 161 | 172 |
########### #################### |
| 162 | 173 |
|
| 163 |
pred_temp_s <- raster("/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990101_reg4_1999_m_gam_CAI_dailyTmax_19990101_reg4_1999.tif")
|
|
| 174 |
#start_date <- day_to_mosaic_range[1] |
|
| 175 |
#end_date <- day_to_mosaic_range[2] |
|
| 176 |
start_date <- day_start #PARAM 12 arg 12 |
|
| 177 |
end_date <- day_end #PARAM 13 arg 13 |
|
| 178 |
|
|
| 179 |
date_to_plot <- seq(as.Date(strptime(start_date,"%Y%m%d")), as.Date(strptime(end_date,"%Y%m%d")), 'day') |
|
| 180 |
l_dates <- format(date_to_plot,"%Y%m%d") #format back to the relevant date format for files |
|
| 181 |
|
|
| 182 |
raster_name_lf <- c("/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990101_reg4_1999_m_gam_CAI_dailyTmax_19990101_reg4_1999.tif",
|
|
| 183 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990102_reg4_1999_m_gam_CAI_dailyTmax_19990102_reg4_1999.tif", |
|
| 184 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990103_reg4_1999_m_gam_CAI_dailyTmax_19990103_reg4_1999.tif") |
|
| 185 |
|
|
| 186 |
pred_temp_s <- raster(raster_name_lf[1]) |
|
| 187 |
|
|
| 188 |
lf_files <- c(raster_name_in) #match to mask |
|
| 189 |
rast_ref <- infile_mask |
|
| 190 |
##Maching resolution is probably only necessary for the r mosaic function |
|
| 191 |
#Modify later to take into account option R or python... |
|
| 192 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,python_bin,out_suffix,out_dir) |
|
| 193 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","python_bin","out_suffix","out_dir_str")
|
|
| 194 |
r_pred_matched <- raster_match(1,list_param_raster_match) |
|
| 195 |
raster_name_in <- c(r_pred_matched) |
|
| 164 | 196 |
|
| 165 | 197 |
if(mask_pred==TRUE){
|
| 166 | 198 |
|
| 167 | 199 |
r_mask <- raster(infile_mask) |
| 168 |
extension_str <- extension(filename(pred_temp_s )) |
|
| 169 |
raster_name_tmp <- |
|
| 170 |
gsub(extension_str,"",basename(filename(pred_temp_s ))) |
|
| 200 |
extension_str <- extension(raster_name_in) |
|
| 201 |
raster_name_tmp <- gsub(extension_str,"",basename(raster_name_in)) |
|
| 171 | 202 |
raster_name <- file.path(out_dir,paste(raster_name_tmp,"_masked.tif",sep = "")) |
| 172 |
r_pred <- mask(pred_temp_s,r_mask,filename = raster_name,overwrite = TRUE) |
|
| 173 |
|
|
| 203 |
r_pred <- mask(raster(raster_name_in),r_mask,filename = raster_name,overwrite = TRUE) |
|
| 174 | 204 |
} |
| 175 | 205 |
|
| 176 |
r_pred <- mask(pred_temp_s*1/scaling,r_mask,filename=)
|
|
| 177 |
#predictions<-mask(predictions,mask_rast)
|
|
| 178 |
|
|
| 206 |
extension_str <- extension(filename(r_pred))
|
|
| 207 |
raster_name_tmp <- gsub(extension_str,"",basename(filename(r_pred)))
|
|
| 208 |
raster_name <- file.path(out_dir,paste(raster_name_tmp,"_rescaled.tif",sep = "")) |
|
| 179 | 209 |
|
| 180 |
NAvalue(pred_temp_s) <- NA_flag_val |
|
| 181 |
pred_temp_s <- setMinMax(pred_temp_s) |
|
| 210 |
r_pred <- overlay(r_pred, fun=function(x){x*1/scaling},filename=raster_name)
|
|
| 211 |
NAvalue(r_pred) <- NA_flag_val |
|
| 212 |
r_pred <- setMinMax(r_pred) |
|
| 182 | 213 |
|
| 214 |
month_name <- month.name() |
|
| 215 |
l_dates <- as.Date(strptime(date_proc,"%Y%m%d")) |
|
| 183 | 216 |
|
| 184 | 217 |
#s.range <- c(min(minValue(pred_temp_s)), max(maxValue(pred_temp_s))) |
| 185 | 218 |
#s.range <- s.range+c(5,-5) |
| ... | ... | |
| 191 | 224 |
#max_val<- -10 |
| 192 | 225 |
min_val <- 0 |
| 193 | 226 |
|
| 194 |
layout_m<-c(1,3) #one row two columns |
|
| 227 |
#layout_m<-c(1,3) #one row two columns |
|
| 228 |
date_proc <- l_dates[i] |
|
| 195 | 229 |
|
| 196 | 230 |
#png(paste("Figure7a__spatial_pattern_tmax_prediction_levelplot_",date_selected,out_prefix,".png", sep=""),
|
| 197 | 231 |
# height=480*layout_m[1],width=480*layout_m[2]) |
| 198 |
plot(pred_temp_s,col=temp.colors(255),zlim=c(-5000,5000)) |
|
| 199 |
plot(pred_temp_s,col=heat.colors(255),zlim=c(-5000,5000)) |
|
| 200 |
plot(pred_temp_s,zlim=c(-5000,5000)) |
|
| 201 |
date_proc <- l_dates[i] |
|
| 232 |
#plot(r_pred,col=temp.colors(255),zlim=c(-3500,4500)) |
|
| 233 |
plot(r_pred,col=matlab.like(255),zlim=c(-35,45)) |
|
| 202 | 234 |
#paste(raster_name[1:7],collapse="_") |
| 203 | 235 |
#add filename option later |
| 204 | 236 |
|
| ... | ... | |
| 207 | 239 |
|
| 208 | 240 |
col_mfrow <- 1 |
| 209 | 241 |
row_mfrow <- 1 |
| 242 |
date_proc <- |
|
| 243 |
png_filename <- file.path(out_dir_str,paste("Figure4_clim_mosaics_day_","_",date_proc,"_",reg_name,"_",out_suffix,".png",sep =""))
|
|
| 244 |
|
|
| 245 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 246 |
|
|
| 247 |
plot(r_pred,main = paste("Predicted ",variable_name, " on ",date_proc , " ", ,sep = ""),cex.main =1.5)
|
|
| 210 | 248 |
|
| 211 |
png( |
|
| 212 |
filename = file.path( |
|
| 213 |
out_dir_str, |
|
| 214 |
paste( |
|
| 215 |
"Figure9_clim_mosaics_day_test","_",date_proc,"_",reg_name,"_",out_suffix,".png",sep = |
|
| 216 |
"" |
|
| 217 |
) |
|
| 218 |
), |
|
| 219 |
width = col_mfrow * res_pix,height = row_mfrow * res_pix |
|
| 220 |
) |
|
| 221 |
|
|
| 222 |
plot( |
|
| 223 |
r_pred,main = paste("Predicted on ",date_proc , " ", reg_name,sep = ""),cex.main =
|
|
| 224 |
1.5 |
|
| 225 |
) |
|
| 226 | 249 |
dev.off() |
| 227 | 250 |
|
| 228 | 251 |
#col.regions=temp.colors(25)) |
| 229 | 252 |
#dev.off() |
| 230 | 253 |
|
| 254 |
|
|
| 255 |
|
|
| 231 | 256 |
############################ END OF SCRIPT ################################## |
Also available in: Unified diff
adding function to screen mosaic using mask and assigning NA as well as scaling