Revision 2835e195
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/19/2014
|
|
| 8 |
#MODIFIED ON: 07/02/2014
|
|
| 9 | 9 |
#Version: 3 |
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 | 11 |
#COMMENTS: analyses for run 3 global using 2 specific tiles |
| ... | ... | |
| 133 | 133 |
|
| 134 | 134 |
### First get background map to display where study area is located |
| 135 | 135 |
#can make this more general later on.. |
| 136 |
if region_name=="USA"{
|
|
| 137 |
usa_map <- getData('GADM', country='USA', level=1) #Get US map
|
|
| 136 |
if region_name=="reg1"{
|
|
| 137 |
usa_map <- getData('GADM', country=c('USA'), level=1) #Get US map
|
|
| 138 |
can_map <- getData('GADM', country=c('CAN'), level=1) #Get Canada map
|
|
| 139 |
mex_map <- getData('GADM', country=c('MEX'), level=1) #Get Mexico map
|
|
| 140 |
|
|
| 138 | 141 |
#usa_map <- getData('GADM', country=region_name,level=1) #Get US map, this is not working right now
|
| 139 | 142 |
usa_map <- usa_map[usa_map$NAME_1!="Alaska",] #remove Alaska |
| 140 |
reg_layer <- usa_map[usa_map$NAME_1!="Hawaii",] #remove Hawai |
|
| 143 |
usa_map <- usa_map[usa_map$NAME_1!="Hawaii",] #remove Hawai |
|
| 144 |
#reg_layer <- combine all three? |
|
| 141 | 145 |
} |
| 142 | 146 |
|
| 143 | 147 |
if region_name=="world"{
|
| ... | ... | |
| 375 | 379 |
#plot(ac_mod1,cex=sqrt(ac_mod1$rmse),pch=1,add=T) |
| 376 | 380 |
plot(ac_mod,cex=(ac_mod$rmse^2)/10,pch=1,add=T) |
| 377 | 381 |
#plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T) |
| 378 |
title(paste("Averrage RMSE per tile and by ",model_name[i]))
|
|
| 382 |
title(paste("Average RMSE per tile and by ",model_name[i]))
|
|
| 379 | 383 |
|
| 380 | 384 |
dev.off() |
| 381 | 385 |
|
| ... | ... | |
| 395 | 399 |
###################### |
| 396 | 400 |
### Figure 9: Plot the number of stations in a processing tile |
| 397 | 401 |
|
| 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 |
|
| 402 |
## Get ID from tile number... |
|
| 403 |
ID_str <- unlist(lapply(1:nrow(df_tile_processed),function(i){unlist(strsplit(as.character(df_tile_processed$tile_id[i]),"_"))[2]}))
|
|
| 404 |
|
|
| 405 |
assign_FID_spatial_polygons_df <-function(list_spdf,ID_str=NULL){
|
|
| 406 |
list_spdf_tmp <- vector("list",length(list_spdf))
|
|
| 407 |
if(is.null(ID_str)){
|
|
| 408 |
nf <- 0 #number of features |
|
| 409 |
#for(i in 1:length(spdf)){
|
|
| 410 |
# shp1 <- list_spdf[[i]] |
|
| 411 |
# f <- nrow(shp1) |
|
| 412 |
# nf <- nf + f |
|
| 413 |
#} |
|
| 414 |
#This assumes that the list has one feature per item list |
|
| 415 |
nf <- length(list_spdf) |
|
| 416 |
ID_str <- as.character(1:nf) |
|
| 417 |
} |
|
| 418 |
for(i in 1:length(list_spdf)){
|
|
| 419 |
#test=spRbind(shps_tiles[[1]],shps_tiles[[2]]) |
|
| 420 |
shp1 <- list_spdf[[i]] |
|
| 421 |
shp1$FID <- ID_str |
|
| 422 |
shp1<- spChFIDs(shp1, as.character(shp1$FID)) #assign ID |
|
| 423 |
list_spdf_tmp[[i]] <-shp1 |
|
| 424 |
} |
|
| 425 |
return(list_spdf_tmp) |
|
| 411 | 426 |
} |
| 412 | 427 |
|
| 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")) |
|
| 428 |
combine_spatial_polygons_df_fun <- function(list_spdf_tmp,ID_str=NULL){
|
|
| 429 |
if(is.null(ID_str)){
|
|
| 430 |
#call function |
|
| 431 |
list_spdf_tmp <- assign_FID_spatial_polygons_df |
|
| 432 |
} |
|
| 433 |
combined_spdf <- list_spdf_tmp[[1]] |
|
| 434 |
for(i in 2:length(list_spdf_tmp)){
|
|
| 435 |
combined_spdf <- rbind(combined_spdf,list_spdf_tmp[[i]]) |
|
| 436 |
#sapply(slot(shps_tiles[[2]], "polygons"), function(x) slot(x, "ID")) |
|
| 437 |
#rownames(as(alaska.tract, "data.frame")) |
|
| 438 |
} |
|
| 439 |
return(combined_spdf) |
|
| 418 | 440 |
} |
| 419 | 441 |
|
| 420 | 442 |
combined_shp$tile_id <- df_tile_processed$tile_id |
| 421 | 443 |
|
| 422 |
test <- combined_shp |
|
| 423 |
test2 <- merge(test,pred_data_month_info, by="tile_id") |
|
| 424 |
|
|
| 425 | 444 |
r <- raster(lf_pred_list[i]) |
| 426 | 445 |
|
| 427 | 446 |
plot(combined_shp) |
| 428 |
# polygons |
|
| 429 |
plot(combined_shp, col = fillColour, border = outlineColour) |
|
| 430 | 447 |
|
| 431 | 448 |
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 |
|
| 449 |
p1 <- spplot(usa_map,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color? |
|
| 450 |
p2 <- spplot(can_map,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color? |
|
| 451 |
p3 <- spplot(mex_map,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color? |
|
| 452 |
|
|
| 453 |
p0 +p1+p2+p3 |
|
| 434 | 454 |
|
| 435 | 455 |
### Now plot number of training for monthly data |
| 436 | 456 |
|
| ... | ... | |
| 446 | 466 |
spp <- SpatialPolygonsDataFrame(pol,data=shp_dat) |
| 447 | 467 |
|
| 448 | 468 |
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 |
|
| 469 |
#p1 <- spplot(reg_layer,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color? |
|
| 470 |
p0 + p1 + p2 + p3 |
|
| 471 |
|
|
| 472 |
###################### |
|
| 473 |
### Figure 10: Plot the number of stations in a processing tile |
|
| 474 |
|
|
| 475 |
boxplot(n_mod~tile_id,pred_data_month_info) |
|
| 476 |
boxplot(n_mod~date,pred_data_month_info) |
|
| 451 | 477 |
|
| 452 | 478 |
##################### END OF SCRIPT ###################### |
Also available in: Unified diff
assessment NEX run 3 part2: adding function to combine shapefiles to create plots of number of stations per tiles and accuracy