1 |
5e7d95a7
|
Benoit Parmentier
|
##################################### METHODS COMPARISON ##########################################
|
2 |
3b090bfa
|
Benoit Parmentier
|
#################################### Spatial Analysis ########################################
|
3 |
5e7d95a7
|
Benoit Parmentier
|
#This script is not aimed at producing new interpolation surfaces. It utilizes the R ojbects created
|
4 |
|
|
# during the interpolation phase. #
|
5 |
|
|
# At this stage the script produces figures of various accuracy metrics and compare methods: #
|
6 |
3b090bfa
|
Benoit Parmentier
|
#- multisampling plots #
|
7 |
|
|
#- spatial accuracy in terms of distance to closest station #
|
8 |
5e7d95a7
|
Benoit Parmentier
|
#- spatial density of station network and accuracy metric
|
9 |
|
|
#- visualization of maps of prediction and difference for comparison
|
10 |
3b090bfa
|
Benoit Parmentier
|
#AUTHOR: Benoit Parmentier #
|
11 |
5e7d95a7
|
Benoit Parmentier
|
#DATE: 10/30/2012 #
|
12 |
|
|
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491 -- #
|
13 |
3b090bfa
|
Benoit Parmentier
|
###################################################################################################
|
14 |
|
|
|
15 |
|
|
###Loading R library and packages
|
16 |
|
|
#library(gtools) # loading some useful tools
|
17 |
|
|
library(mgcv) # GAM package by Wood 2006 (version 2012)
|
18 |
|
|
library(sp) # Spatial pacakge with class definition by Bivand et al. 2008
|
19 |
|
|
library(spdep) # Spatial package with methods and spatial stat. by Bivand et al. 2012
|
20 |
|
|
library(rgdal) # GDAL wrapper for R, spatial utilities (Keitt et al. 2012)
|
21 |
|
|
library(gstat) # Kriging and co-kriging by Pebesma et al. 2004
|
22 |
|
|
library(automap) # Automated Kriging based on gstat module by Hiemstra et al. 2008
|
23 |
|
|
library(spgwr)
|
24 |
|
|
library(gpclib)
|
25 |
|
|
library(maptools)
|
26 |
|
|
library(graphics)
|
27 |
|
|
library(parallel) # Urbanek S. and Ripley B., package for multi cores & parralel processing
|
28 |
|
|
library(raster)
|
29 |
|
|
library(rasterVis)
|
30 |
|
|
library(plotrix) #Draw circle on graph
|
31 |
|
|
library(reshape)
|
32 |
|
|
## Functions
|
33 |
|
|
#loading R objects that might have similar names
|
34 |
|
|
load_obj <- function(f)
|
35 |
|
|
{
|
36 |
|
|
env <- new.env()
|
37 |
|
|
nm <- load(f, env)[1]
|
38 |
|
|
env[[nm]]
|
39 |
|
|
}
|
40 |
|
|
|
41 |
|
|
###Parameters and arguments
|
42 |
|
|
|
43 |
|
|
infile1<- "ghcn_or_tmax_covariates_06262012_OR83M.shp" #GHCN shapefile containing variables for modeling 2010
|
44 |
|
|
#infile2<-"list_10_dates_04212012.txt" #List of 10 dates for the regression
|
45 |
|
|
infile2<-"list_365_dates_04212012.txt" #list of dates
|
46 |
|
|
infile3<-"LST_dates_var_names.txt" #LST dates name
|
47 |
|
|
infile4<-"models_interpolation_05142012.txt" #Interpolation model names
|
48 |
|
|
infile5<-"mean_day244_rescaled.rst" #mean LST for day 244
|
49 |
|
|
inlistf<-"list_files_05032012.txt" #list of raster images containing the Covariates
|
50 |
|
|
infile6<-"OR83M_state_outline.shp"
|
51 |
|
|
#stat_loc<-read.table(paste(path,"/","location_study_area_OR_0602012.txt",sep=""),sep=",", header=TRUE)
|
52 |
|
|
|
53 |
|
|
|
54 |
|
|
obj_list<-"list_obj_08262012.txt" #Results of fusion from the run on ATLAS
|
55 |
|
|
path<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison" #Jupiter LOCATION on Atlas for kriging #Jupiter Location on XANDERS
|
56 |
|
|
#path<-"/Users/benoitparmentier/Dropbox/Data/NCEAS/Oregon_covariates/" #Local dropbox folder on Benoit's laptop
|
57 |
|
|
setwd(path)
|
58 |
|
|
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";
|
59 |
|
|
out_prefix<-"methods_09262012_" #User defined output prefix
|
60 |
|
|
|
61 |
|
|
sampling_CAI<-load_obj("results2_CAI_sampling_obj_09132012_365d_GAM_CAI2_multisampling2.RData")
|
62 |
|
|
sampling_fus<-load_obj("results2_fusion_sampling_obj_10d_GAM_fusion_multisamp4_09192012.RData")
|
63 |
|
|
fus_CAI_mod<-load_obj("results2_CAI_Assessment_measure_all_09132012_365d_GAM_CAI2_multisampling2.RData")
|
64 |
|
|
gam_fus_mod1<-load_obj("results2_fusion_Assessment_measure_all_10d_GAM_fusion_multisamp4_09192012.RData")
|
65 |
|
|
|
66 |
|
|
filename<-sub(".shp","",infile1) #Removing the extension from file.
|
67 |
|
|
ghcn<-readOGR(".", filename) #reading shapefile
|
68 |
|
|
|
69 |
|
|
#CRS<-proj4string(ghcn) #Storing projection information (ellipsoid, datum,etc.)
|
70 |
|
|
lines<-read.table(paste(path,"/",inlistf,sep=""), sep="") #Column 1 contains the names of raster files
|
71 |
|
|
inlistvar<-lines[,1]
|
72 |
|
|
inlistvar<-paste(path,"/",as.character(inlistvar),sep="")
|
73 |
|
|
covar_names<-as.character(lines[,2]) #Column two contains short names for covaraites
|
74 |
|
|
|
75 |
|
|
s_raster<- stack(inlistvar) #Creating a stack of raster images from the list of variables.
|
76 |
|
|
layerNames(s_raster)<-covar_names #Assigning names to the raster layers
|
77 |
|
|
projection(s_raster)<-proj_str
|
78 |
|
|
|
79 |
|
|
#Create mask using land cover data
|
80 |
|
|
pos<-match("LC10",layerNames(s_raster)) #Find the layer which contains water bodies
|
81 |
|
|
LC10<-subset(s_raster,pos)
|
82 |
|
|
LC10[is.na(LC10)]<-0 #Since NA values are 0, we assign all zero to NA
|
83 |
|
|
mask_land<-LC10<100 #All values below 100% water are assigned the value 1, value 0 is "water"
|
84 |
|
|
mask_land_NA<-mask_land
|
85 |
|
|
mask_land_NA[mask_land_NA==0]<-NA #Water bodies are assigned value 1
|
86 |
|
|
|
87 |
|
|
data_name<-"mask_land_OR"
|
88 |
|
|
raster_name<-paste(data_name,".rst", sep="")
|
89 |
|
|
writeRaster(mask_land, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI)
|
90 |
|
|
#writeRaster(r2, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI)
|
91 |
|
|
|
92 |
|
|
pos<-match("mm_01",layerNames(s_raster))
|
93 |
|
|
mm_01<-subset(s_raster,pos)
|
94 |
|
|
mm_01<-mm_01-273.15
|
95 |
|
|
mm_01<-mask(mm_01,mask_land_NA) #Raster image used as backround
|
96 |
|
|
#mention this is the last... files
|
97 |
|
|
|
98 |
|
|
### RESULTS COMPARISON
|
99 |
|
|
|
100 |
|
|
### CODE BEGIN #####
|
101 |
|
|
|
102 |
|
|
### PART I MULTISAMPLING COMPARISON ####
|
103 |
|
|
|
104 |
|
|
tb_diagnostic2<-sampling_CAI$tb #Extracting the accuracy metric information...
|
105 |
|
|
tb_diagnostic<-sampling_fus$tb
|
106 |
|
|
|
107 |
|
|
tb_diagnostic[["prop"]]<-as.factor(tb_diagnostic[["prop"]])
|
108 |
|
|
tb_diagnostic2[["prop"]]<-as.factor(tb_diagnostic2[["prop"]])
|
109 |
|
|
|
110 |
|
|
#Preparing the data for the plot
|
111 |
|
|
#fus data
|
112 |
|
|
t<-melt(tb_diagnostic,
|
113 |
|
|
measure=c("mod1","mod2","mod3","mod4", "mod5", "mod6", "mod7", "mod8","mod9"),
|
114 |
|
|
id=c("dates","metric","prop"),
|
115 |
|
|
na.rm=F)
|
116 |
|
|
avg_tb<-cast(t,metric+prop~variable,mean)
|
117 |
|
|
sd_tb<-cast(t,metric+prop~variable,sd)
|
118 |
|
|
n_tb<-cast(t,metric+prop~variable,length)
|
119 |
|
|
avg_tb[["prop"]]<-as.numeric(as.character(avg_tb[["prop"]]))
|
120 |
|
|
avg_RMSE<-subset(avg_tb,metric=="RMSE")
|
121 |
|
|
|
122 |
|
|
#CAI data
|
123 |
|
|
t2<-melt(tb_diagnostic2,
|
124 |
|
|
measure=c("mod1","mod2","mod3","mod4", "mod5", "mod6", "mod7", "mod8","mod9"),
|
125 |
|
|
id=c("dates","metric","prop"),
|
126 |
|
|
na.rm=F)
|
127 |
|
|
avg_tb2<-cast(t2,metric+prop~variable,mean)
|
128 |
|
|
sd_tb2<-cast(t2,metric+prop~variable,sd)
|
129 |
|
|
n_tb2<-cast(t2,metric+prop~variable,length)
|
130 |
|
|
avg_tb2[["prop"]]<-as.numeric(as.character(avg_tb2[["prop"]]))
|
131 |
|
|
avg_RMSE2<-subset(avg_tb2,metric=="RMSE")
|
132 |
|
|
|
133 |
|
|
#Select only information related to FUSION
|
134 |
|
|
|
135 |
|
|
x<-avg_RMSE[["prop"]]
|
136 |
|
|
i=9
|
137 |
|
|
mod_name<-paste("mod",i,sep="")
|
138 |
|
|
y<-avg_RMSE[[mod_name]]
|
139 |
|
|
|
140 |
|
|
sd_tb_RMSE <- subset(sd_tb, metric=="RMSE")
|
141 |
|
|
x_sd<-sd_tb_RMSE[["prop"]]
|
142 |
|
|
i=9
|
143 |
|
|
mod_name<-paste("mod",i,sep="")
|
144 |
|
|
y_sd<-sd_tb_RMSE[[mod_name]]
|
145 |
|
|
|
146 |
|
|
#Select only information related to CAI
|
147 |
|
|
|
148 |
|
|
x2<-avg_RMSE2[["prop"]]
|
149 |
|
|
i=9
|
150 |
|
|
mod_name<-paste("mod",i,sep="")
|
151 |
|
|
y2<-avg_RMSE2[[mod_name]]
|
152 |
|
|
|
153 |
|
|
sd_tb_RMSE2 <- subset(sd_tb2, metric=="RMSE")
|
154 |
|
|
x_sd2<-sd_tb_RMSE2[["prop"]]
|
155 |
|
|
i=9
|
156 |
|
|
mod_name<-paste("mod",i,sep="")
|
157 |
|
|
y_sd2<-sd_tb_RMSE2[[mod_name]]
|
158 |
|
|
|
159 |
|
|
n=150
|
160 |
|
|
ciw <- qt(0.975, n) * y_sd / sqrt(n)
|
161 |
|
|
ciw2 <- qt(0.975, n) * y_sd2 / sqrt(n)
|
162 |
|
|
|
163 |
|
|
#Comparison of MAE for different proportions for FUSION and CAI using CI
|
164 |
|
|
X11()
|
165 |
|
|
plotCI(y=y, x=x, uiw=ciw, col="red", main=" FUS: RMSE proportion of validation hold out", barcol="blue", lwd=1,
|
166 |
|
|
ylab="RMSE (C)", xlab="Proportions of validation hold out (in %)")
|
167 |
|
|
lines(x,y,col="red")
|
168 |
|
|
legend("bottomright",legend=c("fus"), cex=1.2, col=c("red"),
|
169 |
|
|
lty=1, title="RMSE")
|
170 |
|
|
savePlot(paste("Comparison_multisampling_fus_RMSE_CI",out_prefix,".png", sep=""), type="png")
|
171 |
|
|
|
172 |
|
|
plotCI(y=y2, x=x2, uiw=ciw2, col="black", main=" CAI: RMSE proportion of validation hold out", barcol="blue", lwd=1,
|
173 |
|
|
ylab="RMSE (C)", xlab="Proportions of validation hold out (in %)")
|
174 |
|
|
lines(x2,y2,col="grey")
|
175 |
|
|
legend("bottomright",legend=c("CAI"), cex=1.2, col=c("grey"),
|
176 |
|
|
lty=1, title="RMSE")
|
177 |
|
|
savePlot(paste("Comparison_multisampling_CAI_RMSE_CI",out_prefix,".png", sep=""), type="png")
|
178 |
|
|
dev.off()
|
179 |
|
|
|
180 |
|
|
#Comparison of MAE for different proportions for FUSION and CAI
|
181 |
|
|
X11()
|
182 |
|
|
plot(x,y,col="red",type="b", ylab="RMSE (C)", xlab="Proportions of validation hold out (in %)")
|
183 |
|
|
lines(x2,y2,col="grey")
|
184 |
|
|
points(x2,y2,col="grey")
|
185 |
|
|
title("MAE in terms of proportions and random sampling")
|
186 |
|
|
legend("bottomright",legend=c("fus","CAI"), cex=1.2, col=c("red","grey"),
|
187 |
|
|
lty=1, title="RMSE")
|
188 |
|
|
savePlot(paste("Comparison_multisampling_fus_CAI_RMSE_averages",out_prefix,".png", sep=""), type="png")
|
189 |
|
|
dev.off()
|
190 |
|
|
|
191 |
5e7d95a7
|
Benoit Parmentier
|
############################################################################
|
192 |
|
|
#### PART II EXAMINIG PREDICTIONS AND RESIDUALS TEMPORAL PROFILES ##########
|
193 |
3b090bfa
|
Benoit Parmentier
|
|
194 |
|
|
l_f<-list.files(pattern="*tmax_predicted.*fusion5.rst$") #Search for files in relation to fusion
|
195 |
|
|
l_f2<-list.files(pattern="CAI_tmax_predicted.*_GAM_CAI2.rst$")
|
196 |
|
|
inlistpred<-paste(path,"/",as.character(l_f),sep="")
|
197 |
|
|
inlistpred2<-paste(path,"/",as.character(l_f2),sep="")
|
198 |
|
|
|
199 |
|
|
fus_rast<- stack(inlistpred) #Creating a stack of raster images from the list of variables.
|
200 |
|
|
cai_rast<- stack(inlistpred2) #Creating a stack of raster images from the list of variables.
|
201 |
|
|
|
202 |
|
|
id<-unique(ghcn$station)
|
203 |
|
|
ghcn_id<-as.data.frame(subset(ghcn,select=c("station","x_OR83M","y_OR83M")))
|
204 |
|
|
|
205 |
|
|
ghcn_melt<-melt(ghcn_id,
|
206 |
|
|
measure=c("x_OR83M","y_OR83M"),
|
207 |
|
|
id=c("station"),
|
208 |
|
|
na.rm=F)
|
209 |
|
|
|
210 |
|
|
ghcn_cast<-cast(ghcn_melt,station~variable,mean)
|
211 |
|
|
ghcn_locs<-as.data.frame(ghcn_cast)
|
212 |
|
|
|
213 |
|
|
coords<- ghcn_locs[,c('x_OR83M','y_OR83M')]
|
214 |
|
|
coordinates(ghcn_locs)<-coords
|
215 |
|
|
proj4string(ghcn_locs)<-proj_str #Need to assign coordinates...
|
216 |
|
|
|
217 |
|
|
tmp<-extract(fus_rast,ghcn_locs)
|
218 |
|
|
tmp2<-extract(cai_rast,ghcn_locs)
|
219 |
|
|
|
220 |
|
|
tmp_names<-paste("fusd",seq(1,365),sep="")
|
221 |
|
|
colnames(tmp)<-tmp_names
|
222 |
|
|
tmp_names<-paste("caid",seq(1,365),sep="")
|
223 |
|
|
colnames(tmp2)<-tmp_names
|
224 |
|
|
ghcn_fus_pred<-cbind(as.data.frame(ghcn_locs),as.data.frame(tmp))
|
225 |
|
|
ghcn_cai_pred<-cbind(as.data.frame(ghcn_locs),as.data.frame(tmp2))
|
226 |
|
|
|
227 |
|
|
write.table(ghcn_fus_pred,file="extract3_fus_y2010.txt",sep=",")
|
228 |
|
|
write.table(ghcn_cai_pred,file="extract3_cai_y2010.txt",sep=",")
|
229 |
|
|
|
230 |
|
|
ghcn$value[ghcn$value< -150 | ghcn$value>400]<-NA #screenout values out of range
|
231 |
|
|
ghcn$value<-ghcn$value/10
|
232 |
|
|
ghcn_m<-melt(as.data.frame(ghcn),
|
233 |
|
|
measure=c("value"),
|
234 |
|
|
id=c("station","date"),
|
235 |
|
|
na.rm=F)
|
236 |
|
|
|
237 |
|
|
ghcn_mc<-cast(ghcn_m,station~date~variable,mean) #This creates an array of dimension 186,366,1
|
238 |
|
|
|
239 |
|
|
ghcn_value<-as.data.frame(ghcn_mc[,,1])
|
240 |
5e7d95a7
|
Benoit Parmentier
|
ghcn_value<-cbind(ghcn_locs,ghcn_value[,1:365]) #This data frame contains values for 365 days
|
241 |
|
|
#for 186 stations of year 2010...
|
242 |
3b090bfa
|
Benoit Parmentier
|
write.table(ghcn_value,na="",file="extract3_dailyTmax_y2010.txt",sep=",")
|
243 |
|
|
|
244 |
|
|
id<-c("USW00094261","USW00004141","USC00356252","USC00357208")
|
245 |
5e7d95a7
|
Benoit Parmentier
|
#id<-c("USW00024284","USC00354126","USC00358536","USC00354835",
|
246 |
3b090bfa
|
Benoit Parmentier
|
"USC00356252","USC00359316","USC00358246","USC00350694",
|
247 |
|
|
"USC00350699","USW00024230","USC00353542")
|
248 |
|
|
|
249 |
|
|
m<-match(id,ghcn_locs$station)
|
250 |
|
|
dat_id<-ghcn_locs[m,] #creating new subset
|
251 |
|
|
#dat_id<-subset(ghcn_locs[gj])
|
252 |
|
|
|
253 |
|
|
filename<-sub(".shp","",infile6) #Removing the extension from file.
|
254 |
|
|
reg_outline<-readOGR(".", filename) #reading shapefile
|
255 |
|
|
X11()
|
256 |
|
|
s.range <- c(min(minValue(mm_01)), max(maxValue(mm_01)))
|
257 |
|
|
col.breaks <- pretty(s.range, n=50)
|
258 |
|
|
lab.breaks <- pretty(s.range, n=5)
|
259 |
|
|
temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
|
260 |
|
|
plot(mm_01, breaks=col.breaks, col=temp.colors(length(col.breaks)-1),
|
261 |
|
|
axis=list(at=lab.breaks, labels=lab.breaks))
|
262 |
|
|
plot(reg_outline, add=TRUE)
|
263 |
|
|
plot(dat_id,cex=1.5,add=TRUE)
|
264 |
|
|
title("Selected stations for comparison",line=3)
|
265 |
|
|
title("(Background: mean January LST)", cex=0.5, line=2)
|
266 |
|
|
coords<-coordinates(dat_id)
|
267 |
|
|
text(x=coords[,1],y=coords[,2],labels=id,cex=0.8, adj=c(0,1),offset=2) #c(0,1) for lower right corner!
|
268 |
|
|
savePlot(paste("temporal_profile_station_locations_map",out_prefix,".png", sep=""), type="png")
|
269 |
|
|
dev.off()
|
270 |
|
|
|
271 |
5e7d95a7
|
Benoit Parmentier
|
stat_list<-vector("list",3 ) #list containing the selected stations...
|
272 |
3b090bfa
|
Benoit Parmentier
|
stat_list[[1]]<-ghcn_fus_pred
|
273 |
|
|
stat_list[[2]]<-ghcn_cai_pred
|
274 |
|
|
stat_list[[3]]<-ghcn_value
|
275 |
|
|
ac_temp<-matrix(NA,length(id),2)
|
276 |
|
|
|
277 |
|
|
#id<-ghcn_value$station #if runinng on all the station...
|
278 |
|
|
for (i in 1:length(id)){
|
279 |
|
|
m1<-match(id[i],ghcn_fus_pred$station)
|
280 |
|
|
m2<-match(id[i],ghcn_cai_pred$station)
|
281 |
|
|
m3<-match(id[i],ghcn_value$station)
|
282 |
5e7d95a7
|
Benoit Parmentier
|
y1<-as.numeric(ghcn_fus_pred[m1,6:ncol(ghcn_fus_pred)]) #vector containing fusion time series of predictecd tmax
|
283 |
|
|
y2<-as.numeric(ghcn_cai_pred[m2,6:ncol(ghcn_cai_pred)]) #vector containing CAI time series of predictecd
|
284 |
|
|
y3<-as.numeric(ghcn_value[m3,6:ncol(ghcn_value)]) #vector containing observed time series of predictecd
|
285 |
|
|
res2<-y2-y3 #CAI time series residuals
|
286 |
|
|
res1<-y1-y3 #fusion time series residuals
|
287 |
3b090bfa
|
Benoit Parmentier
|
x<-1:365
|
288 |
|
|
X11(6,15)
|
289 |
5e7d95a7
|
Benoit Parmentier
|
plot(x,y1,type="l",col="red",ylab="tmax (C)",xlab="Day of year")
|
290 |
3b090bfa
|
Benoit Parmentier
|
lines(x,y2,col="blue")
|
291 |
5e7d95a7
|
Benoit Parmentier
|
lines(x,y3,col="green")
|
292 |
3b090bfa
|
Benoit Parmentier
|
title(paste("temporal profile for station ", id[i],sep=""))
|
293 |
|
|
# add a legend
|
294 |
5e7d95a7
|
Benoit Parmentier
|
legend("topright",legend=c("fus","CAI","OBS"), cex=1.2, col=c("red","blue","green"),
|
295 |
3b090bfa
|
Benoit Parmentier
|
lty=1, title="tmax")
|
296 |
|
|
savePlot(paste("Temporal_profile_",id[i],out_prefix,".png", sep=""), type="png")
|
297 |
5e7d95a7
|
Benoit Parmentier
|
|
298 |
|
|
### RESIDUALS PLOT
|
299 |
3b090bfa
|
Benoit Parmentier
|
zero<-rep(0,365)
|
300 |
5e7d95a7
|
Benoit Parmentier
|
plot(x,res2,type="l",col="blue", ylab="tmax (C)",xlab="Day of year") #res2 contains residuals from cai
|
301 |
|
|
lines(x,res1,col="red") #res1 contains fus
|
302 |
|
|
lines(x,zero,col="green")
|
303 |
|
|
legend("topright",legend=c("fus","CAI"), cex=1.2, col=c("red","blue"),
|
304 |
|
|
lty=1)
|
305 |
|
|
title(paste("temporal profile for station ", id[i],sep=""))
|
306 |
|
|
|
307 |
3b090bfa
|
Benoit Parmentier
|
savePlot(paste("Temporal_profile_res",id[i],out_prefix,".png", sep=""), type="png")
|
308 |
|
|
|
309 |
|
|
ac_temp[i,1]<-mean(abs(res1),na.rm=T)
|
310 |
|
|
ac_temp[i,2]<-mean(abs(res2),na.rm=T)
|
311 |
|
|
dev.off()
|
312 |
|
|
}
|
313 |
|
|
ac_temp<-as.data.frame(ac_temp)
|
314 |
|
|
ac_temp$station<-id
|
315 |
|
|
names(ac_temp)<-c("fus","CAI","station") #ac_temp contains the MAE per station
|
316 |
|
|
|
317 |
5e7d95a7
|
Benoit Parmentier
|
### RESIDUALS FOR EVERY STATION...############
|
318 |
|
|
|
319 |
3b090bfa
|
Benoit Parmentier
|
id<-ghcn_value$station #if runinng on all the station...
|
320 |
|
|
ac_temp2<-matrix(NA,length(id),2)
|
321 |
|
|
|
322 |
|
|
for (i in 1:length(id)){
|
323 |
|
|
m1<-match(id[i],ghcn_fus_pred$station)
|
324 |
|
|
m2<-match(id[i],ghcn_cai_pred$station)
|
325 |
|
|
m3<-match(id[i],ghcn_value$station)
|
326 |
|
|
y1<-as.numeric(ghcn_fus_pred[m1,6:ncol(ghcn_fus_pred)])
|
327 |
|
|
y2<-as.numeric(ghcn_cai_pred[m2,6:ncol(ghcn_cai_pred)])
|
328 |
|
|
y3<-as.numeric(ghcn_value[m3,6:ncol(ghcn_value)])
|
329 |
|
|
res2<-y2-y3
|
330 |
|
|
res1<-y1-y3
|
331 |
|
|
ac_temp2[i,1]<-mean(abs(res1),na.rm=T)
|
332 |
|
|
ac_temp2[i,2]<-mean(abs(res2),na.rm=T)
|
333 |
|
|
}
|
334 |
|
|
|
335 |
5e7d95a7
|
Benoit Parmentier
|
|
336 |
3b090bfa
|
Benoit Parmentier
|
ac_temp2<-as.data.frame(ac_temp2)
|
337 |
|
|
ac_temp2$station<-id
|
338 |
|
|
names(ac_temp2)<-c("fus","CAI","station")
|
339 |
|
|
|
340 |
|
|
ac_temp2<-ac_temp2[order(ac_temp2$fus,ac_temp2$CAI), ]
|
341 |
|
|
ghcn_MAE<-merge(ghcn_locs,ac_temp2,by.x=station,by.y=station)
|
342 |
|
|
|
343 |
5e7d95a7
|
Benoit Parmentier
|
########### TRANSECT-- DAY OF YEAR PLOT...#########
|
344 |
|
|
|
345 |
|
|
id<-c("USW00024284","USC00354126","USC00358536","USC00354835",
|
346 |
|
|
"USC00356252","USC00359316","USC00358246","USC00350694",
|
347 |
|
|
"USC00350699","USW00024230","USC00353542")
|
348 |
|
|
id_order<-1:11
|
349 |
|
|
m<-match(id,ghcn_locs$station)
|
350 |
|
|
dat_id<-ghcn_locs[m,] #creating new subset
|
351 |
|
|
#dat_id<-subset(ghcn_locs[gj])
|
352 |
|
|
|
353 |
|
|
filename<-sub(".shp","",infile6) #Removing the extension from file.
|
354 |
|
|
reg_outline<-readOGR(".", filename) #reading shapefile
|
355 |
|
|
X11()
|
356 |
|
|
s.range <- c(min(minValue(mm_01)), max(maxValue(mm_01)))
|
357 |
|
|
col.breaks <- pretty(s.range, n=50)
|
358 |
|
|
lab.breaks <- pretty(s.range, n=5)
|
359 |
|
|
temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
|
360 |
|
|
plot(mm_01, breaks=col.breaks, col=temp.colors(length(col.breaks)-1),
|
361 |
|
|
axis=list(at=lab.breaks, labels=lab.breaks))
|
362 |
|
|
plot(reg_outline, add=TRUE)
|
363 |
|
|
plot(dat_id,cex=1.5,add=TRUE)
|
364 |
|
|
title("Selected stations for comparison",line=3)
|
365 |
|
|
title("(Background: mean January LST)", cex=0.5, line=2)
|
366 |
|
|
coords<-coordinates(dat_id)
|
367 |
|
|
text(x=coords[,1],y=coords[,2],labels=as.character(id_order),cex=1.5, adj=c(0,1),offset=2) #c(0,1) for lower right corner!
|
368 |
|
|
savePlot(paste("temporal_profile_station_locations_map",out_prefix,".png", sep=""), type="png")
|
369 |
|
|
dev.off()
|
370 |
|
|
|
371 |
3b090bfa
|
Benoit Parmentier
|
m1<-match(id,ghcn_fus_pred$station)
|
372 |
|
|
m2<-match(id,ghcn_cai_pred$station)
|
373 |
|
|
m3<-match(id,ghcn_value$station)
|
374 |
|
|
ghcn_dsub<-subset(ghcn,ghcn$station==id)
|
375 |
|
|
all.equal(m1,m2,m3) #OK order is the same
|
376 |
5e7d95a7
|
Benoit Parmentier
|
date_selection<-c("01-01-2010","01-09-2010")
|
377 |
3b090bfa
|
Benoit Parmentier
|
#date_str<-gsub(date_selection,"-","")
|
378 |
5e7d95a7
|
Benoit Parmentier
|
date_str<-c("20100101","20100901")
|
379 |
3b090bfa
|
Benoit Parmentier
|
covar_dsub<-subset(ghcn_dsub,ghcn_dsub$date==date_str[i],select=c("station","ELEV_SRTM","LC1"))
|
380 |
|
|
|
381 |
|
|
date_pred<-as.Date(date_selection)
|
382 |
|
|
#mo<-as.integer(strftime(date_pred, "%m")) # current month of the date being processed
|
383 |
|
|
doy_pred<-(strptime(date_pred, "%d-%m-%Y")$yday+1)
|
384 |
|
|
|
385 |
|
|
for (i in 1:length(date_pred)){
|
386 |
|
|
doy<-doy_pred[i]+5 #column label
|
387 |
5e7d95a7
|
Benoit Parmentier
|
doy<-243+5
|
388 |
3b090bfa
|
Benoit Parmentier
|
stat_subset<-cbind(id,ghcn_fus_pred[m1,doy],ghcn_cai_pred[m1,doy],ghcn_value[m1,doy])
|
389 |
|
|
colnames(stat_subset)<-c("station","fus","cai","value")
|
390 |
|
|
stat_subset<-as.data.frame(stat_subset)
|
391 |
|
|
for(j in 2:4){ # start of the for loop #1
|
392 |
|
|
stat_subset[,j]<-as.numeric(as.character(stat_subset[,j]))
|
393 |
|
|
}
|
394 |
|
|
X11()
|
395 |
5e7d95a7
|
Benoit Parmentier
|
plot(1:11,stat_subset$value,type="b",col="green",ylab="tmax",xlab="station transtect number")
|
396 |
3b090bfa
|
Benoit Parmentier
|
# xlabels())
|
397 |
5e7d95a7
|
Benoit Parmentier
|
lines(1:11,stat_subset$fus,type="b",col="red")
|
398 |
|
|
lines(1:11,stat_subset$cai,type="b",col="blue")
|
399 |
3b090bfa
|
Benoit Parmentier
|
|
400 |
5e7d95a7
|
Benoit Parmentier
|
legend("bottomright",legend=c("obs","fus","cai"), cex=1.2, col=c("green","red","blue"),
|
401 |
3b090bfa
|
Benoit Parmentier
|
lty=1, title="tmax")
|
402 |
|
|
title(paste("Daily tmax prediction ",date_selection[i],sep=" "))
|
403 |
|
|
savePlot(paste("transect_profile_tmax_",date_str[i],out_prefix,".png", sep=""), type="png")
|
404 |
|
|
dev.off()
|
405 |
|
|
}
|
406 |
|
|
|
407 |
5e7d95a7
|
Benoit Parmentier
|
##############################################
|
408 |
|
|
########## USING TEMPORAL IMAGES...############
|
409 |
3b090bfa
|
Benoit Parmentier
|
|
410 |
|
|
date_list<-vector("list", length(l_f))
|
411 |
|
|
for (k in 1:length(l_f)){
|
412 |
|
|
tmp<-(unlist(strsplit(l_f[k],"_"))) #spliting file name to obtain the prediction date
|
413 |
|
|
date_list[k]<-tmp[4]
|
414 |
|
|
}
|
415 |
|
|
|
416 |
|
|
date_list2<-vector("list", length(l_f2))
|
417 |
|
|
for (k in 1:length(l_f2)){
|
418 |
|
|
tmp<-(unlist(strsplit(l_f2[k],"_"))) #spliting file name to obtain the prediction date
|
419 |
|
|
date_list2[k]<-tmp[4]
|
420 |
|
|
}
|
421 |
|
|
|
422 |
|
|
setdiff(date_list,date_list2)
|
423 |
|
|
all.equal(date_list,date_list2) #This checks that both lists are equals
|
424 |
|
|
|
425 |
|
|
nel<-length(gam_fus_mod1)
|
426 |
|
|
list_fus_data_s<-vector("list", nel)
|
427 |
|
|
list_cai_data_s<-vector("list", nel)
|
428 |
|
|
list_fus_data_v<-vector("list", nel)
|
429 |
|
|
list_cai_data_v<-vector("list", nel)
|
430 |
|
|
|
431 |
|
|
list_fus_data<-vector("list", nel)
|
432 |
|
|
list_cai_data<-vector("list", nel)
|
433 |
|
|
|
434 |
|
|
list_dstspat_er<-vector("list", nel)
|
435 |
|
|
list_dstspat_er2<-vector("list", nel)
|
436 |
|
|
k=1
|
437 |
|
|
|
438 |
|
|
for (k in 1:nel){
|
439 |
|
|
#for (k in 1:365){
|
440 |
|
|
|
441 |
|
|
#Start loop over the full year!!!
|
442 |
|
|
names(gam_fus_mod1[[k]])
|
443 |
|
|
data_s<-gam_fus_mod1[[k]]$data_s
|
444 |
|
|
data_v<-gam_fus_mod1[[k]]$data_v
|
445 |
|
|
|
446 |
|
|
date_proc<-unique(data_s$date)
|
447 |
|
|
index<-match(as.character(date_proc),unlist(date_list)) #find the correct date..
|
448 |
|
|
#raster_pred<-raster(rp_raster,index)
|
449 |
|
|
|
450 |
|
|
#####second series added
|
451 |
|
|
data_v2<-fus_CAI_mod[[k]]$data_v
|
452 |
|
|
data_s2<-fus_CAI_mod[[k]]$data_s
|
453 |
|
|
|
454 |
|
|
date_proc<-unique(data_s$date)
|
455 |
|
|
index<-match(as.character(date_proc),unlist(date_list)) #find the correct date..
|
456 |
|
|
#raster_pred<-raster(rp_raster,index)
|
457 |
|
|
|
458 |
|
|
###Checking if training and validation have the same columns
|
459 |
|
|
nd<-setdiff(names(data_s),names(data_v))
|
460 |
|
|
nd2<-setdiff(names(data_s2),names(data_v2))
|
461 |
|
|
|
462 |
|
|
data_v[[nd]]<-NA #daily_delta is not the same
|
463 |
|
|
|
464 |
|
|
data_v$training<-rep(0,nrow(data_v))
|
465 |
|
|
data_v2$training<-rep(0,nrow(data_v2))
|
466 |
|
|
data_s$training<-rep(1,nrow(data_s))
|
467 |
|
|
data_s2$training<-rep(1,nrow(data_s2))
|
468 |
|
|
|
469 |
|
|
list_fus_data_s[[k]]<-data_s
|
470 |
|
|
list_cai_data_s[[k]]<-data_s2
|
471 |
|
|
list_fus_data_v[[k]]<-data_v
|
472 |
|
|
list_cai_data_v[[k]]<-data_v2
|
473 |
|
|
list_fus_data[[k]]<-rbind(data_v,data_s)
|
474 |
|
|
list_cai_data[[k]]<-rbind(data_v2,data_s2)
|
475 |
|
|
|
476 |
|
|
d_s_v<-matrix(0,nrow(data_v),nrow(data_s))
|
477 |
|
|
for(i in 1:nrow(data_s)){
|
478 |
|
|
pt<-data_s[i,]
|
479 |
|
|
d_pt<-(spDistsN1(data_v,pt,longlat=FALSE))/1000 #Distance to stataion i in km
|
480 |
|
|
d_s_v[,i]<-d_pt
|
481 |
|
|
}
|
482 |
|
|
|
483 |
|
|
d_s_v2<-matrix(0,nrow(data_v2),nrow(data_s2))
|
484 |
|
|
for(i in 1:nrow(data_s2)){
|
485 |
|
|
pt2<-data_s2[i,]
|
486 |
|
|
d_pt2<-(spDistsN1(data_v2,pt2,longlat=FALSE))/1000 #Distance to stataion i in km
|
487 |
|
|
d_s_v2[,i]<-d_pt2
|
488 |
|
|
}
|
489 |
|
|
|
490 |
|
|
#Create data.frame with position, ID, dst and residuals...YOU HAVE TO DO IT SEPARATELY FOR EACH MODEL!!!
|
491 |
|
|
#Do first fusion and then CAI
|
492 |
|
|
pos<-vector("numeric",nrow(data_v))
|
493 |
|
|
y<-vector("numeric",nrow(data_v))
|
494 |
|
|
dst<-vector("numeric",nrow(data_v))
|
495 |
|
|
|
496 |
|
|
pos2<-vector("numeric",nrow(data_v))
|
497 |
|
|
y2<-vector("numeric",nrow(data_v))
|
498 |
|
|
dst2<-vector("numeric",nrow(data_v))
|
499 |
|
|
|
500 |
|
|
for (i in 1:nrow(data_v)){
|
501 |
|
|
pos[i]<-match(min(d_s_v[i,]),d_s_v[i,])
|
502 |
|
|
dst[i]<-min(d_s_v[i,])
|
503 |
|
|
}
|
504 |
|
|
|
505 |
|
|
for (i in 1:nrow(data_v2)){
|
506 |
|
|
pos2[i]<-match(min(d_s_v2[i,]),d_s_v2[i,])
|
507 |
|
|
dst2[i]<-min(d_s_v2[i,])
|
508 |
|
|
}
|
509 |
|
|
|
510 |
|
|
res_fus<-data_v$res_mod9
|
511 |
|
|
res_CAI<-data_v2$res_mod9
|
512 |
|
|
|
513 |
|
|
dstspat_er<-as.data.frame(cbind(as.vector(data_v$id),as.vector(data_s$id[pos]),pos, data_v$lat, data_v$lon, data_v$x_OR83M,data_v$y_OR83M,
|
514 |
|
|
dst,
|
515 |
|
|
res_fus))
|
516 |
|
|
dstspat_er2<-as.data.frame(cbind(as.vector(data_v2$id),as.vector(data_s2$id[pos2]),pos2, data_v2$lat, data_v2$lon, data_v2$x_OR83M,data_v2$y_OR83M,
|
517 |
|
|
dst2,
|
518 |
|
|
res_CAI))
|
519 |
|
|
names(dstspat_er2)[1:7]<-c("v_id","s_id","pos","lat","lon","x_OR83M","y_OR83M")
|
520 |
|
|
|
521 |
|
|
names(dstspat_er)[1:7]<-c("v_id","s_id","pos","lat","lon","x_OR83M","y_OR83M")
|
522 |
|
|
# names(dstspat_er)[10:15]<-c("res_mod1","res_mod2","res_mod3","res_mod4","res_mod5","res_CAI")
|
523 |
|
|
#names(dstspat_er)[10:15]<-c("res_mod1","res_mod2","res_mod3","res_mod4","res_mod5","res_CAI")
|
524 |
|
|
list_dstspat_er[[k]]<-dstspat_er
|
525 |
|
|
list_dstspat_er2[[k]]<-dstspat_er2
|
526 |
|
|
|
527 |
|
|
}
|
528 |
|
|
save(list_dstspat_er,file="spat1_ac5.RData")
|
529 |
|
|
save(list_dstspat_er2,file="spat2_ac5.RData")
|
530 |
|
|
|
531 |
|
|
#obj_tmp2<-load_obj("spat_ac4.RData")
|
532 |
|
|
save(list_fus_data,file="list_fus_data_combined.RData")
|
533 |
|
|
save(list_cai_data,file="list_cai_data_combined.RData")
|
534 |
|
|
|
535 |
|
|
save(list_fus_data_s,file="list_fus_data_s_combined.RData")
|
536 |
|
|
save(list_cai_data_s,file="list_cai_data_s_combined.RData")
|
537 |
|
|
save(list_fus_data_v,file="list_fus_data_v_combined.RData")
|
538 |
|
|
save(list_cai_data_v,file="list_cai_data_v_combined.RData")
|
539 |
|
|
|
540 |
|
|
for (k in 1:nel){
|
541 |
|
|
data_s<-as.data.frame(list_fus_data_s[[k]])
|
542 |
|
|
data_v<-as.data.frame(list_fus_data_v[[k]])
|
543 |
|
|
list_fus_data[[k]]<-rbind(data_s,data_v)
|
544 |
|
|
data_s2<-as.data.frame(list_cai_data_s[[k]])
|
545 |
|
|
data_v2<-as.data.frame(list_cai_data_v[[k]])
|
546 |
|
|
list_cai_data[[k]]<-rbind(data_s2,data_v2)
|
547 |
|
|
}
|
548 |
|
|
data_fus<-do.call(rbind.fill,list_fus_data)
|
549 |
|
|
data_cai<-do.call(rbind.fill,list_cai_data)
|
550 |
|
|
|
551 |
|
|
data_fus_melt<-melt(data_fus,
|
552 |
|
|
measure=c("x_OR83M","y_OR83M","res_fus","res_mod1","res_mod2","res_mod3","res_mod4","res_mod5","pred_fus","dailyTmax","TMax","LST","training"),
|
553 |
|
|
id=c("id","date"),
|
554 |
|
|
na.rm=F)
|
555 |
|
|
data_fus_cast<-cast(data_fus_melt,id+date~variable,mean)
|
556 |
|
|
|
557 |
|
|
test_dst<-list_dstspat_er
|
558 |
|
|
test<-do.call(rbind,list_dstspat_er)
|
559 |
|
|
|
560 |
|
|
test_dst2<-list_dstspat_er2
|
561 |
|
|
test2<-do.call(rbind,list_dstspat_er2)
|
562 |
|
|
|
563 |
|
|
for(i in 4:ncol(test)){ # start of the for loop #1
|
564 |
|
|
test[,i]<-as.numeric(as.character(test[,i]))
|
565 |
|
|
}
|
566 |
|
|
|
567 |
|
|
for(i in 4:ncol(test2)){ # start of the for loop #1
|
568 |
|
|
test2[,i]<-as.numeric(as.character(test2[,i]))
|
569 |
|
|
}
|
570 |
|
|
|
571 |
|
|
# Plot results
|
572 |
|
|
plot(test$dst,abs(test$res_fus))
|
573 |
|
|
limit<-seq(0,150, by=10)
|
574 |
|
|
tmp<-cut(test$dst,breaks=limit)
|
575 |
|
|
tmp2<-cut(test2$dst,breaks=limit)
|
576 |
|
|
|
577 |
|
|
erd1<-tapply(test$res_fus,tmp, mean)
|
578 |
|
|
erd2<-as.numeric(tapply(abs(test$res_fus),tmp, mean))
|
579 |
|
|
plot(erd2)
|
580 |
|
|
|
581 |
|
|
erd2_CAI<-tapply(abs(test2$res_CAI),tmp2, mean)
|
582 |
|
|
n<-tapply(abs(test$res_fus),tmp, length)
|
583 |
|
|
n2<-tapply(abs(test2$res_CAI),tmp2, length)
|
584 |
|
|
|
585 |
|
|
distance<-seq(5,145,by=10)
|
586 |
|
|
|
587 |
|
|
X11()
|
588 |
|
|
plot(distance,erd2,ylim=c(1,3), type="b", col="red",ylab=" Average MAE",
|
589 |
|
|
xlab="distance to closest training station (km)")
|
590 |
|
|
lines(distance,erd2_CAI,col="grey")
|
591 |
|
|
title("MAE in terms of distance to closest station GAM and FUSION")
|
592 |
|
|
legend("bottomright",legend=c("fus","CAI"), cex=1.2, col=c("red","grey"),
|
593 |
|
|
lty=1, title="MAE")
|
594 |
|
|
savePlot(paste("Comparison_models_er_spat",out_prefix,".png", sep=""), type="png")
|
595 |
|
|
dev.off()
|
596 |
|
|
|
597 |
|
|
means <- erd2_CAI
|
598 |
|
|
means2<- erd2
|
599 |
|
|
stdev <-tapply(abs(test2$res_CAI),tmp2, sd)
|
600 |
|
|
stdev2 <-tapply(abs(test$res_fus),tmp, sd)
|
601 |
|
|
|
602 |
|
|
ciw <- qt(0.975, n) * stdev / sqrt(n)
|
603 |
|
|
ciw2 <- qt(0.975, n) * stdev2 / sqrt(n)
|
604 |
|
|
|
605 |
|
|
X11()
|
606 |
|
|
plotCI(y=means, x=distance, uiw=ciw, col="black", main=" CAI: MAE and distance to clostest training station", barcol="blue", lwd=1)
|
607 |
|
|
lines(distance,erd2_CAI,col="grey")
|
608 |
|
|
points(distance,erd2_CAI,col="grey")
|
609 |
|
|
savePlot(paste("CI_CAI_er_spat_",out_prefix,".png", sep=""), type="png")
|
610 |
|
|
dev.off()
|
611 |
|
|
|
612 |
|
|
X11()
|
613 |
|
|
plotCI(y=means2, x=distance, uiw=ciw2, col="black", main=" FUSION: MAE and distance to clostest training station", barcol="blue", lwd=1)
|
614 |
|
|
lines(distance,erd2,col="black")
|
615 |
|
|
savePlot(paste("CI_fusion_er_spat_",out_prefix,".png", sep=""), type="png")
|
616 |
|
|
dev.off()
|
617 |
|
|
|
618 |
|
|
X11()
|
619 |
|
|
barplot(n,names.arg=as.character(distance))
|
620 |
|
|
savePlot(paste("Barplot_freq_er_spat_",out_prefix,".png", sep=""), type="png")
|
621 |
|
|
dev.off()
|
622 |
|
|
|
623 |
5e7d95a7
|
Benoit Parmentier
|
############################################################
|
624 |
3b090bfa
|
Benoit Parmentier
|
############## PART III #############
|
625 |
|
|
### Average MAE per station and coarse grid box (0.5 deg)
|
626 |
|
|
|
627 |
5e7d95a7
|
Benoit Parmentier
|
#For all stations
|
628 |
|
|
|
629 |
|
|
ghcn$station
|
630 |
|
|
|
631 |
|
|
# For validation and training stations...
|
632 |
|
|
|
633 |
3b090bfa
|
Benoit Parmentier
|
test$abs_res_fus<-abs(test$res_fus)
|
634 |
|
|
test2$abs_res_CAI<-abs(test2$res_CAI)
|
635 |
|
|
|
636 |
|
|
station_melt<-melt(test,
|
637 |
|
|
measure=c("x_OR83M","y_OR83M","res_mod_v","res_mod1","res_mod2","res_mod3","res_mod4","res_mod5","abs_res_fus","abs_res_CAI"),
|
638 |
|
|
id=c("v_id"),
|
639 |
|
|
na.rm=F)
|
640 |
|
|
station_v_er<-cast(station_melt,v_id~variable,mean)
|
641 |
|
|
#station_v_er2<-as.data.frame(station_v_er)
|
642 |
|
|
station_v_er<-as.data.frame(station_v_er)
|
643 |
|
|
oc<-vector("numeric",nrow(station_v_er))
|
644 |
|
|
oc<-oc+1
|
645 |
|
|
station_v_er$oc<-oc
|
646 |
|
|
|
647 |
5e7d95a7
|
Benoit Parmentier
|
unique(ghcn$station)
|
648 |
3b090bfa
|
Benoit Parmentier
|
|
649 |
|
|
coords<- station_v_er[,c('x_OR83M','y_OR83M')]
|
650 |
|
|
coordinates(station_v_er)<-coords
|
651 |
|
|
proj4string(station_v_er)<-CRS #Need to assign coordinates...
|
652 |
|
|
|
653 |
|
|
bubble(station_v_er,"abs_res_fus")
|
654 |
|
|
list_agg_MAE<-vector("list",nel)
|
655 |
|
|
list_agg_RMSE<-vector("list",nel)
|
656 |
|
|
list_density_training<-vector("list",nel)
|
657 |
|
|
list_density_station<-vector("list",nel)
|
658 |
|
|
|
659 |
|
|
for (k in 1:nel){
|
660 |
|
|
data_s<-as.data.frame(list_fus_data_s[[k]])
|
661 |
|
|
data_v<-as.data.frame(list_fus_data_v[[k]])
|
662 |
|
|
list_fus_data[[k]]<-rbind(data_s,data_v)
|
663 |
|
|
data_s2<-as.data.frame(list_cai_data_s[[k]])
|
664 |
|
|
data_v2<-as.data.frame(list_cai_data_v[[k]])
|
665 |
|
|
list_cai_data[[k]]<-rbind(data_s2,data_v2)
|
666 |
|
|
}
|
667 |
|
|
|
668 |
5e7d95a7
|
Benoit Parmentier
|
############ GRID BOX AVERAGING ####################
|
669 |
|
|
####### Create the averaged grid box...##############
|
670 |
3b090bfa
|
Benoit Parmentier
|
|
671 |
|
|
rast_agg<-aggregate(raster_pred,fact=50,fun=mean,na.rm=TRUE) #Changing the raster resolution by aggregation factor
|
672 |
5e7d95a7
|
Benoit Parmentier
|
|
673 |
|
|
ghcn_sub<-as.data.frame(subset(ghcn, select=c("station","x_OR83M","y_OR83M")))
|
674 |
|
|
ghcn_sub_melt<-melt(ghcn_sub,
|
675 |
|
|
measure=c("x_OR83M","y_OR83M"),
|
676 |
|
|
id=c("station"),
|
677 |
|
|
na.rm=F)
|
678 |
|
|
ghcn_stations<-as.data.frame(cast(ghcn_sub_melt,station~variable,mean))
|
679 |
|
|
coords<- ghcn_stations[,c('x_OR83M','y_OR83M')]
|
680 |
|
|
coordinates(ghcn_stations)<-coords
|
681 |
|
|
proj4string(ghcn_stations)<-proj_str #Need to assign coordinates...
|
682 |
|
|
oc_all<-vector("numeric",nrow(ghcn_stations))
|
683 |
|
|
oc_all<-oc_all+1
|
684 |
|
|
|
685 |
|
|
ghcn_stations$oc_all<-oc_all
|
686 |
|
|
rast_oc_all<-rasterize(ghcn_stations,rast_agg,"oc_all",na.rm=TRUE,fun=sum)
|
687 |
|
|
ac_agg50$oc_all<-values(rast_oc_all)
|
688 |
|
|
|
689 |
|
|
plot(rast_oc_all, main="Number of stations in coarsened 50km grid")
|
690 |
|
|
plot(reg_outline, add=TRUE)
|
691 |
|
|
fdens_all50<-as.data.frame(freq(rast_oc_all))
|
692 |
|
|
tot50<-sum(fdens_all50$count[1:(nrow(fdens_all50)-1)])
|
693 |
|
|
percent<-(fdens_all50$count/tot50)*100
|
694 |
|
|
percent[length(percent)]<-NA
|
695 |
|
|
fdens_all50$percent<-percent
|
696 |
3b090bfa
|
Benoit Parmentier
|
#list_agg_MAE<-vector("list",nel)
|
697 |
|
|
#list_agg_RMSE<-vector("list",nel)
|
698 |
|
|
#list_density_training<-vector("list",nel)
|
699 |
|
|
#list_density_station<-vector("list",nel)
|
700 |
|
|
#start loop here for grid box aggregation
|
701 |
|
|
list_density_fus<-vector("list",nel)
|
702 |
|
|
list_density_cai<-vector("list",nel)
|
703 |
|
|
|
704 |
|
|
nel<-365
|
705 |
|
|
for (k in 1:nel){
|
706 |
|
|
|
707 |
|
|
data_s<-list_fus_data_s[[k]] #Extracting the relevant spdf from the list: this is 1050?
|
708 |
|
|
data_s2<-list_cai_data_s[[k]]
|
709 |
|
|
data_v<-list_fus_data_v[[k]]
|
710 |
|
|
data_v2<-list_cai_data_v[[k]]
|
711 |
|
|
|
712 |
|
|
data_v$abs_res_fus<-abs(data_v$res_mod9)
|
713 |
|
|
data_v2$abs_res_CAI<-abs(data_v2$res_mod9)
|
714 |
|
|
data_s$abs_res_fus<-abs(data_s$res_mod9)
|
715 |
|
|
data_s2$abs_res_CAI<-abs(data_s2$res_mod9)
|
716 |
|
|
data_s$oc<-rep(1,nrow(data_s))
|
717 |
|
|
data_s2$oc<-rep(1,nrow(data_s2))
|
718 |
|
|
|
719 |
|
|
#Computing MAE per grid box
|
720 |
|
|
rast_MAE_fus<-rasterize(data_v,rast_agg,"abs_res_fus",na.rm=TRUE,fun=mean)
|
721 |
|
|
rast_MAE_cai<-rasterize(data_v2,rast_agg,"abs_res_CAI",na.rm=TRUE,fun=mean)
|
722 |
|
|
#Computing density of station
|
723 |
|
|
rast_oc<-rasterize(data_s,rast_agg,"oc",na.rm=TRUE,fun=sum)
|
724 |
|
|
rast_oc2<-rasterize(data_s2,rast_agg,"oc",na.rm=TRUE,fun=sum)
|
725 |
|
|
|
726 |
|
|
#Creating plots adding to the list...to get a data frame...
|
727 |
|
|
ac_agg50<-as.data.frame(values(rast_oc))
|
728 |
|
|
ac_agg50$MAE_fus<-as.numeric(values(rast_MAE_fus))
|
729 |
|
|
ac_agg50$MAE_cai<-as.numeric(values(rast_MAE_cai))
|
730 |
|
|
names(ac_agg50)<-c("oc","MAE_fus","MAE_cai")
|
731 |
|
|
|
732 |
5e7d95a7
|
Benoit Parmentier
|
#ghcn_sub<-as.data.frame(subset(ghcn, select=c("station","x_OR83M","y_OR83M")))
|
733 |
|
|
#ghcn_sub_melt<-melt(ghcn_sub,
|
734 |
|
|
# measure=c("x_OR83M","y_OR83M"),
|
735 |
|
|
# id=c("station"),
|
736 |
|
|
# na.rm=F)
|
737 |
|
|
#ghcn_stations<-as.data.frame(cast(ghcn_sub_melt,station~variable,mean))
|
738 |
|
|
#coords<- ghcn_stations[,c('x_OR83M','y_OR83M')]
|
739 |
|
|
#coordinates(ghcn_stations)<-coords
|
740 |
|
|
#proj4string(ghcn_stations)<-proj_str #Need to assign coordinates...
|
741 |
|
|
#oc_all<-vector("numeric",nrow(ghcn_stations))
|
742 |
|
|
#oc_all<-oc_all+1
|
743 |
|
|
|
744 |
|
|
#ghcn_stations$oc_all<-oc_all
|
745 |
|
|
#rast_oc_all<-rasterize(ghcn_stations,rast_agg,"oc_all",na.rm=TRUE,fun=sum)
|
746 |
|
|
#ac_agg50$oc_all<-values(rast_oc_all)
|
747 |
3b090bfa
|
Benoit Parmentier
|
|
748 |
|
|
td1<-aggregate(MAE_fus~oc,data=ac_agg50,mean)
|
749 |
|
|
td2<-aggregate(MAE_cai~oc,data=ac_agg50,mean)
|
750 |
|
|
td<-merge(td1,td2,by="oc")
|
751 |
|
|
|
752 |
|
|
td1_all<-aggregate(MAE_fus~oc_all,data=ac_agg50,mean)
|
753 |
|
|
td2_all<-aggregate(MAE_cai~oc_all,data=ac_agg50,mean)
|
754 |
|
|
td_all<-merge(td1_all,td2_all,by="oc_all")
|
755 |
|
|
|
756 |
|
|
plot(MAE_fus~oc,data=td,type="b")
|
757 |
|
|
lines(td$oc,td$MAE_cai, type="b", lwd=1.5,co="red")
|
758 |
|
|
plot(MAE_fus~oc_all,data=td_all,type="b")
|
759 |
|
|
lines(td_all$oc_all,td_all$MAE_cai, type="b", lwd=1.5,co="red")
|
760 |
|
|
|
761 |
|
|
filename<-sub(".shp","",infile6) #Removing the extension from file.
|
762 |
|
|
reg_outline<-readOGR(".", filename) #reading shapefile
|
763 |
|
|
plot(rast_MAE_fus, main="Fusion MAE in coarsened 50km grid")
|
764 |
|
|
plot(reg_outline, add=TRUE)
|
765 |
|
|
|
766 |
|
|
plot(rast_MAE_cai, main="CAI MAE in coarsened 50km grid")
|
767 |
|
|
plot(reg_outline, add=TRUE)
|
768 |
|
|
|
769 |
|
|
plot(rast_oc, main="Number of val stations in coarsened 50km grid")
|
770 |
|
|
plot(reg_outline, add=TRUE)
|
771 |
5e7d95a7
|
Benoit Parmentier
|
plot(rast_oc_t, main="Number of training stations in coarsened 50km grid")
|
772 |
3b090bfa
|
Benoit Parmentier
|
plot(reg_outline, add=TRUE)
|
773 |
|
|
|
774 |
|
|
#MAKE IT AN OBJECT for future function return...
|
775 |
|
|
|
776 |
|
|
#list(rast_MAE,rast_RMSE,rast_oc,rast_all)
|
777 |
|
|
list_density_fus[[k]]<-list(rast_MAE_fus,rast_oc,rast_oc_all,ac_agg50)
|
778 |
|
|
list_density_cai[[k]]<-list(rast_MAE_cai,rast_oc,rast_oc_all,ac_agg50)
|
779 |
|
|
}
|
780 |
|
|
|
781 |
|
|
#meean over stack oc
|
782 |
|
|
#mean over stack MAE
|
783 |
|
|
do.call(rbind,list_density_)
|
784 |
|
|
list_var_stat<-vector("list", 365)
|
785 |
|
|
#list_var_stat<-vector("list", 2)
|
786 |
|
|
#k=2
|
787 |
|
|
|
788 |
|
|
for (k in 1:length(l_f)){
|
789 |
|
|
|
790 |
|
|
raster_pred<-raster(l_f[[k]])
|
791 |
|
|
layerNames(raster_pred)<-"fus"
|
792 |
|
|
projection(raster_pred)<-proj_str
|
793 |
|
|
|
794 |
|
|
raster_pred2<-raster(l_f2[[k]])
|
795 |
|
|
layerNames(raster_pred2)<-"fus"
|
796 |
|
|
projection(raster_pred2)<-proj_str
|
797 |
|
|
|
798 |
|
|
tmp_rast<-mask(raster_pred2,raster_pred)
|
799 |
|
|
raster_pred2<-tmp_rast
|
800 |
|
|
|
801 |
|
|
t1<-cellStats(raster_pred,na.rm=TRUE,stat=sd) #Calculating the standard deviation for the
|
802 |
|
|
t2<-cellStats(raster_pred2,na.rm=TRUE,stat=sd)
|
803 |
|
|
|
804 |
|
|
m1<-Moran(raster_pred,w=3) #Calculating Moran's I with window of 3 an default Queen's case
|
805 |
|
|
m2<-Moran(tmp_rast,w=3) #Calculating Moran's I with window of 3 an default Queen's case for prediction 2 (CAI)
|
806 |
|
|
stat<-as.data.frame(t(c(m1,m2,t1,t2)))
|
807 |
|
|
names(stat)<-c("moran_fus","moran_CAI","sd_fus","sd_CAI")
|
808 |
|
|
list_var_stat[[k]]<-stat
|
809 |
|
|
}
|
810 |
|
|
|
811 |
|
|
var_stat<-do.call(rbind,list_var_stat)
|
812 |
|
|
|
813 |
|
|
|
814 |
|
|
pos<-match("ELEV_SRTM",layerNames(s_raster)) #Find column with name "value"
|
815 |
|
|
elev<-raster(s_raster,layer=pos) #Select layer from stack
|
816 |
|
|
elev<-mask(elev,raster_pred)
|
817 |
|
|
te<-cellStats(elev,na.rm=TRUE,stat=sd)
|
818 |
|
|
|
819 |
|
|
pos<-match("mm_12",layerNames(s_raster)) #Find column with name "value"
|
820 |
|
|
m_12<-raster(s_raster,layer=pos) #Select layer from stack
|
821 |
|
|
m_LST<-Moran(m_12,w=3)
|
822 |
|
|
m_e<-Moran(elev,w=3)
|
823 |
|
|
m_12<-m_12-273.15
|
824 |
|
|
plot(MAE_fus~oc,data=td,type="b")
|
825 |
|
|
lines(td$oc,td$MAE_CAI, type="b", lwd=1.5,co="red")
|
826 |
|
|
|
827 |
|
|
data_dist<-as.data.frame(cbind(distance,erd2,erd2_mod1,erd2_mod2,erd2_mod3,erd2_mod4,erd2_mod5,erd2_CAI,n))
|
828 |
|
|
rownames(data_dist)<-NULL
|
829 |
|
|
|
830 |
|
|
############# PART IV COMPARISON OF SPATIAL PATTERN BY EXAMING MAPS OF PREDICTION
|
831 |
|
|
#PLOTING CAI AND FUSION TO COMPARE
|
832 |
|
|
|
833 |
|
|
infile2<-"list_10_dates_04212012.txt" #List of 10 dates for the regression
|
834 |
|
|
dates2<-read.table(paste(path,"/",infile2,sep=""), sep="") #Column 1 contains the names of raster files
|
835 |
|
|
date_list2<-as.list(as.character(dates2[,1]))
|
836 |
5e7d95a7
|
Benoit Parmentier
|
names_statistics<-c("mean","sd","min","max")
|
837 |
|
|
stat_val_m<-matrix(NA,nrow=length(date_list2),ncol=length(names_statistics))
|
838 |
|
|
colnames(stat_val_m)<-names_statistics
|
839 |
|
|
rownames(stat_val_m)<-date_list2
|
840 |
|
|
stat_val_m<-as.data.frame(stat_val_m)
|
841 |
3b090bfa
|
Benoit Parmentier
|
|
842 |
|
|
for (k in 1:length(date_list2)){
|
843 |
|
|
|
844 |
|
|
date_proc2<-date_list2[[k]]
|
845 |
|
|
#date_proc<-date_list[[k]]
|
846 |
|
|
index<-match(as.character(date_proc2),unlist(date_list)) #find the correct date... in the 365 stack
|
847 |
|
|
#raster_pred<-raster(rp_raster,index)
|
848 |
5e7d95a7
|
Benoit Parmentier
|
raster_pred1<-raster(l_f[[index]]) # Fusion image
|
849 |
3b090bfa
|
Benoit Parmentier
|
projection(raster_pred1)<-proj_str
|
850 |
|
|
raster_pred1<-mask(raster_pred1,mask_land_NA)
|
851 |
|
|
|
852 |
5e7d95a7
|
Benoit Parmentier
|
raster_pred2<-raster(l_f2[[index]]) # CAI image
|
853 |
3b090bfa
|
Benoit Parmentier
|
projection(raster_pred2)<-proj_str
|
854 |
|
|
raster_pred2<-mask(raster_pred2,mask_land_NA)
|
855 |
|
|
|
856 |
|
|
predictions <- stack(raster_pred1,raster_pred2)
|
857 |
|
|
layerNames(predictions)<-c(paste('fusion',date_list2[[k]],sep=" "),paste('CAI',date_list2[[k]],sep=" "))
|
858 |
|
|
# use overall min and max values to generate an nice, consistent set
|
859 |
|
|
# of breaks for both colors (50 values) and legend labels (5 values)
|
860 |
|
|
s.range <- c(min(minValue(predictions)), max(maxValue(predictions)))
|
861 |
|
|
col.breaks <- pretty(s.range, n=50)
|
862 |
|
|
lab.breaks <- pretty(s.range, n=5)
|
863 |
|
|
temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
|
864 |
|
|
|
865 |
|
|
# plot using these (common) breaks; note use of _reverse_ heat.colors,
|
866 |
|
|
# making it so that larger numbers are redder
|
867 |
|
|
X11(6,12)
|
868 |
|
|
#plot(predictions, breaks=col.breaks, col=rev(heat.colors(length(col.breaks)-1)),
|
869 |
|
|
# axis=list(at=lab.breaks, labels=lab.breaks))
|
870 |
|
|
plot(predictions, breaks=col.breaks, col=temp.colors(length(col.breaks)-1),
|
871 |
|
|
axis=list(at=lab.breaks, labels=lab.breaks))
|
872 |
|
|
#plot(reg_outline, add=TRUE)
|
873 |
|
|
savePlot(paste("comparison_raster1_CAI_fusion_tmax_prediction_",date_list2[[k]],out_prefix,".png", sep=""), type="png")
|
874 |
5e7d95a7
|
Benoit Parmentier
|
|
875 |
|
|
|
876 |
|
|
stat_val_m$mean[i]<-cellStats(raster_pred1,na.rm=TRUE,stat=mean) #Calculating the standard deviation for the
|
877 |
|
|
t1<-cellStats(raster_pred1,na.rm=TRUE,stat=mean) #Calculating the standard deviation for the
|
878 |
|
|
t2<-cellStats(raster_pred2,na.rm=TRUE,stat=mean) #Calculating the standard deviation for the
|
879 |
|
|
t1<-cellStats(raster_pred1,na.rm=TRUE,stat=sd) #Calculating the standard deviation for the
|
880 |
|
|
t2<-cellStats(raster_pred2,na.rm=TRUE,stat=sd) #Calculating the standard deviation for the
|
881 |
|
|
t1<-cellStats(raster_pred1,na.rm=TRUE,stat=min) #Calculating the standard deviation for the
|
882 |
|
|
t2<-cellStats(raster_pred2,na.rm=TRUE,stat=min) #Calculating the standard deviation for the
|
883 |
|
|
t1<-cellStats(raster_pred1,na.rm=TRUE,stat=max) #Calculating the standard deviation for the
|
884 |
|
|
t2<-cellStats(raster_pred2,na.rm=TRUE,stat=max) #Calculating the standard deviation for the
|
885 |
|
|
|
886 |
|
|
hist(predictions,freq=FALSE,maxpixels=ncells(predictions),xlabel="tmax (C)")
|
887 |
3b090bfa
|
Benoit Parmentier
|
savePlot(paste("comparison_histo_CAI_fusion_tmax_prediction_",date_list2[[k]],out_prefix,".png", sep=""), type="png")
|
888 |
|
|
#plot(predictions)
|
889 |
|
|
dev.off()
|
890 |
|
|
|
891 |
5e7d95a7
|
Benoit Parmentier
|
X11(6,12)
|
892 |
3b090bfa
|
Benoit Parmentier
|
diff<-raster_pred2-raster_pred1
|
893 |
|
|
s.range <- c(min(minValue(diff)), max(maxValue(diff)))
|
894 |
|
|
col.breaks <- pretty(s.range, n=50)
|
895 |
|
|
lab.breaks <- pretty(s.range, n=5)
|
896 |
|
|
temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
|
897 |
|
|
plot(diff, breaks=col.breaks, col=temp.colors(length(col.breaks)-1),
|
898 |
|
|
axis=list(at=lab.breaks, labels=lab.breaks))
|
899 |
5e7d95a7
|
Benoit Parmentier
|
title(paste("Difference between CAI and fusion for ",date_list2[[k]],sep=""))
|
900 |
|
|
savePlot(paste("comparison_diff_CAI_fusion_tmax_prediction_",date_list2[[k]],out_prefix,".png", sep=""), type="png")
|
901 |
|
|
dev.off()
|
902 |
3b090bfa
|
Benoit Parmentier
|
|
903 |
5e7d95a7
|
Benoit Parmentier
|
}
|
904 |
3b090bfa
|
Benoit Parmentier
|
|
905 |
|
|
write.table(data_dist,file=paste("data_dist_",out_prefix,".txt",sep=""),sep=",")
|
906 |
|
|
write.table(test,file=paste("ac_spat_dist",out_prefix,".txt",sep=""),sep=",")
|
907 |
|
|
write.table(var_stat,file=paste("moran_var_stat_",out_prefix,".txt",sep=""),sep=",")
|
908 |
|
|
write.table(td,file=paste("MAE_density_station_",out_prefix,".txt",sep=""),sep=",")
|
909 |
|
|
write.table(td_all,file=paste("MAE_density_station_all_",out_prefix,".txt",sep=""),sep=",")
|
910 |
|
|
|
911 |
5e7d95a7
|
Benoit Parmentier
|
symbols(c(2e5, 4e5), c(2e5, 3e5), circles=rep(2e4, 2), inches=FALSE, add=TRUE)
|
912 |
|
|
text(c(2e5, 4e5), c(2e5, 3e5), labels=1:2,cex=0.5)
|
913 |
|
|
points(coordinates(pts), type="c")
|
914 |
|
|
text(coordinates(pts), labels=9:11, cex=0.8)
|
915 |
|
|
points(coordinates(pts), type="b", pch=as.character(1:length(pts))
|
916 |
|
|
points(coordinates(pts), type="b", pch=as.character(9:11)
|
917 |
|
|
|
918 |
|
|
########### END OF THE SCRIPT #############
|