Revision c5732f39
Added by Benoit Parmentier almost 12 years ago
climate/research/oregon/interpolation/methods_comparison_assessment_part4.R | ||
---|---|---|
1 |
##################################### METHODS COMPARISON part 3 ##########################################
|
|
1 |
##################################### METHODS COMPARISON part 4 ##########################################
|
|
2 | 2 |
#################################### Spatial Analysis ############################################ |
3 | 3 |
#This script utilizes the R ojbects created during the interpolation phase. # |
4 | 4 |
#At this stage the script produces figures of various accuracy metrics and compare methods: # |
... | ... | |
7 | 7 |
#- spatial density of station network and accuracy metric # |
8 | 8 |
#- visualization of maps of prediction and difference for comparison # |
9 | 9 |
#AUTHOR: Benoit Parmentier # |
10 |
#DATE: 10/23/2012 #
|
|
10 |
#DATE: 11/23/2012 #
|
|
11 | 11 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491 -- # |
12 | 12 |
################################################################################################### |
13 | 13 |
|
... | ... | |
60 | 60 |
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"; |
61 | 61 |
#User defined output prefix |
62 | 62 |
|
63 |
#CRS<-proj4string(ghcn) #Storing projection information (ellipsoid, datum,etc.) |
|
64 |
lines<-read.table(paste(path,"/",inlistf,sep=""), sep="") #Column 1 contains the names of raster files |
|
65 |
inlistvar<-lines[,1] |
|
66 |
inlistvar<-paste(path,"/",as.character(inlistvar),sep="") |
|
67 |
covar_names<-as.character(lines[,2]) #Column two contains short names for covaraites |
|
68 |
|
|
69 |
s_raster<- stack(inlistvar) #Creating a stack of raster images from the list of variables. |
|
70 |
layerNames(s_raster)<-covar_names #Assigning names to the raster layers |
|
71 |
projection(s_raster)<-proj_str |
|
72 |
|
|
73 |
#Create mask using land cover data |
|
74 |
pos<-match("LC10",layerNames(s_raster)) #Find the layer which contains water bodies |
|
75 |
LC10<-subset(s_raster,pos) |
|
76 |
LC10[is.na(LC10)]<-0 #Since NA values are 0, we assign all zero to NA |
|
77 |
mask_land<-LC10<100 #All values below 100% water are assigned the value 1, value 0 is "water" |
|
78 |
mask_land_NA<-mask_land |
|
79 |
mask_land_NA[mask_land_NA==0]<-NA #Water bodies are assigned value 1 |
|
80 |
|
|
81 |
data_name<-"mask_land_OR" |
|
82 |
raster_name<-paste(data_name,".rst", sep="") |
|
83 |
writeRaster(mask_land, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
84 |
#writeRaster(r2, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
85 |
|
|
86 |
pos<-match("ELEV_SRTM",layerNames(s_raster)) #Find column with name "ELEV_SRTM" |
|
87 |
ELEV_SRTM<-raster(s_raster,layer=pos) #Select layer from stack on 10/30 |
|
88 |
s_raster<-dropLayer(s_raster,pos) |
|
89 |
ELEV_SRTM[ELEV_SRTM <0]<-NA |
|
90 |
mask_ELEV_SRTM<-ELEV_SRTM>0 |
|
91 |
|
|
92 |
#Change this a in loop... |
|
93 |
pos<-match("LC1",layerNames(s_raster)) #Find column with name "value" |
|
94 |
LC1<-raster(s_raster,layer=pos) #Select layer from stack |
|
95 |
s_raster<-dropLayer(s_raster,pos) |
|
96 |
LC1[is.na(LC1)]<-0 |
|
97 |
pos<-match("LC2",layerNames(s_raster)) #Find column with name "value" |
|
98 |
LC2<-raster(s_raster,layer=pos) #Select layer from stack |
|
99 |
s_raster<-dropLayer(s_raster,pos) |
|
100 |
LC2[is.na(LC2)]<-0 |
|
101 |
pos<-match("LC3",layerNames(s_raster)) #Find column with name "value" |
|
102 |
LC3<-raster(s_raster,layer=pos) #Select layer from stack |
|
103 |
s_raster<-dropLayer(s_raster,pos) |
|
104 |
LC3[is.na(LC3)]<-0 |
|
105 |
pos<-match("LC4",layerNames(s_raster)) #Find column with name "value" |
|
106 |
LC4<-raster(s_raster,layer=pos) #Select layer from stack |
|
107 |
s_raster<-dropLayer(s_raster,pos) |
|
108 |
LC4[is.na(LC4)]<-0 |
|
109 |
pos<-match("LC6",layerNames(s_raster)) #Find column with name "value" |
|
110 |
LC6<-raster(s_raster,layer=pos) #Select layer from stack |
|
111 |
s_raster<-dropLayer(s_raster,pos) |
|
112 |
LC6[is.na(LC6)]<-0 |
|
113 |
pos<-match("LC7",layerNames(s_raster)) #Find column with name "value" |
|
114 |
LC7<-raster(s_raster,layer=pos) #Select layer from stack |
|
115 |
s_raster<-dropLayer(s_raster,pos) |
|
116 |
LC7[is.na(LC7)]<-0 |
|
117 |
pos<-match("LC9",layerNames(s_raster)) #Find column with name "LC9", this is wetland... |
|
118 |
LC9<-raster(s_raster,layer=pos) #Select layer from stack |
|
119 |
s_raster<-dropLayer(s_raster,pos) |
|
120 |
LC9[is.na(LC9)]<-0 |
|
121 |
|
|
122 |
LC_s<-stack(LC1,LC2,LC3,LC4,LC6,LC7) |
|
123 |
layerNames(LC_s)<-c("LC1_forest","LC2_shrub","LC3_grass","LC4_crop","LC6_urban","LC7_barren") |
|
124 |
LC_s <-mask(LC_s,mask_ELEV_SRTM) |
|
125 |
plot(LC_s) |
|
126 |
|
|
127 |
s_raster<-addLayer(s_raster, LC_s) |
|
128 |
|
|
129 |
#mention this is the last... files |
|
130 |
|
|
131 |
#Read region outline... |
|
132 |
filename<-sub(".shp","",infile6) #Removing the extension from file. |
|
133 |
reg_outline<-readOGR(".", filename) #reading shapefile |
|
134 |
|
|
63 | 135 |
############ PART 4: RESIDUALS ANALYSIS: ranking, plots, focus regions ################## |
64 | 136 |
############## EXAMINING STATION RESIDUALS ########### |
65 | 137 |
########### CONSTANT OVER 365 AND SAMPLING OVER 365 |
... | ... | |
93 | 165 |
rast_cai_pred<-raster(rast_cai2c,1) |
94 | 166 |
|
95 | 167 |
#list files that contain model objects and ratingin-testing information for CAI and Fusion |
96 |
fus1_c<-load_obj("results2_fusion_Assessment_measure_all_365d_GAM_fusion_const_10172012.RData") |
|
97 | 168 |
|
98 | 169 |
gam_fus<-load_obj(file.path(path_data_fus, |
99 | 170 |
"results_mod_obj__365d_GAM_fusion_const_all_lstd_11022012.RData")) |
... | ... | |
102 | 173 |
sampling_date_list<-gam_fus$sampling_obj$sampling_dat$date |
103 | 174 |
|
104 | 175 |
#Create a residual table... |
105 |
res_mod9_list<-vector("list",365) |
|
106 |
res_mod2_list<-vector("list",365) |
|
107 |
|
|
108 |
tab_nv<-matrix(NA,365,1) |
|
109 |
for(k in 1:365){ |
|
110 |
tab_nv[k]<-length(fus1_c[[k]]$data_v$res_mod9) |
|
111 |
} |
|
112 |
#note that there might be some variation in the number!!! |
|
176 |
res_mod9f_list<-vector("list",365) |
|
177 |
res_mod9c_list<-vector("list",365) |
|
178 |
data_vf_list<-vector("list",365) |
|
179 |
data_vc_list<-vector("list",365) |
|
113 | 180 |
|
114 | 181 |
for(k in 1:365){ |
115 |
res_mod9<-gam_fus$gam_fus_mod[[k]]$data_v$res_mod7
|
|
116 |
res_mod2<-gam_fus$gam_fus_mod[[k]]$data_v$res_mod2
|
|
117 |
res_mod9_list[[k]]<-res_mod9
|
|
118 |
res_mod2_list[[k]]<-res_mod2
|
|
182 |
data_vf<-as.data.frame(c)
|
|
183 |
data_vc<-as.data.frame(gam_cai$gam_CAI_mod[[k]]$data_v)
|
|
184 |
data_vf<-data_vf[,c("id","date","res_mod7","x_OR83M", "y_OR83M")]
|
|
185 |
data_vc<-data_vc[,c("id","date","res_mod9","x_OR83M", "y_OR83M")]
|
|
119 | 186 |
#subset data frame? or rbind them...and reshape?? think about it |
187 |
data_vf_list[[k]]<-data_vf |
|
188 |
data_vc_list[[k]]<-data_vc |
|
120 | 189 |
} |
121 |
tab_resmod9<-do.call(rbind,res_mod9_list) |
|
122 |
|
|
123 |
data_v_list<-vector("list",365) |
|
190 |
tab_resmod9f<-do.call(rbind,data_vf_list) |
|
191 |
tab_resmod9c<-do.call(rbind,data_vc_list) |
|
124 | 192 |
|
125 |
for(k in 1:365){ |
|
126 |
data_v<-as.data.frame(gam_fus$gam_fus_mod[[k]]$data_v) |
|
127 |
data_v<-data_v[,c("id","date","res_mod7","x_OR83M", "y_OR83M")] |
|
128 |
#subset data frame? or rbind them...and reshape?? think about it |
|
129 |
data_v_list[[k]]<-data_v |
|
130 |
} |
|
131 |
tab_resmod9<-do.call(rbind,data_v_list) |
|
132 |
tab_locs<-melt(tab_resmod9, |
|
193 |
#Get the unique location information... |
|
194 |
tab_locsc<-melt(tab_resmod9c, |
|
133 | 195 |
measure=c("x_OR83M","y_OR83M"), |
134 | 196 |
id=c("id"), |
135 | 197 |
na.rm=F) |
136 |
tab_locs_cast<-cast(tab_locs,id~variable,mean) |
|
137 |
tab_locs<-as.data.frame(tab_locs_cast) |
|
138 |
coords<- tab_locs[,c('x_OR83M','y_OR83M')] |
|
139 |
coordinates(tab_locs)<-coords |
|
140 |
proj4string(tab_locs)<-proj_str #Need to assign coordinates... |
|
141 |
|
|
142 |
tab_melt<-melt(tab_resmod9[c("id","date","res_mod7")], |
|
198 |
tab_locsc_cast<-cast(tab_locsc,id~variable,mean) |
|
199 |
tab_locsc<-as.data.frame(tab_locsc_cast) |
|
200 |
coords<- tab_locsc[,c('x_OR83M','y_OR83M')] |
|
201 |
coordinates(tab_locsc)<-coords |
|
202 |
proj4string(tab_locsc)<-proj_str #Need to assign coordinates... |
|
203 |
|
|
204 |
tab_locsf<-melt(tab_resmod9f, |
|
205 |
measure=c("x_OR83M","y_OR83M"), |
|
206 |
id=c("id"), |
|
207 |
na.rm=F) |
|
208 |
tab_locsf_cast<-cast(tab_locsf,id~variable,mean) |
|
209 |
tab_locsf<-as.data.frame(tab_locsf_cast) |
|
210 |
coords<- tab_locsf[,c('x_OR83M','y_OR83M')] |
|
211 |
coordinates(tab_locsf)<-coords |
|
212 |
proj4string(tab_locsf)<-proj_str #Need to assign coordinates.. |
|
213 |
|
|
214 |
all.equal(tab_locsc,tab_locsf) #Checking that stations used over the year are equal!!! |
|
215 |
#since all equal we can use one object, it will be sufficient... |
|
216 |
tab_locs<-tab_locsf |
|
217 |
|
|
218 |
####### Now join and summarize objects together... |
|
219 |
|
|
220 |
tabf_melt<-melt(tab_resmod9f[c("id","date","res_mod7")], |
|
143 | 221 |
measure=c("res_mod7"), |
144 | 222 |
id=c("id","date"), |
145 | 223 |
na.rm=F) |
224 |
tabc_melt<-melt(tab_resmod9c[c("id","date","res_mod9")], |
|
225 |
measure=c("res_mod9"), |
|
226 |
id=c("id","date"), |
|
227 |
na.rm=F) |
|
228 |
tabf_cast<-cast(tabf_melt,id~date~variable,mean) |
|
229 |
tabf_resmod9_locs<-as.data.frame(tabf_cast[,,1]) |
|
230 |
tabc_cast<-cast(tabc_melt,id~date~variable,mean) |
|
231 |
tabc_resmod9_locs<-as.data.frame(tabc_cast[,,1]) |
|
146 | 232 |
|
147 |
tab_cast<-cast(tab_melt,id~date~variable,mean) |
|
148 |
tab_resmod9_locs<-as.data.frame(tab_cast[,,1]) |
|
149 |
sd_v<-sapply(tab_resmod9_locs,sd,na.rm=T) |
|
150 |
mean_v<-sapply(tab_resmod9_locs,mean,na.rm=T) |
|
151 |
tab_res_rec<-tab_resmod9_locs |
|
152 |
|
|
153 |
for (k in 1:365){ |
|
154 |
mean<-mean_v[k] |
|
155 |
sd<-sd_v[k] |
|
156 |
y<-tab_resmod9_locs[,k] |
|
157 |
breaks<-c("-inf",mean-2*sd,mean+2*sd,"+inf") |
|
158 |
rec<-cut(y,breaks,labels=c(-1,0,1)) |
|
159 |
rec<-as.numeric(as.character(rec)) |
|
160 |
tab_res_rec[,k]<-rec |
|
161 |
} |
|
162 |
|
|
163 |
|
|
164 |
tab_res_rec<-tab_res_rec[,1:365] #mae_fun<-function (x){mean(abs(x),na.rm=T)} |
|
165 |
for (i in 1:nrow(tab_resmod9_locs)){ |
|
166 |
rec<-as.numeric(tab_res_rec[i,]) |
|
167 |
tmp<-table(rec) |
|
168 |
tab_res_rec$c1[i]<-as.numeric(tmp[1]) |
|
169 |
tab_res_rec$c2[i]<-as.numeric(tmp[2]) |
|
170 |
tab_res_rec$c3[i]<-as.numeric(tmp[3]) |
|
171 |
} |
|
172 | 233 |
|
173 |
tab_res_rec$id<-as.character(rownames(tab_res_rec)) |
|
174 |
#mae_fun<-function (x){mean(abs(x),na.rm=T)} |
|
175 |
for (i in 1:nrow(tab_resmod9_locs)){ |
|
176 |
x<-tab_resmod9_locs[i,] |
|
177 |
tab_locs$mae[i]<-mean(abs(x),na.rm=T) |
|
178 |
#tab_locs$mae2[i]<-mae_fun(x) |
|
179 |
rec<-as.numeric(tab_res_rec[i,]) |
|
180 |
tmp<-table(rec) |
|
181 |
tab_res_rec$c1<-as.numeric(tmp[1]) |
|
182 |
tab_res_rec$c2<-as.numeric(tmp[2]) |
|
183 |
tab_res_rec$c3<-as.numeric(tmp[3]) |
|
234 |
#Now reclass the data into outliers: Class with -1 (negative exteme),class 0 (not extreme), class 1 (positive extre) |
|
235 |
calculate_extremes_stat <-function (tab_resmod9_locs,tab_locs){ |
|
236 |
|
|
237 |
#tab_resmod9_locs<-tabf_resmod9_locs #comment out when running function |
|
238 |
#tab_locs<-tab_locsf |
|
239 |
#tab_res_rec<-tabf_resmod9_locs |
|
240 |
|
|
241 |
#start here |
|
242 |
sd_v<-sapply(tab_resmod9_locs,sd,na.rm=T) #This the sandard deviation of residuals for every day for the fusion prediction |
|
243 |
mean_v<-sapply(tab_resmod9_locs,mean,na.rm=T) #This the mean of residuals for every day for the fusion prediction |
|
244 |
|
|
245 |
for (k in 1:365){ |
|
246 |
mean_val<-mean_v[k] |
|
247 |
sd_val<-sd_v[k] |
|
248 |
y<-tab_resmod9_locs[,k] |
|
249 |
breaks<-c("-inf",mean_val-2*sd_val,mean_val+2*sd_val,"+inf") |
|
250 |
rec<-cut(y,breaks,labels=c(-1,0,1)) |
|
251 |
rec<-as.numeric(as.character(rec)) |
|
252 |
tab_res_rec[,k]<-rec |
|
253 |
} |
|
184 | 254 |
|
255 |
tab_res_rec<-tab_res_rec[,1:365] #mae_fun<-function (x){mean(abs(x),na.rm=T)} |
|
256 |
|
|
257 |
tab_res_rec$id<-as.character(rownames(tab_res_rec)) |
|
258 |
#mae_fun<-function (x){mean(abs(x),na.rm=T)} |
|
259 |
for (i in 1:nrow(tab_resmod9_locs)){ |
|
260 |
x<-tab_resmod9_locs[i,] |
|
261 |
tab_locs$mae[i]<-mean(abs(x),na.rm=T) |
|
262 |
#tab_locs$mae2[i]<-mae_fun(x) |
|
263 |
rec<-as.numeric(tab_res_rec[i,]) |
|
264 |
tmp<-table(rec) |
|
265 |
tab_res_rec$c1[i]<-as.numeric(tmp[1]) |
|
266 |
tab_res_rec$c2[i]<-as.numeric(tmp[2]) |
|
267 |
tab_res_rec$c3[i]<-as.numeric(tmp[3]) |
|
268 |
} |
|
269 |
|
|
270 |
projstr_info<-proj4string(tab_locs) |
|
271 |
tab_resmod9_locs$id<-rownames(tab_resmod9_locs) |
|
272 |
tab_resmod9_locs[is.na(tab_resmod9_locs)]<-NA #replace NaN by NA |
|
273 |
tab_locs$id<-as.character(tab_locs$id) |
|
274 |
data_v_res<-merge(tab_locs,tab_resmod9_locs,by="id") |
|
275 |
data_v_res<-merge(tab_res_rec[,c("id","c1","c2","c3")],data_v_res,by="id") |
|
276 |
|
|
277 |
coords<- data_v_res[,c('x_OR83M','y_OR83M')] |
|
278 |
coordinates(data_v_res)<-coords |
|
279 |
proj4string(data_v_res)<-projstr_info #Need to assign coordinates... |
|
280 |
|
|
281 |
#tmp<-subset(data_v_res,select=date_selected) |
|
282 |
data_v_res$idx<-1:length(data_v_res) |
|
283 |
#Plotting results... |
|
284 |
return (data_v_res) |
|
185 | 285 |
} |
186 |
tab_resmod9_locs$id<-rownames(tab_resmod9_locs) |
|
187 |
tab_resmod9_locs[is.na(tab_resmod9_locs)]<-NA #replace NaN by NA |
|
188 |
tab_locs$id<-as.character(tab_locs$id) |
|
189 |
data_v_res<-merge(tab_locs,tab_resmod9_locs,by="id") |
|
190 |
data_v_res<-merge(tab_res_rec[,c("id","c1","c2","c3")],data_v_res,by="id") |
|
191 |
coords<- data_v_res[,c('x_OR83M','y_OR83M')] |
|
192 |
coordinates(data_v_res)<-coords |
|
193 |
proj4string(data_v_res)<-proj_str #Need to assign coordinates... |
|
194 |
tmp<-subset(data_v_res,select=date_selected) |
|
195 |
|
|
196 |
data_v_res$idx<-1:length(data_v_res) |
|
197 |
#Plotting results... |
|
198 |
|
|
199 |
data_v_plot<-subset(data_v_res,!is.na(mae),"mae") |
|
286 |
|
|
287 |
### usef the function... |
|
288 |
data_v_resf<-calculate_extremes_stat(tabf_resmod9_locs,tab_locs) # call function |
|
289 |
data_v_resc<-calculate_extremes_stat(tabc_resmod9_locs,tab_locs) # call functio |
|
290 |
|
|
291 |
#Some plotting of results... |
|
292 |
data_v_plot<-subset(data_v_resf,!is.na(data_v_resf$mae),"mae") |
|
200 | 293 |
bubble(data_v_plot,"mae",add=TRUE) |
201 |
spplot(data_v_res,"mae") |
|
294 |
spplot(data_v_resf,"mae")
|
|
202 | 295 |
|
203 | 296 |
data_v_plot<-subset(data_v_res,!is.na(c1),"c1") |
204 |
bubble(data_v_res,"c3")
|
|
297 |
bubble(data_v_res,"c1")
|
|
205 | 298 |
plot(data_v_res[5,]) |
206 | 299 |
plot(reg_outline,add=TRUE) |
207 | 300 |
#Figure on 11/07 |
208 |
s_range<-c(minValue(rast_fus_pred),maxValue(rast_fus_pred)) #stack min and max |
|
209 |
s_range<-c(min(s_range),max(s_range)) |
|
210 |
col_breaks <- pretty(s_range, n=50) |
|
211 |
lab_breaks <- pretty(s_range, n=5) |
|
212 |
temp_colors <- colorRampPalette(c('blue', 'white', 'red')) |
|
213 |
plot(rast_fus_pred, breaks=col_breaks, col=temp_colors(49), |
|
214 |
axis=list(at=lab_breaks, labels=lab_breaks)) |
|
301 |
s.range <- c(min(minValue(rast_fus_pred)), max(maxValue(rast_fus_pred))) |
|
302 |
col.breaks <- pretty(s.range, n=50) |
|
303 |
lab.breaks <- pretty(s.range, n=5) |
|
304 |
temp.colors <- colorRampPalette(c('blue', 'white', 'red')) |
|
305 |
|
|
306 |
#Training plot |
|
307 |
plot(rast_fus_pred, breaks=col.breaks, col=temp.colors(length(col.breaks)-1), |
|
308 |
axis=list(at=lab.breaks, labels=lab.breaks)) |
|
215 | 309 |
plot(reg_outline,add=TRUE) |
216 | 310 |
|
217 | 311 |
plot(data_v_res, pch=1,cex=sqrt(data_v_res$c1)/5,add=TRUE) |
... | ... | |
220 | 314 |
# Calculate RMSE and MAE for each station over 365 dates and plot it on a map as well as on a scatterplot |
221 | 315 |
#Look at the number of observation. |
222 | 316 |
|
317 |
#bubble(data_v_res,"X20101230",na.rm=T,do.sqrt=TRUE) |
|
223 | 318 |
|
224 | 319 |
#1) Digitize features from high resolution imagery and see if you can see differences in the way it is detected in CAI and fusion method... |
225 | 320 |
#in particular how much variation there is in such polygons... |
226 | 321 |
#Polygon to digitize: valley, crop area, mountain...Show that CAI does not capture crop as well as Fusion |
227 |
#2) Select transect with slope changes and examining variation in temperatures...: calculate average MAE on the transect for examle river
|
|
228 |
#3) Land cover: examine MAE per land cover...for CAI and Fusion, LST bias and ecoregions...
|
|
229 |
#4) Look at differences...regions |
|
230 |
#5) Summarize by season/month |
|
322 |
#2) Transect with slope changes and examining variation in temperatures...: calculate average MAE on the transect for examle river
|
|
323 |
#3) X Land cover: examine MAE per land cover...for CAI and Fusion, LST bias and ecoregions...: boxplots at
|
|
324 |
#4) Look at differences...regions with box plot of MAE inside and outside regions of differences...
|
|
325 |
#5) X Summarize by season/month
|
|
231 | 326 |
#6) PCA on differences and CAI-Fusion |
232 |
#7) PCA on residuals... |
|
233 |
#8) Plot residuals at station by buble plot and kriging... |
|
327 |
#7) X PCA on residuals...
|
|
328 |
#8) X Plot residuals at station by buble plot and kriging...
|
|
234 | 329 |
#9) PCA MAE and other variables ... |
235 | 330 |
|
236 |
### RESIDUALS FOR SPECIFIC DATE
|
|
331 |
### PLOTS OF RESIDUALS AND COVARIATES FOR SPECIFIC DATE
|
|
237 | 332 |
|
333 |
#This should be in a loop... |
|
334 |
setwd(path) |
|
238 | 335 |
date_selected<-"20100103" |
336 |
dates<-c("20100103","20100901") |
|
337 |
|
|
338 |
for (i in 1:length(dates)){ |
|
339 |
date_selected<-dates[i] |
|
340 |
k<-match(date_selected,sampling_date_list) |
|
341 |
names(gam_fus$gam_fus_mod[[k]]) #Show the name structure of the object/list |
|
342 |
|
|
343 |
#Extract the training and testing information for the given date... |
|
344 |
data_sf<-gam_fus$gam_fus_mod[[k]]$data_s #object for the first date...20100103 |
|
345 |
data_vf<-gam_fus$gam_fus_mod[[k]]$data_v #object for the first date...20100103 |
|
346 |
data_sc<-gam_cai$gam_CAI_mod[[k]]$data_s #object for the first date...20100103 |
|
347 |
data_vc<-gam_cai$gam_CAI_mod[[k]]$data_v #object for the first date...20100103 |
|
348 |
|
|
349 |
#Join information extracted using the function previously |
|
350 |
id_selected<-intersect(data_v_resf$id,data_vf$id) #Match station using their official IDSs |
|
351 |
pos<-match(id_selected,data_v_resf$id) |
|
352 |
tmp<-as.data.frame(data_v_resf[pos,c("id","idx","mae","c1","c2","c3")]) |
|
353 |
tmp<-tmp[,c("id","idx","mae","c1","c2","c3")] |
|
354 |
data_vf<-merge(data_vf,tmp,by="id") |
|
355 |
id_selected<-intersect(data_v_resc$id,data_vc$id) #Match station using their official IDSs |
|
356 |
pos<-match(id_selected,data_v_resc$id) |
|
357 |
tmp<-as.data.frame(data_v_resc[pos,c("id","idx","mae","c1","c2","c3")]) |
|
358 |
tmp<-tmp[,c("id","idx","mae","c1","c2","c3")] |
|
359 |
data_vc<-merge(data_vc,tmp,by="id") |
|
360 |
#Turn the data frame back into a spdf |
|
361 |
coords<- data_vf[,c('x_OR83M','y_OR83M')] |
|
362 |
coordinates(data_vf)<-coords |
|
363 |
proj4string(data_vf)<-proj_str #Need to assign coordinates... |
|
364 |
coords<- data_vc[,c('x_OR83M','y_OR83M')] |
|
365 |
coordinates(data_vc)<-coords |
|
366 |
proj4string(data_vc)<-proj_str #Need to assign coordinates... |
|
367 |
|
|
368 |
#Prepare for plotting |
|
369 |
X11(width=16,height=9) |
|
370 |
par(mfrow=c(1,2)) |
|
371 |
|
|
372 |
#Select background image for plotting and Plot validation residuals for fusion on 20100103 |
|
373 |
s.range <- c(min(minValue(rast_fus_pred)), max(maxValue(rast_fus_pred))) |
|
374 |
col.breaks <- pretty(s.range, n=50) |
|
375 |
lab.breaks <- pretty(s.range, n=5) |
|
376 |
temp.colors <- colorRampPalette(c('blue', 'white', 'red')) |
|
377 |
|
|
378 |
#Training plot |
|
379 |
plot(rast_fus_pred, breaks=col.breaks, col=temp.colors(length(col.breaks)-1), |
|
380 |
axis=list(at=lab.breaks, labels=lab.breaks)) |
|
381 |
plot(data_sf,cex=1, add=TRUE) |
|
382 |
plot(reg_outline,add=TRUE) |
|
383 |
title(paste("Training stations",date_selected,sep=" ")) |
|
384 |
#Testing plot |
|
385 |
plot(rast_fus_pred, breaks=col.breaks, col=temp.colors(length(col.breaks)-1), |
|
386 |
axis=list(at=lab.breaks, labels=lab.breaks)) |
|
387 |
plot(data_vf,cex=1,pch=1,add=TRUE) |
|
388 |
plot(reg_outline,add=TRUE) |
|
389 |
text(data_vf,labels=data_vf$idx,pos=3,cex=1.3) |
|
390 |
title(paste("Testing stations",date_selected,sep=" ")) |
|
391 |
|
|
392 |
#savePlot here... |
|
393 |
savePlot(paste("fig1_training_testing_",date_selected,out_prefix,".png", sep=""), type="png") |
|
394 |
|
|
395 |
#plot(rast_fus_pred, breaks=col.breaks, col=temp.colors(length(col.breaks)-1), |
|
396 |
# axis=list(at=lab.breakts, labels=lab.breaks)) |
|
397 |
#plot(data_vf,cex=sqrt(data_vf$res_mod7),pch=1,add=TRUE) |
|
398 |
#plot(reg_outline,add=TRUE) |
|
399 |
#text(data_vf,labels=data_vf$idx,pos=3) |
|
400 |
|
|
401 |
#X11(width=16,height=9) |
|
402 |
#par(mfrow=c(1,2)) |
|
403 |
|
|
404 |
# CREATE A FUNCTION TO AUTOMATE THE PLOT: improve code here |
|
405 |
#Plot residuals |
|
406 |
y_range<-range(c(data_vf$res_mod7,data_vc$res_mod9)) |
|
407 |
plot(data_vf$idx,data_vf$res_mod7, ylab="Residuals", xlab="Index", ylim=y_range) |
|
408 |
text(data_vf$idx,data_vf$res_mod7,labels=data_vf$idx,pos=3) |
|
409 |
grid(lwd=0.5,col="black") |
|
410 |
title(paste("Testing stations residuals fusion",date_selected,sep=" ")) |
|
411 |
plot(data_vc$idx,data_vc$res_mod9,ylab="Residuals", xlab="Index", ylim=y_range) |
|
412 |
text(data_vc$idx,data_vc$res_mod9,data_vc$idx,labels=data_vc$idx,pos=3) |
|
413 |
grid(lwd=0.5, col="black") |
|
414 |
title(paste("Testing stations residuals CAI",date_selected,sep=" ")) |
|
415 |
#savePlot here... |
|
416 |
savePlot(paste("fig2_testing_residuals_fusion_CAI_",date_selected,out_prefix,".png", sep=""), type="png") |
|
417 |
|
|
418 |
#Plot predicted vs observed |
|
419 |
y_range<-range(c(data_vf$pred_mod7,data_vc$pred_mod9)) |
|
420 |
x_range<-range(c(data_vf$dailyTmax,data_vc$dailyTamx)) |
|
421 |
plot(data_vf$pred_mod7,data_vf$dailyTmax, ylab="Observed daily tmax (C)", xlab="Predicted daily tmax (C)", |
|
422 |
ylim=y_range,xlim=x_range) |
|
423 |
text(data_vf$pred_mod7,data_vf$dailyTmax,labels=data_vf$idx,pos=3) |
|
424 |
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
425 |
grid(lwd=0.5,col="black") |
|
426 |
title(paste("Testing stations tmax fusion vs daily tmax",date_selected,sep=" ")) |
|
427 |
plot(data_vc$pred_mod9,data_vc$dailyTmax,ylab="Observed daily tmax (C)", xlab="Predicted daily tmax (C)", |
|
428 |
ylim=y_range,xlim=x_range) |
|
429 |
text(data_vc$pred_mod9,data_vc$dailyTmax,labels=data_vc$idx,pos=3) |
|
430 |
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
431 |
grid(lwd=0.5, col="black") |
|
432 |
title(paste("Testing stations tmax CAI vs daily tmax",date_selected,sep=" ")) |
|
433 |
#savePlot here... |
|
434 |
savePlot(paste("fig3_testing_predicted_observed_fusion_CAI_",date_selected,out_prefix,".png", sep=""), type="png") |
|
435 |
|
|
436 |
###########Plot residuals and covariates |
|
437 |
##Elevation and residulas |
|
438 |
y_range<-range(c(data_vf$res_mod7,data_vc$res_mod9)) |
|
439 |
x_range<-range(c(data_vf$ELEV_SRTM,data_vc$ELEV_SRTM)) |
|
440 |
plot(data_vf$ELEV_SRTM,data_vf$res_mod7, ylab="Residuals", xlab="ELEV_SRTM (m) ", |
|
441 |
ylim=y_range, xlim=x_range) |
|
442 |
text(data_vf$ELEV_SRTM,data_vf$res_mod7,labels=data_vf$idx,pos=3) |
|
443 |
grid(lwd=0.5,col="black") |
|
444 |
title(paste("Testing stations residuals fusion vs Elevation",date_selected,sep=" ")) |
|
445 |
plot(data_vc$ELEV_SRTM,data_vc$res_mod9, ylab="Residuals", xlab="ELEV_SRTM (m) ", |
|
446 |
ylim=y_range, xlim=x_range) |
|
447 |
text(data_vc$ELEV_SRTM,data_vc$res_mod9,labels=data_vc$idx,pos=3) |
|
448 |
grid(lwd=0.5,col="black") |
|
449 |
title(paste("Testing stations residuals CAI vs Elevation",date_selected,sep=" ")) |
|
450 |
savePlot(paste("fig4_testing_residuals_fusion_CAI_Elev_",date_selected,out_prefix,".png", sep=""), type="png") |
|
451 |
|
|
452 |
##Forest (LC1) and residuals |
|
453 |
y_range<-range(c(data_vf$res_mod7,data_vc$res_mod9)) |
|
454 |
x_range<-range(c(data_vf$LC1,data_vc$LC1)) |
|
455 |
plot(data_vf$LC1,data_vf$res_mod7,ylab="Residuals", xlab="LC1 (%)", |
|
456 |
ylim=y_range, xlim=x_range) |
|
457 |
text(data_vf$LC1,data_vf$res_mod7,labels=data_vf$idx,pos=3) |
|
458 |
grid(lwd=0.5,col="black") |
|
459 |
title(paste("Testing stations residuals CAI vs LC1 (forest)",date_selected,sep=" ")) |
|
460 |
plot(data_vc$LC1,data_vc$res_mod9,ylab="Residuals", xlab="LC1 (%)", |
|
461 |
ylim=y_range, xlim=x_range) |
|
462 |
text(data_vc$LC1,data_vc$res_mod9,labels=data_vc$idx,pos=3) |
|
463 |
grid(lwd=0.5,col="black") |
|
464 |
title(paste("Testing stations residuals CAI vs LC1(forest)",date_selected,sep=" ")) |
|
465 |
savePlot(paste("fig5_testing_residuals_fusion_CAI_LC1_",date_selected,out_prefix,".png", sep=""), type="png") |
|
466 |
|
|
467 |
##Grass (LC3) and residuals |
|
468 |
y_range<-range(c(data_vf$res_mod7,data_vc$res_mod9)) |
|
469 |
x_range<-range(c(data_vf$LC3,data_vc$LC3)) |
|
470 |
plot(data_vf$LC3,data_vf$res_mod7,ylab="Residuals", xlab="LC3 (%)", |
|
471 |
ylim=y_range, xlim=x_range) |
|
472 |
text(data_vf$LC3,data_vf$res_mod7,labels=data_vf$idx,pos=3) |
|
473 |
grid(lwd=0.5,col="black") |
|
474 |
title(paste("Testing stations residuals CAI vs LC3 (grass)",date_selected,sep=" ")) |
|
475 |
plot(data_vc$LC3,data_vc$res_mod9,ylab="Residuals", xlab="LC3 (%)", |
|
476 |
ylim=y_range, xlim=x_range) |
|
477 |
text(data_vc$LC3,data_vc$res_mod9,labels=data_vc$idx,pos=3) |
|
478 |
grid(lwd=0.5,col="black") |
|
479 |
title(paste("Testing stations residuals CAI vs LC3(grass)",date_selected,sep=" ")) |
|
480 |
savePlot(paste("fig6_testing_residuals_fusion_CAI_LC3_",date_selected,out_prefix,".png", sep=""), type="png") |
|
481 |
|
|
482 |
#LSTD_bias and residuals |
|
483 |
y_range<-range(c(data_vf$res_mod7,data_vc$res_mod9)) |
|
484 |
x_range<-range(c(data_vf$LSTD_bias,data_vc$LSTD_bias)) |
|
485 |
plot(data_vf$LSTD_bias,data_vf$res_mod7) |
|
486 |
text(data_vf$LSTD_bias,data_vf$res_mod7,labels=data_vf$idx,pos=3) |
|
487 |
grid(lwd=0.5,col="black") |
|
488 |
title(paste("Testing stations LST bias vs residuals",date_selected,sep=" ")) |
|
489 |
plot(data_sf$LSTD_bias,data_sf$res_mod7) |
|
490 |
#text(data_sf$LSTD_bias,data_sf$res_mod7,labels=data_vf$idx,pos=3) # Also labels for training? |
|
491 |
grid(lwd=0.5,col="black") |
|
492 |
title(paste("Training stations LST bias vs residuals",date_selected,sep=" ")) |
|
493 |
savePlot(paste("fig7_testing_training_residuals_fusion_LSTD_bias_",date_selected,out_prefix,".png", sep=""), type="png") |
|
494 |
|
|
495 |
#LSTD vs TMax and residuals |
|
496 |
y_range<-range(c(data_vf$TMax,data_vc$TMax)) |
|
497 |
x_range<-range(c(data_vf$LSTD_bias,data_vc$LSTD_bias)) |
|
498 |
plot(data_vf$LSTD_bias,data_vf$res_mod7,xlab="LST bias", ylab="Residuals") |
|
499 |
text(data_vf$LSTD_bias,data_vf$res_mod7,labels=data_vf$idx,pos=3) |
|
500 |
grid(lwd=0.5,col="black") |
|
501 |
title(paste("Testing stations LST bias vs fusion residuals",date_selected,sep=" ")) |
|
502 |
plot((data_vf$LST-273.16),data_vf$TMax,xlab="LST", ylab="TMax (monthly tmax in C)") |
|
503 |
text((data_vf$LST-273.16),data_vf$TMax,labels=data_vf$idx,pos=3) |
|
504 |
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
505 |
grid(lwd=0.5,col="black") |
|
506 |
title(paste("Testing stations LST vs TMax ",date_selected,sep=" ")) |
|
507 |
savePlot(paste("fig8_testing_TMax_fusion_",date_selected,out_prefix,".png", sep=""), type="png") |
|
508 |
dev.off() |
|
509 |
} |
|
510 |
|
|
511 |
##Check...difference at diff stations... |
|
239 | 512 |
|
240 |
k<-match(date_selected,sampling_date_list) |
|
241 |
names(gam_fus$gam_fus_mod[[k]]) #Show the name structure of the object/list |
|
242 | 513 |
|
514 |
#Extract the training and testing information for the given date... |
|
243 | 515 |
data_sf<-gam_fus$gam_fus_mod[[k]]$data_s #object for the first date...20100103 |
244 | 516 |
data_vf<-gam_fus$gam_fus_mod[[k]]$data_v #object for the first date...20100103 |
245 | 517 |
data_sc<-gam_cai$gam_CAI_mod[[k]]$data_s #object for the first date...20100103 |
246 | 518 |
data_vc<-gam_cai$gam_CAI_mod[[k]]$data_v #object for the first date...20100103 |
247 | 519 |
|
248 |
id_selected<-intersect(data_v_res$id,data_vf$id) |
|
249 |
pos<-match(id_selected,data_v_res$id) |
|
250 |
tmp<-as.data.frame(data_v_res[pos,c("id","idx","mae","c1","c2","c3")]) |
|
520 |
#Join information extracted using the function previously |
|
521 |
id_selected<-intersect(data_v_resf$id,data_vf$id) #Match station using their official IDSs |
|
522 |
pos<-match(id_selected,data_v_resf$id) |
|
523 |
tmp<-as.data.frame(data_v_resf[pos,c("id","idx","mae","c1","c2","c3")]) |
|
251 | 524 |
tmp<-tmp[,c("id","idx","mae","c1","c2","c3")] |
252 | 525 |
data_vf<-merge(data_vf,tmp,by="id") |
253 |
|
|
526 |
id_selected<-intersect(data_v_resc$id,data_vc$id) #Match station using their official IDSs |
|
527 |
pos<-match(id_selected,data_v_resc$id) |
|
528 |
tmp<-as.data.frame(data_v_resc[pos,c("id","idx","mae","c1","c2","c3")]) |
|
529 |
tmp<-tmp[,c("id","idx","mae","c1","c2","c3")] |
|
530 |
data_vc<-merge(data_vc,tmp,by="id") |
|
531 |
#Turn the data frame back into a spdf |
|
254 | 532 |
coords<- data_vf[,c('x_OR83M','y_OR83M')] |
255 | 533 |
coordinates(data_vf)<-coords |
256 | 534 |
proj4string(data_vf)<-proj_str #Need to assign coordinates... |
535 |
coords<- data_vc[,c('x_OR83M','y_OR83M')] |
|
536 |
coordinates(data_vc)<-coords |
|
537 |
proj4string(data_vc)<-proj_str #Need to assign coordinates... |
|
538 |
|
|
539 |
################################################### |
|
540 |
############## RESIDUALS BY TIME!!! ############## |
|
541 |
|
|
542 |
#tab_resmod9 |
|
543 |
for (i in 1:nrow(tab_resmod9f)){ |
|
544 |
date<-strptime(tab_resmod9f$date[i], "%Y%m%d") # interpolation date being processed, converting the string using specific format |
|
545 |
month<-as.integer(strftime(date, "%m")) |
|
546 |
tab_resmod9f$month[i]<-month |
|
547 |
} |
|
257 | 548 |
|
258 |
#Select background image for plotting and Plot validation residuals for fusion on 20100103 |
|
259 |
s.range <- c(min(minValue(rast_fus_pred)), max(maxValue(rast_fus_pred))) |
|
260 |
col.breaks <- pretty(s.range, n=50) |
|
261 |
lab.breaks <- pretty(s.range, n=5) |
|
262 |
temp.colors <- colorRampPalette(c('blue', 'white', 'red')) |
|
263 |
|
|
264 |
#Training plot |
|
265 |
plot(rast_fus_pred, breaks=col.breaks, col=temp.colors(length(col.breaks)-1), |
|
266 |
axis=list(at=lab.breaks, labels=lab.breaks)) |
|
267 |
plot(data_sf,cex=1,main="Training stations", add=TRUE) |
|
268 |
plot(reg_outline,add=TRUE) |
|
269 |
|
|
270 |
#savePlot... |
|
271 |
|
|
272 |
#Testing plot |
|
273 |
plot(rast_fus_pred, breaks=col.breaks, col=temp.colors(length(col.breaks)-1), |
|
274 |
axis=list(at=lab.breaks, labels=lab.breaks)) |
|
275 |
plot(data_vf,cex=1,pch=1,add=TRUE) |
|
276 |
plot(reg_outline,add=TRUE) |
|
277 |
text(data_vf,labels=data_vf$idx,pos=3) |
|
278 |
#savePlot... |
|
279 |
|
|
280 |
#Plot residuals |
|
281 |
res_modvf<-res_mod7 |
|
282 |
|
|
283 |
plot(data_vf$res_mod7) |
|
284 |
text(data_vf$res_mod7,labels=data_vf$idx,pos=3) |
|
285 |
|
|
286 |
#Plot residuals |
|
287 |
res_modvf<-res_mod7 |
|
288 |
|
|
289 |
plot(data_vf$res_mod7) |
|
290 |
text(data_vf$res_mod7,labels=data_vf$idx,pos=3) |
|
291 |
|
|
292 |
|
|
293 |
|
|
294 |
|
|
295 |
|
|
296 |
#NOW USE RESHAPE TO CREATE TABLE.... |
|
549 |
mae_function<-function(res){mae<-mean(abs(res),na.rm=T)} |
|
550 |
tab_test<-melt(tab_resmod9f, |
|
551 |
measure=c("res_mod7","x_OR83M","y_OR83M"), |
|
552 |
id=c("id","month"), |
|
553 |
na.rm=F) |
|
554 |
tab_test_cast<-cast(data=tab_test,formula=id~month~variable,mean) |
|
555 |
table_id_month<-tab_test_cast[,,1] |
|
556 |
tab_month_cast<-cast(data=tab_test,formula=month~variable,mean) |
|
557 |
tab_month_cast<-cast(data=tab_test,formula=month~variable,mae_function) |
|
558 |
X11() |
|
559 |
barplot(tab_month_cast$res_mod7, |
|
560 |
names.arg=1:12,cex.names=0.8, #names of the teleconnections indices and size of fonts of axis labes |
|
561 |
xlab="Month", # font.lab is 2 to make the font bold |
|
562 |
ylab="MAE",font.lab=2) |
|
563 |
|
|
564 |
savePlot(paste("fig9_mae_fusion_per_season_",date_selected,out_prefix,".png", sep=""), type="png") |
|
565 |
|
|
566 |
plot(table_id_month[55,]) |
|
567 |
plot(table_id_month[5,]) |
|
568 |
|
|
569 |
#tab_resmod9c |
|
570 |
|
|
571 |
for (i in 1:nrow(tab_resmod9c)){ |
|
572 |
date<-strptime(tab_resmod9c$date[i], "%Y%m%d") # interpolation date being processed, converting the string using specific format |
|
573 |
month<-as.integer(strftime(date, "%m")) |
|
574 |
tab_resmod9c$month[i]<-month |
|
575 |
} |
|
297 | 576 |
|
298 |
#updated the analysis |
|
577 |
mae_function<-function(res){mae<-mean(abs(res),na.rm=T)} |
|
578 |
tab_test<-melt(tab_resmod9c, |
|
579 |
measure=c("res_mod9","x_OR83M","y_OR83M"), |
|
580 |
id=c("id","month"), |
|
581 |
na.rm=F) |
|
582 |
tab_test_cast<-cast(data=tab_test,formula=id~month~variable,mean) |
|
583 |
table_id_month_cai<-tab_test_cast[,,1] |
|
584 |
tab_month_cast_cai<-cast(data=tab_test,formula=month~variable,mean) |
|
585 |
tab_month_cast_cai<-cast(data=tab_test,formula=month~variable,mae_function) |
|
299 | 586 |
|
300 |
"tmax_predicted*20100103_30_1_365d_GAM_fusion_all_lstd_10242012.rst" |
|
301 |
oldpath<-getwd() |
|
302 |
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_GAM" |
|
303 |
setwd(path) |
|
304 |
date_selected<-"20100103" |
|
305 |
#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 |
|
306 |
file_pat<-glob2rx(paste("*tmax_pred*",date_selected,"*_365d_GAM_fusion_all_lstd_10242012.rst",sep="")) #Search for files in relation to fusion |
|
307 |
lf_fus1a<-list.files(pattern=file_pat) #Search for files in relation to fusion |
|
587 |
barplot(tab_month_cast_cai$res_mod9, |
|
588 |
names.arg=1:12,cex.names=0.8, #names of the teleconnections indices and size of fonts of axis labes |
|
589 |
xlab="Month", # font.lab is 2 to make the font bold |
|
590 |
ylab="MAE",font.lab=2) |
|
591 |
grid(nx=12,ny=10) |
|
308 | 592 |
|
309 |
rast_fus1a<-stack(lf_fus1a) |
|
310 |
rast_fus1a<-mask(rast_fus1a,mask_ELEV_SRTM) |
|
593 |
savePlot(paste("fig10_mae_cai_per_season_",date_selected,out_prefix,".png", sep=""), type="png") |
|
311 | 594 |
|
312 |
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_07192012_GAM" |
|
313 |
setwd(path) |
|
314 |
date_selected<-"20100103" |
|
315 |
#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 |
|
316 |
#lf_fus1<-list.files(pattern=paste("*.tmax_predicted.*.",date_selected,".*._365d_GAM_fusion_lstd_10062012.rst$",sep="")) |
|
317 |
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_fusion_lstd_10062012.rst",sep="")) #Search for files in relation to fusion |
|
318 |
lf_fus1s<-list.files(pattern=file_pat) #Search for files in relation to fusion |
|
319 |
|
|
320 |
rast_fus1s<-stack(lf_fus1s) |
|
321 |
rast_fus1s<-mask(rast_fus1s,mask_ELEV_SRTM) |
|
322 |
|
|
323 |
s_range<-c(minValue(rast_fusa1),maxValue(rast_fusa1)) #stack min and max |
|
324 |
s_range<-c(min(s_range),max(s_range)) |
|
325 |
col_breaks <- pretty(s_range, n=50) |
|
326 |
lab_breaks <- pretty(s_range, n=5) |
|
327 |
temp_colors <- colorRampPalette(c('blue', 'white', 'red')) |
|
328 |
X11(width=18,height=12) |
|
329 |
par(mfrow=c(3,3)) |
|
330 |
for (k in 1:length(lf_fus1a)){ |
|
331 |
fus1a_r<-raster(rast_fus1a,k) |
|
332 |
plot(fus1a_r, breaks=col_breaks, col=temp_colors(length(col_breaks)-1), |
|
333 |
axis=list(at=lab_breaks, labels=lab_breaks)) |
|
334 |
} |
|
595 |
plot(table_id_month_cai[55,]) |
|
596 |
plot(table_id_month_cai[5,]) |
|
335 | 597 |
|
598 |
###### |
|
599 |
overlay(data_v_resf,diff_r) |
|
600 |
plot(avg_LC_rec) |
|
336 | 601 |
#Wihtout setting range |
337 | 602 |
s_range<-c(-12,18) |
338 | 603 |
#col_breaks <- pretty(s_range, n=50) |
... | ... | |
350 | 615 |
axis=list(at=lab_breaks, labels=lab_breaks),main=rname[k]) |
351 | 616 |
} |
352 | 617 |
|
353 |
#Wihtout setting range |
|
618 |
### Wihtout setting range
|
|
354 | 619 |
s_range<-c(-12,23) |
355 | 620 |
#col_breaks <- pretty(s_range, n=50) |
356 | 621 |
#lab_breaks <- pretty(s_range, n=5) |
... | ... | |
366 | 631 |
plot(fus1s_r, breaks=col_breaks, col=temp_colors(col_breaks-1), |
367 | 632 |
axis=list(at=lab_breaks, labels=lab_breaks),main=rname[k]) |
368 | 633 |
} |
369 |
|
|
634 |
|
|
635 |
############### NOW USE RESHAPE TO CREATE TABLE....########## |
|
636 |
# DO A PCA HERE...also look at the difference between prediction at stations... |
|
637 |
|
|
638 |
var_pat<-glob2rx("*2010*") #Search for files in relation to fusion |
|
639 |
var<-match("*2010*",names(data_v_resf)) #Search for files in relation to fusion |
|
640 |
xp<-data_v_resf[,names(data_v_resf)[10:374]] # subsetting the original training data, note that that this is a "SpatialPointsDataFrame" |
|
641 |
|
|
642 |
#xp<-data_s[,var] # subsetting the original training data, note that that this is a "SpatialPointsDataFrame" |
|
643 |
xp<-as.data.frame(xp) |
|
644 |
dropsc<-c("x_OR83M","y_OR83M") #columns to be removed |
|
645 |
xp<-xp[,!(names(xp) %in% dropsc)] # dropping columns |
|
646 |
|
|
647 |
X<-as.matrix(na.omit(xp)) |
|
648 |
|
|
649 |
PCA9var<-prcomp(~.,data=xp, retx = TRUE, center= TRUE, scale = TRUE, na.action=na.omit) |
|
650 |
E<-matrix() |
|
651 |
E<-PCA9var$rotation #E are the eigenvectors of the standardized PCA base on xp data.frame |
|
652 |
Xs<-scale(X) #Use scale to standardize the original data matrix (center on mean and divide by standard deviation) |
|
653 |
Xpc<- Xs %*% E #Rotate the original standardized matrix Xs by the eigenvectors from the standardized PCA |
|
654 |
plot(PCA9var) |
|
655 |
#plot(PCA9var$scores) |
|
656 |
loadings<-cor(X, Xpc) #Correlation between the original variables and the principal components... |
|
657 |
loadings<-as.data.frame(loadings) |
|
658 |
#quartz(width=6, height=6) |
|
659 |
plot(PC2~PC1, xlab= "PC1", ylab= "PC2", xlim=c(-1, 1), ylim=c(-1,1), pch =16, width=1, |
|
660 |
main= paste("This is PCA for ",date_selected, sep=""), height=1,asp=1,data=loadings) #Add aspect options to get |
|
661 |
axis(1,pos=0) |
|
662 |
axis(2,pos=0) |
|
663 |
text(loadings$PC1,loadings$PC2,rownames(loadings), adj = c(0,0),offset=2, col="red",cex=1.2) |
|
664 |
draw.circle(0,0,radius=1) |
|
665 |
|
|
370 | 666 |
|
371 | 667 |
|
Also available in: Unified diff
Methods comp part4- task#491-, residuals analses, cleaning, added functions, residuals by land cover and temporal profiles