1
|
##################################### METHODS COMPARISON part 3 ##########################################
|
2
|
#################################### Spatial Analysis ############################################
|
3
|
#This script utilizes the R ojbects created during the interpolation phase. #
|
4
|
#At this stage the script produces figures of various accuracy metrics and compare methods: #
|
5
|
#- rank station in terms of residuals and metrics (MAE and RMSE) #
|
6
|
#- calculate accuary metrics for climtology predictions... #
|
7
|
#- spatial density of station network and accuracy metric #
|
8
|
#- visualization of maps of prediction and difference for comparison #
|
9
|
#AUTHOR: Benoit Parmentier #
|
10
|
#DATE: 10/23/2012 #
|
11
|
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491 -- #
|
12
|
###################################################################################################
|
13
|
|
14
|
###Loading R library and packages
|
15
|
library(gtools) # loading some useful tools such as mixedsort
|
16
|
library(mgcv) # GAM package by Wood 2006 (version 2012)
|
17
|
library(sp) # Spatial pacakge with class definition by Bivand et al. 2008
|
18
|
library(spdep) # Spatial package with methods and spatial stat. by Bivand et al. 2012
|
19
|
library(rgdal) # GDAL wrapper for R, spatial utilities (Keitt et al. 2012)
|
20
|
library(gstat) # Kriging and co-kriging by Pebesma et al. 2004
|
21
|
library(automap) # Automated Kriging based on gstat module by Hiemstra et al. 2008
|
22
|
library(spgwr)
|
23
|
library(gpclib)
|
24
|
library(maptools)
|
25
|
library(graphics)
|
26
|
library(parallel) # Urbanek S. and Ripley B., package for multi cores & parralel processing
|
27
|
library(raster)
|
28
|
library(rasterVis)
|
29
|
library(plotrix) #Draw circle on graph
|
30
|
library(reshape)
|
31
|
## Functions
|
32
|
#loading R objects that might have similar names
|
33
|
load_obj <- function(f)
|
34
|
{
|
35
|
env <- new.env()
|
36
|
nm <- load(f, env)[1]
|
37
|
env[[nm]]
|
38
|
}
|
39
|
|
40
|
###Parameters and arguments
|
41
|
|
42
|
infile1<- "ghcn_or_tmax_covariates_06262012_OR83M.shp" #GHCN shapefile containing variables for modeling 2010
|
43
|
#infile2<-"list_10_dates_04212012.txt" #List of 10 dates for the regression
|
44
|
infile2<-"list_365_dates_04212012.txt" #list of dates
|
45
|
infile3<-"LST_dates_var_names.txt" #LST dates name
|
46
|
infile4<-"models_interpolation_05142012.txt" #Interpolation model names
|
47
|
infile5<-"mean_day244_rescaled.rst" #mean LST for day 244
|
48
|
inlistf<-"list_files_05032012.txt" #list of raster images containing the Covariates
|
49
|
infile6<-"OR83M_state_outline.shp"
|
50
|
#stat_loc<-read.table(paste(path,"/","location_study_area_OR_0602012.txt",sep=""),sep=",", header=TRUE)
|
51
|
|
52
|
out_prefix<-"methods_11072012_"
|
53
|
#output path
|
54
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_07192012_GAM"
|
55
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_07192012_GAM"
|
56
|
|
57
|
setwd(path)
|
58
|
|
59
|
##### LOAD USEFUL DATA
|
60
|
|
61
|
obj_list<-"list_obj_08262012.txt" #Results of fusion from the run on ATLAS
|
62
|
path<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison_10242012" #Jupiter LOCATION on Atlas for kriging #Jupiter Location on XANDERS
|
63
|
#path<-"/Users/benoitparmentier/Dropbox/Data/NCEAS/Oregon_covariates/" #Local dropbox folder on Benoit's laptop
|
64
|
setwd(path)
|
65
|
proj_str="+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";
|
66
|
#User defined output prefix
|
67
|
|
68
|
sampling_CAI<-load_obj("results2_CAI_sampling_obj_09132012_365d_GAM_CAI2_multisampling2.RData")
|
69
|
sampling_fus<-load_obj("results2_fusion_sampling_obj_10d_GAM_fusion_multisamp4_09192012.RData")
|
70
|
fus_CAI_mod<-load_obj("results2_CAI_Assessment_measure_all_09132012_365d_GAM_CAI2_multisampling2.RData")
|
71
|
gam_fus_mod1<-load_obj("results2_fusion_Assessment_measure_all_10d_GAM_fusion_multisamp4_09192012.RData")
|
72
|
|
73
|
filename<-sub(".shp","",infile1) #Removing the extension from file.
|
74
|
ghcn<-readOGR(".", filename) #reading shapefile
|
75
|
filename<-sub(".shp","",infile6) #Removing the extension from file.
|
76
|
#reg_outline<-readOGR(".", filename) #reading shapefile
|
77
|
reg_outline<-readOGR(dsn=".", layer=filename) #reading shapefile
|
78
|
|
79
|
#CRS<-proj4string(ghcn) #Storing projection information (ellipsoid, datum,etc.)
|
80
|
lines<-read.table(paste(path,"/",inlistf,sep=""), sep="") #Column 1 contains the names of raster files
|
81
|
inlistvar<-lines[,1]
|
82
|
inlistvar<-paste(path,"/",as.character(inlistvar),sep="")
|
83
|
covar_names<-as.character(lines[,2]) #Column two contains short names for covaraites
|
84
|
|
85
|
s_raster<- stack(inlistvar) #Creating a stack of raster images from the list of variables.
|
86
|
layerNames(s_raster)<-covar_names #Assigning names to the raster layers
|
87
|
projection(s_raster)<-proj_str
|
88
|
|
89
|
#Create mask using land cover data
|
90
|
pos<-match("LC10",layerNames(s_raster)) #Find the layer which contains water bodies
|
91
|
LC10<-subset(s_raster,pos)
|
92
|
LC10[is.na(LC10)]<-0 #Since NA values are 0, we assign all zero to NA
|
93
|
#LC10_mask[LC10_mask<100]<-1
|
94
|
mask_land<-LC10==100 #All values below 100% water are assigned the value 1, value 0 is "water"
|
95
|
mask_land_NA<-mask_land
|
96
|
mask_land_NA[mask_land_NA==0]<-NA #Water bodies are assigned value 1
|
97
|
|
98
|
data_name<-"mask_land_OR"
|
99
|
raster_name<-paste(data_name,".rst", sep="")
|
100
|
writeRaster(mask_land, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI)
|
101
|
#writeRaster(r2, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI)
|
102
|
|
103
|
pos<-match("mm_01",layerNames(s_raster))
|
104
|
mm_01<-subset(s_raster,pos)
|
105
|
mm_01<-mm_01-273.15
|
106
|
mm_01<-mask(mm_01,mask_land_NA) #Raster image used as backround
|
107
|
#mention this is the last... files
|
108
|
|
109
|
################# PART 1: COMPARING RESULTS USING ALL MONTHNLY DATA AND SAMPLED FROM DAILY ###################
|
110
|
######## Visualize data: USING ALL MONTHLY DATA...
|
111
|
|
112
|
file_pat<-glob2rx("GAMCAI_tmax_pred_predicted_mod*_365d_GAM_CAI2_all_lstd_10262012.rst")
|
113
|
lf_cai_all<-list.files(pattern=file_pat)
|
114
|
image_file<-grep(file_pat,lf_ck,value=TRUE) # using grep with "value" extracts the matching names
|
115
|
raster_pred_k<-raster(image_file)
|
116
|
|
117
|
#### CHECKING RESULTS ON 10/28 BY VISUALIZING PREDICTIONS... DO THAT FOR ALL MODELS AND TWO DATES...
|
118
|
|
119
|
date_selected<-"20100103"
|
120
|
lf_fus<-list.files(pattern=paste("GAM_bias_tmax_predicted_mod8_",date_selected,"_30_1_365d_GAM_fusion_all_lstd_10242012.rst$",sep="")) #Search for files in relation to fusion
|
121
|
lf_cai<-list.files(pattern=paste("*CAI_tmax_pred.*",date_selected,"*.08072012_365d_GAM_CAI2.rst$",sep="")) #Search for files in relation to fusion
|
122
|
|
123
|
lf_cai<-list.files(pattern=paste("GAMCAI_tmax_pred_predicted_mod8_20100103_30_1_365d_GAM_CAI2_all_lstd_10262012.rst",sep=""))
|
124
|
|
125
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_GAM"
|
126
|
setwd(path)
|
127
|
tmax_predicted_mod8<-raster("GAM_bias_tmax_predicted_mod8_20100103_30_1_365d_GAM_fusion_all_lstd_10242012.rst")
|
128
|
tmax_bias_mod8<-raster("GAM_bias_predicted_mod8_20100103_30_1_365d_GAM_fusion_all_lstd_10242012.rst")
|
129
|
daily_delta_tmax<-raster("fusion_daily_delta_20100103_30_1_365d_GAM_fusion_all_lstd_10242012.rst")
|
130
|
fusion<-raster("fusion_tmax_predicted_20100103_30_1_365d_GAM_fusion_all_lstd_10242012.rst")
|
131
|
LST<-(raster("mean_month1_rescaled.rst"))-273.16
|
132
|
|
133
|
tmp_rast<-LST+daily_delta_tmax-tmax_bias_mod8
|
134
|
dif_rast<-tmp_rast-tmax_predicted_mod8
|
135
|
eq<-tmp_rast==tmax_predicted_mod8
|
136
|
plot(tmax_predicted_mod8)
|
137
|
|
138
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_CAI"
|
139
|
setwd(path)
|
140
|
cai_tmax_mod8_rast<-raster("GAMCAI_tmax_pred_predicted_mod8_20100103_30_1_365d_GAM_CAI2_all_lstd_10262012.rst")
|
141
|
cai_delta_rast<-raster("CAI_deltaclim_20100103_30_1_365d_GAM_CAI2_all_lstd_10262012.rst")
|
142
|
cai_clim_rast<-raster("CAI_clim_20100103_30_1_365d_GAM_CAI2_all_lstd_10262012.rst")
|
143
|
cai_climod8_rast<-raster("GAMCAI_predicted_mod8_20100103_30_1_365d_GAM_CAI2_all_lstd_10262012.rst")
|
144
|
|
145
|
# s.range <- c(min(minValue(diff)), max(maxValue(diff)))
|
146
|
# col.breaks <- pretty(s.range, n=50)
|
147
|
# lab.breaks <- pretty(s.range, n=5)
|
148
|
temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
|
149
|
# plot(diff, breaks=col.breaks, col=temp.colors(length(col.breaks)-1)
|
150
|
#
|
151
|
pred_cai_mod8<-stack(cai_climod8_rast,cai_delta_rast,cai_tmax_mod8_rast)
|
152
|
layerNames(pred_cai_mod8)<-paste(c("cai_climod8_rast","cai_delta_rast","cai_tmax_mod8_rast"),"20100103",sep=" ")
|
153
|
pred_fus_mod8<-stack(tmax_bias_mod8,daily_delta_tmax,LST,tmax_predicted_mod8)
|
154
|
layerNames(pred_fus_mod8)<-paste(c("tmax_bias_mod8","daily_delta_tmax","LST","tmax_predicted_mod8"),"20100103",sep=" ")
|
155
|
pred_fus_cai_comp<-stack(pred_fus_mod8,pred_cai_mod8)
|
156
|
|
157
|
X11(12,12)
|
158
|
#for (i in 1:)
|
159
|
#plot(pred_cai_mod8,col=temp.colors)
|
160
|
plot(pred_cai_mod8,col=temp.colors(25))
|
161
|
#title("prediction for 20100103")
|
162
|
plot(pred_fus_mod8,col=temp.colors(25))
|
163
|
plot(pred_fus_cai_comp,col=temp.colors(25))
|
164
|
dev.off
|
165
|
tmp_rast2<-cai_delta_rast + cai_climod8_rast
|
166
|
dif2_rast<-tmp_rast2-cai_tmax_mod8_rast
|
167
|
|
168
|
dif2<-cai_rast-tmax_predicted_mod8 #This seems to be the same!!!
|
169
|
plot(cai_climod8_rast) #THe bias surface is different when compared to the climatology based on mod8
|
170
|
plot(tmax_bias_mod8)
|
171
|
plot(tmax_predicted_mod8)
|
172
|
plot(cai_tmax_rast)
|
173
|
|
174
|
cellStats(pred_fus_cai_comp)
|
175
|
|
176
|
#t<-"results2_CAI_Assessment_measure_all_08072012_365d_GAM_CAI2.RData"
|
177
|
|
178
|
## Examining data with simplified CAI3 models in comparison to other: 11/03/2012
|
179
|
|
180
|
oldpath<-getwd()
|
181
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_CAI"
|
182
|
setwd(path)
|
183
|
|
184
|
date_selected<-"20100103"
|
185
|
#lf_fus<-list.files(pattern=paste("GAM_bias_tmax_predicted_mod8_",date_selected,"_30_1_365d_GAM_fusion_all_lstd_10242012.rst$",sep="")) #Search for files in relation to fusion
|
186
|
file_pat<-glob2rx(paste("*tmax_pred*",date_selected,"*_365d_GAM_CAI3_all_10302012.rst",sep="")) #Search for files in relation to fusion
|
187
|
lf_cai3<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
188
|
|
189
|
cai3_mod<-load_obj("results_mod_obj__365d_GAM_CAI3_all_10302012.RData")
|
190
|
names(cai3_mod$gam_CAI_mod[[1]]) #check names within object
|
191
|
formulas<-as.character(cai3_mod$gam_CAI_mod[[1]]$formulas)
|
192
|
rmse<-(cai3_mod$gam_CAI_mod[[1]]$tb_metrics1)
|
193
|
nmodels<-6
|
194
|
for(i in 4:(nmodels+3)){ # start of the for loop #1
|
195
|
rmse[,i]<-as.numeric(as.character(rmse[,i]))
|
196
|
}
|
197
|
rmse_val<-as.numeric(c(rmse[1,9],rmse[1,4:8]))
|
198
|
rmse_val<-format(rmse_val,digits=3)
|
199
|
rmse_val<-paste("(",rmse_val,")",sep="")
|
200
|
rname<-c("CAI_Kr",formulas)
|
201
|
rname<-paste(rname,rmse_val,sep=" ")
|
202
|
|
203
|
rast_cai3<-stack(lf_cai3)
|
204
|
rast_cai3<-mask(rast_cai3,mask_ELEV_SRTM)
|
205
|
#s_range <- c(min(minValue(diff)), max(maxValue(diff)))
|
206
|
s_range<-c(minValue(rast_cai3),maxValue(rast_cai3)) #stack min and max
|
207
|
s_range<-c(min(s_range),max(s_range))
|
208
|
col_breaks <- pretty(s_range, n=50)
|
209
|
lab_breaks <- pretty(s_range, n=5)
|
210
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
211
|
par(mfrow=c(2,3))
|
212
|
for (k in 1:length(lf_cai3)){
|
213
|
cai3_r<-raster(rast_cai3,k)
|
214
|
plot(cai3_r, breaks=col_breaks, col=temp_colors(length(col_breaks)-1),
|
215
|
axis=list(at=lab_breaks, labels=lab_breaks),main=rname[k])
|
216
|
}
|
217
|
|
218
|
#Wihtout setting range
|
219
|
s_range<-c(-12,18)
|
220
|
#col_breaks <- pretty(s_range, n=50)
|
221
|
#lab_breaks <- pretty(s_range, n=5)
|
222
|
col_breaks<-49
|
223
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
224
|
par(mfrow=c(2,3))
|
225
|
for (k in 1:length(lf_cai3)){
|
226
|
cai3_r<-raster(rast_cai3,k)
|
227
|
plot(cai3_r, col=temp_colors(col_breaks),
|
228
|
,main=rname[k])
|
229
|
}
|
230
|
|
231
|
###### PART 2: ################
|
232
|
#### EXAMINING LST AND ELEVATION DATASETS TO CHECK FOR EXTREMES AND SELECT SCREENING VALUES....
|
233
|
|
234
|
# Examining ELEV_SRTM and screening out...
|
235
|
|
236
|
pos<-match("ELEV_SRTM",layerNames(s_raster))
|
237
|
ELEV_SRTM<-raster(s_raster,pos)
|
238
|
head(freq(ELEV_SRTM)) #Show how many pixels are below 0
|
239
|
|
240
|
# mask out values below zero
|
241
|
elev<-ELEV_SRTM
|
242
|
elev[0<elev]<-NA #Remove all negative elevation
|
243
|
#Show how many were NA...
|
244
|
|
245
|
mask_ELEV_SRTM<-elev >0
|
246
|
quartz()
|
247
|
plot(elev, main="Elevation in Oregon")
|
248
|
savePlot(paste("elevation_Oregon",out_prefix,".png", sep=""), type="png")
|
249
|
|
250
|
# Examining LST data: extremes values and number of valid observation
|
251
|
|
252
|
l<-list.files(pattern="mean_month.*rescaled.rst")
|
253
|
ln<-list.files(pattern="mean_month_valid_obs_.*_Sum.rst") #number of observation missing
|
254
|
l<-mixedsort(l)
|
255
|
ln<-mixedsort(ln)
|
256
|
|
257
|
#l <-readLines(paste(path,"/",infile6, sep=""))
|
258
|
molst<-stack(l) #Creating a raster stack...
|
259
|
nobslst<-stack(ln) #Creating a raster stack...
|
260
|
|
261
|
#setwd(old)
|
262
|
molst<-molst-273.16 #K->C #LST stack of monthly average...
|
263
|
idx <- seq(as.Date('2010-01-15'), as.Date('2010-12-15'), 'month')
|
264
|
molst <- setZ(molst, idx)
|
265
|
layerNames(molst) <- month.abb
|
266
|
layerNames(nobslst)<-month.abb
|
267
|
|
268
|
X11(width=18,height=12)
|
269
|
#PLOT THE STACK AND SAVE IMAGE
|
270
|
plot(molst,col=temp.colors(25))
|
271
|
savePlot("lst_climatology_OR.png",type="png")
|
272
|
levelplot(molst,col.regions=temp.colors(25))
|
273
|
plot(nobslst)
|
274
|
savePlot("lst_climatology_OR_nobs.png",type="png")
|
275
|
dev.off)
|
276
|
|
277
|
#note differences in patternin agricultural areas and
|
278
|
extremeValues(molst) #not yet implemented...
|
279
|
min_values<-cellStats(molst,"min")
|
280
|
max_values<-cellStats(molst,"max")
|
281
|
mean_values<-cellStats(molst,"mean")
|
282
|
sd_values<-cellStats(molst,"sd")
|
283
|
#median_values<-cellStats(molst,"median") Does not extist
|
284
|
statistics_LST_s<-cbind(min_values,max_values,mean_values,sd_values) #This shows that some values are extremes...especially in October
|
285
|
LST_stat_data<-as.data.frame(statistics_LST_s)
|
286
|
names(LST_stat_data)<-c("min","max","mean","sd")
|
287
|
# Statistics for number of valid observation stack
|
288
|
min_values<-cellStats(nobslst,"min")
|
289
|
max_values<-cellStats(nobslst,"max")
|
290
|
mean_values<-cellStats(nobslst,"mean")
|
291
|
sd_values<-cellStats(nobslst,"sd")
|
292
|
statistics_LSTnobs_s<-cbind(min_values,max_values,mean_values,sd_values) #This shows that some values are extremes...especially in October
|
293
|
LSTnobs_stat_data<-as.data.frame(statistics_LSTnobs_s)
|
294
|
|
295
|
X11(width=12,height=12)
|
296
|
#Plot statiscs (mean,min,max) for monthly LST images
|
297
|
plot(1:12,LST_stat_data$mean,type="b",ylim=c(-15,60),col="black",xlab="month",ylab="tmax (degree C)")
|
298
|
lines(1:12,LST_stat_data$min,type="b",col="blue")
|
299
|
lines(1:12,LST_stat_data$max,type="b",col="red")
|
300
|
text(1:12,LST_stat_data$mean,rownames(LST_stat_data),cex=1,pos=2)
|
301
|
|
302
|
legend("topleft",legend=c("min","mean","max"), cex=1.5, col=c("blue","black","red"),
|
303
|
lty=1)
|
304
|
title(paste("LST statistics for Oregon", "2010",sep=" "))
|
305
|
savePlot("lst_statistics_OR.png",type="png")
|
306
|
|
307
|
#Plot number of valid observations for LST
|
308
|
plot(1:12,LSTnobs_stat_data$mean,type="b",ylim=c(0,280),col="black",xlab="month",ylab="tmax (degree C)")
|
309
|
lines(1:12,LSTnobs_stat_data$min,type="b",col="blue")
|
310
|
lines(1:12,LSTnobs_stat_data$max,type="b",col="red")
|
311
|
text(1:12,LSTnobs_stat_data$mean,rownames(LSTnobs_stat_data),cex=1,pos=2)
|
312
|
|
313
|
legend("topleft",legend=c("min","mean","max"), cex=1.5, col=c("blue","black","red"),
|
314
|
lty=1)
|
315
|
title(paste("LST number of valid observations for Oregon", "2010",sep=" "))
|
316
|
savePlot("lst_nobs_OR.png",type="png")
|
317
|
|
318
|
mm_10<-raster(molst,10)
|
319
|
tab_mm_10<-freq(mm_10) #This shows that there are only 6 pixels with values between -86 and -89.. and 8 pixels with values between 38 and 31
|
320
|
head(tab_mm_10)
|
321
|
#screening maybe good...?
|
322
|
|
323
|
mm_01<-raster(molst,1)
|
324
|
tab_mm_10<-freq(mm_01) #This shows that there are only 6 pixels with values between -86 and -89.. and 8 pixels with values between 38 and 31
|
325
|
#head(tab_mm_10)
|
326
|
|
327
|
mm_07<-raster(molst,7)
|
328
|
tab_mm_10<-freq(mm_07) #This shows that there are only 6 pixels with values between -86 and -89.. and 8 pixels with values between 38 and 31
|
329
|
#head(tab_mm_10)
|
330
|
|
331
|
LST_comp<-stack(mm_01,mm_07)
|
332
|
|
333
|
X11(width=18,height=10)
|
334
|
#PLOT THE STACK AND SAVE IMAGE
|
335
|
plot(LST_comp,col=temp.colors(25))
|
336
|
#title("LST spatial pattern in January and August")
|
337
|
savePlot("lst_climatology_OR.png",type="png")
|
338
|
|
339
|
|
340
|
###### PART 3: EXAMANING MODEL OUTPUTS USING MOD OBJECTS TO CHECK AND COMPARE RESULTS ############
|
341
|
############ LOOKING AT SPECIFIC STATIONS...
|
342
|
path_data_cai<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_CAI"
|
343
|
path_data_fus<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_GAM"
|
344
|
setwd(path_data) #output data
|
345
|
gam_fus<-load_obj(file.path(path_data_fus,
|
346
|
"results2_fusion_Assessment_measure_all_365d_GAM_fusion_all_lstd_10242012.RData"))
|
347
|
gam_cai<-load_obj(file.path(path_dat_cai,
|
348
|
"results2_CAI_Assessment_measure_all_365d_GAM_CAI2_all_lstd_10262012.RData") #This contains all the info
|
349
|
sampling_obj_cai<-load_obj(file.path(path_data_cai,
|
350
|
"results2_CAI_sampling_obj_365d_GAM_CAI2_all_lstd_10262012.RData"))
|
351
|
sampling_date_list<-sampling_obj_cai$sampling_dat$date
|
352
|
date_selected="20100103"
|
353
|
|
354
|
k<-match(date_selected,sampling_date_list)
|
355
|
|
356
|
fus_d<-gam_fus[[k]] #object for the first date...20100103
|
357
|
fus_d$tb_metrics1 #This is a summary for validation metrics...(RMSE and MAE)
|
358
|
|
359
|
cai_d<-gam_cai[[k]] #object for the first date...20100103
|
360
|
cai_d$tb_metrics1 #This is a summary for validation metrics...(RMSE and MAE)
|
361
|
|
362
|
names(fus_d$mod_obj) #list models, nine models: mod1, mod2,...,mod9a, mod9b
|
363
|
names(cai_d$mod_obj)
|
364
|
|
365
|
mod1_f<-fus_d$mod_obj$mod1 # extract model 1 for fusion on
|
366
|
mod4_f<-fus_d$mod_obj$mod4
|
367
|
mod7_f<-fus_d$mod_obj$mod7
|
368
|
data_sf<-fus_d$data_s
|
369
|
data_vf<-fus_d$data_v
|
370
|
|
371
|
plot(mod1_f) #Examining the fit
|
372
|
plot(mod1) #Examing the fit
|
373
|
#data_monthf<-fus_d1$data_month
|
374
|
|
375
|
mod1<-cai_d1$mod_obj$mod1
|
376
|
mod4<-cai_d1$mod_obj$mod4
|
377
|
mod7<-cai_d1$mod_obj$mod7
|
378
|
data_s<-cai_d1$data_s
|
379
|
data_v<-cai_d1$data_v
|
380
|
data_month<-cai_d1$data_month
|
381
|
|
382
|
#### COMPARE LST AND TMAX averages over the year... ###############
|
383
|
#### Screening necessary?
|
384
|
|
385
|
LSTD_bias_month<-vector("list",length(gam_cai)) # list to store bias information
|
386
|
data_month_list<-vector("list",length(gam_cai)) # list to store data_month used for training
|
387
|
|
388
|
for (i in 1:length(gam_cai)){
|
389
|
data_month_list[[i]]<-gam_cai[[i]]$data_month
|
390
|
#LSTD_bias_month[[i]]<-aggregate(LSTD_bias~month,data=data_month_list[[i]],mean)
|
391
|
}
|
392
|
|
393
|
#tmp<-do.call(rbind,LSTD_bias_month)
|
394
|
data_month_all<-do.call(rbind,data_month_list)
|
395
|
|
396
|
data_month_all$LST<-data_month_all$LST-273.16
|
397
|
LSTD_bias_avgm<-aggregate(LSTD_bias~month,data=data_month_all,mean)
|
398
|
LSTD_bias_sdm<-aggregate(LSTD_bias~month,data=data_month_all,sd)
|
399
|
|
400
|
LST_avgm<-aggregate(LST~month,data=data_month_all,mean)
|
401
|
TMax_avgm<-aggregate(TMax~month,data=data_month_all,mean)
|
402
|
|
403
|
#plot(LST_avgm,type="b")
|
404
|
#lines(TMax_avgm,col="blue",add=T)
|
405
|
plot(data_month_all$month,data_month_all$LST, xlab="month",ylab="LST at station")
|
406
|
title(paste("Monthly LST at stations in Oregon", "2010",sep=" "))
|
407
|
savePlot(paste("LST_at_stations_per_month_",out_prefix,".png", sep=""), type="png")
|
408
|
|
409
|
tab_freq<-cbind(hist_LST$breaks,hist_LST$counts)
|
410
|
median(data_month_all$LST,na.rm=T)
|
411
|
median(data_month_all$TMax,na.rm=T)
|
412
|
|
413
|
statistics_LST_s<-as.data.frame(statistics_LST_s)
|
414
|
plot(TMax_avgm,type="b",ylim=c(0,35),col="red",xlab="month",ylab="tmax (degree C)")
|
415
|
lines(1:12,statistics_LST_s$mean_values,type="b",col="blue")
|
416
|
text(TMax_avgm[,1],TMax_avgm[,2],rownames(statistics_LST_s),cex=0.8,pos=2)
|
417
|
|
418
|
legend("topleft",legend=c("TMax","LST"), cex=1, col=c("red","blue"),
|
419
|
lty=1)
|
420
|
title(paste("Monthly mean tmax and LST at stations in Oregon", "2010",sep=" "))
|
421
|
savePlot(paste("Temporal_profile_res",id[i],out_prefix,".png", sep=""), type="png")
|
422
|
|
423
|
#Account for forest
|
424
|
data_month_all_forest<-data_month_all[data_month_all$LC1>=50,]
|
425
|
data_month_all_grass<-data_month_all[data_month_all$LC3>=10,]
|
426
|
LST_avgm_forest<-aggregate(LST~month,data=data_month_all_forest,mean)
|
427
|
LST_avgm_grass<-aggregate(LST~month,data=data_month_all_grass,mean)
|
428
|
|
429
|
plot(TMax_avgm,type="b",ylim=c(-7,42))
|
430
|
lines(LST_avgm,col="blue",add=T)
|
431
|
lines(LST_avgm_forest,col="green",add=T)
|
432
|
lines(LST_avgm_grass,col="red",add=T)
|
433
|
legend("topleft",legend=c("TMax","LST","LST_forest","LST_grass"), cex=0.8, col=c("black","blue","green","red"),
|
434
|
lty=1)
|
435
|
title(paste("Monthly average tmax for stations in Oregon ", "2010",sep=""))
|
436
|
|
437
|
savePlot(paste("Temporal_profile_res",id[i],out_prefix,".png", sep=""), type="png")
|
438
|
|
439
|
#Check where forest >50 and grass>10
|
440
|
#Create mask using land cover data
|
441
|
pos<-match("LC1",layerNames(s_raster)) #Find the layer which contains water bodies
|
442
|
LC1<-subset(s_raster,pos)
|
443
|
LC1[is.na(LC1)]<-0 #Since NA values are 0, we assign all zero to NA
|
444
|
LC1_50<-LC1>= 50
|
445
|
plot(LC1_50)
|
446
|
pos<-match("LC3",layerNames(s_raster)) #Find the layer which contains water bodies
|
447
|
LC3<-subset(s_raster,pos)
|
448
|
LC3[is.na(LC3)]<-0 #Since NA values are 0, we assign all zero to NA
|
449
|
LC3_10<-LC3>= 10
|
450
|
|
451
|
LC1_mean <- cellStats(LC3, mean)
|
452
|
LC1_mean <- cellStats(LC3, mean)
|
453
|
#
|
454
|
avl<-c(0,10,1,10,20,2,20,30,3,30,40,4,40,50,5,50,60,6,60,70,7,70,80,8,80,90,9,90,100,10)
|
455
|
rclmat<-matrix(avl,ncol=3,byrow=TRUE)
|
456
|
LC3_rec<-reclass(LC3,rclmat) #Loss of layer names when using reclass
|
457
|
tab_freq<-freq(LC3_rec)
|
458
|
|
459
|
crosstab(LC1_50,LC3_10) #Very little overalp between classes...
|
460
|
X11(12,12)
|
461
|
LC_rec<-stack(LC1_50,LC3_10)
|
462
|
layerNames(LC_rec)<-c("forest: LC1>50%","grass: LC3>10%")
|
463
|
plot(LC_rec,col=c("black","red"))
|
464
|
dev.off()
|
465
|
#legend("topleft",legend=c("fo","LST","LST_forest","LST_grass"), cex=0.8, col=c("black","blue","green","red"),
|
466
|
#lty=1)
|
467
|
|
468
|
###################### Visualization of results ######################
|
469
|
|
470
|
"tmax_predicted*20100103_30_1_365d_GAM_fusion_all_lstd_10242012.rst"
|
471
|
oldpath<-getwd()
|
472
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_GAM"
|
473
|
setwd(path)
|
474
|
date_selected<-"20100103"
|
475
|
#lf_fus<-list.files(pattern=paste("GAM_bias_tmax_predicted_mod8_",date_selected,"_30_1_365d_GAM_fusion_all_lstd_10242012.rst$",sep="")) #Search for files in relation to fusion
|
476
|
file_pat<-glob2rx(paste("*tmax_pred*",date_selected,"*_365d_GAM_fusion_all_lstd_10242012.rst",sep="")) #Search for files in relation to fusion
|
477
|
lf_fus1a<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
478
|
|
479
|
rast_fus1a<-stack(lf_fus1a)
|
480
|
rast_fus1a<-mask(rast_fus1a,mask_ELEV_SRTM)
|
481
|
|
482
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_07192012_GAM"
|
483
|
setwd(path)
|
484
|
date_selected<-"20100103"
|
485
|
#lf_fus<-list.files(pattern=paste("GAM_bias_tmax_predicted_mod8_",date_selected,"_30_1_365d_GAM_fusion_all_lstd_10242012.rst$",sep="")) #Search for files in relation to fusion
|
486
|
#lf_fus1<-list.files(pattern=paste("*.tmax_predicted.*.",date_selected,".*._365d_GAM_fusion_lstd_10062012.rst$",sep=""))
|
487
|
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_fusion_lstd_10062012.rst",sep="")) #Search for files in relation to fusion
|
488
|
lf_fus1s<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
489
|
|
490
|
rast_fus1s<-stack(lf_fus1s)
|
491
|
rast_fus1s<-mask(rast_fus1s,mask_ELEV_SRTM)
|
492
|
|
493
|
s_range<-c(minValue(rast_fusa1),maxValue(rast_fusa1)) #stack min and max
|
494
|
s_range<-c(min(s_range),max(s_range))
|
495
|
col_breaks <- pretty(s_range, n=50)
|
496
|
lab_breaks <- pretty(s_range, n=5)
|
497
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
498
|
X11(width=18,height=12)
|
499
|
par(mfrow=c(3,3))
|
500
|
for (k in 1:length(lf_fus1a)){
|
501
|
fus1a_r<-raster(rast_fus1a,k)
|
502
|
plot(fus1a_r, breaks=col_breaks, col=temp_colors(length(col_breaks)-1),
|
503
|
axis=list(at=lab_breaks, labels=lab_breaks))
|
504
|
}
|
505
|
|
506
|
#Wihtout setting range
|
507
|
s_range<-c(-12,18)
|
508
|
#col_breaks <- pretty(s_range, n=50)
|
509
|
#lab_breaks <- pretty(s_range, n=5)
|
510
|
col_breaks<-49
|
511
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
512
|
lab_breaks <- pretty(s_range, n=5)
|
513
|
X11(width=18,height=12)
|
514
|
par(mfrow=c(3,3))
|
515
|
mod_name<-paste("mod",1:8, sep="")
|
516
|
rname<-c("FUS_kr",mod_name)
|
517
|
for (k in 1:length(lf_fus1a)){
|
518
|
fus1a_r<-raster(rast_fus1a,k)
|
519
|
plot(fus1a_r, breaks=col_breaks, col=temp_colors(col_breaks-1),
|
520
|
axis=list(at=lab_breaks, labels=lab_breaks),main=rname[k])
|
521
|
}
|
522
|
|
523
|
#Wihtout setting range
|
524
|
s_range<-c(-12,23)
|
525
|
#col_breaks <- pretty(s_range, n=50)
|
526
|
#lab_breaks <- pretty(s_range, n=5)
|
527
|
col_breaks<-49
|
528
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
529
|
lab_breaks <- pretty(s_range, n=5)
|
530
|
X11(width=18,height=12)
|
531
|
par(mfrow=c(3,3))
|
532
|
mod_name<-paste("mod",1:8, sep="")
|
533
|
rname<-c("FUS_kr",mod_name)
|
534
|
for (k in 1:length(lf_fus1s)){
|
535
|
fus1s_r<-raster(rast_fus1s,k)
|
536
|
plot(fus1s_r, breaks=col_breaks, col=temp_colors(col_breaks-1),
|
537
|
axis=list(at=lab_breaks, labels=lab_breaks),main=rname[k])
|
538
|
}
|
539
|
|
540
|
|
541
|
|
542
|
|