Revision 95bb7a22
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part1.R | ||
|---|---|---|
| 4 | 4 |
#Combining tables and figures for individual runs for years and tiles. |
| 5 | 5 |
#AUTHOR: Benoit Parmentier |
| 6 | 6 |
#CREATED ON: 05/15/2016 |
| 7 |
#MODIFIED ON: 05/24/2016
|
|
| 7 |
#MODIFIED ON: 07/30/2016
|
|
| 8 | 8 |
#Version: 1 |
| 9 | 9 |
#PROJECT: Environmental Layers project |
| 10 | 10 |
#COMMENTS: Initial commit, script based on part NASA biodiversity conferenc |
| ... | ... | |
| 118 | 118 |
|
| 119 | 119 |
#reg1 (North Am), reg2(Europe),reg3(Asia), reg4 (South Am), reg5 (Africa), reg6 (Australia-Asia) |
| 120 | 120 |
#master directory containing the definition of tile size and tiles predicted |
| 121 |
in_dir <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/assessment"
|
|
| 122 |
in_dir_mosaic <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/mosaic"
|
|
| 121 |
in_dir <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg5/assessment"
|
|
| 122 |
in_dir_mosaic <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg5/mosaic/mosaic"
|
|
| 123 | 123 |
|
| 124 |
region_name <- c("reg4") #param 6, arg 3
|
|
| 124 |
region_name <- c("reg5") #param 6, arg 3
|
|
| 125 |
out_suffix <- "_global_assessment_reg5_07292016" |
|
| 125 | 126 |
|
| 126 | 127 |
create_out_dir_param <- TRUE #param 9, arg 6 |
| 127 |
out_suffix <- "_global_assessment_reg4_05152016" |
|
| 128 | 128 |
|
| 129 |
out_dir <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/assessment" |
|
| 129 |
|
|
| 130 |
out_dir <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg5/assessment" |
|
| 130 | 131 |
|
| 131 | 132 |
#run_figure_by_year <- TRUE # param 10, arg 7 |
| 132 | 133 |
|
| ... | ... | |
| 143 | 144 |
day_end <- "19981231" #PARAM 13 arg 13 |
| 144 | 145 |
|
| 145 | 146 |
#infile_mask <- "/nobackupp8/bparmen1/NEX_data/regions_input_files/r_mask_LST_reg4.tif" |
| 146 |
infile_mask <- "/data/project/layers/commons/NEX_data/regions_input_files/r_mask_LST_reg4.tif"
|
|
| 147 |
infile_mask <- "/data/project/layers/commons/NEX_data/regions_input_files/r_mask_LST_reg5.tif"
|
|
| 147 | 148 |
|
| 148 | 149 |
#run_figure_by_year <- TRUE # param 10, arg 7 |
| 149 | 150 |
list_year_predicted <- "1984,2014" |
| 150 | 151 |
scaling <- 0.01 #was scaled on 100 |
| 151 | 152 |
#if scaling is null then perform no scaling!! |
| 152 | 153 |
|
| 153 |
df_centroids_fname <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/output_reg4_1999/df_centroids_19990701_reg4_1999.txt"
|
|
| 154 |
df_centroids_fname <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/output_reg5_1999/df_centroids_19990701_reg5_1999.txt"
|
|
| 154 | 155 |
|
| 155 |
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_19920101_reg4_1992_m_gam_CAI_dailyTmax_19920101_reg4_1992.tif",
|
|
| 156 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920102_reg4_1992_m_gam_CAI_dailyTmax_19920102_reg4_1992.tif",
|
|
| 157 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920103_reg4_1992_m_gam_CAI_dailyTmax_19920103_reg4_1992.tif",
|
|
| 158 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920701_reg4_1992_m_gam_CAI_dailyTmax_19920701_reg4_1992.tif",
|
|
| 159 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920702_reg4_1992_m_gam_CAI_dailyTmax_19920702_reg4_1992.tif",
|
|
| 160 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920703_reg4_1992_m_gam_CAI_dailyTmax_19920703_reg4_1992.tif")
|
|
| 156 |
raster_name_lf <- c("/data/project/layers/commons/NEX_data/climateLayers/out/reg5/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920101_reg4_1992_m_gam_CAI_dailyTmax_19920101_reg4_1992.tif",
|
|
| 157 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg5/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920102_reg4_1992_m_gam_CAI_dailyTmax_19920102_reg4_1992.tif",
|
|
| 158 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg5/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920103_reg4_1992_m_gam_CAI_dailyTmax_19920103_reg4_1992.tif",
|
|
| 159 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg5/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920701_reg4_1992_m_gam_CAI_dailyTmax_19920701_reg4_1992.tif",
|
|
| 160 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg5/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920702_reg4_1992_m_gam_CAI_dailyTmax_19920702_reg4_1992.tif",
|
|
| 161 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg5/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920703_reg4_1992_m_gam_CAI_dailyTmax_19920703_reg4_1992.tif")
|
|
| 161 | 162 |
|
| 162 | 163 |
#l_dates <- c("19990101","19990102","19990103","19990701","19990702","19990703")
|
| 163 |
l_dates <- c("19920101","19920102","19920103","19920701","19920702","19920703") #dates to plot and analze
|
|
| 164 |
l_dates <- c("19910110","19910111","19910112","19910113","19910114") #dates to plot and analze
|
|
| 164 | 165 |
|
| 165 |
df_points_extracted_fname <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/data_points_extracted.txt"
|
|
| 166 |
df_points_extracted_fname <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg5/mosaic/int_mosaics/data_points_extracted.txt"
|
|
| 166 | 167 |
df_points_extracted_fname <- NULL #if null compute on the fly |
| 167 | 168 |
#r_mosaic_fname <- "r_mosaic.RData" |
| 168 | 169 |
r_mosaic_fname <- NULL #if null create a stack from input dir |
| ... | ... | |
| 170 | 171 |
#NA_flag_val_mosaic <- -3399999901438340239948148078125514752.000 |
| 171 | 172 |
NA_flag_val_mosaic <- -32768 |
| 172 | 173 |
in_dir_list_filename <- NULL #if NULL, use the in_dir directory to search for info |
| 174 |
countries_shp <-"/data/project/layers/commons/NEX_data/countries.shp" #Atlas |
|
| 175 |
|
|
| 173 | 176 |
|
| 174 | 177 |
##################### START SCRIPT ################# |
| 175 | 178 |
|
| ... | ... | |
| 210 | 213 |
|
| 211 | 214 |
### Get station data information: create a data.frame with stations info |
| 212 | 215 |
data_date <- unlist(lapply(list_data_s_fname, function(x){unlist(strsplit(basename(x),"_"))[5]}))
|
| 213 |
df_data_s_fname <- data.frame(data_date,region_name,dataset="data_s",file=list_data_s_fname) |
|
| 216 |
#df_data_s_fname <- data.frame(data_date,region_name,dataset="data_s",file=list_data_s_fname) |
|
| 217 |
year_str <- as.character(unlist(lapply(list_data_s_fname, function(x){unlist(strsplit(basename(x),"_"))[5]})))
|
|
| 218 |
df_data_s_fname <- data.frame(year_str,region_name,dataset="data_s",file=list_data_s_fname) |
|
| 214 | 219 |
|
| 215 | 220 |
year_str <- as.character(unlist(lapply(list_data_v_fname, function(x){unlist(strsplit(basename(x),"_"))[5]})))
|
| 216 | 221 |
df_data_v_fname <- data.frame(year_str,region_name,dataset="data_v",file=list_data_s_fname) |
| ... | ... | |
| 219 | 224 |
df_data_v_fname$year <- as.character(df_data_v_fname$year_str) |
| 220 | 225 |
df_data_v_fname$file <- as.character(df_data_v_fname$file) |
| 221 | 226 |
df_data_s_fname$file <- as.character(df_data_s_fname$file) |
| 227 |
df_data_s_fname$year <- as.character(df_data_s_fname$year_str) |
|
| 222 | 228 |
|
| 223 | 229 |
############### PART1: Select stations for specific list of dates ############# |
| 224 |
#### Extract corresponding stationsfor given dates and plot stations used |
|
| 230 |
#### Extract corresponding stationsfor given dates and/or plot stations used
|
|
| 225 | 231 |
## Use station from specific year and date? |
| 226 | 232 |
|
| 227 | 233 |
list_dates_val <- as.Date(strptime(l_dates,"%Y%m%d")) |
| ... | ... | |
| 246 | 252 |
list_df_s_stations[[i]] <- read.table(filename_tmp,stringsAsFactors=F,sep=",") |
| 247 | 253 |
} |
| 248 | 254 |
|
| 255 |
df_combined <- do.call(rbind,list_df_v_stations) |
|
| 256 |
|
|
| 257 |
#### Get df points for specific dates |
|
| 258 |
#lapply(1:length(l_dates)list_df_v_stations,function(x){x$date==l_dates})
|
|
| 259 |
#dim(list_df_v_stations[[1]]) |
|
| 260 |
list_df_points_dates <- vector("list",length=length(l_dates))
|
|
| 261 |
for(i in 1:length(l_dates)){
|
|
| 262 |
#year_str <- list_year_str[[i]] |
|
| 263 |
list_df_points_dates[[i]] <- subset(df_combined,df_combined$date==l_dates[i]) |
|
| 264 |
} |
|
| 265 |
|
|
| 266 |
df_combined <- do.call(rbind,list_df_points_dates) |
|
| 267 |
|
|
| 249 | 268 |
#function(x){x$date==}
|
| 250 |
df_points <- subset(df_stations,df_stations_tmp$date==l_dates[1]) |
|
| 251 |
table(df_points$tile_id) |
|
| 252 |
plot(df_points,add=T) |
|
| 253 |
coordinates(df_points) <- cbind(df_points$x,df_points$y) |
|
| 254 |
proj4string(df_points) <- CRS_locs_WGS84 |
|
| 255 |
## No spatial duplicates |
|
| 256 |
df_points_day <- remove.duplicates(df_points) #remove duplicates... |
|
| 269 |
#df_points <- subset(df_stations,df_stations_tmp$date==l_dates[1]) |
|
| 270 |
reg_layer <- readOGR(dsn=dirname(countries_shp),sub(".shp","",basename(countries_shp)))
|
|
| 271 |
|
|
| 272 |
## Now plot by dates: |
|
| 273 |
|
|
| 274 |
model_name <- "res_mod1" |
|
| 275 |
var_selected <- "res_mod1" |
|
| 276 |
|
|
| 277 |
for(i in 1:length(l_dates)){
|
|
| 278 |
#d |
|
| 279 |
date_processed <- l_dates[i] |
|
| 280 |
|
|
| 281 |
df_points <- list_df_points_dates[[i]] |
|
| 282 |
#plot(df_points) |
|
| 283 |
table(df_points$tile_id) |
|
| 284 |
#plot(df_points,add=T) |
|
| 285 |
coordinates(df_points) <- cbind(df_points$x,df_points$y) |
|
| 286 |
proj4string(df_points) <- CRS_locs_WGS84 |
|
| 287 |
# # No spatial duplicates |
|
| 288 |
df_points_day <- remove.duplicates(df_points) #remove duplicates... |
|
| 289 |
plot(df_points_day) |
|
| 290 |
dim(df_points_day) |
|
| 291 |
dim(df_points) |
|
| 292 |
#plot(df_points) |
|
| 293 |
|
|
| 294 |
#data_stations_temp <- aggregate(id ~ date, data = data_stations, min) |
|
| 295 |
#data_stations_temp <- aggregate(id ~ x + y + date + dailyTmax + mod1 + res_mod1 , data = data_stations, min) |
|
| 296 |
data_stations_temp <- aggregate(id ~ x + y + date + dailyTmax + res_mod1,data = df_points, mean ) #+ mod1 + res_mod1 , data = data_stations, min) |
|
| 297 |
dim(data_stations_temp) |
|
| 298 |
#> dim(data_stations_temp) |
|
| 299 |
#[1] 11171 5 |
|
| 300 |
|
|
| 301 |
data_stations_temp$date_str <- data_stations_temp$date |
|
| 302 |
data_stations_temp$date <- as.Date(strptime(data_stations_temp$date_str,"%Y%m%d")) |
|
| 303 |
coordinates(data_stations_temp) <- cbind(data_stations_temp$x,data_stations_temp$y) |
|
| 304 |
data_stations_temp$constant <- c(1000,2000) |
|
| 305 |
#plot(data_stations_temp) |
|
| 306 |
#plot(reg_layer) |
|
| 307 |
#res_pix <- 1200 |
|
| 308 |
res_pix <- 400 |
|
| 309 |
col_mfrow <- 1 |
|
| 310 |
row_mfrow <- 1 |
|
| 311 |
|
|
| 312 |
png_filename <- paste("Figure_ac_metrics_map_stations_locations_validation_",model_name,"_",date_processed,
|
|
| 313 |
"_",out_suffix,".png",sep="") |
|
| 314 |
png(filename=png_filename, |
|
| 315 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 316 |
#plot(data_stations_temp |
|
| 317 |
p_shp <- spplot(reg_layer,"ISO3" ,col.regions=NA, col="black") #ok problem solved!! |
|
| 318 |
#title("(a) Mean for 1 January")
|
|
| 319 |
#p <- bubble(data_stations_temp,"constant",main=paste0("Average Residuals by validation stations ",
|
|
| 320 |
# date_processed, |
|
| 321 |
# "for ",var_selected)) |
|
| 322 |
p <- spplot(data_stations_temp,"constant",col.regions=NA, col="black", |
|
| 323 |
main=paste0("Average Residuals by validation stations ",pch=3,cex=10,
|
|
| 324 |
date_processed, "for ",var_selected)) |
|
| 325 |
|
|
| 326 |
p1 <- p+p_shp |
|
| 327 |
try(print(p1)) #error raised if number of missing values below a threshold does not exist |
|
| 328 |
|
|
| 329 |
dev.off() |
|
| 330 |
|
|
| 331 |
res_pix <- 800 |
|
| 332 |
col_mfrow <- 1 |
|
| 333 |
row_mfrow <- 1 |
|
| 334 |
png_filename <- paste("Figure_ac_metrics_map_stations_validation_averaged_",model_name,"_",date_processed,
|
|
| 335 |
"_",out_suffix,".png",sep="") |
|
| 336 |
png(filename=png_filename, |
|
| 337 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 338 |
|
|
| 339 |
#model_name[j] |
|
| 340 |
|
|
| 341 |
#p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black')) |
|
| 342 |
p_shp <- spplot(reg_layer,"ISO3" ,col.regions=NA, col="black") #ok problem solved!! |
|
| 343 |
#title("(a) Mean for 1 January")
|
|
| 344 |
p <- bubble(data_stations_temp,"res_mod1",main=paste0("Average Residuals by validation stations ",date_processed,
|
|
| 345 |
"for ",var_selected)) |
|
| 346 |
p1 <- p+p_shp |
|
| 347 |
try(print(p1)) #error raised if number of missing values below a threshold does not exist |
|
| 348 |
|
|
| 349 |
dev.off() |
|
| 350 |
|
|
| 351 |
} |
|
| 352 |
|
|
| 353 |
|
|
| 257 | 354 |
|
| 258 | 355 |
############### PART2: Select stations by ID to build a time series ############# |
| 259 | 356 |
#### Extract corresponding stationsfor given dates and plot stations used |
Also available in: Unified diff
global product part1, combined figure of stations per day for Africa