Revision 5415478e
Added by Benoit Parmentier almost 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 |
################################################################################################## |
|
17 |
|
|
18 |
###Loading R library and packages |
|
19 |
#library(gtools) # loading some useful tools |
|
20 |
library(mgcv) # GAM package by Wood 2006 (version 2012) |
|
21 |
library(sp) # Spatial pacakge with class definition by Bivand et al. 2008 |
|
22 |
library(spdep) # Spatial pacakge with methods and spatial stat. by Bivand et al. 2012 |
|
23 |
library(rgdal) # GDAL wrapper for R, spatial utilities (Keitt et al. 2012) |
|
24 |
library(gstat) # Kriging and co-kriging by Pebesma et al. 2004 |
|
25 |
library(automap) # Automated Kriging based on gstat module by Hiemstra et al. 2008 |
|
26 |
library(spgwr) |
|
27 |
library(maptools) |
|
28 |
library(graphics) |
|
29 |
library(parallel) # Urbanek S. and Ripley B., package for multi cores & parallel processing |
|
30 |
library(raster) |
|
31 |
library(rasterVis) |
|
32 |
library(fields) # May be used later... |
|
33 |
|
|
34 |
### Parameters and argument |
|
35 |
|
|
36 |
infile1<- "ghcn_or_tmax_covariates_06262012_OR83M.shp" #GHCN shapefile containing variables for modeling 2010 |
|
37 |
#infile2<-"list_10_dates_04212012.txt" #List of 10 dates for the regression |
|
38 |
#infile2<-"list_2_dates_04212012.txt" |
|
39 |
infile2<-"list_365_dates_04212012.txt" |
|
40 |
infile3<-"LST_dates_var_names.txt" #LST dates name |
|
41 |
infile4<-"models_interpolation_05142012.txt" #Interpolation model names |
|
42 |
infile5<-"mean_day244_rescaled.rst" #Raster or grid for the locations of predictions |
|
43 |
#infile6<-"lst_climatology.txt" |
|
44 |
infile6<-"LST_files_monthly_climatology.txt" |
|
45 |
inlistf<-"list_files_05032012.txt" #Stack of images containing the Covariates |
|
46 |
|
|
47 |
#path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_07192012_GAM" |
|
48 |
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_07152012_GWR" #Jupiter LOCATION on Atlas for kriging |
|
49 |
|
|
50 |
#Station location of the study area |
|
51 |
#stat_loc<-read.table(paste(path,"/","location_study_area_OR_0602012.txt",sep=""),sep=",", header=TRUE) |
|
52 |
#GHCN Database for 1980-2010 for study area (OR) |
|
53 |
#data3<-read.table(paste(path,"/","ghcn_data_TMAXy1980_2010_OR_0602012.txt",sep=""),sep=",", header=TRUE) #Not needing at this stage... |
|
54 |
|
|
55 |
nmodels<-8 #number of models running |
|
56 |
y_var_name<-"dailyTmax" #variable value being modeled...("value" in the GHCND database) |
|
57 |
predval<-1 # if set to 1, full interpolation raster produced for the study area |
|
58 |
prederr<-0 # if set to 0, no uncertain error (e.g. standard error or kriging std dev) is produced |
|
59 |
prop<-0.3 #Proportion of testing retained for validation |
|
60 |
#prop<-0.25 |
|
61 |
seed_number<- 100 #Seed number for random sampling |
|
62 |
out_prefix<-"test2_07312012_365d_gwr" #User defined output prefix |
|
63 |
setwd(path) |
|
64 |
|
|
65 |
#source("fusion_function_07192012.R") |
|
66 |
#source("KrigingUK_function_07262012.R") |
|
67 |
source("GWR_function_07312012.R") |
|
68 |
############ START OF THE SCRIPT ################## |
|
69 |
|
|
70 |
###Reading the station data and setting up for models' comparison |
|
71 |
filename<-sub(".shp","",infile1) #Removing the extension from file. |
|
72 |
ghcn<-readOGR(".", filename) #reading shapefile |
|
73 |
|
|
74 |
CRS<-proj4string(ghcn) #Storing projection information (ellipsoid, datum,etc.) |
|
75 |
|
|
76 |
mean_LST<- readGDAL(infile5) #Reading the whole raster in memory. This provides a grid for kriging |
|
77 |
proj4string(mean_LST)<-CRS #Assigning coordinate information to prediction grid. |
|
78 |
|
|
79 |
ghcn <- transform(ghcn,Northness = cos(ASPECT*pi/180)) #Adding a variable to the dataframe |
|
80 |
ghcn <- transform(ghcn,Eastness = sin(ASPECT*pi/180)) #adding variable to the dataframe. |
|
81 |
ghcn <- transform(ghcn,Northness_w = sin(slope*pi/180)*cos(ASPECT*pi/180)) #Adding a variable to the dataframe |
|
82 |
ghcn <- transform(ghcn,Eastness_w = sin(slope*pi/180)*sin(ASPECT*pi/180)) #adding variable to the dataframe. |
|
83 |
|
|
84 |
#Remove NA for LC and CANHEIGHT |
|
85 |
ghcn$LC1[is.na(ghcn$LC1)]<-0 |
|
86 |
ghcn$LC3[is.na(ghcn$LC3)]<-0 |
|
87 |
ghcn$CANHEIGHT[is.na(ghcn$CANHEIGHT)]<-0 |
|
88 |
|
|
89 |
dates <-readLines(paste(path,"/",infile2, sep="")) |
|
90 |
LST_dates <-readLines(paste(path,"/",infile3, sep="")) |
|
91 |
models <-readLines(paste(path,"/",infile4, sep="")) |
|
92 |
|
|
93 |
##Extracting the variables values from the raster files |
|
94 |
|
|
95 |
lines<-read.table(paste(path,"/",inlistf,sep=""), sep=" ") #Column 1 contains the names of raster files |
|
96 |
inlistvar<-lines[,1] |
|
97 |
inlistvar<-paste(path,"/",as.character(inlistvar),sep="") |
|
98 |
covar_names<-as.character(lines[,2]) #Column two contains short names for covaraites |
|
99 |
|
|
100 |
s_raster<- stack(inlistvar) #Creating a stack of raster images from the list of variables. |
|
101 |
layerNames(s_raster)<-covar_names #Assigning names to the raster layers |
|
102 |
projection(s_raster)<-CRS |
|
103 |
|
|
104 |
#stat_val<- extract(s_raster, ghcn3) #Extracting values from the raster stack for every point location in coords data frame. |
|
105 |
pos<-match("ASPECT",layerNames(s_raster)) #Find column with name "value" |
|
106 |
r1<-raster(s_raster,layer=pos) #Select layer from stack |
|
107 |
pos<-match("slope",layerNames(s_raster)) #Find column with name "value" |
|
108 |
r2<-raster(s_raster,layer=pos) #Select layer from stack |
|
109 |
N<-cos(r1*pi/180) |
|
110 |
E<-sin(r1*pi/180) |
|
111 |
Nw<-sin(r2*pi/180)*cos(r1*pi/180) #Adding a variable to the dataframe |
|
112 |
Ew<-sin(r2*pi/180)*sin(r1*pi/180) #Adding variable to the dataframe. |
|
113 |
|
|
114 |
pos<-match("LC1",layerNames(s_raster)) #Find column with name "value" |
|
115 |
LC1<-raster(s_raster,layer=pos) #Select layer from stack |
|
116 |
s_raster<-dropLayer(s_raster,pos) |
|
117 |
LC1[is.na(LC1)]<-0 |
|
118 |
pos<-match("LC3",layerNames(s_raster)) #Find column with name "value" |
|
119 |
LC3<-raster(s_raster,layer=pos) #Select layer from stack |
|
120 |
s_raster<-dropLayer(s_raster,pos) |
|
121 |
LC3[is.na(LC3)]<-0 |
|
122 |
|
|
123 |
xy<-coordinates(r1) #get x and y projected coordinates... |
|
124 |
xy_latlon<-project(xy, CRS, inv=TRUE) # find lat long for projected coordinats (or pixels...) |
|
125 |
lon<-raster(xy_latlon) #Transform a matrix into a raster object ncol=ncol(r1), nrow=nrow(r1)) |
|
126 |
ncol(lon)<-ncol(r1) |
|
127 |
nrow(lon)<-nrow(r1) |
|
128 |
extent(lon)<-extent(r1) |
|
129 |
projection(lon)<-CRS #At this stage this is still an empty raster with 536 nrow and 745 ncell |
|
130 |
lat<-lon |
|
131 |
values(lon)<-xy_latlon[,1] |
|
132 |
values(lat)<-xy_latlon[,2] |
|
133 |
|
|
134 |
r<-stack(N,E,Nw,Ew,lon,lat,LC1,LC3) |
|
135 |
rnames<-c("Northness","Eastness","Northness_w","Eastness_w", "lon","lat","LC1","LC3") |
|
136 |
layerNames(r)<-rnames |
|
137 |
s_raster<-addLayer(s_raster, r) |
|
138 |
|
|
139 |
#s_sgdf<-as(s_raster,"SpatialGridDataFrame") #Conversion to spatial grid data frame |
|
140 |
|
|
141 |
####### Preparing LST stack of climatology... |
|
142 |
|
|
143 |
#l=list.files(pattern="mean_month.*rescaled.rst") |
|
144 |
l <-readLines(paste(path,"/",infile6, sep="")) |
|
145 |
molst<-stack(l) #Creating a raster stack... |
|
146 |
#setwd(old) |
|
147 |
molst<-molst-273.16 #K->C #LST stack of monthly average... |
|
148 |
idx <- seq(as.Date('2010-01-15'), as.Date('2010-12-15'), 'month') |
|
149 |
molst <- setZ(molst, idx) |
|
150 |
layerNames(molst) <- month.abb |
|
151 |
|
|
152 |
####### CREATE A MASK FOR WATER AND |
|
153 |
|
|
154 |
pos<-match("LC10",layerNames(s_raster)) #Find column with the current month for instance mm12 |
|
155 |
LC10<-raster(s_raster,pos) |
|
156 |
|
|
157 |
LC10_mask<-LC10 |
|
158 |
LC10_mask[is.na(LC10_mask)]<-0 |
|
159 |
LC10_mask[LC10==100]<-NA |
|
160 |
LC10_mask[LC10_mask<100]<-1 |
|
161 |
LC10_mask[is.na(LC10_mask)]<-0 |
|
162 |
mask_land_NA<-LC10_mask |
|
163 |
mask_land_NA[mask_land_NA==0]<-NA |
|
164 |
|
|
165 |
breakpoints<-c(0,0.99) |
|
166 |
colors<-c("black","red") |
|
167 |
#plot(LC10_mask, breaks=breakpoints,col=colors) |
|
168 |
#quartz() |
|
169 |
plot(LC10_mask, col=colors) |
|
170 |
|
|
171 |
# ELEV_SRTM[ELEV_SRTM==-9999]<-NA |
|
172 |
# avl<-c(0,60000,1) |
|
173 |
# rclmat<-matrix(avl,ncol=3,byrow=TRUE) |
|
174 |
# ELEV_rc<-reclass(ELEV_SRTM,rclmat) |
|
175 |
# ELEV_SRTM_m <-mask(ELEV_STRM,mask_land_NA) |
|
176 |
# mask2<-ELEV_SRTM_m |
|
177 |
# mask2[is.na(mask)]<-1 |
|
178 |
# s_rst_m <-mask(s_raster,mask_land_NA) #This assigns NA to all values from LC1 that are NA in the mask_land_NA |
|
179 |
# |
|
180 |
####### Preparing tables for model assessment: specific diagnostic/metrics |
|
181 |
|
|
182 |
#Model assessment: specific diagnostics/metrics |
|
183 |
results_m1<- matrix(1,1,nmodels+3) #Diagnostic metrics specific to the modeleling framework |
|
184 |
results_m2<- matrix(1,1,nmodels+3) |
|
185 |
results_m3<- matrix(1,1,nmodels+3) |
|
186 |
#results_RMSE_f<- matrix(1,length(models)+3) |
|
187 |
|
|
188 |
#Model assessment: general diagnostic/metrics |
|
189 |
results_RMSE <- matrix(1,1,nmodels+3) |
|
190 |
results_MAE <- matrix(1,1,nmodels+3) |
|
191 |
results_ME <- matrix(1,1,nmodels+3) #There are 8 models for kriging!!! |
|
192 |
results_R2 <- matrix(1,1,nmodels+3) #Coef. of determination for the validation dataset |
|
193 |
|
|
194 |
results_RMSE_f<- matrix(1,1,nmodels+3) #RMSE fit, RMSE for the training dataset |
|
195 |
results_MAE_f <- matrix(1,1,nmodels+3) |
|
196 |
results_R2_f <- matrix(1,1,nmodels+3) |
|
197 |
|
|
198 |
######### Preparing daily values for training and testing |
|
199 |
|
|
200 |
#Screening for bad values: value is tmax in this case |
|
201 |
#ghcn$value<-as.numeric(ghcn$value) |
|
202 |
ghcn_all<-ghcn |
|
203 |
ghcn_test<-subset(ghcn,ghcn$value>-150 & ghcn$value<400) |
|
204 |
ghcn_test2<-subset(ghcn_test,ghcn_test$ELEV_SRTM>0) |
|
205 |
ghcn<-ghcn_test2 |
|
206 |
#coords<- ghcn[,c('x_OR83M','y_OR83M')] |
|
207 |
|
|
208 |
##Sampling: training and testing sites... |
|
209 |
|
|
210 |
set.seed(seed_number) #Using a seed number allow results based on random number to be compared... |
|
211 |
ghcn.subsets <-lapply(dates, function(d) subset(ghcn, date==d)) #this creates a list of 10 or 365 subsets dataset based on dates |
|
212 |
sampling<-vector("list",length(dates)) |
|
213 |
|
|
214 |
for(i in 1:length(dates)){ |
|
215 |
n<-nrow(ghcn.subsets[[i]]) |
|
216 |
ns<-n-round(n*prop) #Create a sample from the data frame with 70% of the rows |
|
217 |
nv<-n-ns #create a sample for validation with prop of the rows |
|
218 |
ind.training <- sample(nrow(ghcn.subsets[[i]]), size=ns, replace=FALSE) #This selects the index position for 70% of the rows taken randomly |
|
219 |
ind.testing <- setdiff(1:nrow(ghcn.subsets[[i]]), ind.training) |
|
220 |
sampling[[i]]<-ind.training |
|
221 |
} |
|
222 |
|
|
223 |
######## Prediction for the range of dates |
|
224 |
|
|
225 |
#gwr_mod<-mclapply(1:length(dates), runGWR,mc.preschedule=FALSE,mc.cores = 8) #This is the end bracket from mclapply(...) statement |
|
226 |
#fusion_mod357<-mclapply(357:365,runFusion, mc.cores=8)# for debugging |
|
227 |
#test<-runKriging(1) |
|
228 |
#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 |
|
230 |
|
|
231 |
#test<-mclapply(357,runFusion, mc.cores=1)# for debugging |
|
232 |
|
|
233 |
tb<-gwr_mod[[1]][[3]][0,] #empty data frame with metric table structure that can be used in rbinding... |
|
234 |
tb_tmp<-gwr_mod #copy |
|
235 |
|
|
236 |
for (i in 1:length(tb_tmp)){ |
|
237 |
tmp<-tb_tmp[[i]][[3]] |
|
238 |
tb<-rbind(tb,tmp) |
|
239 |
} |
|
240 |
rm(tb_tmp) |
|
241 |
|
|
242 |
for(i in 4:nmodels+3){ # start of the for loop #1 |
|
243 |
tb[,i]<-as.numeric(as.character(tb[,i])) |
|
244 |
} |
|
245 |
|
|
246 |
tb_RMSE<-subset(tb, metric=="RMSE") |
|
247 |
tb_MAE<-subset(tb,metric=="MAE") |
|
248 |
tb_ME<-subset(tb,metric=="ME") |
|
249 |
tb_R2<-subset(tb,metric=="R2") |
|
250 |
tb_RMSE_f<-subset(tb, metric=="RMSE_f") |
|
251 |
tb_MAE_f<-subset(tb,metric=="MAE_f") |
|
252 |
|
|
253 |
tb_diagnostic1<-rbind(tb_RMSE,tb_MAE,tb_ME,tb_R2) |
|
254 |
#tb_diagnostic2<-rbind(tb_,tb_MAE,tb_ME,tb_R2) |
|
255 |
|
|
256 |
mean_RMSE<-sapply(tb_RMSE[,4:(nmodels+3)],mean) |
|
257 |
mean_MAE<-sapply(tb_MAE[,4:(nmodels+3)],mean) |
|
258 |
mean_R2<-sapply(tb_R2[,4:(nmodels+3)],mean) |
|
259 |
mean_ME<-sapply(tb_ME[,4:(nmodels+3)],mean) |
|
260 |
mean_MAE_f<-sapply(tb_MAE[,4:(nmodels+3)],mean) |
|
261 |
mean_RMSE_f<-sapply(tb_RMSE_f[,4:(nmodels+3)],mean) |
|
262 |
|
|
263 |
#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 |
|
|
268 |
#### END OF SCRIPT |
Also available in: Unified diff
GWR, raster prediction full year main script