Revision e9998467
Added by Benoit Parmentier over 11 years ago
| climate/research/oregon/interpolation/methods_comparison_assessment_part5.R | ||
|---|---|---|
| 1 | 1 |
##################################### METHODS COMPARISON part 5 ########################################## |
| 2 | 2 |
#################################### Spatial Analysis ############################################ |
| 3 | 3 |
#This script utilizes the R ojbects created during the interpolation phase. # |
| 4 |
#At this stage the script produces figures of various accuracy metrics and compare methods: # |
|
| 5 |
#This scripts focuses on a detailed studay of differences in the predictions of CAI_kr and FUsion_Kr # |
|
| 4 |
#This scripts focuses on a detailed study of differences in the predictions of CAI_kr and FUsion_Kr |
|
| 5 |
#Differences are examined through: |
|
| 6 |
#1) per land cover classes |
|
| 7 |
#2) per elevation classes |
|
| 8 |
#3) through spiatial transects |
|
| 9 |
# |
|
| 10 |
#Note this code is for exploratory analyses so some sections are not succinct and |
|
| 11 |
#can be improve for repeatability and clarity. |
|
| 6 | 12 |
#AUTHOR: Benoit Parmentier # |
| 7 |
#DATE: 11/23/2012 #
|
|
| 13 |
#DATE: 12/04/2012 #
|
|
| 8 | 14 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491 -- # |
| 9 | 15 |
################################################################################################### |
| 10 | 16 |
|
| ... | ... | |
| 182 | 188 |
infile6<-"OR83M_state_outline.shp" |
| 183 | 189 |
#stat_loc<-read.table(paste(path,"/","location_study_area_OR_0602012.txt",sep=""),sep=",", header=TRUE) |
| 184 | 190 |
|
| 185 |
out_prefix<-"methods_11292012_"
|
|
| 191 |
out_prefix<-"methods_comp5_12042012_"
|
|
| 186 | 192 |
nb_transect<-4 |
| 187 | 193 |
##### LOAD USEFUL DATA |
| 188 | 194 |
|
| 189 | 195 |
#obj_list<-"list_obj_08262012.txt" #Results of fusion from the run on ATLAS |
| 190 |
path<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison_10242012" #Jupiter LOCATION on Atlas for kriging #Jupiter Location on XANDERS |
|
| 196 |
path_wd<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison_10242012" #Jupiter LOCATION on Atlas for kriging #Jupiter Location on XANDERS
|
|
| 191 | 197 |
#path<-"/Users/benoitparmentier/Dropbox/Data/NCEAS/Oregon_covariates/" #Local dropbox folder on Benoit's laptop |
| 192 |
setwd(path) |
|
| 198 |
setwd(path_wd) |
|
| 199 |
path_data_cai<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_CAI" #Change to constant |
|
| 200 |
path_data_fus<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_GAM" |
|
| 201 |
#list files that contain model objects and ratingin-testing information for CAI and Fusion |
|
| 202 |
obj_mod_fus_name<-"results_mod_obj__365d_GAM_fusion_const_all_lstd_11022012.RData" |
|
| 203 |
obj_mod_cai_name<-"results_mod_obj__365d_GAM_CAI2_const_all_10312012.RData" |
|
| 204 |
|
|
| 205 |
gam_fus<-load_obj(file.path(path_data_fus,obj_mod_fus_name)) |
|
| 206 |
gam_cai<-load_obj(file.path(path_data_cai,obj_mod_cai_name)) #This contains all the info |
|
| 207 |
sampling_date_list<-gam_fus$sampling_obj$sampling_dat$date |
|
| 208 |
|
|
| 209 |
### Projection for the current region |
|
| 193 | 210 |
proj_str="+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs"; |
| 194 | 211 |
#User defined output prefix |
| 195 | 212 |
|
| 213 |
### MAKE THIS A FUNCTION TO LOAD STACK AND DEFINE VALID RANGE... |
|
| 196 | 214 |
#CRS<-proj4string(ghcn) #Storing projection information (ellipsoid, datum,etc.) |
| 197 | 215 |
lines<-read.table(paste(path,"/",inlistf,sep=""), sep="") #Column 1 contains the names of raster files |
| 198 | 216 |
inlistvar<-lines[,1] |
| ... | ... | |
| 516 | 534 |
avg_elev_rec_forest<-zonal(rast_diff,zones=elev_rec_forest,stat="mean",na.rm=TRUE) |
| 517 | 535 |
std_elev_rec_forest<-zonal(rast_diff,zones=elev_rec_forest,stat="sd",na.rm=TRUE) |
| 518 | 536 |
|
| 519 |
|
|
| 520 |
|
|
| 521 | 537 |
## CREATE plots |
| 522 | 538 |
X11() |
| 523 | 539 |
plot(avg_elev_rec[,1],avg_elev_rec[,2],type="b",ylim=c(-10,1), |
| ... | ... | |
| 534 | 550 |
} |
| 535 | 551 |
|
| 536 | 552 |
################################################################### |
| 537 |
################ TRANSECT THROUGH THE IMAGE: #################### |
|
| 553 |
################ SPATIAL TRANSECT THROUGH THE IMAGE: ####################
|
|
| 538 | 554 |
|
| 539 | 555 |
#select date |
| 540 | 556 |
dates<-c("20100103","20100901")
|
| 541 |
j=1
|
|
| 557 |
#j=2
|
|
| 542 | 558 |
|
| 543 | 559 |
for (j in 1:length(dates)){
|
| 544 | 560 |
|
| 545 | 561 |
#Read predicted tmax raster surface and modeling information |
| 546 | 562 |
date_selected<-dates[j] |
| 563 |
|
|
| 547 | 564 |
oldpath<-getwd() |
| 548 | 565 |
setwd(path_data_cai) |
| 549 | 566 |
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_CAI2_const_all_10312012.rst",sep="")) #Search for files in relation to fusion
|
| ... | ... | |
| 562 | 579 |
rast_fus_pred<-raster(rast_fus1c,1) |
| 563 | 580 |
rast_cai_pred<-raster(rast_cai2c,1) |
| 564 | 581 |
rast_diff_fc<-rast_fus_pred-rast_cai_pred |
| 582 |
|
|
| 565 | 583 |
#Read in data_s and data_v |
| 584 |
|
|
| 585 |
k<-match(date_selected,sampling_date_list) |
|
| 586 |
names(gam_fus$gam_fus_mod[[k]]) #Show the name structure of the object/list |
|
| 566 | 587 |
|
| 588 |
#Extract the training and testing information for the given date... |
|
| 589 |
data_sf<-gam_fus$gam_fus_mod[[k]]$data_s #object for the first date...20100103 #Make this a function?? |
|
| 590 |
data_vf<-gam_fus$gam_fus_mod[[k]]$data_v #object for the first date...20100103 |
|
| 591 |
data_sc<-gam_cai$gam_CAI_mod[[k]]$data_s #object for the first date...20100103 |
|
| 592 |
data_vc<-gam_cai$gam_CAI_mod[[k]]$data_v #object for the first date...20100103 |
|
| 567 | 593 |
|
| 568 | 594 |
### CREATE A NEW TRANSECT BASED ON LOCATION OF SPECIFIED STATIONS |
| 569 | 595 |
|
| ... | ... | |
| 574 | 600 |
data_vf$training<-rep(0,nrow(data_vf)) |
| 575 | 601 |
data_sf$training<-rep(1,nrow(data_sf)) |
| 576 | 602 |
|
| 577 |
data_stat<-rbind(data_vf[,c("id","training")],data_sf[,c("id","training")])
|
|
| 603 |
data_stat<-rbind(data_vf[,c("id","training")],data_sf[,c("id","training")]) #bringing together data_v and data_s
|
|
| 578 | 604 |
m<-match(selected_stations,data_stat$id) |
| 579 |
|
|
| 605 |
m<-as.integer(na.omit(m)) |
|
| 580 | 606 |
trans4_stations<-transect_from_spdf(data_stat,m) |
| 607 |
point4_stations<-data_stat[m,] |
|
| 581 | 608 |
#tmp<-as.data.frame(data_stat[1,]) |
| 582 | 609 |
#row.names(tmp)<-rep("X",1)
|
| 583 | 610 |
#test<-SpatialLinesDataFrame(trans4_stations,data=tmp) |
| ... | ... | |
| 597 | 624 |
list_transect2[[1]]<-c("t1_line.shp",paste("figure_13_tmax_elevation_transect1_OR",date_selected,out_prefix,sep="_"))
|
| 598 | 625 |
list_transect2[[2]]<-c("t2_line.shp",paste("figure_14_tmax_elevation_transect2_OR",date_selected,out_prefix,sep="_"))
|
| 599 | 626 |
list_transect2[[3]]<-c("t3_line.shp",paste("figure_15_tmax_elevation_transect3_OR",date_selected,out_prefix,sep="_"))
|
| 600 |
list_transect2[[4]]<-c("t4_line.shp",paste("figure_16_tmax_elevation_transect3_OR",date_selected,out_prefix,sep="_"))
|
|
| 627 |
list_transect2[[4]]<-c("t4_line.shp",paste("figure_16_tmax_elevation_transect4_OR",date_selected,out_prefix,sep="_"))
|
|
| 601 | 628 |
|
| 602 | 629 |
names(list_transect2)<-c("transect_OR1","transect_OR2","transect_OR3","transect_OR4")
|
| 603 | 630 |
|
| ... | ... | |
| 634 | 661 |
trans_data2<-plot_transect_m(list_transect2,rast_pred2,title_plot2,disp=TRUE,m_layers_sc) |
| 635 | 662 |
dev.off() |
| 636 | 663 |
|
| 664 |
### PLOT LOCATIONS OF STATION ON FIGURES |
|
| 637 | 665 |
|
| 638 |
X11(width=18,height=9) |
|
| 639 |
trans_elev<-vector("list",nb_transect)
|
|
| 640 |
for (k in 1:nb_transect){
|
|
| 641 |
|
|
| 642 |
trans_file<-list_transect[[k]] |
|
| 643 |
filename<-sub(".shp","",trans_file) #Removing the extension from file.
|
|
| 644 |
transect<-readOGR(".", filename) #reading shapefile
|
|
| 645 |
trans_elev[[k]]<-extract(ELEV_SRTM,transect) |
|
| 646 |
y<-as.numeric(trans_elev[[k]][[1]]) |
|
| 647 |
elev_y<-y |
|
| 648 |
x<-1:length(y) |
|
| 649 |
plot(x,y,type="l", ylab="Elevation (in meters)",xlab="Transect position (in km)") |
|
| 650 |
data_y<-(trans_data[[k]][[1]]) # data for the first transect |
|
| 651 |
#as.data.frame(data_y) |
|
| 652 |
par(new=TRUE) # key: ask for new plot without erasing old |
|
| 653 |
y<-data_y[,1] |
|
| 654 |
x <- 1:length(y) |
|
| 655 |
fus_y<-y |
|
| 656 |
plot(x,y,type="l",col="red",axes=F) #plotting fusion profile |
|
| 657 |
axis(4,xlab="",ylab="tmax (in degree C)") |
|
| 658 |
y<-data_y[,2] |
|
| 659 |
cai_y<-y |
|
| 660 |
lines(x,y,col="green") |
|
| 661 |
|
|
| 662 |
#title(title_plot[i])) |
|
| 663 |
legend("topleft",legend=c("elev","fus","CAI"),
|
|
| 664 |
cex=1.2, col=c("black","red","green"),
|
|
| 666 |
data_stat<-rbind(data_vf[,c("id","training")],data_sf[,c("id","training")]) #bringing together data_v and data_s
|
|
| 667 |
m<-match(selected_stations,data_stat$id) |
|
| 668 |
m<-as.integer(na.omit(m)) |
|
| 669 |
trans4_stations<-transect_from_spdf(data_stat,m) |
|
| 670 |
point4_stations<-data_stat[m,] |
|
| 671 |
|
|
| 672 |
pos<-match(c("x_OR83M","y_OR83M"),layerNames(s_raster)) #Find column with name "value"
|
|
| 673 |
xy_stack<-subset(s_raster,pos) #Select multiple layers from the stack |
|
| 674 |
r_stack<-stack(xy_stack, rast_pred2) |
|
| 675 |
trans4_data<-extract(r_stack,trans4_stations,cellnumbers=TRUE) #This extracts a list |
|
| 676 |
trans4_data<-as.data.frame(trans4_data[[1]]) |
|
| 677 |
point4_cellID<-cellFromXY(r_stack,coordinates(point4_stations)) #This contains the cell ID the points |
|
| 678 |
pos<-match(point4_cellID,trans4_data$cell) |
|
| 679 |
|
|
| 680 |
#Plots lines where there are stations... |
|
| 681 |
X11(width=18,height=9) |
|
| 682 |
y<-trans4_data$fus |
|
| 683 |
x <- 1:length(y) |
|
| 684 |
plot(x,y,type="l",col="red", #plotting fusion profile |
|
| 685 |
,xlab="",ylab="tmax (in degree C)") |
|
| 686 |
y<-trans4_data$CAI |
|
| 687 |
lines(x,y,col="green") |
|
| 688 |
abline(v=pos)#addlines whtere the stations area... |
|
| 689 |
#plot(elev) |
|
| 690 |
#title(title_plot[i])) |
|
| 691 |
legend("topleft",legend=c("fus","CAI"),
|
|
| 692 |
cex=1.2, col=c("red","green"),
|
|
| 693 |
lty=1) |
|
| 694 |
savePlot(paste("fig17_transect_path_tmax_diff_CAI_fusion_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 695 |
|
|
| 696 |
|
|
| 697 |
y<-trans4_data$fus[1:150] |
|
| 698 |
x <- 1:150 |
|
| 699 |
plot(x,y,type="l",col="red", #plotting fusion profile |
|
| 700 |
,xlab="",ylab="tmax (in degree C)") |
|
| 701 |
y<-trans4_data$CAI[1:150] |
|
| 702 |
lines(x,y,col="green") |
|
| 703 |
abline(v=pos)#addlines whtere the stations area... |
|
| 704 |
#plot(elev) |
|
| 705 |
#title(title_plot[i])) |
|
| 706 |
legend("topleft",legend=c("fus","CAI"),
|
|
| 707 |
cex=1.2, col=c("red","green"),
|
|
| 665 | 708 |
lty=1) |
| 666 |
savePlot(file=paste(list_transect[[k]][2],".png",sep=""),type="png") |
|
| 667 |
|
|
| 668 |
cor(fus_y,elev_y) |
|
| 669 |
cor(cai_y,elev_y) |
|
| 670 |
cor(fus_y,cai_y) |
|
| 671 |
} |
|
| 672 |
dev.off |
|
| 709 |
savePlot(paste("fig18a_transect_path_tmax_diff_CAI_fusion_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 710 |
|
|
| 711 |
|
|
| 712 |
y<-trans4_data$fus[151:300] |
|
| 713 |
x <- 151:300 |
|
| 714 |
plot(x,y,type="l",col="red", #plotting fusion profile |
|
| 715 |
,xlab="",ylab="tmax (in degree C)") |
|
| 716 |
y<-trans4_data$CAI[151:300] |
|
| 717 |
lines(x,y,col="green") |
|
| 718 |
abline(v=pos)#addlines whtere the stations area... |
|
| 719 |
#plot(elev) |
|
| 720 |
#title(title_plot[i])) |
|
| 721 |
legend("topleft",legend=c("fus","CAI"),
|
|
| 722 |
cex=1.2, col=c("red","green"),
|
|
| 723 |
lty=1) |
|
| 724 |
savePlot(paste("fig18b_transect_path_tmax_diff_CAI_fusion_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 725 |
|
|
| 726 |
dev.off() |
|
| 727 |
|
|
| 728 |
#cor(fus_y,elev_y) |
|
| 729 |
#cor(cai_y,elev_y) |
|
| 730 |
#cor(fus_y,cai_y) |
|
| 673 | 731 |
|
| 674 | 732 |
} |
Also available in: Unified diff
Methods comp part5-task#491- residuals analyses, spatial transect through stations with diff and elevation