Revision 03aa4873
Added by Benoit Parmentier almost 12 years ago
| climate/research/oregon/interpolation/results_interpolation_date_output_analyses.R | ||
|---|---|---|
| 53 | 53 |
validation_obj<-load_obj("gam_fus_validation_mod_365d_GAM_fus5_all_lstd_02202013.RData")
|
| 54 | 54 |
clim_obj<-load_obj("gamclim_fus_mod_365d_GAM_fus5_all_lstd_02202013.RData")
|
| 55 | 55 |
|
| 56 |
## Read covariate stack... |
|
| 57 |
rnames <-c("x","y","lon","lat","N","E","N_w","E_w","elev","slope","aspect","CANHEIGHT","DISTOC")
|
|
| 58 |
lc_names<-c("LC1","LC2","LC3","LC4","LC5","LC6","LC7","LC8","LC9","LC10","LC11","LC12")
|
|
| 59 |
lst_names<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12",
|
|
| 60 |
"nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08", |
|
| 61 |
"nobs_09","nobs_10","nobs_11","nobs_12") |
|
| 62 |
|
|
| 63 |
covar_names<-c(rnames,lc_names,lst_names) |
|
| 64 |
|
|
| 65 |
s_raster<-stack(infile3) #read in the data stack |
|
| 66 |
names(s_raster)<-covar_names #Assigning names to the raster layers: making sure it is included in the extraction |
|
| 67 |
|
|
| 68 |
## Figure 0: study area based on LC12 (water) mask |
|
| 69 |
|
|
| 70 |
LC_mask<-subset(s_raster,"LC12") |
|
| 71 |
LC_mask[LC_mask==100]<-NA |
|
| 72 |
LC_mask <- LC_mask < 100 |
|
| 73 |
LC_mask_rec<-LC_mask |
|
| 74 |
LC_mask_rec[is.na(LC_mask_rec)]<-0 |
|
| 75 |
|
|
| 76 |
png(paste("Study_area_",
|
|
| 77 |
out_prefix,".png", sep="")) |
|
| 78 |
plot(LC_mask_rec,legend=FALSE,col=c("black","red"))
|
|
| 79 |
legend("topright",legend=c("Outside","Inside"),title="Study area",
|
|
| 80 |
pt.cex=0.9,fill=c("black","red"),bty="n")
|
|
| 81 |
title("Study area")
|
|
| 82 |
dev.off() |
|
| 83 |
|
|
| 56 | 84 |
#determine index position matching date selected |
| 57 | 85 |
|
| 58 | 86 |
i_dates<-vector("list",length(date_selected))
|
| ... | ... | |
| 82 | 110 |
data_v<-validation_obj[[index]]$data_v |
| 83 | 111 |
data_s<-validation_obj[[index]]$data_s |
| 84 | 112 |
data_month<-clim_obj[[index]]$data_month |
| 113 |
formulas<-clim_obj[[index]]$formulas |
|
| 85 | 114 |
|
| 86 | 115 |
#Adding layer LST to the raster stack of covariates |
| 87 | 116 |
#The names of covariates can be changed... |
| 88 |
rnames <-c("x","y","lon","lat","N","E","N_w","E_w","elev","slope","aspect","CANHEIGHT","DISTOC")
|
|
| 89 |
lc_names<-c("LC1","LC2","LC3","LC4","LC5","LC6","LC7","LC8","LC9","LC10","LC11","LC12")
|
|
| 90 |
lst_names<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12",
|
|
| 91 |
"nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08", |
|
| 92 |
"nobs_09","nobs_10","nobs_11","nobs_12") |
|
| 93 |
|
|
| 94 |
covar_names<-c(rnames,lc_names,lst_names) |
|
| 95 |
|
|
| 96 |
s_raster<-stack(infile3) #read in the data stack |
|
| 97 |
names(s_raster)<-covar_names #Assigning names to the raster layers: making sure it is included in the extraction |
|
| 98 | 117 |
|
| 99 | 118 |
LST_month<-paste("mm_",month,sep="") # name of LST month to be matched
|
| 100 | 119 |
pos<-match("LST",layerNames(s_raster)) #Find the position of the layer with name "LST", if not present pos=NA
|
| ... | ... | |
| 143 | 162 |
#This is for mod_kr!! add other models later... |
| 144 | 163 |
png(paste("Predicted_tmax_versus_observed_scatterplot_",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,
|
| 145 | 164 |
out_prefix,".png", sep="")) |
| 146 |
plot(data_s$mod_kr~data_s[[y_var_name]],xlab=paste("Actual daily for",datelabel),ylab="Pred daily")
|
|
| 165 |
#plot(data_s$mod_kr~data_s[[y_var_name]],xlab=paste("Actual daily for",datelabel),ylab="Pred daily")
|
|
| 166 |
|
|
| 167 |
y_range<-range(c(data_s$mod_kr,data_v$mod_kr),na.rm=T) |
|
| 168 |
x_range<-range(c(data_s[[y_var_name]],data_v[[y_var_name]]),na.rm=T) |
|
| 169 |
col_t<- c("black","red")
|
|
| 170 |
pch_t<- c(1,2) |
|
| 171 |
plot(data_s$mod_kr,data_s[[y_var_name]], |
|
| 172 |
xlab=paste("Actual daily for",datelabel),ylab="Pred daily",
|
|
| 173 |
ylim=y_range,xlim=x_range,col=col_t[1],pch=pch_t[1]) |
|
| 174 |
points(data_v$mod_kr,data_v[[y_var_name]],col=col_t[2],pch=pch_t[2]) |
|
| 175 |
grid(lwd=0.5, col="black") |
|
| 147 | 176 |
#plot(data_v$mod_kr~data_v[[y_var_name]],xlab=paste("Actual daily for",datelabel),ylab="Pred daily")
|
| 148 | 177 |
abline(0,1) |
| 178 |
legend("topleft",legend=c("training","testing"),pch=pch_t,col=col_t,bty="n",cex=0.8)
|
|
| 149 | 179 |
title(paste("Predicted_tmax_versus_observed_scatterplot for",datelabel,sep=" "))
|
| 150 | 180 |
nb_point1<-paste("ns_obs=",length(data_s$TMax)-sum(is.na(data_s[[y_var_name]])),sep="")
|
| 151 | 181 |
rmse_str1<-paste("RMSE= ",format(rmse,digits=3),sep="")
|
| 152 | 182 |
rmse_str2<-paste("RMSE_f= ",format(rmse_f,digits=3),sep="")
|
| 153 | 183 |
|
| 154 | 184 |
#Add the number of data points on the plot |
| 155 |
legend("topleft",legend=c(nb_point1,rmse_str1,rmse_str2),bty="n",cex=0.8)
|
|
| 185 |
legend("bottomright",legend=c(nb_point1,rmse_str1,rmse_str2),bty="n",cex=0.8)
|
|
| 156 | 186 |
dev.off() |
| 157 | 187 |
|
| 158 |
## Figure 4: delta surface and bias |
|
| 159 |
|
|
| 160 |
#Plot bias,delta and prediction? |
|
| 161 |
|
|
| 162 |
#To do |
|
| 163 |
#Delta surface |
|
| 164 |
#png(paste("Delta_surface_LST_TMax_",sampling_dat$date[i],"_",sampling_dat$prop[i],
|
|
| 165 |
# "_",sampling_dat$run_samp[i],out_prefix,".png", sep="")) |
|
| 166 |
#surface(fitdelta,col=rev(terrain.colors(100)),asp=1,main=paste("Interpolated delta for",datelabel,sep=" "))
|
|
| 167 |
#dev.off() |
|
| 168 |
# |
|
| 169 |
#bias_d_rast<-raster("fusion_bias_LST_20100103_30_1_10d_GAM_fus5_all_lstd_02082013.rst")
|
|
| 170 |
#plot(bias_d_rast) |
|
| 171 |
|
|
| 172 | 188 |
## Figure 5: prediction raster images |
| 173 | 189 |
png(paste("Raster_prediction_",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,
|
| 174 | 190 |
out_prefix,".png", sep="")) |
| 175 | 191 |
#paste(metrics_v$pred_mod,format(metrics_v$rmse,digits=3),sep=":") |
| 176 | 192 |
names(rast_pred_temp)<-paste(metrics_v$pred_mod,format(metrics_v$rmse,digits=3),sep=":") |
| 177 |
plot(rast_pred_temp) |
|
| 193 |
#plot(rast_pred_temp) |
|
| 194 |
levelplot(rast_pred_temp) |
|
| 178 | 195 |
dev.off() |
| 179 | 196 |
|
| 180 | 197 |
## Figure 6: training and testing stations used |
| 181 |
|
|
| 198 |
png(paste("Training_testing_stations_map_",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,
|
|
| 199 |
out_prefix,".png", sep="")) |
|
| 182 | 200 |
plot(raster(rast_pred_temp,layer=5)) |
| 183 |
plot(data_s,col="black",cex=1.2,pch=4,add=TRUE) |
|
| 184 |
plot(data_v,col="red",cex=1.2,pch=2,add=TRUE) |
|
| 201 |
plot(data_s,col="black",cex=1.2,pch=2,add=TRUE) |
|
| 202 |
plot(data_v,col="red",cex=1.2,pch=1,add=TRUE) |
|
| 203 |
legend("topleft",legend=c("training stations", "testing stations"),
|
|
| 204 |
cex=1, col=c("black","red"),
|
|
| 205 |
pch=c(2,1),bty="n") |
|
| 206 |
dev.off() |
|
| 185 | 207 |
|
| 186 |
## Figure 7: monthly stations used
|
|
| 208 |
## Figure 7: monthly stations used |
|
| 187 | 209 |
|
| 210 |
png(paste("Monthly_data_study_area",
|
|
| 211 |
out_prefix,".png", sep="")) |
|
| 188 | 212 |
plot(raster(rast_pred_temp,layer=5)) |
| 189 | 213 |
plot(data_month,col="black",cex=1.2,pch=4,add=TRUE) |
| 190 | 214 |
title("Monthly ghcn station in Venezuela for January")
|
| 215 |
dev.off() |
|
| 216 |
|
|
| 217 |
## Figure 8: delta surface and bias |
|
| 218 |
|
|
| 219 |
png(paste("Bias_delta_surface_",sampling_dat$date[i],"_",sampling_dat$prop[i],
|
|
| 220 |
"_",sampling_dat$run_samp[i],out_prefix,".png", sep="")) |
|
| 191 | 221 |
|
| 192 |
## Summarize information for the day: |
|
| 222 |
bias_rast<-stack(clim_obj[[index]]$bias) |
|
| 223 |
delta_rast<-raster(gam_fus_mod[[index]]$delta) #only one delta image!!! |
|
| 224 |
names(delta_rast)<-"delta" |
|
| 225 |
rast_temp_date<-stack(bias_rast,delta_rast) |
|
| 226 |
rast_temp_date<-mask(rast_temp_date,LC_mask,file="test.tif",overwrite=TRUE) |
|
| 227 |
#bias_d_rast<-raster("fusion_bias_LST_20100103_30_1_10d_GAM_fus5_all_lstd_02082013.rst")
|
|
| 228 |
plot(rast_temp_date) |
|
| 229 |
|
|
| 230 |
dev.off() |
|
| 231 |
|
|
| 232 |
#Figure 9: histogram for all images... |
|
| 233 |
|
|
| 234 |
histogram(rast_pred_temp) |
|
| 235 |
|
|
| 236 |
## Summarize information for the day: write out textfile... |
|
| 237 |
|
|
| 238 |
#Number of station per month |
|
| 239 |
#Number of station per day (training, testing,NA) |
|
| 240 |
#metrics_v,metrics_s |
|
| 241 |
# |
|
| 193 | 242 |
|
| 194 | 243 |
# ################ |
| 195 | 244 |
# #PART 2: Region Covariate analyses ### |
| ... | ... | |
| 300 | 349 |
# plot(data_month,col="black",cex=1.2,pch=4,add=TRUE) |
| 301 | 350 |
# title("Monthly ghcn station in Venezuela for 2000-2010")
|
| 302 | 351 |
# |
| 352 |
#png...output? |
|
| 353 |
# plot(interp_area, axes =TRUE) |
|
| 354 |
# plot(stat_reg, pch=1, col="red", cex= 0.7, add=TRUE) |
|
| 355 |
# plot(data_reg,pch=2,col="blue",cex=2,add=TRUE) |
|
| 303 | 356 |
|
| 304 | 357 |
#### End of script #### |
Also available in: Unified diff
Ouptut analyses, 9 figures generated from prediction step for specific date