Revision be7d0e21
Added by Benoit Parmentier over 10 years ago
| climate/research/oregon/interpolation/global_run_scalingup_assessment_part2.R | ||
|---|---|---|
| 84 | 84 |
#in_dir_shp <- "/nobackupp4/aguzman4/climateLayers/output4/subset/shapefiles/" |
| 85 | 85 |
|
| 86 | 86 |
df_tile_processed <- read.table(file=file.path(out_dir,paste("df_tile_processed_",out_prefix,".txt",sep="")),sep=",")
|
| 87 |
in_dir_list <- file.path(in_dir1,read.table(file.path(in_dir1,"processed.txt"))$V1) |
|
| 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 | 90 |
out_prefix<-"run3_global_analyses_06192014" |
| ... | ... | |
| 111 | 111 |
#lf_tables <- list.files(out_dir,) |
| 112 | 112 |
summary_metrics_v <- read.table(file=file.path(out_dir,paste("summary_metrics_v2_NA_",out_prefix,".txt",sep="")),sep=",")
|
| 113 | 113 |
tb <- read.table(file=file.path(out_dir,paste("tb_diagnostic_v_NA","_",out_prefix,".txt",sep="")),sep=",")
|
| 114 |
df_tile_processed <- read.table(file=file.path(out_dir,paste("df_tile_processed_",out_prefix,".txt",sep="")),sep=",")
|
|
| 114 |
df_tile_processed <- read.table(file=file.path(out_dir,paste("df_tile_processed_",out_prefix,".txt",sep="")),sep=","),paste("pred_data_month_info_",out_prefix,".txt",sep="")) <- read.table(file=file.path(out_dir,paste("pred_data_month_info_",out_prefix,".txt",sep="")),sep=",")
|
|
| 115 |
pred_data_month_info <- read.table(file=file.path(out_dir,paste("pred_data_month_info_",out_prefix,".txt",sep="")),sep=",")
|
|
| 116 |
pred_data_day_info <- read.table(file=file.path(out_dir,paste("pred_data_day_info_",out_prefix,".txt",sep="")),sep=",")
|
|
| 115 | 117 |
|
| 116 | 118 |
########################## START SCRIPT ############################## |
| 117 | 119 |
|
| ... | ... | |
| 123 | 125 |
#matching_index <- match(l_shp,list_shp_reg_files) |
| 124 | 126 |
matching_index <- match(list_shp_reg_files,l_shp) |
| 125 | 127 |
df_tile_processed$shp_files <-list_shp_world[matching_index] |
| 128 |
list_shp_reg_files <- df_tile_processed$shp_files |
|
| 126 | 129 |
|
| 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 |
|
| 130 |
#df_tiled_processed <- na.omit(df_tile_processed) #remove other list of folders irrelevant
|
|
| 131 |
#list_shp_reg_files <- df_tiled_processed$shp_files
|
|
| 129 | 132 |
list_shp_reg_files<- file.path(in_dir_shp,list_shp_reg_files) |
| 130 | 133 |
|
| 131 | 134 |
### First get background map to display where study area is located |
| ... | ... | |
| 133 | 136 |
if region_name=="USA"{
|
| 134 | 137 |
usa_map <- getData('GADM', country='USA', level=1) #Get US map
|
| 135 | 138 |
#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
|
|
| 139 |
usa_map <- usa_map[usa_map$NAME_1!="Alaska",] #remove Alaska
|
|
| 137 | 140 |
reg_layer <- usa_map[usa_map$NAME_1!="Hawaii",] #remove Hawai |
| 138 | 141 |
} |
| 139 | 142 |
|
| ... | ... | |
| 162 | 165 |
shps_tiles[[i]] <- shp1 |
| 163 | 166 |
} |
| 164 | 167 |
|
| 168 |
|
|
| 165 | 169 |
#plot info: with labels |
| 166 | 170 |
res_pix <- 480 |
| 167 | 171 |
col_mfrow <- 1 |
| ... | ... | |
| 389 | 393 |
## plotting of delta and clim for later scripts... |
| 390 | 394 |
|
| 391 | 395 |
###################### |
| 392 |
### Figure 8: Histograms... |
|
| 396 |
### Figure 9: Plot the number of stations in a processing tile |
|
| 397 |
|
|
| 398 |
dd <- merge(df_tile_processed,pred_data_month_info,"tile_id") |
|
| 399 |
coordinates(dd) <- c(dd$x,dd$y) |
|
| 400 |
|
|
| 401 |
## Make this a function later... |
|
| 402 |
list_shp_tmp <- vector("list",length(shps_tiles))
|
|
| 403 |
for(i in 1:length(shps_tiles)){
|
|
| 404 |
#test=spRbind(shps_tiles[[1]],shps_tiles[[2]]) |
|
| 405 |
shp1 <- shps_tiles[[i]] |
|
| 406 |
|
|
| 407 |
ID_str <- unlist(strsplit(as.character(df_tile_processed$tile_id[i]),"_"))[2] |
|
| 408 |
shp1$FID <- ID_str |
|
| 409 |
shp1<- spChFIDs(shp1, as.character(shp1$FID)) #assign ID |
|
| 410 |
list_shp_tmp[[i]] <-shp1 |
|
| 411 |
} |
|
| 412 |
|
|
| 413 |
combined_shp <- list_shp_tmp[[1]] |
|
| 414 |
for(i in 2:length(list_shp_tmp)){
|
|
| 415 |
combined_shp <- rbind(combined_shp,list_shp_tmp[[i]]) |
|
| 416 |
#sapply(slot(shps_tiles[[2]], "polygons"), function(x) slot(x, "ID")) |
|
| 417 |
#rownames(as(alaska.tract, "data.frame")) |
|
| 418 |
} |
|
| 419 |
|
|
| 420 |
combined_shp$tile_id <- df_tile_processed$tile_id |
|
| 421 |
|
|
| 422 |
test <- combined_shp |
|
| 423 |
test2 <- merge(test,pred_data_month_info, by="tile_id") |
|
| 424 |
|
|
| 425 |
r <- raster(lf_pred_list[i]) |
|
| 426 |
|
|
| 427 |
plot(combined_shp) |
|
| 428 |
# polygons |
|
| 429 |
plot(combined_shp, col = fillColour, border = outlineColour) |
|
| 430 |
|
|
| 431 |
p0 <- spplot(combined_shp, "Stations",col.regions=matlab.like(100)) |
|
| 432 |
p1 <- spplot(reg_layer,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color? |
|
| 433 |
p0 +p1 |
|
| 434 |
|
|
| 435 |
### Now plot number of training for monthly data |
|
| 436 |
|
|
| 437 |
df_dat <- subset(pred_data_month_info, pred_mod == "mod1" & date =="20100115") |
|
| 438 |
#shp_dat <-merge(combined_shp,df_dat,by="tile_id") |
|
| 439 |
shp_dat <-merge(x=combined_shp,y=df_dat,by="tile_id",all.x=T) #if tile is missing then add rows with NA |
|
| 440 |
#shp_dat <- merge(shp_dat,df_tile_processed,by="tile_id") |
|
| 441 |
#coordinates(shp_dat) <- cbind(shp_dat$lon,shp_dat$lat) |
|
| 442 |
#proj4string(shp_dat) <- CRS_WGS84 |
|
| 443 |
|
|
| 444 |
#test <- overlay(combined_shp,shp_dat) |
|
| 445 |
pol <- SpatialPolygons(combined_shp@polygons,proj4string=CRS(CRS_WGS84)) |
|
| 446 |
spp <- SpatialPolygonsDataFrame(pol,data=shp_dat) |
|
| 393 | 447 |
|
| 448 |
p0 <- spplot(spp, "n_mod",col.regions=matlab.like(100)) |
|
| 449 |
p1 <- spplot(reg_layer,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color? |
|
| 450 |
p0 +p1 |
|
| 394 | 451 |
|
| 395 | 452 |
##################### END OF SCRIPT ###################### |
Also available in: Unified diff
assessment NEX run part2: adding plots of number of stations per tiles