Revision a408f6cb
Added by Benoit Parmentier over 10 years ago
| climate/research/oregon/interpolation/global_run_scalingup_assessment_part2.R | ||
|---|---|---|
| 5 | 5 |
#Analyses, figures, tables and data are also produced in the script. |
| 6 | 6 |
#AUTHOR: Benoit Parmentier |
| 7 | 7 |
#CREATED ON: 03/23/2014 |
| 8 |
#MODIFIED ON: 06/01/2014
|
|
| 8 |
#MODIFIED ON: 06/19/2014
|
|
| 9 | 9 |
#Version: 3 |
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 | 11 |
#COMMENTS: analyses for run 3 global using 2 specific tiles |
| ... | ... | |
| 69 | 69 |
#on ATLAS |
| 70 | 70 |
#in_dir1 <- "/data/project/layers/commons/NEX_data/test_run1_03232014/output" #On Atlas |
| 71 | 71 |
#parent output dir : contains subset of the data produced on NEX |
| 72 |
in_dir1 <- "/data/project/layers/commons/NEX_data/output_run3_global_analyses_05292014/output"
|
|
| 72 |
in_dir1 <- "/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/output"
|
|
| 73 | 73 |
# parent output dir for the curent script analyes |
| 74 |
out_dir <-"/data/project/layers/commons/NEX_data/output_run3_global_analyses_05292014/" #On NCEAS Atlas
|
|
| 75 |
out_dir <-"/data/project/layers/commons/NEX_data/output_run3_global_analyses_05292014/"
|
|
| 74 |
#out_dir <-"/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/" #On NCEAS Atlas
|
|
| 75 |
out_dir <-"/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/"
|
|
| 76 | 76 |
# input dir containing shapefiles defining tiles |
| 77 |
in_dir_shp <- "/data/project/layers/commons/NEX_data/output_run3_global_analyses_05292014/output/subset/shapefiles"
|
|
| 77 |
in_dir_shp <- "/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/output/subset/shapefiles"
|
|
| 78 | 78 |
|
| 79 | 79 |
#On NEX |
| 80 | 80 |
#contains all data from the run by Alberto |
| ... | ... | |
| 87 | 87 |
in_dir_list <- file.path(in_dir1,read.table(file.path(in_dir1,"processed.txt"))$V1) |
| 88 | 88 |
y_var_name <- "dailyTmax" |
| 89 | 89 |
interpolation_method <- c("gam_CAI")
|
| 90 |
out_prefix<-"run3_global_analyses_05292014"
|
|
| 90 |
out_prefix<-"run3_global_analyses_06192014"
|
|
| 91 | 91 |
|
| 92 | 92 |
#out_dir <-paste(out_dir,"_",out_prefix,sep="") |
| 93 | 93 |
create_out_dir_param <- FALSE |
| ... | ... | |
| 102 | 102 |
setwd(out_dir) |
| 103 | 103 |
|
| 104 | 104 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 105 |
CRS_WGS84<-c("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
|
| 105 | 106 |
|
| 106 |
region_name <- "USA"
|
|
| 107 |
region_name <- "world"
|
|
| 107 | 108 |
|
| 108 | 109 |
###Table 1: Average accuracy metrics |
| 109 | 110 |
###Table 2: daily accuracy metrics for all tiles |
| ... | ... | |
| 116 | 117 |
|
| 117 | 118 |
############### |
| 118 | 119 |
### Figure 1: plot location of the study area with tiles processed |
| 120 |
list_shp_world <- list.files(in_dir_shp,".shp") |
|
| 121 |
l_shp <- unlist(lapply(1:length(list_shp_world),FUN=function(i){paste(strsplit(list_shp_world[i],"_")[[1]][2:3],collapse="_")}))
|
|
| 122 |
list_shp_reg_files <- as.character(df_tile_processed$tile_coord) |
|
| 123 |
#matching_index <- match(l_shp,list_shp_reg_files) |
|
| 124 |
matching_index <- match(list_shp_reg_files,l_shp) |
|
| 125 |
df_tile_processed$shp_files <-list_shp_world[matching_index] |
|
| 126 |
|
|
| 127 |
df_tiled_processed <- na.omit(df_tile_processed) #remove other list of folders irrelevant |
|
| 128 |
list_shp_reg_files <- df_tiled_processed$shp_files |
|
| 129 |
list_shp_reg_files<- file.path(in_dir_shp,list_shp_reg_files) |
|
| 119 | 130 |
|
| 120 |
list_shp_reg_files <- file.path(in_dir_shp,df_tile_processed$shp_files) |
|
| 121 |
|
|
| 122 | 131 |
### First get background map to display where study area is located |
| 123 |
#can make this more general later on.. |
|
| 124 |
#usa_map <- getData('GADM', country='USA', level=1) #Get US map
|
|
| 125 |
usa_map <- getData('GADM', country=region_name,level=1) #Get US map, this is not working right now
|
|
| 126 |
usa_map_2 <- usa_map[usa_map$NAME_1!="Alaska",] #remove Alaska |
|
| 127 |
usa_map_2 <- usa_map_2[usa_map_2$NAME_1!="Hawaii",] #remove Hawai |
|
| 132 |
#can make this more general later on.. |
|
| 133 |
if region_name=="USA"{
|
|
| 134 |
usa_map <- getData('GADM', country='USA', level=1) #Get US map
|
|
| 135 |
#usa_map <- getData('GADM', country=region_name,level=1) #Get US map, this is not working right now
|
|
| 136 |
#usa_map_2 <- usa_map[usa_map$NAME_1!="Alaska",] #remove Alaska |
|
| 137 |
reg_layer <- usa_map[usa_map$NAME_1!="Hawaii",] #remove Hawai |
|
| 138 |
} |
|
| 139 |
|
|
| 140 |
if region_name=="world"{
|
|
| 141 |
#http://www.diva-gis.org/Data |
|
| 142 |
countries_shp <-"/data/project/layers/commons/NEX_data/countries.shp" |
|
| 143 |
path_to_shp <- dirname(countries_shp) |
|
| 144 |
layer_name <- sub(".shp","",basename(countries_shp))
|
|
| 145 |
reg_layer <- readOGR(path_to_shp, layer_name) |
|
| 146 |
proj4string(reg_layer) <- CRS_locs_WGS84 |
|
| 147 |
#reg_shp<-readOGR(dirname(list_shp_reg_files[[i]]),sub(".shp","",basename(list_shp_reg_files[[i]])))
|
|
| 148 |
|
|
| 149 |
} |
|
| 128 | 150 |
|
| 129 | 151 |
centroids_pts <- vector("list",length(list_shp_reg_files))
|
| 130 | 152 |
shps_tiles <- vector("list",length(list_shp_reg_files))
|
| 131 | 153 |
#collect info: read in all shapfiles |
| 154 |
#This is slow...make a function and use mclapply?? |
|
| 132 | 155 |
for(i in 1:length(list_shp_reg_files)){
|
| 133 | 156 |
path_to_shp <- dirname(list_shp_reg_files[[i]]) |
| 134 |
layer_name <- basename(list_shp_reg_files[[i]])
|
|
| 157 |
layer_name <- sub(".shp","",basename(list_shp_reg_files[[i]]))
|
|
| 135 | 158 |
shp1 <- readOGR(path_to_shp, layer_name) |
| 136 |
shp1<-readOGR(dirname(list_shp_reg_files[[i]]),sub(".shp","",basename(list_shp_reg_files[[i]])))
|
|
| 159 |
#shp1<-readOGR(dirname(list_shp_reg_files[[i]]),sub(".shp","",basename(list_shp_reg_files[[i]])))
|
|
| 137 | 160 |
pt <- gCentroid(shp1) |
| 138 | 161 |
centroids_pts[[i]] <-pt |
| 139 | 162 |
shps_tiles[[i]] <- shp1 |
| ... | ... | |
| 144 | 167 |
col_mfrow <- 1 |
| 145 | 168 |
row_mfrow <- 1 |
| 146 | 169 |
|
| 147 |
png(filename=paste("Figure1_tile_processed_region_",out_prefix,".png",sep=""),
|
|
| 170 |
png(filename=paste("Figure1_tile_processed_region_",region_name,"_",out_prefix,".png",sep=""),
|
|
| 148 | 171 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 149 | 172 |
|
| 150 |
plot(usa_map_2)
|
|
| 173 |
plot(reg_layer)
|
|
| 151 | 174 |
#Add polygon tiles... |
| 152 | 175 |
for(i in 1:length(list_shp_reg_files)){
|
| 153 | 176 |
shp1 <- shps_tiles[[i]] |
| ... | ... | |
| 237 | 260 |
#y_var_name <-"dailyTmax" |
| 238 | 261 |
#index <-244 #index corresponding to Sept 1 |
| 239 | 262 |
index <- 1 #index corresponding to Jan 1 |
| 240 |
date_selected <- "20100101"
|
|
| 263 |
date_selected <- "20100901"
|
|
| 241 | 264 |
name_method_var <- paste(interpolation_method,"_",y_var_name,"_",sep="") |
| 242 | 265 |
|
| 243 | 266 |
pattern_str <- paste("mosaiced","_",name_method_var,"predicted",".*.",date_selected,".*.tif",sep="")
|
| ... | ... | |
| 251 | 274 |
col_mfrow <- 1 |
| 252 | 275 |
row_mfrow <- 1 |
| 253 | 276 |
|
| 254 |
png(filename=paste("Figure4_models_predicted_surfaces_",model_name[i],"_",name_method_var,out_prefix,".png",sep=""),
|
|
| 277 |
png(filename=paste("Figure4_models_predicted_surfaces_",model_name[i],"_",name_method_var,"_",data_selected,"_",out_prefix,".png",sep=""),
|
|
| 255 | 278 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 256 | 279 |
|
| 257 | 280 |
plot(r_pred) |
| ... | ... | |
| 260 | 283 |
|
| 261 | 284 |
} |
| 262 | 285 |
|
| 263 |
## plotting of delta and clim for later scripts... |
|
| 286 |
####### Figure 5... |
|
| 287 |
### Adding tiles do a plot of mod1 with tiles |
|
| 288 |
|
|
| 289 |
r_pred <- raster(lf_pred_list[i]) |
|
| 290 |
|
|
| 291 |
res_pix <- 480 |
|
| 292 |
col_mfrow <- 1 |
|
| 293 |
row_mfrow <- 1 |
|
| 294 |
|
|
| 295 |
png(filename=paste("Figure5_tiles_with_models_predicted_surfaces_",model_name[i],"_",name_method_var,out_prefix,".png",sep=""),
|
|
| 296 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 297 |
|
|
| 298 |
plot(r_pred) |
|
| 299 |
title(paste("Mosaiced",model_name[i],name_method_var,date_selected,"with tiles",sep=" "))
|
|
| 300 |
|
|
| 301 |
#Add polygon tiles... |
|
| 302 |
for(i in 1:length(list_shp_reg_files)){
|
|
| 303 |
shp1 <- shps_tiles[[i]] |
|
| 304 |
pt <- centroids_pts[[i]] |
|
| 305 |
plot(shp1,add=T,border="blue") |
|
| 306 |
#plot(pt,add=T,cex=2,pch=5) |
|
| 307 |
label_id <- df_tile_processed$tile_id[i] |
|
| 308 |
text(coordinates(pt)[1],coordinates(pt)[2],labels=i,cex=1,col=c("red"))
|
|
| 309 |
} |
|
| 310 |
#title(paste("Prediction with tiles location 10x10 degrees for ", region_name,sep=""))
|
|
| 311 |
dev.off() |
|
| 312 |
|
|
| 313 |
### |
|
| 264 | 314 |
|
| 265 | 315 |
#### Now combined plot... |
| 316 |
#Use the function to match extent... |
|
| 266 | 317 |
|
| 267 | 318 |
#pred_s <- stack(lf_list) #problem different extent!! |
| 268 | 319 |
#methods_names <-c("gam","kriging","gwr")
|
| ... | ... | |
| 295 | 346 |
#dev.off() |
| 296 | 347 |
|
| 297 | 348 |
###################### |
| 298 |
### Figure 5: plot map of average RMSE per tile at centroids
|
|
| 349 |
### Figure 6: plot map of average RMSE per tile at centroids
|
|
| 299 | 350 |
|
| 300 | 351 |
#Turn summary table to a point shp |
| 301 | 352 |
|
| 302 | 353 |
coordinates(summary_metrics_v) <- cbind(summary_metrics_v$lon,summary_metrics_v$lat) |
| 303 | 354 |
proj4string(summary_metrics_v) <- CRS_WGS84 |
| 355 |
lf_list <- lf_pred_list |
|
| 304 | 356 |
list_df_ac_mod <- vector("list",length=length(lf_pred_list))
|
| 305 | 357 |
for (i in 1:length(lf_list)){
|
| 306 | 358 |
|
| ... | ... | |
| 311 | 363 |
col_mfrow <- 1 |
| 312 | 364 |
row_mfrow <- 1 |
| 313 | 365 |
|
| 314 |
png(filename=paste("Figure5_ac_metrics_map_centroids_tile_",model_name[i],"_",out_prefix,".png",sep=""),
|
|
| 366 |
png(filename=paste("Figure6_ac_metrics_map_centroids_tile_",model_name[i],"_",out_prefix,".png",sep=""),
|
|
| 315 | 367 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 316 | 368 |
|
| 317 | 369 |
plot(r_pred) |
| ... | ... | |
| 332 | 384 |
#autokrige(rmse~1,r2,) |
| 333 | 385 |
|
| 334 | 386 |
###################### |
| 335 |
### Figure 6: Histograms... |
|
| 387 |
### Figure 7: Delta and clim... |
|
| 388 |
|
|
| 389 |
## plotting of delta and clim for later scripts... |
|
| 390 |
|
|
| 391 |
###################### |
|
| 392 |
### Figure 8: Histograms... |
|
| 336 | 393 |
|
| 337 | 394 |
|
| 338 | 395 |
##################### END OF SCRIPT ###################### |
Also available in: Unified diff
assessment NEX run 3: modifications to analyze results with figures