Project

General

Profile

« Previous | Next » 

Revision e9998467

Added by Benoit Parmentier about 12 years ago

Methods comp part5-task#491- residuals analyses, spatial transect through stations with diff and elevation

View differences:

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