Revision 5f28f8d6
Added by Benoit Parmentier almost 12 years ago
| climate/research/oregon/interpolation/methods_comparison_assessment_part7.R | ||
|---|---|---|
| 27 | 27 |
library(reshape) |
| 28 | 28 |
library(RCurl) |
| 29 | 29 |
######### Functions used in the script |
| 30 |
# |
|
| 30 | 31 |
|
| 31 | 32 |
load_obj <- function(f) |
| 32 | 33 |
{
|
| ... | ... | |
| 153 | 154 |
######### |
| 154 | 155 |
#loading R objects that might have similar names |
| 155 | 156 |
|
| 156 |
out_prefix<-"_method_comp7_12042012_"
|
|
| 157 |
out_prefix<-"_method_comp7_12102012b_"
|
|
| 157 | 158 |
infile2<-"list_365_dates_04212012.txt" |
| 159 |
infile1<- "ghcn_or_tmax_covariates_06262012_OR83M.shp" #GHCN shapefile containing variables for modeling 2010 |
|
| 160 |
#infile2<-"list_10_dates_04212012.txt" #List of 10 dates for the regression |
|
| 161 |
infile2<-"list_365_dates_04212012.txt" #list of dates |
|
| 162 |
infile3<-"LST_dates_var_names.txt" #LST dates name |
|
| 163 |
infile4<-"models_interpolation_05142012.txt" #Interpolation model names |
|
| 164 |
infile5<-"mean_day244_rescaled.rst" #mean LST for day 244 |
|
| 165 |
inlistf<-"list_files_05032012.txt" #list of raster images containing the Covariates |
|
| 166 |
infile6<-"OR83M_state_outline.shp" |
|
| 167 |
#stat_loc<-read.table(paste(path,"/","location_study_area_OR_0602012.txt",sep=""),sep=",", header=TRUE) |
|
| 158 | 168 |
|
| 159 | 169 |
i=2 |
| 160 | 170 |
##### LOAD USEFUL DATA |
| 161 | 171 |
|
| 162 | 172 |
#obj_list<-"list_obj_08262012.txt" #Results of fusion from the run on ATLAS |
| 163 |
path_wd<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison_10242012" #Jupiter LOCATION on Atlas for kriging #Jupiter Location on XANDERS |
|
| 173 |
path<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison_10242012" #Jupiter LOCATION on Atlas for kriging #Jupiter Location on XANDERS |
|
| 174 |
path_wd<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison_10242012" #Jupiter LOCATION on Atlas for kriging |
|
| 164 | 175 |
#path<-"/Users/benoitparmentier/Dropbox/Data/NCEAS/Oregon_covariates/" #Local dropbox folder on Benoit's laptop |
| 165 |
setwd(path_wd)
|
|
| 176 |
setwd(path) |
|
| 166 | 177 |
path_data_cai<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_CAI" #Change to constant |
| 167 | 178 |
path_data_fus<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_GAM" |
| 168 | 179 |
#list files that contain model objects and ratingin-testing information for CAI and Fusion |
| 169 | 180 |
obj_mod_fus_name<-"results_mod_obj__365d_GAM_fusion_const_all_lstd_11022012.RData" |
| 170 | 181 |
obj_mod_cai_name<-"results_mod_obj__365d_GAM_CAI2_const_all_10312012.RData" |
| 171 | 182 |
|
| 183 |
#external function |
|
| 184 |
source("function_methods_comparison_assessment_part7_12102012.R")
|
|
| 172 | 185 |
|
| 173 | 186 |
### Projection for the current region |
| 174 | 187 |
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"; |
| 175 | 188 |
#User defined output prefix |
| 176 | 189 |
|
| 190 |
### MAKE THIS A FUNCTION TO LOAD STACK AND DEFINE VALID RANGE... |
|
| 177 | 191 |
#CRS<-proj4string(ghcn) #Storing projection information (ellipsoid, datum,etc.) |
| 178 | 192 |
lines<-read.table(paste(path,"/",inlistf,sep=""), sep="") #Column 1 contains the names of raster files |
| 179 | 193 |
inlistvar<-lines[,1] |
| ... | ... | |
| 184 | 198 |
layerNames(s_raster)<-covar_names #Assigning names to the raster layers |
| 185 | 199 |
projection(s_raster)<-proj_str |
| 186 | 200 |
|
| 201 |
#Create mask using land cover data |
|
| 202 |
pos<-match("LC10",layerNames(s_raster)) #Find the layer which contains water bodies
|
|
| 203 |
LC10<-subset(s_raster,pos) |
|
| 204 |
LC10[is.na(LC10)]<-0 #Since NA values are 0, we assign all zero to NA |
|
| 205 |
mask_land<-LC10<100 #All values below 100% water are assigned the value 1, value 0 is "water" |
|
| 206 |
mask_land_NA<-mask_land |
|
| 207 |
mask_land_NA[mask_land_NA==0]<-NA #Water bodies are assigned value 1 |
|
| 208 |
|
|
| 209 |
data_name<-"mask_land_OR" |
|
| 210 |
raster_name<-paste(data_name,".rst", sep="") |
|
| 211 |
writeRaster(mask_land, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
| 212 |
#writeRaster(r2, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
| 213 |
|
|
| 214 |
pos<-match("ELEV_SRTM",layerNames(s_raster)) #Find column with name "ELEV_SRTM"
|
|
| 215 |
ELEV_SRTM<-raster(s_raster,layer=pos) #Select layer from stack on 10/30 |
|
| 216 |
s_raster<-dropLayer(s_raster,pos) |
|
| 217 |
ELEV_SRTM[ELEV_SRTM <0]<-NA |
|
| 218 |
mask_ELEV_SRTM<-ELEV_SRTM>0 |
|
| 219 |
|
|
| 220 |
#Change this a in loop... |
|
| 221 |
pos<-match("LC1",layerNames(s_raster)) #Find column with name "value"
|
|
| 222 |
LC1<-raster(s_raster,layer=pos) #Select layer from stack |
|
| 223 |
s_raster<-dropLayer(s_raster,pos) |
|
| 224 |
LC1[is.na(LC1)]<-0 |
|
| 225 |
pos<-match("LC2",layerNames(s_raster)) #Find column with name "value"
|
|
| 226 |
LC2<-raster(s_raster,layer=pos) #Select layer from stack |
|
| 227 |
s_raster<-dropLayer(s_raster,pos) |
|
| 228 |
LC2[is.na(LC2)]<-0 |
|
| 229 |
pos<-match("LC3",layerNames(s_raster)) #Find column with name "value"
|
|
| 230 |
LC3<-raster(s_raster,layer=pos) #Select layer from stack |
|
| 231 |
s_raster<-dropLayer(s_raster,pos) |
|
| 232 |
LC3[is.na(LC3)]<-0 |
|
| 233 |
pos<-match("LC4",layerNames(s_raster)) #Find column with name "value"
|
|
| 234 |
LC4<-raster(s_raster,layer=pos) #Select layer from stack |
|
| 235 |
s_raster<-dropLayer(s_raster,pos) |
|
| 236 |
LC4[is.na(LC4)]<-0 |
|
| 237 |
pos<-match("LC6",layerNames(s_raster)) #Find column with name "value"
|
|
| 238 |
LC6<-raster(s_raster,layer=pos) #Select layer from stack |
|
| 239 |
s_raster<-dropLayer(s_raster,pos) |
|
| 240 |
LC6[is.na(LC6)]<-0 |
|
| 241 |
pos<-match("LC7",layerNames(s_raster)) #Find column with name "value"
|
|
| 242 |
LC7<-raster(s_raster,layer=pos) #Select layer from stack |
|
| 243 |
s_raster<-dropLayer(s_raster,pos) |
|
| 244 |
LC7[is.na(LC7)]<-0 |
|
| 245 |
pos<-match("LC9",layerNames(s_raster)) #Find column with name "LC9", this is wetland...
|
|
| 246 |
LC9<-raster(s_raster,layer=pos) #Select layer from stack |
|
| 247 |
s_raster<-dropLayer(s_raster,pos) |
|
| 248 |
LC9[is.na(LC9)]<-0 |
|
| 249 |
|
|
| 250 |
LC_s<-stack(LC1,LC2,LC3,LC4,LC6,LC7) |
|
| 251 |
layerNames(LC_s)<-c("LC1_forest","LC2_shrub","LC3_grass","LC4_crop","LC6_urban","LC7_barren")
|
|
| 252 |
LC_s <-mask(LC_s,mask_ELEV_SRTM) |
|
| 253 |
plot(LC_s) |
|
| 254 |
|
|
| 255 |
s_raster<-addLayer(s_raster, LC_s) |
|
| 256 |
|
|
| 257 |
#mention this is the last... files |
|
| 258 |
|
|
| 259 |
#Read region outline... |
|
| 260 |
filename<-sub(".shp","",infile6) #Removing the extension from file.
|
|
| 261 |
reg_outline<-readOGR(".", filename) #reading shapefile
|
|
| 262 |
|
|
| 187 | 263 |
########## Load Snotel data |
| 188 | 264 |
infile_snotname<-"snot_OR_2010_sp2_methods_11012012_.shp" #load Snotel data |
| 189 | 265 |
snot_OR_2010_sp<-readOGR(".",sub(".shp","",infile_snotname))
|
| ... | ... | |
| 206 | 282 |
#Load GHCN data used in modeling: training and validation site |
| 207 | 283 |
|
| 208 | 284 |
### load specific date...and plot: make a function to extract the diff and prediction... |
| 209 |
rast_diff_fc<-rast_fus_pred-rast_cai_pred |
|
| 210 |
layerNames(rast_diff)<-paste("diff",date_selected,sep="_")
|
|
| 285 |
#rast_diff_fc<-rast_fus_pred-rast_cai_pred
|
|
| 286 |
#layerNames(rast_diff)<-paste("diff",date_selected,sep="_")
|
|
| 211 | 287 |
|
| 212 | 288 |
####COMPARE WITH LOCATION OF GHCN and SNOTEL NETWORK |
| 213 | 289 |
|
| ... | ... | |
| 215 | 291 |
i=1 |
| 216 | 292 |
date_selected<-dates[i] |
| 217 | 293 |
|
| 218 |
X11(width=16,height=9) |
|
| 219 |
par(mfrow=c(1,2)) |
|
| 294 |
X11(12,12) |
|
| 295 |
# #plot(rast_diff_fc) |
|
| 296 |
# plot(snot_OR_2010_sp,pch=2,col="red",add=T) |
|
| 297 |
# plot(data_stat,add=T) #This is the GHCN network |
|
| 298 |
# legend("bottom",legend=c("SNOTEL", "GHCN"),
|
|
| 299 |
# cex=0.8, col=c("red","black"),
|
|
| 300 |
# pch=c(2,1)) |
|
| 301 |
# title(paste("SNOTEL and GHCN networks on ", date_selected, sep=""))
|
|
| 220 | 302 |
|
| 221 |
plot(rast_diff_fc)
|
|
| 303 |
plot(ELEV_SRTM)
|
|
| 222 | 304 |
plot(snot_OR_2010_sp,pch=2,col="red",add=T) |
| 223 |
plot(data_stat,add=T) #This is the GHCN network
|
|
| 305 |
#plot(data_stat,add=T)
|
|
| 224 | 306 |
legend("bottom",legend=c("SNOTEL", "GHCN"),
|
| 225 | 307 |
cex=0.8, col=c("red","black"),
|
| 226 | 308 |
pch=c(2,1)) |
| 227 |
title(paste("SNOTEL and GHCN networks on ", date_selected, sep=""))
|
|
| 228 |
|
|
| 229 |
plot(ELEV_SRTM) |
|
| 230 |
plot(snot_OR_2010_sp,pch=2,col="red",add=T) |
|
| 231 |
plot(data_stat,add=T) |
|
| 232 |
legend("bottom",legend=c("SNOTEL", "GHCN"),
|
|
| 233 |
cex=0.8, col=c("red","black"),
|
|
| 234 |
pch=c(2,1)) |
|
| 235 | 309 |
title(paste("SNOTEL and GHCN networks", sep=""))
|
| 236 | 310 |
savePlot(paste("fig1_map_SNOT_GHCN_network_diff_elev_bckgd",date_selected,out_prefix,".png", sep=""), type="png")
|
| 311 |
dev.off() |
|
| 237 | 312 |
|
| 238 | 313 |
#add histogram of elev for SNOT and GHCN |
| 239 |
hist(snot_data_selected$ELEV_SRTM,main="") |
|
| 240 |
title(paste("SNOT stations and Elevation",date_selected,sep=" "))
|
|
| 241 |
hist(data_vc$ELEV_SRTM,main="") |
|
| 242 |
title(paste("GHCN stations and Elevation",date_selected,sep=" "))
|
|
| 243 |
savePlot(paste("fig2_hist_elev_SNOT_GHCN_",out_prefix,".png", sep=""), type="png")
|
|
| 244 |
dev.off() |
|
| 314 |
#X11(width=16,height=9) |
|
| 315 |
#par(mfrow=c(1,2)) |
|
| 316 |
#hist(snot_data_selected$ELEV_SRTM,main="") |
|
| 317 |
#title(paste("SNOT stations and Elevation",date_selected,sep=" "))
|
|
| 318 |
#hist(data_vc$ELEV_SRTM,main="") |
|
| 319 |
#title(paste("GHCN stations and Elevation",date_selected,sep=" "))
|
|
| 320 |
#savePlot(paste("fig2_hist_elev_SNOT_GHCN_",out_prefix,".png", sep=""), type="png")
|
|
| 321 |
#dev.off() |
|
| 245 | 322 |
## Select date from SNOT |
| 246 | 323 |
#not_selected<-subset(snot_OR_2010_sp, date=="90110" ) |
| 247 | 324 |
list_ac_tab <-vector("list", length(dates)) #storing the accuracy metric data.frame in a list...
|
| 248 | 325 |
names(list_ac_tab)<-paste("date",1:length(dates),sep="")
|
| 249 |
X11(width=16,height=9) |
|
| 250 |
par(mfrow=c(1,2)) |
|
| 251 |
#for(i in 1:length(dates)){
|
|
| 252 |
for(i in 163:length(dates)){
|
|
| 253 |
date_selected<-dates[i] |
|
| 254 |
|
|
| 255 |
## Get the relevant raster layers with prediction for fusion and CAI |
|
| 256 |
oldpath<-getwd() |
|
| 257 |
setwd(path_data_cai) |
|
| 258 |
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_CAI2_const_all_10312012.rst",sep="")) #Search for files in relation to fusion
|
|
| 259 |
lf_cai2c<-list.files(pattern=file_pat) #Search for files in relation to fusion |
|
| 260 |
rast_cai2c<-stack(lf_cai2c) #lf_cai2c CAI results with constant sampling over 365 dates |
|
| 261 |
rast_cai2c<-mask(rast_cai2c,mask_ELEV_SRTM) |
|
| 262 |
|
|
| 263 |
oldpath<-getwd() |
|
| 264 |
setwd(path_data_fus) |
|
| 265 |
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_fusion_const_all_lstd_11022012.rst",sep="")) #Search for files in relation to fusion
|
|
| 266 |
lf_fus1c<-list.files(pattern=file_pat) #Search for files in relation to fusion |
|
| 267 |
rast_fus1c<-stack(lf_fus1c) |
|
| 268 |
rast_fus1c<-mask(rast_fus1c,mask_ELEV_SRTM) |
|
| 269 |
|
|
| 270 |
#PLOT ALL MODELS |
|
| 271 |
#Prepare for plotting |
|
| 272 |
|
|
| 273 |
setwd(path) #set path to the output path |
|
| 274 |
|
|
| 275 |
rast_fus_pred<-raster(rast_fus1c,1) # Select the first model from the stack i.e fusion with kriging for both steps |
|
| 276 |
rast_cai_pred<-raster(rast_cai2c,1) |
|
| 277 |
layerNames(rast_cai_pred)<-paste("cai",date_selected,sep="_")
|
|
| 278 |
layerNames(rast_fus_pred)<-paste("fus",date_selected,sep="_")
|
|
| 279 |
rast_pred2<-stack(rast_fus_pred,rast_cai_pred) |
|
| 280 |
#function to extract training and test from object from object models created earlier during interpolation... |
|
| 281 |
|
|
| 282 |
#load training and testing date for the specified date for fusion and CAI |
|
| 283 |
data_vf<-station_data_interp(date_selected,file.path(path_data_fus,obj_mod_fus_name),training=FALSE,testing=TRUE) |
|
| 284 |
#data_sf<-station_data_interp(date_selected,file.path(path_data_fus,obj_mod_fus_name),training=TRUE,testing=FALSE) |
|
| 285 |
data_vc<-station_data_interp(date_selected,file.path(path_data_cai,obj_mod_cai_name),training=FALSE,testing=TRUE) |
|
| 286 |
#data_sc<-station_data_interp(date_selected,file.path(path_data_cai,obj_mod_cai_name),training=TRUE,testing=FALSE) |
|
| 287 |
|
|
| 288 |
date_selected_snot<-strptime(date_selected,"%Y%m%d") |
|
| 289 |
snot_selected <-snot_OR_2010_sp[snot_OR_2010_sp$date_formatted==date_selected_snot,] |
|
| 290 |
#snot_selected<-na.omit(as.data.frame(snot_OR_2010_sp[snot_OR_2010_sp$date==90110,])) |
|
| 291 |
rast_diff_fc<-rast_fus_pred-rast_cai_pred |
|
| 292 |
LC_stack<-stack(LC1,LC2,LC3,LC4,LC6,LC7) |
|
| 293 |
rast_pred3<-stack(rast_diff_fc,rast_pred2,ELEV_SRTM,LC_stack) |
|
| 294 |
layerNames(rast_pred3)<-c("diff_fc","fus","CAI","ELEV_SRTM","LC1","LC2","LC3","LC4","LC6","LC7") #extract amount of veg...
|
|
| 295 |
|
|
| 296 |
#extract predicted tmax corresponding to |
|
| 297 |
extract_snot<-extract(rast_pred3,snot_selected) #return value from extract is a matrix (with input SPDF) |
|
| 298 |
snot_data_selected<-cbind(as.data.frame(snot_selected),extract_snot) #bind data together |
|
| 299 |
snot_data_selected$res_f<-snot_data_selected$fus-snot_data_selected$tmax #calculate the residuals for Fusion |
|
| 300 |
snot_data_selected$res_c<-snot_data_selected$CAI-snot_data_selected$tmax #calculate the residuals for CAI |
|
| 301 |
#snot_data_selected<-(na.omit(as.data.frame(snot_data_selected))) #remove rows containing NA, this may need to be modified later. |
|
| 302 |
|
|
| 303 |
###fig3: Plot predicted vs observed tmax |
|
| 304 |
#fig3a: FUS |
|
| 305 |
x_range<-range(c(data_vf$pred_mod7,snot_data_selected$fus,data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T) |
|
| 306 |
y_range<-range(c(data_vf$dailyTmax,snot_data_selected$tmax,data_vc$dailyTmax,snot_data_selected$tmax),na.rm=T) |
|
| 307 |
plot(data_vf$pred_mod7,data_vf$dailyTmax, ylab="Observed daily tmax (C)", xlab="Fusion predicted daily tmax (C)", |
|
| 308 |
ylim=y_range,xlim=x_range) |
|
| 309 |
#text(data_vf$pred_mod7,data_vf$dailyTmax,labels=data_vf$idx,pos=3) |
|
| 310 |
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
| 311 |
grid(lwd=0.5,col="black") |
|
| 312 |
points(snot_data_selected$fus,snot_data_selected$tmax,pch=2,co="red") |
|
| 313 |
title(paste("Testing stations tmax fusion vs daily tmax",date_selected,sep=" "))
|
|
| 314 |
legend("topleft",legend=c("GHCN", "SNOT"),
|
|
| 315 |
cex=1.2, col=c("black","red"),
|
|
| 316 |
pch=c(1,2)) |
|
| 317 |
#fig 3b: CAI |
|
| 318 |
#x_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI)) |
|
| 319 |
#y_range<-range(c(data_vc$dailyTmax,snot_data_selected$tmax)) |
|
| 320 |
plot(data_vc$pred_mod9,data_vc$dailyTmax, ylab="Observed daily tmax (C)", xlab="CAI predicted daily tmax (C)", |
|
| 321 |
ylim=y_range,xlim=x_range) |
|
| 322 |
#text(data_vc$pred_mod9,data_vc$dailyTmax,labels=data_vf$idx,pos=3) |
|
| 323 |
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
| 324 |
grid(lwd=0.5,col="black") |
|
| 325 |
points(snot_data_selected$CAI,snot_data_selected$tmax,pch=2,co="red") |
|
| 326 |
#text(snot_data_selected$CAI,snot_data_selected$tmax,labels=1:nrow(snot_data_selected),pos=3) |
|
| 327 |
#title(paste("Testing stations tmax CAI vs daily tmax",date_selected,sep=" "))
|
|
| 328 |
legend("topleft",legend=c("GHCN", "SNOT"),
|
|
| 329 |
cex=1.2, col=c("black","red"),
|
|
| 330 |
pch=c(1,2)) |
|
| 331 |
savePlot(paste("fig3_testing_scatterplot_pred_fus_CAI_observed_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 332 |
|
|
| 333 |
##### Fig4a: ELEV-CAI |
|
| 334 |
y_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T) |
|
| 335 |
#y_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T) |
|
| 336 |
x_range<-range(c(data_vc$ELEV_SRTM,snot_data_selected$ELEV_SRTM),na.rm=T) |
|
| 337 |
lm_mod1<-lm(data_vc$pred_mod9~data_vc$ELEV_SRTM) |
|
| 338 |
lm_mod2<-lm(snot_data_selected$CAI~snot_data_selected$ELEV_SRTM) |
|
| 339 |
plot(data_vc$ELEV_SRTM,data_vc$pred_mod9,ylab="Observed daily tmax (C)", xlab="Elevation (m)", |
|
| 340 |
ylim=y_range,xlim=x_range) |
|
| 341 |
#text(data_vc$ELEV_SRTM,data_vc$pred_mod9,labels=data_vc$idx,pos=3) |
|
| 342 |
abline(lm_mod1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
| 343 |
abline(lm_mod2,col="red") #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
| 344 |
grid(lwd=0.5, col="black") |
|
| 345 |
points(snot_data_selected$ELEV_SRTM,snot_data_selected$CAI,pch=2,co="red") |
|
| 346 |
title(paste("Testing stations tmax CAI vs elevation",date_selected,sep=" "))
|
|
| 347 |
legend("topleft",legend=c("GHCN", "SNOT"),
|
|
| 348 |
cex=1.2, col=c("black","red"),
|
|
| 349 |
pch=c(1,2)) |
|
| 350 |
|
|
| 351 |
#Fig4bELEV-FUS |
|
| 352 |
y_range<-range(c(data_vf$pred_mod7,snot_data_selected$fus),na.rm=T) |
|
| 353 |
x_range<-range(c(data_vf$ELEV_SRTM,snot_data_selected$ELEV_SRTM),na.rm=T) |
|
| 354 |
lm_mod1<-lm(data_vf$pred_mod7~data_vf$ELEV_SRTM) |
|
| 355 |
lm_mod2<-lm(snot_data_selected$fus~snot_data_selected$ELEV_SRTM) |
|
| 356 |
plot(data_vf$ELEV_SRTM,data_vf$pred_mod7,ylab="Observed daily tmax (C)", xlab="Elevation (m)", |
|
| 357 |
ylim=y_range,xlim=x_range) |
|
| 358 |
#text(data_vc$ELEV_SRTM,data_vc$pred_mod9,labels=data_vc$idx,pos=3) |
|
| 359 |
abline(lm_mod1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
| 360 |
abline(lm_mod2,col="red") #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
| 361 |
grid(lwd=0.5, col="black") |
|
| 362 |
points(snot_data_selected$ELEV_SRTM,snot_data_selected$fus,pch=2,co="red") |
|
| 363 |
title(paste("Testing stations tmax vs elevation",date_selected,sep=" "))
|
|
| 364 |
legend("topleft",legend=c("GHCN", "SNOT"),
|
|
| 365 |
cex=1.2, col=c("black","red"),
|
|
| 366 |
pch=c(1,2)) |
|
| 367 |
savePlot(paste("fig4_testing_scatterplot_pred_fus_CIA_elev_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 368 |
|
|
| 369 |
############ ACCURACY METRICS AND RESIDUALS ############# |
|
| 370 |
|
|
| 371 |
#START FIG 5 |
|
| 372 |
#####Fig5a: CAI vs FUSION: difference by plotting on in terms of the other |
|
| 373 |
lm_mod<-lm(snot_data_selected$CAI~snot_data_selected$fus) |
|
| 374 |
y_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T) |
|
| 375 |
x_range<-range(c(data_vf$pred_mod7,snot_data_selected$fus),na.rm=T) |
|
| 376 |
|
|
| 377 |
plot(data_vf$pred_mod7,data_vc$pred_mod9,ylab="Predicted CAI daily tmax (C)", xlab="Predicted fusion daily tmax (C)", |
|
| 378 |
ylim=y_range,xlim=x_range) |
|
| 379 |
#text(data_vc$ELEV_SRTM,data_vc$dailyTmax,labels=data_vc$idx,pos=3) |
|
| 380 |
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
| 381 |
abline(lm_mod,col="red") |
|
| 382 |
grid(lwd=0.5, col="black") |
|
| 383 |
points(snot_data_selected$fus,snot_data_selected$CAI,pch=2,co="red") |
|
| 384 |
title(paste("Testing stations predicted tmax fusion vs CAI tmax",date_selected,sep=" "))
|
|
| 385 |
legend("topleft",legend=c("GHCN", "SNOT"),
|
|
| 386 |
cex=1.2, col=c("black","red"),
|
|
| 387 |
pch=c(1,2)) |
|
| 388 |
####Fig5b: diff vs elev: difference by plotting on in terms of elev |
|
| 389 |
diff_fc<-data_vf$pred_mod7-data_vc$pred_mod9 |
|
| 390 |
plot(snot_data_selected$ELEV_SRTM,snot_data_selected$diff_fc,pch=2,col="red") |
|
| 391 |
lm_mod<-lm(snot_data_selected$diff_fc~snot_data_selected$ELEV_SRTM) |
|
| 392 |
abline(lm_mod,col="red") |
|
| 393 |
points(data_vf$ELEV_SRTM,diff_fc) |
|
| 394 |
lm_mod<-lm(diff_fc~data_vf$ELEV_SRTM) |
|
| 395 |
abline(lm_mod) |
|
| 396 |
legend("topleft",legend=c("GHCN", "SNOT"),
|
|
| 397 |
cex=1.2, col=c("black","red"),
|
|
| 398 |
pch=c(1,2)) |
|
| 399 |
title(paste("Prediction tmax difference and elevation ",sep=""))
|
|
| 400 |
savePlot(paste("fig5_testing_scatterplot_pred_fus_CAI_observed_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 401 |
|
|
| 402 |
#DO diff IN TERM OF ELEVATION CLASSES as well as diff.. |
|
| 403 |
|
|
| 404 |
#### START FIG 6: difference fc vs elev |
|
| 405 |
#fig6a |
|
| 406 |
brks<-c(0,500,1000,1500,2000,2500,4000) |
|
| 407 |
lab_brks<-1:6 |
|
| 408 |
elev_rcstat<-cut(snot_data_selected$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F) |
|
| 409 |
snot_data_selected$elev_rec<-elev_rcstat |
|
| 410 |
y_range<-range(c(snot_data_selected$diff_fc),na.rm=T) |
|
| 411 |
x_range<-range(c(elev_rcstat),na.rm=T) |
|
| 412 |
plot(elev_rcstat,snot_data_selected$diff_fc, ylab="diff_fc", xlab="ELEV_SRTM (m) ", |
|
| 413 |
ylim=y_range, xlim=x_range) |
|
| 414 |
#text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3) |
|
| 415 |
grid(lwd=0.5,col="black") |
|
| 416 |
title(paste("SNOT stations diff f vs Elevation",date_selected,sep=" "))
|
|
| 417 |
|
|
| 418 |
###With fewer classes...fig6b |
|
| 419 |
brks<-c(0,1000,2000,3000,4000) |
|
| 420 |
lab_brks<-1:4 |
|
| 421 |
elev_rcstat<-cut(snot_data_selected$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F) |
|
| 422 |
snot_data_selected$elev_rec<-elev_rcstat |
|
| 423 |
y_range<-range(c(snot_data_selected$diff_fc),na.rm=T) |
|
| 424 |
x_range<-range(c(elev_rcstat),na.rm=T) |
|
| 425 |
plot(elev_rcstat,snot_data_selected$diff_fc, ylab="diff_fc", xlab="ELEV_SRTM (m) ", |
|
| 426 |
ylim=y_range, xlim=x_range) |
|
| 427 |
#text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3) |
|
| 428 |
grid(lwd=0.5,col="black") |
|
| 429 |
title(paste("SNOT stations diff f vs Elevation",date_selected,sep=" "))
|
|
| 430 |
savePlot(paste("fig6_elevation_classes_diff_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 431 |
|
|
| 432 |
#START FIG 7 with residuals |
|
| 433 |
#fig 7a |
|
| 434 |
brks<-c(0,1000,2000,3000,4000) |
|
| 435 |
lab_brks<-1:4 |
|
| 436 |
elev_rcstat<-cut(snot_data_selected$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F) |
|
| 437 |
snot_data_selected$elev_rec<-elev_rcstat |
|
| 438 |
y_range<-range(c(snot_data_selected$res_f,snot_data_selected$res_c),na.rm=T) |
|
| 439 |
x_range<-range(c(elev_rcstat),na.rm=T) |
|
| 440 |
plot(elev_rcstat,snot_data_selected$res_f, ylab="res_f", xlab="ELEV_SRTM (m) ", |
|
| 441 |
ylim=y_range, xlim=x_range) |
|
| 442 |
#text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3) |
|
| 443 |
grid(lwd=0.5,col="black") |
|
| 444 |
title(paste("SNOT stations residuals fusion vs Elevation",date_selected,sep=" "))
|
|
| 445 |
#fig 7b |
|
| 446 |
elev_rcstat<-cut(snot_data_selected$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F) |
|
| 447 |
y_range<-range(c(snot_data_selected$res_c,snot_data_selected$res_f),na.rm=T) |
|
| 448 |
x_range<-range(c(elev_rcstat)) |
|
| 449 |
plot(elev_rcstat,snot_data_selected$res_c, ylab="res_c", xlab="ELEV_SRTM (m) ", |
|
| 450 |
ylim=y_range, xlim=x_range) |
|
| 451 |
#text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3) |
|
| 452 |
grid(lwd=0.5,col="black") |
|
| 453 |
title(paste("SNOT stations residuals CAI vs Elevation",date_selected,sep=" "))
|
|
| 454 |
savePlot(paste("fig7_elevation_classes_residuals_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 455 |
|
|
| 456 |
####### COMPARE CAI FUSION USING SNOTEL DATA WITH ACCURACY METRICS############### |
|
| 457 |
################ RESIDUALS and MAE etc. ##################### |
|
| 458 |
|
|
| 459 |
### Run for full list of date? --365 |
|
| 460 |
ac_tab_snot_fus<-calc_accuracy_metrics(snot_data_selected$tmax,snot_data_selected$fus) |
|
| 461 |
ac_tab_snot_cai<-calc_accuracy_metrics(snot_data_selected$tmax,snot_data_selected$CAI) |
|
| 462 |
ac_tab_ghcn_fus<-calc_accuracy_metrics(data_vf$dailyTmax,data_vf$pred_mod7) |
|
| 463 |
ac_tab_ghcn_cai<-calc_accuracy_metrics(data_vc$dailyTmax,data_vc$pred_mod9) |
|
| 464 |
|
|
| 465 |
ac_tab<-do.call(rbind,list(ac_tab_snot_fus,ac_tab_snot_cai,ac_tab_ghcn_fus,ac_tab_ghcn_cai)) |
|
| 466 |
rownames(ac_tab)<-c("snot_fus","snot_cai","ghcn_fus","ghcn_cai")
|
|
| 467 |
ac_tab$date<-date_selected |
|
| 468 |
list_ac_tab[[i]]<-ac_tab #storing the accuracy metric data.frame in a list... |
|
| 469 |
#save(list_ac_tab,) |
|
| 470 |
save(list_ac_tab,file= paste("list_ac_tab_", date_selected,out_prefix,".RData",sep=""))
|
|
| 471 |
|
|
| 472 |
#FIG8: boxplot of residuals for methods (fus, cai) using SNOT and GHCN data |
|
| 473 |
#fig8a |
|
| 474 |
y_range<-range(c(snot_data_selected$res_f,snot_data_selected$res_c,data_vf$res_mod7,data_vc$res_mod9),na.rm=T) |
|
| 475 |
boxplot(snot_data_selected$res_f,snot_data_selected$res_c,names=c("FUS","CAI"),ylim=y_range,ylab="Residuals tmax degree C")
|
|
| 476 |
title(paste("Residuals for fusion and CAI methods for SNOT data ",date_selected,sep=" "))
|
|
| 477 |
#fig8b |
|
| 478 |
boxplot(data_vf$res_mod7,data_vc$res_mod9,names=c("FUS","CAI"),ylim=y_range,ylab="Residuals tmax degree C")
|
|
| 479 |
title(paste("Residuals for fusion and CAI methods for GHCN data ",date_selected,sep=" "))
|
|
| 480 |
savePlot(paste("fig8_residuals_boxplot_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 481 |
|
|
| 482 |
mae_fun<-function(residuals){
|
|
| 483 |
mean(abs(residuals),na.rm=T) |
|
| 484 |
} |
|
| 485 |
|
|
| 486 |
mean_diff_fc<-aggregate(diff_fc~elev_rec,data=snot_data_selected,mean) |
|
| 487 |
mean_mae_c<-aggregate(res_c~elev_rec,data=snot_data_selected,mae_fun) |
|
| 488 |
mean_mae_f<-aggregate(res_f~elev_rec,data=snot_data_selected,mae_fun) |
|
| 489 |
|
|
| 490 |
####FIG 9: plot MAE for fusion and CAI as well as boxplots of both thechnique |
|
| 491 |
#fig 9a: boxplot of residuals for MAE and CAI |
|
| 492 |
height<-cbind(snot_data_selected$res_f,snot_data_selected$res_c) |
|
| 493 |
boxplot(height,names=c("FUS","CAI"),ylab="Residuals tmax degree C")
|
|
| 494 |
title(paste("Residuals for fusion and CAI methods for SNOT data ",date_selected,sep=" "))
|
|
| 495 |
#par(new=TRUE) |
|
| 496 |
#abline(h=ac_tab[1,1],col="red") |
|
| 497 |
points(1,ac_tab[1,1],pch=5,col="red") |
|
| 498 |
points(2,ac_tab[2,1],pch=5,col="black") |
|
| 499 |
legend("bottom",legend=c("FUS_MAE", "CAI_MAE"),
|
|
| 500 |
cex=0.8, col=c("red","black"),
|
|
| 501 |
pch=c(2,1)) |
|
| 502 |
#fig 9b: MAE per 3 elevation classes:0-1000,1000-2000,2000-3000,3000-4000 |
|
| 503 |
y_range<-c(0,max(c(mean_mae_c[,2],mean_mae_f[,2]),na.rm=T)) |
|
| 504 |
plot(1:3,mean_mae_c[,2],ylim=y_range,type="n",ylab="MAE in degree C",xlab="elevation classes") |
|
| 505 |
points(mean_mae_c,ylim=y_range) |
|
| 506 |
lines(1:3,mean_mae_c[,2],col="black") |
|
| 507 |
par(new=TRUE) # key: ask for new plot without erasing old |
|
| 508 |
points(mean_mae_f,ylim=y_range) |
|
| 509 |
lines(1:3,mean_mae_f[,2],col="red") |
|
| 510 |
legend("bottom",legend=c("FUS_MAE", "CAI_MAE"),
|
|
| 511 |
cex=0.8, col=c("red","black"),
|
|
| 512 |
pch=c(2,1)) |
|
| 513 |
title(paste("MAE per elevation classes for SNOT data ",date_selected,sep=" "))
|
|
| 514 |
savePlot(paste("fig9_residuals_boxplot_MAE_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 515 |
|
|
| 516 |
### LM MODELS for difference and elevation categories |
|
| 517 |
## Are the differences plotted on fig 9 significant?? |
|
| 518 |
diffelev_mod<-lm(diff_fc~elev_rec,data=snot_data_selected) |
|
| 519 |
summary(diffelev_mod) |
|
| 520 |
##LM MODEL MAE PER ELEVATION CLASS: residuals for CAI |
|
| 521 |
diffelev_mod<-lm(res_c~elev_rec,data=snot_data_selected) |
|
| 522 |
summary(diffelev_mod) |
|
| 523 |
##LM MODEL MAE PER ELEVATION CLASS: residuals for Fusions |
|
| 524 |
diffelev_mod<-lm(res_f~elev_rec,data=snot_data_selected) |
|
| 525 |
summary(diffelev_mod) |
|
| 526 |
|
|
| 527 |
### LM MODELS for RESIDUALS BETWEEN CAI AND FUSION |
|
| 528 |
## Are the differences plotted on fig 9 significant?? |
|
| 529 |
## STORE THE p values...?? overall and per cat? |
|
| 530 |
|
|
| 531 |
#diffelev_mod<-lm(res_f~elev_rec,data=snot_data_selected) |
|
| 532 |
#table(snot_data_selected$elev_rec) #Number of observation per class |
|
| 533 |
#max(snot_data_selected$E_STRM) |
|
| 534 |
|
|
| 535 |
#res |
|
| 536 |
|
|
| 537 |
############################################# |
|
| 538 |
#USING BOTH validation and training |
|
| 539 |
#This part is exploratory.... |
|
| 540 |
################## EXAMINING RESIDUALS AND DIFFERENCES IN LAND COVER......############ |
|
| 541 |
###### |
|
| 542 |
|
|
| 543 |
#LC_names<-c("LC1_rec","LC2_rec","LC3_rec","LC4_rec","LC6_rec")
|
|
| 544 |
suf_name<-c("rec1")
|
|
| 545 |
sum_var<-c("diff_fc")
|
|
| 546 |
LC_names<-c("LC1","LC2","LC3","LC4","LC6")
|
|
| 547 |
brks<-c(-1,20,40,60,80,101) |
|
| 548 |
lab_brks<-seq(1,5,1) |
|
| 549 |
#var_name<-LC_names; suffix<-"rec1"; s_function<-"mean";df<-snot_data_selected;summary_var<-"diff_fc" |
|
| 550 |
#reclassify_df(snot_data_selected,LC_names,var_name,brks,lab_brks,suffix,summary_var) |
|
| 551 |
|
|
| 552 |
#Calculate mean per land cover percentage |
|
| 553 |
data_agg<-reclassify_df(snot_data_selected,LC_names,brks,lab_brks,suf_name,sum_var) |
|
| 554 |
data_lc<-data_agg[[1]] |
|
| 555 |
snot_data_selected<-data_agg[[2]] |
|
| 556 |
|
|
| 557 |
by_name<-"rec1" |
|
| 558 |
df_lc_diff_fc<-merge_multiple_df(data_lc,by_name) |
|
| 559 |
|
|
| 560 |
###### FIG10: PLOT LAND COVER |
|
| 561 |
zones_stat<-df_lc_diff_fc #first land cover |
|
| 562 |
#names(zones_stat)<-c("lab_brks","LC")
|
|
| 563 |
y_range<-range(as.vector(t(zones_stat[,-1])),na.rm=T) |
|
| 564 |
lab_brks_mid<-c(10,30,50,70,90) |
|
| 565 |
plot(lab_brks_mid,zones_stat[,2],type="b",ylim=y_range,col="black", lwd=2, |
|
| 566 |
ylab="difference between fusion and CAI",xlab="land cover percent classes") |
|
| 567 |
lines(lab_brks_mid,zones_stat[,3],col="red",type="b") |
|
| 568 |
lines(lab_brks_mid,zones_stat[,4],col="blue",type="b") |
|
| 569 |
lines(lab_brks_mid,zones_stat[,5],col="darkgreen",type="b") |
|
| 570 |
lines(lab_brks_mid,zones_stat[,6],col="purple",type="b") |
|
| 571 |
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
|
| 572 |
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
|
| 573 |
lty=1,lwd=1.8) |
|
| 574 |
title(paste("Prediction tmax difference and land cover ",date_selected,sep=""))
|
|
| 575 |
|
|
| 576 |
###NOW USE RESIDUALS FOR FUSION |
|
| 577 |
sum_var<-"res_f" |
|
| 578 |
suf_name<-"rec2" |
|
| 579 |
data_agg2<-reclassify_df(snot_data_selected,LC_names,brks,lab_brks,suf_name,sum_var) |
|
| 580 |
data_resf_lc<-data_agg2[[1]] |
|
| 581 |
#snot_data_selected<-data_agg[[2]] |
|
| 582 |
|
|
| 583 |
by_name<-"rec2" |
|
| 584 |
df_lc_resf<-merge_multiple_df(data_resf_lc,by_name) |
|
| 585 |
|
|
| 586 |
zones_stat<-df_lc_resf #first land cover |
|
| 587 |
#names(zones_stat)<-c("lab_brks","LC")
|
|
| 588 |
lab_brks_mid<-c(10,30,50,70,90) |
|
| 589 |
plot(lab_brks_mid,zones_stat[,2],type="b",ylim=y_range,col="black",lwd=2, |
|
| 590 |
ylab="tmax residuals fusion ",xlab="land cover percent classes") |
|
| 591 |
lines(lab_brks_mid,zones_stat[,3],col="red",type="b") |
|
| 592 |
lines(lab_brks_mid,zones_stat[,4],col="blue",type="b") |
|
| 593 |
lines(lab_brks_mid,zones_stat[,5],col="darkgreen",type="b") |
|
| 594 |
lines(lab_brks_mid,zones_stat[,6],col="purple",type="b") |
|
| 595 |
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
|
| 596 |
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
|
| 597 |
lty=1,lwd=1.2) |
|
| 598 |
title(paste("Prediction tmax residuals and land cover ",date_selected,sep=""))
|
|
| 599 |
savePlot(paste("fig10_diff_prediction_tmax_diff_res_f_land cover",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 600 |
|
|
| 601 |
#### FIGURE11: res_f and res_c per land cover |
|
| 602 |
|
|
| 603 |
sum_var<-"res_c" |
|
| 604 |
suf_name<-"rec3" |
|
| 605 |
data_agg3<-reclassify_df(snot_data_selected,LC_names,brks,lab_brks,suf_name,sum_var) |
|
| 606 |
data_resc_lc<-data_agg3[[1]] |
|
| 607 |
snot_data_selected<-data_agg3[[2]] |
|
| 608 |
|
|
| 609 |
by_name<-"rec3" |
|
| 610 |
df_lc_resc<-merge_multiple_df(data_resc_lc,by_name) |
|
| 611 |
|
|
| 612 |
zones_stat<-df_lc_resc #first land cover |
|
| 613 |
#names(zones_stat)<-c("lab_brks","LC")
|
|
| 614 |
y_range<-range(as.vector(t(zones_stat[,-1])),na.rm=T) |
|
| 615 |
lab_brks_mid<-c(10,30,50,70,90) |
|
| 616 |
plot(lab_brks_mid,zones_stat[,2],type="b",ylim=y_range,col="black",lwd=2, |
|
| 617 |
ylab="tmax residuals CAI",xlab="land cover percent classes") |
|
| 618 |
lines(lab_brks_mid,zones_stat[,3],col="red",type="b") |
|
| 619 |
lines(lab_brks_mid,zones_stat[,4],col="blue",type="b") |
|
| 620 |
lines(lab_brks_mid,zones_stat[,5],col="darkgreen",type="b") |
|
| 621 |
lines(lab_brks_mid,zones_stat[,6],col="purple",type="b") |
|
| 622 |
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
|
| 623 |
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
|
| 624 |
lty=1,lwd=1.2) |
|
| 625 |
title(paste("Prediction tmax residuals CAI and land cover ",date_selected,sep=""))
|
|
| 626 |
|
|
| 627 |
#fig11b |
|
| 628 |
zones_stat<-df_lc_resf #first land cover |
|
| 629 |
#names(zones_stat)<-c("lab_brks","LC")
|
|
| 630 |
y_range<-range(as.vector(t(zones_stat[,-1])),na.rm=T) |
|
| 631 |
lab_brks_mid<-c(10,30,50,70,90) |
|
| 632 |
plot(lab_brks_mid,zones_stat[,2],type="b",ylim=y_range,col="black",lwd=2, |
|
| 633 |
ylab="tmax residuals fusion ",xlab="land cover percent classes") |
|
| 634 |
lines(lab_brks_mid,zones_stat[,3],col="red",type="b") |
|
| 635 |
lines(lab_brks_mid,zones_stat[,4],col="blue",type="b") |
|
| 636 |
lines(lab_brks_mid,zones_stat[,5],col="darkgreen",type="b") |
|
| 637 |
lines(lab_brks_mid,zones_stat[,6],col="purple",type="b") |
|
| 638 |
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
|
| 639 |
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
|
| 640 |
lty=1,lwd=1.2) |
|
| 641 |
title(paste("Prediction tmax residuals and land cover ",date_selected,sep=""))
|
|
| 642 |
#savePlot(paste("fig10_diff_prediction_tmax_diff_res_f_land cover",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 643 |
savePlot(paste("fig11_prediction_tmax_res_f_res_c_land cover",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 644 |
|
|
| 645 |
} |
|
| 646 | 326 |
|
| 647 |
#Collect accuracy information for different dates |
|
| 648 |
ac_data_xdates<-do.call(rbind,list_ac_tab) |
|
| 649 | 327 |
|
| 650 |
ac_data_xdates$mod_id<-rownames(ac_data_xdates) |
|
| 328 |
#ac_mod<-mclapply(1:length(dates), accuracy_comp_CAI_fus_function,mc.preschedule=FALSE,mc.cores = 8) #This is the end bracket from mclapply(...) statement |
|
| 329 |
source("function_methods_comparison_assessment_part7_12102012.R")
|
|
| 330 |
#Use mcMap or mappply for function with multiple arguments... |
|
| 331 |
#ac_mod<-mclapply(1:6, accuracy_comp_CAI_fus_function,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement |
|
| 332 |
ac_mod<-mclapply(1:length(dates), accuracy_comp_CAI_fus_function,mc.preschedule=FALSE,mc.cores = 8) #This is the end bracket from mclapply(...) statement |
|
| 651 | 333 |
|
| 652 |
tmp_rownames<-rownames(ac_data_xdates) |
|
| 653 |
rowstr<-strsplit(tmp_rownames,"\\.") |
|
| 654 |
for (i in 1:length(rowstr)){
|
|
| 655 |
ac_data_xdates$mod_id[i]<-rowstr[[i]][[2]] |
|
| 334 |
tb<-ac_mod[[1]][[4]][0,] #empty data frame with metric table structure that can be used in rbinding... |
|
| 335 |
tb_tmp<-ac_mod #copy |
|
| 336 |
|
|
| 337 |
for (i in 1:length(tb_tmp)){
|
|
| 338 |
tmp<-tb_tmp[[i]][[4]] |
|
| 339 |
tb<-rbind(tb,tmp) |
|
| 656 | 340 |
} |
| 341 |
rm(tb_tmp) |
|
| 342 |
#Collect accuracy information for different dates |
|
| 343 |
#ac_data_xdates<-do.call(rbind,tb) |
|
| 344 |
ac_data_xdates<-tb |
|
| 657 | 345 |
##Now subset for each model... |
| 658 | 346 |
|
| 659 | 347 |
mod_names<-unique(ac_data_xdates$mod_id) |
| ... | ... | |
| 684 | 372 |
mean(data_ac_ghcn_fus) |
| 685 | 373 |
mean(data_ac_ghcn_cai) |
| 686 | 374 |
|
| 375 |
### END OF CODE |
|
| 687 | 376 |
### END OF CODE |
| 688 | 377 |
#Write a part to caculate MAE per date... |
| 689 | 378 |
#ac_table_metrics<-do.call(rbind,ac_tab_list) |
Also available in: Unified diff
Methods comp part7-task#491- SNOT-GHCN data updated script debugged to run through mcapply