Project

General

Profile

« Previous | Next » 

Revision 55056785

Added by Benoit Parmentier almost 12 years ago

GAM fusion function, added GAM models for bias surface and extraction of monthly mean tmax OR

View differences:

climate/research/oregon/interpolation/GAM_fusion_function_multisampling.R
8 8
  
9 9
  #Adding layer LST to the raster stack
10 10
  
11
  pos<-match("LST",layerNames(s_raster)) #Find column with name "LST"
12
  s_raster<-dropLayer(s_raster,pos) # If it exists drop layer
11 13
  pos<-match(LST_month,layerNames(s_raster)) #Find column with the current month for instance mm12
12 14
  r1<-raster(s_raster,layer=pos)             #Select layer from stack
13 15
  layerNames(r1)<-"LST"
14
  s_raster<-addLayer(s_raster,r1)            #Adding current month
16
  s_raster<-addLayer(s_raster,r1)            #Adding current month as "LST"
15 17
  
16 18
  ###Regression part 1: Creating a validation dataset by creating training and testing datasets
17 19
  
18 20
  mod_LST <-ghcn.subsets[[i]][,match(LST_month, names(ghcn.subsets[[i]]))]  #Match interpolation date and monthly LST average
19 21
  ghcn.subsets[[i]] = transform(ghcn.subsets[[i]],LST = mod_LST)            #Add the variable LST to the subset dataset
22
  dst$LST<-dst[[LST_month]]
20 23
  #n<-nrow(ghcn.subsets[[i]])
21 24
  #ns<-n-round(n*prop)   #Create a sample from the data frame with 70% of the rows
22 25
  #nv<-n-ns              #create a sample for validation with prop of the rows
......
48 51
  # STEP 2 - Weather station means across same days: Monthly mean calculation
49 52
  ###########
50 53
  
51
  modst=dst[dst$month==mo,] #Subsetting dataset for the relevant month of the date being processed
54
  modst<-dst[dst$month==mo,] #Subsetting dataset for the relevant month of the date being processed
52 55
  
53 56
  ##########
54
  # STEP 3 - get LST at stations
57
  # STEP 3 - get LST at stations  
55 58
  ##########
56 59
  
57 60
  sta_lola=modst[,c("lon","lat")] #Extracting locations of stations for the current month..
......
94 97
  names(d)[pos]<-c("id")
95 98
  names(x)[pos]<-c("id")
96 99
  names(modst)[1]<-c("id")       #modst contains the average tmax per month for every stations...
97
  dmoday=merge(modst,d,by="id")  #LOOSING DATA HERE!!! from 113 t0 103
98
  xmoday=merge(modst,x,by="id")  #LOOSING DATA HERE!!! from 48 t0 43
99
  names(dmoday)[4]<-c("lat")
100
  names(dmoday)[5]<-c("lon")     #dmoday contains all the the information: BIAS, monn
101
  names(xmoday)[4]<-c("lat")
102
  names(xmoday)[5]<-c("lon")     #dmoday contains all the the information: BIAS, monn
100
  
101
  dmoday=merge(modst,d,by="id",suffixes=c("",".y2"))  #LOOSING DATA HERE!!! from 113 t0 103
102
  xmoday=merge(modst,x,by="id",suffixes=c("",".y2"))  #LOOSING DATA HERE!!! from 48 t0 43
103
  mod_pat<-glob2rx("*.y2")   
104
  var_pat<-grep(mod_pat,names(dmoday),value=FALSE) # using grep with "value" extracts the matching names
105
  dmoday<-dmoday[,-var_pat]
106
  mod_pat<-glob2rx("*.y2")   
107
  var_pat<-grep(mod_pat,names(xmoday),value=FALSE) # using grep with "value" extracts the matching names
108
  xmoday<-xmoday[,-var_pat] #Removing duplicate columns
103 109
  
104 110
  data_v<-xmoday
105 111
  ###
......
107 113
  #dmoday contains the daily tmax values for training with TMax being the monthly station tmax mean
108 114
  #xmoday contains the daily tmax values for validation with TMax being the monthly station tmax mean
109 115
  
110
  # windows()
111
  #png(paste("LST_TMax_scatterplot_",sampling_dat$date[i],out_prefix,".png", sep=""))
112 116
  png(paste("Daily_tmax_monthly_TMax_scatterplot_",sampling_dat$date[i],"_",sampling_dat$prop[i],
113 117
            "_",sampling_dat$run_samp[i],out_prefix,".png", sep=""))
114 118
  plot(dailyTmax~TMax,data=dmoday,xlab="Mo Tmax",ylab=paste("Daily for",datelabel),main="across stations in OR")
115
  #savePlot(paste("Daily_tmax_monthly_TMax_scatterplot_",sampling_dat$date[i],out_prefix,".png", sep=""), type="png")
116
  #png(paste("LST_TMax_scatterplot_",sampling_dat$date[i],out_prefix,".png", sep=""))
117 119
  dev.off()
118 120
  
119 121
  ########
120 122
  # STEP 5 - interpolate bias
121 123
  ########
122 124
  
123
  # ?? include covariates like elev, distance to coast, cloud frequency, tree height
124
  #library(fields)
125
  #windows()
126
  quilt.plot(sta_lola,sta_bias,main="Bias at stations",asp=1)
127
  US(add=T,col="magenta",lwd=2)
125
  # ?? include covariates like elev, distance to coast, cloud frequency, tree heig
128 126
  #fitbias<-Tps(bias_xy,sta_bias) #use TPS or krige
129 127
  
130
  #Adding options to use only training stations: 07/11/2012
128
  #Adding options to use only training stations : 07/11/2012
131 129
  bias_xy=project(as.matrix(sta_lola),proj_str)
132 130
  #bias_xy2=project(as.matrix(c(dmoday$lon,dmoday$lat),proj_str)
133 131
  if(bias_val==1){
......
138 136
  fitbias<-Krig(bias_xy,sta_bias,theta=1e5) #use TPS or krige 
139 137
  #The output is a krig object using fields
140 138
  mod9a<-fitbias
141
  # Creating plot of bias surface and saving it
142
  #X11()
143
  png(paste("Bias_surface_LST_TMax_",sampling_dat$date[i],"_",sampling_dat$prop[i],
144
            "_",sampling_dat$run_samp[i],out_prefix,".png", sep="")) #Create file to write a plot
145
  datelabel2=format(ISOdate(year,mo,day),"%B ") #added by Benoit, label
146
  surface(fitbias,col=rev(terrain.colors(100)),asp=1,main=paste("Interpolated bias for",datelabel2,sep=" ")) #Plot to file
147
  #savePlot(paste("Bias_surface_LST_TMax_",sampling_dat$date[i],out_prefix,".png", sep=""), type="png")
148
  dev.off()  #Release the hold to the file
149
  
150
  #US(add=T,col="magenta",lwd=2)
151
  
139
    
152 140
  ##########
153 141
  # STEP 7 - interpolate delta across space
154 142
  ##########
......
156 144
  daily_sta_lola=dmoday[,c("lon","lat")] #could be same as before but why assume merge does this - assume not
157 145
  daily_sta_xy=project(as.matrix(daily_sta_lola),proj_str)
158 146
  daily_delta=dmoday$dailyTmax-dmoday$TMax
159
  #windows()
160
  quilt.plot(daily_sta_lola,daily_delta,asp=1,main="Station delta for Jan 15")
161
  US(add=T,col="magenta",lwd=2)
147
  
162 148
  #fitdelta<-Tps(daily_sta_xy,daily_delta) #use TPS or krige
163 149
  fitdelta<-Krig(daily_sta_xy,daily_delta,theta=1e5) #use TPS or krige
164 150
  #Kriging using fields package
165 151
  mod9b<-fitdelta
166
  # Creating plot of bias surface and saving it
167
  #X11()
152

  
168 153
  png(paste("Delta_surface_LST_TMax_",sampling_dat$date[i],"_",sampling_dat$prop[i],
169 154
            "_",sampling_dat$run_samp[i],out_prefix,".png", sep=""))
170 155
  surface(fitdelta,col=rev(terrain.colors(100)),asp=1,main=paste("Interpolated delta for",datelabel,sep=" "))
171
  #savePlot(paste("Delta_surface_LST_TMax_",sampling_dat$date[i],out_prefix,".png", sep=""), type="png")
172 156
  dev.off()
173
  #US(add=T,col="magenta",lwd=2)
174 157
  #
175 158
  
176 159
  #### Added by Benoit on 06/19
......
237 220
  mae<-MAE_fun(sta_pred,tmax)
238 221
  r2<-cor(sta_pred,tmax)^2              #R2, coef. of var
239 222
  me<-mean(sta_pred-tmax)
240
 
241
  #plot(sta_pred~dmoday$dailyTmax,xlab=paste("Actual daily for",datelabel),ylab="Pred daily",main=paste("RMSE=",rmse))
242
  
223
   
243 224
  png(paste("Predicted_tmax_versus_observed_scatterplot_",sampling_dat$date[i],"_",sampling_dat$prop[i],
244 225
            "_",sampling_dat$run_samp[i],out_prefix,".png", sep=""))
245 226
  plot(sta_pred~tmax,xlab=paste("Actual daily for",datelabel),ylab="Pred daily",main=paste("RMSE=",rmse))
246 227
  abline(0,1)
247
  #savePlot(paste("Predicted_tmax_versus_observed_scatterplot_",sampling_dat$date[i],out_prefix,".png", sep=""), type="png")
248 228
  dev.off()
249 229
  #resid=sta_pred-dmoday$dailyTmax
250 230
  resid=sta_pred-tmax
251
  #quilt.plot(daily_sta_lola,resid)
252
  
253
  ### END OF BRIAN's code
254
  
255
  ### Added by benoit
256 231
  
257 232
  ###BEFORE GAM prediction the data object must be transformed to SDF
258 233
  
......
262 237
  coords<- data_s[,c('x_OR83M','y_OR83M')]
263 238
  coordinates(data_s)<-coords
264 239
  proj4string(data_s)<-CRS  #Need to assign coordinates..
240
  coords<- modst[,c('x_OR83M','y_OR83M')]
241
  coordinates(modst)<-coords
242
  proj4string(modst)<-CRS  #Need to assign coordinates..
265 243
  
266 244
  ns<-nrow(data_s) #This is added to because some loss of data might have happened because of the averaging...
267 245
  nv<-nrow(data_v)
268 246
  
269 247
  ###GAM PREDICTION
270 248
  
271
  #data_s$y_var<-data_s$dailyTmax  #This shoudl be changed for any variable!!!
272
  #data_v$y_var<-data_v$dailyTmax
273
  data_v$y_var<-data_v[[y_var_name]]
274
  data_s$y_var<-data_s[[y_var_name]]
249
  if (bias_prediction==1){
250
    data_s$y_var<-data_s$LSTD_bias  #This shoudl be changed for any variable!!!
251
    data_v$y_var<-data_v$LSTD_bias
252
    data_month<-modst
253
    data_month$y_var<-modst$LSTD_bias
254
  }
255

  
256
  if (bias_prediction==0){
257
    data_v$y_var<-data_v[[y_var_name]]
258
    data_s$y_var<-data_s[[y_var_name]]
259
  }
275 260
  
276 261
  #Model and response variable can be changed without affecting the script
277 262
  
......
282 267
  formula5 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)", env=.GlobalEnv)
283 268
  formula6 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC1)", env=.GlobalEnv)
284 269
  formula7 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC3)", env=.GlobalEnv)
285
  formula8 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(LST,LC1)", env=.GlobalEnv)
286
  
287
  #formula1 <- as.formula("y_var ~ s(lat,lon,ELEV_SRTM)", env=.GlobalEnv)
288
  #formula2 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM)  + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(CANHEIGHT)", env=.GlobalEnv)
289
  #formula3 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST,CANHEIGHT)", env=.GlobalEnv)
290
  #formula4 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST,LC1)", env=.GlobalEnv)
291
  #formula5 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST,LC3)", env=.GlobalEnv)
292
  #formula6 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC1)", env=.GlobalEnv)
293
  #formula7 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC3)", env=.GlobalEnv)
294
  #formula8 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST,LC1) + s(LST,LC3)", env=.GlobalEnv)
295
  
296
  mod1<- try(gam(formula1, data=data_s))
297
  mod2<- try(gam(formula2, data=data_s)) #modified nesting....from 3 to 2
298
  mod3<- try(gam(formula3, data=data_s))
299
  mod4<- try(gam(formula4, data=data_s))
300
  mod5<- try(gam(formula5, data=data_s))
301
  mod6<- try(gam(formula6, data=data_s))
302
  mod7<- try(gam(formula7, data=data_s))
303
  mod8<- try(gam(formula8, data=data_s))
304
  
305
#   mod1<- try(gam(formula1, data=data_s))
306
#   mod2<- try(gam(formula2, data=data_s)) #modified nesting....from 3 to 2
307
#   mod3<- try(gam(y_var~ s(lat) + s (lon) + s (ELEV_SRTM) +  s (Northness)+ s (Eastness) + s(DISTOC), data=data_s))
308
#   mod4<- try(gam(y_var~ s(lat) + s (lon) + s(ELEV_SRTM) + s(Northness) + s (Eastness) + s(DISTOC) + s(LST), data=data_s))
309
#   mod5<- try(gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST), data=data_s))
310
#   mod6<- try(gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC1), data=data_s))
311
#   mod7<- try(gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC3), data=data_s))
312
#   mod8<- try(gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(LC1,LC3), data=data_s))
313
#   
270
  formula8 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(LC1,LC3)", env=.GlobalEnv)
271
  
272
  if (bias_prediction==1){
273
    mod1<- try(gam(formula1, data=data_month))
274
    mod2<- try(gam(formula2, data=data_month)) #modified nesting....from 3 to 2
275
    mod3<- try(gam(formula3, data=data_month))
276
    mod4<- try(gam(formula4, data=data_month))
277
    mod5<- try(gam(formula5, data=data_month))
278
    mod6<- try(gam(formula6, data=data_month))
279
    mod7<- try(gam(formula7, data=data_month))
280
    mod8<- try(gam(formula8, data=data_month)) 
281
    
282
  } else if (bias_prediction==0){
283
    
284
    mod1<- try(gam(formula1, data=data_s))
285
    mod2<- try(gam(formula2, data=data_s)) #modified nesting....from 3 to 2
286
    mod3<- try(gam(formula3, data=data_s))
287
    mod4<- try(gam(formula4, data=data_s))
288
    mod5<- try(gam(formula5, data=data_s))
289
    mod6<- try(gam(formula6, data=data_s))
290
    mod7<- try(gam(formula7, data=data_s))
291
    mod8<- try(gam(formula8, data=data_s))
292
  }
293

  
314 294
  #Added
315 295
  #tmax_predicted=themolst+daily_delta_rast-bias_rast #Final surface?? but daily_rst
316 296
  
......
351 331
  #ns<-nrow(data_s) #This is added to because some loss of data might have happened because of the averaging...
352 332
  #nv<-nrow(data_v)
353 333
  
334
  pred_mod<-paste("pred_mod",j,sep="")
335
  #Adding the results back into the original dataframes.
336
  data_s[[pred_mod]]<-sta_pred_data_s
337
  data_v[[pred_mod]]<-sta_pred 
338
  
339
  #Model assessment: RMSE and then krig the residuals....!
340
  
341
  res_mod_s<- data_s$dailyTmax - data_s[[pred_mod]]           #Residuals from kriging training
342
  res_mod_v<- data_v$dailyTmax - data_v[[pred_mod]]           #Residuals from kriging validation
343
  
344
  name2<-paste("res_mod",j,sep="")
345
  data_v[[name2]]<-as.numeric(res_mod_v)
346
  data_s[[name2]]<-as.numeric(res_mod_s)
347
  
354 348
  
355 349
  for (j in 1:nmodels){
356 350
    
......
453 447
        s_sgdf<-as(s_raster,"SpatialGridDataFrame") #Conversion to spatial grid data frame
454 448
        
455 449
        rpred<- predict(mod, newdata=s_sgdf, se.fit = TRUE) #Using the coeff to predict new values.
456
        y_pred<-rpred$fit
450
        y_pred<-rpred$fit #rpred is a list with fit being and array
457 451
        raster_pred<-r1
458 452
        layerNames(raster_pred)<-"y_pred"
459 453
        values(raster_pred)<-as.numeric(y_pred)
460
        data_name<-paste("predicted_mod",j,"_",sampling_dat$date[i],"_",sampling_dat$prop[i],
461
                         "_",sampling_dat$run_samp[i],sep="")
462
        raster_name<-paste("GAM_",data_name,out_prefix,".rst", sep="")
463
        writeRaster(raster_pred, filename=raster_name,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)
464
        #writeRaster(r2, filename=raster_name,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)
454
        
455
        if (bias_prediction==1){
456
          data_name<-paste("predicted_mod",j,"_",sampling_dat$date[i],"_",sampling_dat$prop[i],
457
                           "_",sampling_dat$run_samp[i],sep="")
458
          raster_name<-paste("GAM_bias_",data_name,out_prefix,".rst", sep="")
459
          writeRaster(raster_pred, filename=raster_name,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)
460
          bias_rast<-raster_pred
461
          
462
          raster_pred=themolst+daily_delta_rast-bias_rast #Final surface  as a raster layer...
463
          layerNames(raster_pred)<-"y_pred"
464
          #=themolst+daily_delta_rast-bias_rast #Final surface  as a raster layer...
465
          
466
          data_name<-paste("predicted_mod",j,"_",sampling_dat$date[i],"_",sampling_dat$prop[i],
467
                           "_",sampling_dat$run_samp[i],sep="")
468
          raster_name<-paste("GAM_bias_tmax_",data_name,out_prefix,".rst", sep="")
469
          writeRaster(raster_pred, filename=raster_name,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)
470
          
471
        }
472
        
473
        if (bias_prediction==0){
474
          data_name<-paste("predicted_mod",j,"_",sampling_dat$date[i],"_",sampling_dat$prop[i],
475
                           "_",sampling_dat$run_samp[i],sep="")
476
          raster_name<-paste("GAM_",data_name,out_prefix,".rst", sep="")
477
          writeRaster(raster_pred, filename=raster_name,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)
478
          #writeRaster(r2, filename=raster_name,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)
479
          
480
        }
481
        
465 482
        
466 483
        pred_sgdf<-as(raster_pred,"SpatialGridDataFrame") #Conversion to spatial grid data frame
467 484
        #rpred_val_s <- overlay(raster_pred,data_s)             #This overlays the kriged surface tmax and the location of weather stations
468 485
        
469
        rpred_val_s <- overlay(pred_sgdf,data_s)             #This overlays the kriged surface tmax and the location of weather stations
470
        rpred_val_v <- overlay(pred_sgdf,data_v)             #This overlays the kriged surface tmax and the location of weather stations
486
        rpred_val_s <- overlay(pred_sgdf,data_s)             #This overlays the interpolated surface tmax and the location of weather stations
487
        rpred_val_v <- overlay(pred_sgdf,data_v)             #This overlays the interpolated surface tmax and the location of weather stations
471 488
        
472 489
        pred_mod<-paste("pred_mod",j,sep="")
473 490
        #Adding the results back into the original dataframes.
......
477 494
        
478 495
        #Model assessment: RMSE and then krig the residuals....!
479 496
        
480
        res_mod_s<- data_s$y_var - data_s[[pred_mod]]           #Residuals from kriging training
481
        res_mod_v<- data_v$y_var - data_v[[pred_mod]]           #Residuals from kriging validation
497
        res_mod_s<-data_s[[y_var_name]] - data_s[[pred_mod]] #residuals from modeling training
498
        res_mod_v<-data_v[[y_var_name]] - data_v[[pred_mod]] #residuals from modeling validation
482 499
        
483 500
      }
484 501
      
......
493 510
        
494 511
        #Model assessment: RMSE and then krig the residuals....!
495 512
        
496
        res_mod_s<- data_s$y_var - data_s[[pred_mod]]           #Residuals from kriging training
497
        res_mod_v<- data_v$y_var - data_v[[pred_mod]]           #Residuals from kriging validation
513
        #res_mod_s<- data_s$y_var - data_s[[pred_mod]]           #Residuals from modeling training
514
        #res_mod_v<- data_v$y_var - data_v[[pred_mod]]           #Residuals from modeling validation
515
        res_mod_s<-data_s[[y_var_name]] - data_s[[pred_mod]]
516
        res_mod_v<-data_v[[y_var_name]] - data_v[[pred_mod]]
517
        
498 518
      }
499 519

  
500 520
      ####ADDED ON JULY 20th
......
568 588
  mod_obj<-list(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,mod9a,mod9b)
569 589
  names(mod_obj)<-c("mod1","mod2","mod3","mod4","mod5","mod6","mod7","mod8","mod9a","mod9b")
570 590
  #results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2)
571
  results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj)
572
  names(results_list)<-c("data_s","data_v","tb_metrics1","tb_metrics2","mod_obj")
591
  results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj,sampling_dat[i,],data_month)
592
  names(results_list)<-c("data_s","data_v","tb_metrics1","tb_metrics2","mod_obj","sampling_dat","data_month")
573 593
  save(results_list,file= paste(path,"/","results_list_metrics_objects_",sampling_dat$date[i],"_",sampling_dat$prop[i],
574 594
                                "_",sampling_dat$run_samp[i],out_prefix,".RData",sep=""))
575 595
  return(results_list)

Also available in: Unified diff