Project

General

Profile

« Previous | Next » 

Revision 4180b123

Added by Benoit Parmentier about 12 years ago

  • ID 4180b1232d716793def3ae75a4e38c2c15b5e4e9

Transformed the aspect variable and added GAM model closer to PRISM

View differences:

climate/research/oregon/interpolation/Linear_reg.R
11 11 
library(mgcv)
12 12 

  		




  13
  13
  
    
###Parameters and arguments


  		




  14
  
  
    
infile1<-"ghcn_or_b_02122012_OR83M.csv"


  		




  
  14
  
    
#infile1<-"ghcn_or_b_02122012_OR83M.csv"


  		




  
  15
  
    
infile1<-"ghcn_or_tmax_b_03032012_OR83M.csv"


  		




  15
  16
  
    
path<-"C:/Data/Benoit/NCEAS/window_Oregon_data"


  		




  16
  17
  
    
setwd(path)


  		




  17
  18
  
    
#infile2<-"dates_interpolation_03012012.txt"  # list of 10 dates for the regression


  		




  18
  19
  
    
infile2<-"dates_interpolation_03032012.txt"


  		




  19
  20
  
    
prop<-0.3                                                                            #Proportion of testing retained for validation   


  		




  20
  
  
    
out_prefix<-"_03022012_r3"


  		




  
  21
  
    
out_prefix<-"_03042012_r1"


  		




  21
  22
  
    

  		




  22
  23
  
    
#######START OF THE SCRIPT #############


  		




  23
  24
  
    

  		




  24
  25
  
    
###Reading the station data and setting up for models' comparison


  		




  25
  26
  
    
ghcn<-read.csv(paste(path,"/",infile1, sep=""), header=TRUE)                            #The "paste" function concatenates the path and file name in common string. 


  		




  
  27
  
    
ghcn = transform(ghcn,Northness_w = sin(slope)*cos(ASPECT)) #Adding a variable to the dataframe


  		




  
  28
  
    
ghcn = transform(ghcn,Eastness_w = sin(slope)*sin(ASPECT))  #adding variable to the dataframe.


  		




  
  29
  
    
set.seed(100)


  		




  26
  30
  
    
dates <-readLines(paste(path,"/",infile2, sep=""))


  		




  27
  31
  
    

  		




  28
  32
  
    
results <- matrix(1,length(dates),10)            #This is a matrix containing the diagnostic measures from the GAM models.


  		




  29
  
  
    
ghcn.subsets <-lapply(dates, function(d) subset(ghcn, date_==d)) #this creates a list of 10 subsets data


  		




  
  33
  
    
ghcn.subsets <-lapply(dates, function(d) subset(ghcn, date==d)) #this creates a list of 10 subsets data


  		




  
  34
  
    
#note that compare to the previous version date_ column was changed to date


  		




  30
  35
  
    

  		




  31
  36
  
    
## looping through the dates...


  		




  32
  37
  
    

  		




  ......




  47
  52
  
    

  		




  48
  53
  
    
  GAM_ANUSPLIN1<-gam(tmax~ s(lat) + s (lon) + s (ELEV_SRTM), data=data_s)


  		




  49
  54
  
    
  GAM_PRISM1<-gam(tmax~ s(lat) + s (lon) + s (ELEV_SRTM) + s (ASPECT)+ s(DISTOC), data=data_s)


  		




  
  55
  
    
  GAM_PRISM2<-gam(tmax~ s(lat) + s (lon) + s (ELEV_SRTM) + s (Northness_w)+ s (Eastness_w) + s(DISTOC), data=data_s)


  		




  
  56
  
    
  


  		




  50
  57
  
    
  


  		




  51
  58
  
    
  ####Regression part 3: Calculating and storing diagnostic measures


  		




  52
  59
  
    
  


  		




  ......




  102
  109
  
    
win.graph()


  		




  103
  110
  
    
barplot(results_table$AIC_A1,main="AIC for the A1 models",names.arg=results_table$dates,ylab="Temp (0.1 X deg. C)",xlab="interolated date")


  		




  104
  111
  
    
savePlot(paste("GAM_PRISM1_RMSE",out_prefix,".emf", sep=""), type="emf")


  		




  105
  
  
    

  		




  
  112
  
    
win.graph()


  		




  
  113
  
    
barplot(results_table$Deviance_A1,main="Deviance for the A1 models",names.arg=results_table$dates,ylab="Temp (0.1 X deg. C)",xlab="interolated date")


  		




  
  114
  
    
savePlot(paste("GAM_ANUSPLIN1_Deviance",out_prefix,".emf", sep=""), type="emf")


  		




  
  115
  
    
win.graph()


  		




  
  116
  
    
barplot(results_table$Deviance_P1,main="Deviance for the P1 models",names.arg=results_table$dates,ylab="Temp (0.1 X deg. C)",xlab="interolated date")


  		




  
  117
  
    
savePlot(paste("GAM_PRISM1_Deviance",out_prefix,".emf", sep=""), type="emf")


  		




  106
  118
  
    
write.csv(results_table, file= paste(path,"/","results_GAM_Assessment",out_prefix,".txt",sep=""))


  		




  107
  119
  
    

  		




  
  120
  
    

  		




  108
  121
  
    
# End of script##########


  		




  109
  122
  
    

  		




  110
  123
  
    
# ###############################


  		












Also available in:  Unified diff