Project

General

Profile

« Previous | Next » 

Revision 8d91ca4a

Added by Benoit Parmentier over 12 years ago

GAM LST, added model and modified summary plots, run on 365 dates GAM+Kriging, task #364 and #406

View differences:

climate/research/oregon/interpolation/GAM_LST_reg.R
27 27 

  		




  28
  28
  
    
###Parameters and arguments


  		




  29
  29
  
    

  		




  30
  
  
    
infile1<-"ghcn_or_tmax_b_04142012_OR83M.shp"


  		




  31
  
  
    
#infile2<-"dates_interpolation_03052012_2dates_test.txt"


  		




  32
  
  
    
infile2<-"dates_interpolation_03052012.txt"                                          #List of 10 dates for the regression


  		




  
  30
  
    
infile1<- "ghcn_or_tmax_b_04142012_OR83M.shp"


  		




  
  31
  
    
#infile2<-"dates_interpolation_03052012.txt"                                          #List of 10 dates for the regression


  		




  
  32
  
    
infile2<-"list_365_dates_04212012.txt"


  		




  33
  33
  
    
infile3<-"LST_dates_var_names.txt"


  		




  34
  
  
    
infile4<-"models_interpolation_04032012.txt"


  		




  
  34
  
    
infile4<-"models_interpolation_05142012.txt"


  		




  35
  35
  
    
infile5<-"mean_day244_rescaled.rst" #Raster or grid for the locations of predictions


  		




  36
  36
  
    

  		




  37
  37
  
    
#path<-"/data/computer/parmentier/Data/IPLANT_project/data_Oregon_stations"         #Jupiter LOCATION on EOS


  		




  ......




  39
  39
  
    
setwd(path) 


  		




  40
  40
  
    
prop<-0.3                                                                           #Proportion of testing retained for validation   


  		




  41
  41
  
    
seed_number<-100


  		




  42
  
  
    
out_prefix<-"_05062012m_Kr_LST"


  		




  
  42
  
    
out_prefix<-"_05142012_365d_Kr_LST"


  		




  43
  43
  
    

  		




  44
  44
  
    
#######START OF THE SCRIPT #############


  		




  45
  45
  
    

  		




  ......




  115
  115
  
    
  mod5<- gam(tmax~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST), data=data_s)


  		




  116
  116
  
    
  mod6<- gam(tmax~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST,LC1), data=data_s)


  		




  117
  117
  
    
  mod7<- gam(tmax~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST,LC3), data=data_s)


  		




  
  118
  
    
  mod8<- gam(tmax~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(LC1), data=data_s)


  		




  118
  119
  
    
  


  		




  119
  120
  
    
  ####Regression part 3: Calculating and storing diagnostic measures


  		




  120
  121
  
    
  #listmod can be created and looped over. In this case we loop around the objects..


  		




  ......




  221
  222
  
    
  


  		




  222
  223
  
    
  data_name<-paste("ghcn_v_",out_prefix,"_",dates[[i]],sep="")


  		




  223
  224
  
    
  assign(data_name,data_v)


  		




  224
  
  
    
  write.table(data_v, file= paste(path,"/",data_name,".txt",sep=""), sep=" ")


  		




  
  225
  
    
  #write.table(data_v, file= paste(path,"/",data_name,".txt",sep=""), sep=" ")


  		




  225
  226
  
    
  #write out a new shapefile (including .prj component)


  		




  226
  227
  
    
  #outfile<-sub(".shp","",data_name)   #Removing extension if it is present


  		




  227
  228
  
    
  #writeOGR(data_v,".", outfile, driver ="ESRI Shapefile")


  		




  228
  229
  
    
  


  		




  229
  230
  
    
  data_name<-paste("ghcn_s_",out_prefix,"_",dates[[i]],sep="")


  		




  230
  231
  
    
  assign(data_name,data_s)


  		




  231
  
  
    
  write.table(data_s, file= paste(path,"/",data_name,".txt",sep=""), sep=" ")


  		




  
  232
  
    
  #write.table(data_s, file= paste(path,"/",data_name,".txt",sep=""), sep=" ")


  		




  232
  233
  
    
  #outfile<-sub(".shp","",data_name)   #Removing extension if it is present


  		




  233
  234
  
    
  #writeOGR(data_s,".", outfile, driver ="ESRI Shapefile")


  		




  234
  235
  
    
  


  		




  ......




  280
  281
  
    
  }


  		




  281
  282
  
    
  


  		




  282
  283
  
    
r1<-(results_table_RMSE[,3:10]) #selecting only the columns related to models and method 1


  		




  283
  
  
    
r2<-(results_table_RMSE[,3:10]) #selecting only the columns related to models and method 1


  		




  
  284
  
    
r2<-(results_table_RMSE_kr[,3:10]) #selecting only the columns related to models and method 1


  		




  284
  285
  
    
mean_r1<-mean(r1)


  		




  285
  286
  
    
mean_r2<-mean(r2)


  		




  286
  287
  
    
median_r1<-sapply(r1, median)   #Calulcating the mean for every model (median of columns)


  		




  ......




  288
  289
  
    
sd_r1<-sapply(r1, sd)


  		




  289
  290
  
    
sd_r2<-sapply(r2, sd)


  		




  290
  291
  
    

  		




  291
  
  
    
barplot(mean_r,ylim=c(23,26),ylab="RMSE in tenth deg C")


  		




  292
  
  
    
barplot(median_r,ylim=c(23,26),ylab="RMSE in tenth deg C",add=TRUE,inside=FALSE,beside=TRUE) # put both on the same plot


  		




  293
  
  
    
barplot(sd_r,ylim=c(6,8),ylab="RMSE in tenth deg C") # put both on the same plot


  		




  
  292
  
    
barplot(mean_r1,ylim=c(23,26),ylab="RMSE in tenth deg C")


  		




  
  293
  
    
barplot(mean_r2,ylim=c(23,26),ylab="RMSE in tenth deg C")


  		




  
  294
  
    
barplot(median_r1,ylim=c(23,26),ylab="RMSE in tenth deg C",add=TRUE,inside=FALSE,beside=TRUE) # put both on the same plot


  		




  
  295
  
    
barplot(median_r2,ylim=c(23,26),ylab="RMSE in tenth deg C",add=TRUE,inside=FALSE,beside=TRUE) # put both on the same plot


  		




  294
  296
  
    

  		




  295
  
  
    
height<-rbind(mean_r,median_r)


  		




  296
  
  
    
barplot(height,ylim=c(23,26),ylab="RMSE in tenth deg C",beside=TRUE,legend=rownames(height))


  		




  297
  
  
    
barplot(height,ylim=c(23,26),ylab="RMSE in tenth deg C",beside=TRUE, col=c("darkblue","red"),legend=rownames(height)) # put both on the same plot


  		




  
  297
  
    
barplot(sd_r1,ylim=c(6,8),ylab="RMSE in tenth deg C") # put both on the same plot


  		




  
  298
  
    
barplot(sd_r2,ylim=c(6,8),ylab="RMSE in tenth deg C") # put both on the same plot


  		




  298
  299
  
    

  		




  299
  
  
    
barplot2(mean_r,median_r,ylim=c(23,26),ylab="RMSE in tenth deg C") # put both on the same plot


  		




  
  300
  
    
height<-rbind(mean_r1,mean_r2)


  		




  
  301
  
    
barplot(height,ylim=c(20,26),ylab="RMSE in tenth deg C",beside=TRUE,legend=rownames(height))


  		




  
  302
  
    
barplot(height,ylim=c(20,26),ylab="RMSE in tenth deg C",beside=TRUE, col=c("darkblue","red"),legend=rownames(height)) # put both on the same plot


  		




  
  303
  
    

  		




  
  304
  
    
height<-rbind(median_r1,median_r2)


  		




  
  305
  
    
barplot(height,ylim=c(20,26),ylab="RMSE in tenth deg C",beside=TRUE,legend=rownames(height))


  		




  
  306
  
    
barplot(height,ylim=c(20,26),ylab="RMSE in tenth deg C",beside=TRUE, col=c("darkblue","red"),legend=rownames(height)) # put both on the same plot


  		




  
  307
  
    

  		




  
  308
  
    
height<-rbind(mean_r2,median_r2)


  		




  
  309
  
    
barplot(height,ylim=c(20,26),ylab="RMSE in tenth deg C",beside=TRUE,legend=rownames(height))


  		




  
  310
  
    
barplot(height,ylim=c(20,26),ylab="RMSE in tenth deg C",beside=TRUE, col=c("darkblue","red"),legend=rownames(height)) # put both on the same plot


  		




  
  311
  
    

  		




  
  312
  
    
#barplot2(mean_r,median_r,ylim=c(23,26),ylab="RMSE in tenth deg C") # put both on the same plot


  		




  300
  313
  
    
#Collect var explained and p values for each var...


  		




  
  314
  
    

  		




  301
  315
  
    
### End of script  ##########


  		




  302
  316
  
    

  		




  303
  317
  
    

  		












Also available in:  Unified diff