Revision 449d85fa
Added by Benoit Parmentier over 12 years ago
| climate/research/oregon/interpolation/GWR_prediction_reg.R | ||
|---|---|---|
| 1 |
################## Interpolation of Tmax Using Kriging ####################################### |
|
| 2 |
########################### Kriging and Cokriging ############################################### |
|
| 3 |
#This script interpolates station values for the Oregon case study using Kriging and Cokring. # |
|
| 4 |
#The script uses LST monthly averages as input variables and loads the station data # |
|
| 5 |
#from a shape file with projection information. # |
|
| 6 |
#Note that this program: # |
|
| 7 |
#1)assumes that the shape file is in the current working. # |
|
| 8 |
#2)relevant variables were extracted from raster images before performing the regressions # |
|
| 9 |
# and stored shapefile # |
|
| 10 |
#This scripts predicts tmax using autokrige, gstat and LST derived from MOD11A1. # |
|
| 11 |
#also included and assessed using the RMSE,MAE,ME and R2 from validation dataset. # |
|
| 12 |
#TThe dates must be provided as a textfile. # |
|
| 13 |
#AUTHOR: Benoit Parmentier # |
|
| 14 |
#DATE: 07/26/2012 # |
|
| 15 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#364-- # |
|
| 16 |
################################################################################################## |
|
| 1 |
################## Interpolation of Tmax Using GWR ######################################## |
|
| 2 |
########################### GWR WITH LST ####################################################### |
|
| 3 |
#This script interpolates station values for the Oregon case study using Geographically Weighted. # |
|
| 4 |
#The script uses LST monthly averages as input variables and loads the station data # |
|
| 5 |
#from a shape file with projection information. # |
|
| 6 |
#Note that this program: # |
|
| 7 |
#1)assumes that the shape file is in the current working folder # |
|
| 8 |
#2)relevant variables were extracted from raster images before performing the regressions # |
|
| 9 |
# and stored shapefile # |
|
| 10 |
#3)covariate raster images are present in the current working folder # |
|
| 11 |
#This scripts predicts tmax using autokrige, gstat and LST derived from MOD11A1. # |
|
| 12 |
#also included and assessed using the RMSE,MAE,ME and R2 from validation dataset. # |
|
| 13 |
#TThe dates must be provided as a textfile. # |
|
| 14 |
#AUTHOR: Benoit Parmentier # |
|
| 15 |
#DATE: 08/15/2012 # |
|
| 16 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#364-- # |
|
| 17 |
################################################################################################### |
|
| 17 | 18 |
|
| 18 | 19 |
###Loading R library and packages |
| 19 | 20 |
#library(gtools) # loading some useful tools |
| ... | ... | |
| 59 | 60 |
prop<-0.3 #Proportion of testing retained for validation |
| 60 | 61 |
#prop<-0.25 |
| 61 | 62 |
seed_number<- 100 #Seed number for random sampling |
| 62 |
out_prefix<-"test2_07312012_365d_gwr" #User defined output prefix
|
|
| 63 |
out_prefix<-"_08152012_1d_gwr4" #User defined output prefix
|
|
| 63 | 64 |
setwd(path) |
| 64 | 65 |
|
| 65 | 66 |
#source("fusion_function_07192012.R")
|
| 66 | 67 |
#source("KrigingUK_function_07262012.R")
|
| 67 |
source("GWR_function_07312012.R")
|
|
| 68 |
source("GWR_function_08152012b.R")
|
|
| 68 | 69 |
############ START OF THE SCRIPT ################## |
| 69 | 70 |
|
| 70 | 71 |
###Reading the station data and setting up for models' comparison |
| ... | ... | |
| 119 | 120 |
LC3<-raster(s_raster,layer=pos) #Select layer from stack |
| 120 | 121 |
s_raster<-dropLayer(s_raster,pos) |
| 121 | 122 |
LC3[is.na(LC3)]<-0 |
| 123 |
pos<-match("CANHEIGHT",layerNames(s_raster)) #Find column with name "value"
|
|
| 124 |
CANHEIGHT<-raster(s_raster,layer=pos) #Select layer from stack |
|
| 125 |
s_raster<-dropLayer(s_raster,pos) |
|
| 126 |
CANHEIGHT[is.na(CANHEIGHT)]<-0 |
|
| 122 | 127 |
|
| 123 | 128 |
xy<-coordinates(r1) #get x and y projected coordinates... |
| 124 | 129 |
xy_latlon<-project(xy, CRS, inv=TRUE) # find lat long for projected coordinats (or pixels...) |
| ... | ... | |
| 131 | 136 |
values(lon)<-xy_latlon[,1] |
| 132 | 137 |
values(lat)<-xy_latlon[,2] |
| 133 | 138 |
|
| 134 |
r<-stack(N,E,Nw,Ew,lon,lat,LC1,LC3) |
|
| 135 |
rnames<-c("Northness","Eastness","Northness_w","Eastness_w", "lon","lat","LC1","LC3")
|
|
| 139 |
r<-stack(N,E,Nw,Ew,lon,lat,LC1,LC3,CANHEIGHT)
|
|
| 140 |
rnames<-c("Northness","Eastness","Northness_w","Eastness_w", "lon","lat","LC1","LC3","CANHEIGHT")
|
|
| 136 | 141 |
layerNames(r)<-rnames |
| 137 | 142 |
s_raster<-addLayer(s_raster, r) |
| 138 |
|
|
| 139 | 143 |
#s_sgdf<-as(s_raster,"SpatialGridDataFrame") #Conversion to spatial grid data frame |
| 140 | 144 |
|
| 141 | 145 |
####### Preparing LST stack of climatology... |
| ... | ... | |
| 222 | 226 |
|
| 223 | 227 |
######## Prediction for the range of dates |
| 224 | 228 |
|
| 225 |
#gwr_mod<-mclapply(1:length(dates), runGWR,mc.preschedule=FALSE,mc.cores = 8) #This is the end bracket from mclapply(...) statement
|
|
| 229 |
gwr_mod<-mclapply(1:length(dates), runGWR,mc.preschedule=FALSE,mc.cores = 8) #This is the end bracket from mclapply(...) statement |
|
| 226 | 230 |
#fusion_mod357<-mclapply(357:365,runFusion, mc.cores=8)# for debugging |
| 227 | 231 |
#test<-runKriging(1) |
| 228 | 232 |
#test357<-mclapply(357:365,runFusion, mc.cores=8)# for debugging |
| 229 |
gwr_mod<-mclapply(1:1, runGWR,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement |
|
| 233 |
#gwr_mod<-mclapply(1:1, runGWR,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement
|
|
| 230 | 234 |
|
| 231 | 235 |
#test<-mclapply(357,runFusion, mc.cores=1)# for debugging |
| 232 | 236 |
|
| ... | ... | |
| 239 | 243 |
} |
| 240 | 244 |
rm(tb_tmp) |
| 241 | 245 |
|
| 242 |
for(i in 4:nmodels+3){ # start of the for loop #1
|
|
| 246 |
for(i in 4:(nmodels+3)){ # start of the for loop #1
|
|
| 243 | 247 |
tb[,i]<-as.numeric(as.character(tb[,i])) |
| 244 | 248 |
} |
| 245 |
|
|
| 246 | 249 |
tb_RMSE<-subset(tb, metric=="RMSE") |
| 247 | 250 |
tb_MAE<-subset(tb,metric=="MAE") |
| 248 | 251 |
tb_ME<-subset(tb,metric=="ME") |
| ... | ... | |
| 261 | 264 |
mean_RMSE_f<-sapply(tb_RMSE_f[,4:(nmodels+3)],mean) |
| 262 | 265 |
|
| 263 | 266 |
#Wrting metric results in textfile and model objects in .RData file |
| 264 |
write.table(tb_diagnostic1, file= paste(path,"/","results2_kriging_Assessment_measure1",out_prefix,".txt",sep=""), sep=",")
|
|
| 265 |
write.table(tb, file= paste(path,"/","results2_kriging_Assessment_measure_all",out_prefix,".txt",sep=""), sep=",")
|
|
| 266 |
save(gwr_mod,file= paste(path,"/","results2_kriging_Assessment_measure_all",out_prefix,".RData",sep=""))
|
|
| 267 |
write.table(tb_diagnostic1, file= paste(path,"/","results2_gwr_Assessment_measure1",out_prefix,".txt",sep=""), sep=",")
|
|
| 268 |
write.table(tb, file= paste(path,"/","results2_gwr_Assessment_measure_all",out_prefix,".txt",sep=""), sep=",")
|
|
| 269 |
save(gwr_mod,file= paste(path,"/","results2_gwr_Assessment_measure_all",out_prefix,".RData",sep=""))
|
|
| 267 | 270 |
|
| 268 | 271 |
#### END OF SCRIPT |
Also available in: Unified diff
GWR, modification of code for model 8, raste prediction