Revision bd534a95
Added by Benoit Parmentier about 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part0.R | ||
|---|---|---|
| 1 |
#################################### INTERPOLATION OF TEMPERATURES #######################################
|
|
| 2 |
####################### Assessment of product part 1: mosaic and accuracy ##############################
|
|
| 1 |
############################## INTERPOLATION OF TEMPERATURES ####################################### |
|
| 2 |
####################### Script for assessment of scaling up on NEX : part 0 ##############################
|
|
| 3 | 3 |
#This script uses the worklfow code applied to the globe. Results currently reside on NEX/PLEIADES NASA. |
| 4 |
#This part 2 of the assessment focuses on graphics to explore the spatial patterns of raster times series as figures and movie |
|
| 4 |
#This script checks the number of predictions by tiles and years. |
|
| 5 |
#with the goal of predicting potential gaps or missing predictions in fugure mosaics by region. |
|
| 6 |
#The general logic is to check the number of overlap by shapefile polyon tiles |
|
| 7 |
#along with the predicitons for every day of the year (*.tif) |
|
| 8 |
#Summary tables and data are also produced in the script. |
|
| 9 |
# |
|
| 5 | 10 |
#AUTHOR: Benoit Parmentier |
| 6 |
#CREATED ON: 10/03/2016
|
|
| 7 |
#MODIFIED ON: 10/27/2016
|
|
| 11 |
#CREATED ON: 10/27/2016
|
|
| 12 |
#MODIFIED ON: 10/30/2016
|
|
| 8 | 13 |
#Version: 1 |
| 9 | 14 |
#PROJECT: Environmental Layers project |
| 10 |
#COMMENTS: Initial commit, script based on part NASA biodiversity conferenc
|
|
| 15 |
#COMMENTS: Major update of code with changes to the listing of tiles files
|
|
| 11 | 16 |
#TODO: |
| 12 |
#1) Add plot broken down by year and region |
|
| 13 |
#2) Modify code for overall assessment accross all regions and year |
|
| 14 |
#3) Clean up |
|
| 15 |
|
|
| 17 |
#1) |
|
| 16 | 18 |
#First source these files: |
| 17 | 19 |
#Resolved call issues from R. |
| 18 | 20 |
#source /nobackupp6/aguzman4/climateLayers/sharedModules2/etc/environ.sh |
| ... | ... | |
| 107 | 109 |
metric_name <- "var_pred" #use RMSE if accuracy |
| 108 | 110 |
|
| 109 | 111 |
item_no <- 13 |
| 112 |
day_start <- "2000101" #PARAM 12 arg 12 |
|
| 113 |
day_end <- "20001231" #PARAM 13 arg 13 |
|
| 114 |
#date_start <- day_start |
|
| 115 |
#date_end <- day_end |
|
| 110 | 116 |
date_start <- day_start |
| 111 | 117 |
date_end <- day_end |
| 112 | 118 |
#day_start <- "1984101" #PARAM 12 arg 12 |
| 113 | 119 |
#day_end <- "20141231" #PARAM 13 arg 13 |
| 120 |
day_to_mosaic_range <- NULL |
|
| 114 | 121 |
|
| 115 |
############################################ |
|
| 116 |
#### Parameters and constants |
|
| 117 |
|
|
| 118 |
##Add for precip later... |
|
| 119 |
if (var == "TMAX") {
|
|
| 120 |
y_var_name <- "dailyTmax" |
|
| 121 |
y_var_month <- "TMax" |
|
| 122 |
} |
|
| 123 |
if (var == "TMIN") {
|
|
| 124 |
y_var_name <- "dailyTmin" |
|
| 125 |
y_var_month <- "TMin" |
|
| 126 |
} |
|
| 127 |
|
|
| 128 |
##Add for precip later... |
|
| 129 |
if (var == "TMAX") {
|
|
| 130 |
variable_name <- "maximum temperature" |
|
| 131 |
} |
|
| 132 |
if (var == "TMIN") {
|
|
| 133 |
variable_name <- "minimum temperature" |
|
| 134 |
} |
|
| 135 |
|
|
| 136 |
#on ATLAS |
|
| 137 |
#in_dir1 <- "/data/project/layers/commons/NEX_data/test_run1_03232014/output" #On Atlas |
|
| 138 |
#parent output dir : contains subset of the data produced on NEX |
|
| 139 |
#in_dir1 <- "/data/project/layers/commons/NEX_data/output_run6_global_analyses_09162014/output20Deg2" |
|
| 140 |
# parent output dir for the curent script analyes |
|
| 141 |
#out_dir <-"/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/" #On NCEAS Atlas |
|
| 142 |
# input dir containing shapefiles defining tiles |
|
| 143 |
#in_dir_shp <- "/data/project/layers/commons/NEX_data/output_run5_global_analyses_08252014/output/subset/shapefiles" |
|
| 144 | 122 |
in_dir <- "/nobackupp6/aguzman4/climateLayers/out/reg6/assessment" |
| 145 |
in_dir <- "/nobackupp8/bparmen1/climateLayers/out/reg6/assessment" |
|
| 146 |
|
|
| 123 |
#in_dir <- "/nobackupp8/bparmen1/climateLayers/out/reg6/assessment" |
|
| 147 | 124 |
#in_dir_mosaic <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg6/mosaics/mosaic" #predicted mosaic |
| 148 |
#in_dir_mosaic <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg1/mosaics/mosaic" |
|
| 149 |
#in_dir_mosaic <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg5/mosaics/mosaic" |
|
| 150 |
#in_dir_mosaic <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/mosaic" #note dropped the s in mosaics |
|
| 151 | 125 |
|
| 152 | 126 |
region_name <- c("reg6") #param 6, arg 3
|
| 153 | 127 |
out_suffix <- "global_assessment_reg6_10232016" |
| 154 | 128 |
|
| 155 | 129 |
create_out_dir_param <- TRUE #param 9, arg 6 |
| 156 | 130 |
|
| 157 |
|
|
| 158 | 131 |
out_dir <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg6/assessment" |
| 159 | 132 |
|
| 160 | 133 |
#run_figure_by_year <- TRUE # param 10, arg 7 |
| ... | ... | |
| 168 | 141 |
#python_bin <- "/nobackupp6/aguzman4/climateLayers/sharedModules2/bin" #PARAM 30 |
| 169 | 142 |
python_bin <- "/usr/bin" #PARAM 30 |
| 170 | 143 |
|
| 171 |
day_start <- "2000101" #PARAM 12 arg 12 |
|
| 172 |
day_end <- "20001231" #PARAM 13 arg 13 |
|
| 173 |
#date_start <- day_start |
|
| 174 |
#date_end <- day_end |
|
| 175 |
|
|
| 176 | 144 |
#NA_flag_val_mosaic <- -3399999901438340239948148078125514752.000 |
| 177 | 145 |
NA_flag_val_mosaic <- -32768 |
| 178 | 146 |
in_dir_list_filename <- NULL #if NULL, use the in_dir directory to search for info |
| ... | ... | |
| 209 | 177 |
#/nobackupp6/aguzman4/climateLayers/out/reg6/subset/shapefiles |
| 210 | 178 |
list_year_predicted <- c(2000,2012,2013) #year still on disk for reg6 |
| 211 | 179 |
|
| 180 |
|
|
| 181 |
############################################ |
|
| 182 |
#### Parameters and constants |
|
| 183 |
|
|
| 184 |
##Add for precip later... |
|
| 185 |
if (var == "TMAX") {
|
|
| 186 |
y_var_name <- "dailyTmax" |
|
| 187 |
y_var_month <- "TMax" |
|
| 188 |
} |
|
| 189 |
if (var == "TMIN") {
|
|
| 190 |
y_var_name <- "dailyTmin" |
|
| 191 |
y_var_month <- "TMin" |
|
| 192 |
} |
|
| 193 |
|
|
| 194 |
##Add for precip later... |
|
| 195 |
if (var == "TMAX") {
|
|
| 196 |
variable_name <- "maximum temperature" |
|
| 197 |
} |
|
| 198 |
if (var == "TMIN") {
|
|
| 199 |
variable_name <- "minimum temperature" |
|
| 200 |
} |
|
| 201 |
|
|
| 202 |
i<-1 |
|
| 203 |
|
|
| 204 |
predictions_tiles_missing_fun <- function(list_param,i){
|
|
| 205 |
|
|
| 212 | 206 |
#from script: |
| 213 | 207 |
#interpolation/global_run_scalingup_assessment_part1a.R |
| 214 | 208 |
|
| ... | ... | |
| 246 | 240 |
######################## PART0: Read content of predictions first.... ##### |
| 247 | 241 |
#function looped over i, correspoding to year predicted |
| 248 | 242 |
|
| 249 |
list_outfiles <- vector("list", length=35) #collect names of output files, this should be dynamic?
|
|
| 250 |
list_outfiles_names <- vector("list", length=35) #collect names of output files
|
|
| 243 |
#list_outfiles <- vector("list", length=35) #collect names of output files, this should be dynamic?
|
|
| 244 |
#list_outfiles_names <- vector("list", length=35) #collect names of output files
|
|
| 251 | 245 |
|
| 252 |
|
|
| 253 |
#### STart here |
|
| 254 |
###### This is from assessment 1 |
|
| 255 |
|
|
| 256 |
#for each polygon find you overlap!! |
|
| 257 |
#plot number of overlap |
|
| 258 |
#for specific each find prediction... |
|
| 259 | 246 |
year_predicted <- list_param_run_assessment_prediction$list_year_predicted[i] |
| 260 | 247 |
|
| 261 | 248 |
in_dir1_reg <- file.path(in_dir1,region_name) |
| ... | ... | |
| 295 | 282 |
|
| 296 | 283 |
setwd(out_dir) |
| 297 | 284 |
|
| 298 |
|
|
| 299 | 285 |
##raster_prediction object : contains testing and training stations with RMSE and model object |
| 300 | 286 |
in_dir_list_tmp <- file.path(in_dir_list,year_predicted) |
| 301 | 287 |
list_raster_obj_files <- try(lapply(in_dir_list_tmp,FUN=function(x){list.files(path=x,pattern="^raster_prediction_obj.*.RData",full.names=T)}))
|
| ... | ... | |
| 358 | 344 |
table(df_time_series$missing) |
| 359 | 345 |
table(df_time_series$year) |
| 360 | 346 |
|
| 347 |
######################## |
|
| 348 |
#### Step 2: examine overlap |
|
| 361 | 349 |
|
| 362 | 350 |
#combine polygon |
| 363 | 351 |
#http://gis.stackexchange.com/questions/155328/merging-multiple-spatialpolygondataframes-into-1-spdf-in-r |
| ... | ... | |
| 371 | 359 |
#for each date of year report data in table. |
| 372 | 360 |
|
| 373 | 361 |
#go through table and hsow if there are missing data (no prediction) or report min predictions for tile set? |
| 374 |
|
|
| 362 |
|
|
| 363 |
#for each polygon find you overlap!! |
|
| 364 |
#plot number of overlap |
|
| 365 |
#for specific each find prediction... |
|
| 366 |
|
|
| 367 |
######################## |
|
| 368 |
#### Step 3: combine overlap information and number of predictions by day |
|
| 375 | 369 |
|
| 376 | 370 |
|
| 371 |
|
|
| 372 |
|
|
| 373 |
return() |
|
| 374 |
} |
|
| 375 |
|
|
| 377 | 376 |
############################ END OF SCRIPT ################################## |
Also available in: Unified diff
adding function predictions_tiles_missing_fun