Revision 145bb8a6
Added by Benoit Parmentier about 11 years ago
| climate/research/oregon/interpolation/contribution_of_covariates_paper_interpolation.R | ||
|---|---|---|
| 37 | 37 |
|
| 38 | 38 |
#### FUNCTION USED IN SCRIPT |
| 39 | 39 |
|
| 40 |
function_analyses_paper <-"contribution_of_covariates_paper_interpolation_functions_09092013.R"
|
|
| 40 |
function_analyses_paper <-"contribution_of_covariates_paper_interpolation_functions_09232013.R"
|
|
| 41 | 41 |
|
| 42 | 42 |
############################## |
| 43 | 43 |
#### Parameters and constants |
| ... | ... | |
| 63 | 63 |
met_stations_outfiles_obj_file<-"/data/project/layers/commons/data_workflow/output_data_365d_gam_fus_lst_test_run_07172013/met_stations_outfiles_obj_gam_fusion__365d_gam_fus_lst_test_run_07172013.RData" |
| 64 | 64 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 65 | 65 |
y_var_name <- "dailyTmax" |
| 66 |
out_prefix<-"analyses_09232013"
|
|
| 66 |
out_prefix<-"analyses_10102013"
|
|
| 67 | 67 |
|
| 68 | 68 |
#method_interpolation <- "gam_daily" |
| 69 | 69 |
covar_obj_file_1 <- "covar_obj__365d_gam_day_lst_comb3_08132013.RData" |
| ... | ... | |
| 250 | 250 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 251 | 251 |
par(mfrow=c(1,1)) |
| 252 | 252 |
|
| 253 |
|
|
| 254 |
#plot(elev_WGS84) |
|
| 255 |
plot(interp_area_WGS84) |
|
| 253 |
#plot(elev_WGS84,cex.axis=2.1,cex.z=2.1) |
|
| 254 |
plot(elev_WGS84,cex.axis=1.4,axis.args=list(at=seq(0,3500,500), |
|
| 255 |
labels=seq(0, 3500, 500), |
|
| 256 |
cex.axis=1.4)) |
|
| 257 |
plot(interp_area_WGS84,add=T) |
|
| 256 | 258 |
plot(ghcn_dat_WGS84,add=T) |
| 257 |
title("SStudy area with ")
|
|
| 258 |
#this work on non sp plot too: Scale bar postion |
|
| 259 |
#scale_position<-c(450000, 600000) |
|
| 260 |
#arrow_position<-c(900000, 600000) |
|
| 261 |
|
|
| 262 |
#legend("topright",legend=c(0:7),title="Number of change",
|
|
| 263 |
# pt.cex=1.4,cex=2.1,fill=rev(terrain.colors(8)),bty="n") |
|
| 264 |
#label_scalebar<-c("0","125","250")
|
|
| 265 |
#scalebar(d=250000, xy=scale_position, type = 'bar', |
|
| 266 |
# divs=3,label=label_scalebar,below="kilometers", |
|
| 267 |
# cex=1.8) |
|
| 268 |
|
|
| 269 |
## Northern arrow |
|
| 270 |
#SpatialPolygonsRescale(layout.north.arrow(), offset = arrow_position, |
|
| 271 |
# scale = 150000, fill=c("transparent","black"),plot.grid=FALSE)
|
|
| 272 |
##note that scale in SpatialPolygonRescale sets the size of the north arrow!! |
|
| 273 |
|
|
| 274 |
### PLACE INSET MAP OF THE USA |
|
| 275 |
#opar <- par(fig=c(0.7, 0.95, 0.5, 0.75), new=TRUE) |
|
| 276 |
#opar <- par(fig=c(0, 0, 1, 1), new=TRUE) |
|
| 259 |
title_text <-c("Elevation (m) and meteorological"," stations in Oregon")
|
|
| 260 |
legend("topleft",legend=title_text,cex=2.1,bty="n")
|
|
| 277 | 261 |
|
| 278 | 262 |
par(mar = c(0,0,0,0)) # remove margin |
| 279 | 263 |
#opar <- par(fig=c(0.9,0.95,0.8, 0.85), new=TRUE) |
| 280 |
opar <- par(fig=c(0.85,0.95,0.8, 0.9), new=TRUE) |
|
| 281 |
|
|
| 282 |
#p1<-spplot(ghcn_dat,"station") |
|
| 283 |
#p2<-spplot(usa_map_2,"NAMES_1") |
|
| 284 |
#print(p1,position=c(0,0,1,1),more=T) |
|
| 285 |
#print(p2,position=c(0,0,0.3,0.3),more=T) |
|
| 264 |
#opar <- par(fig=c(0.85,0.95,0.8, 0.9), new=TRUE) |
|
| 265 |
#opar <- par(fig=c(0.8,0.95,0.75, 0.9), new=TRUE) |
|
| 266 |
#opar <- par(fig=c(0.75,0.95,0.7, 0.9), new=TRUE) |
|
| 267 |
opar <- par(fig=c(0.65,0.9,0.8, 0.92), new=TRUE) |
|
| 286 | 268 |
|
| 287 | 269 |
plot(usa_map_2,border="black") #border and lwd are options of graphics package polygon object |
| 288 |
plot(usa_map_OR,col="grey",add=T) |
|
| 270 |
plot(usa_map_OR,col="dark grey",add=T)
|
|
| 289 | 271 |
box() |
| 290 | 272 |
dev.off() |
| 291 | 273 |
|
| 292 |
|
|
| 293 | 274 |
### Figure 2: Method comparison workflow |
| 294 | 275 |
|
| 295 | 276 |
# Workflow figure is not generated in R |
| ... | ... | |
| 615 | 596 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 616 | 597 |
par(mfrow=c(row_mfrow,col_mfrow)) |
| 617 | 598 |
|
| 599 |
col_t<-c("red","blue","black")
|
|
| 600 |
pch_t<- 1:length(col_t) |
|
| 601 |
legend_text <- c("GAM","Kriging","GWR")
|
|
| 602 |
|
|
| 618 | 603 |
y_range<-range(tb1_month$mae,tb3_month$mae,tb4_month$mae) |
| 619 | 604 |
xlab_tick <- unique(tb1$month) |
| 620 | 605 |
xlab_text <-"Month" |
| 621 |
|
|
| 622 |
plot(1:12,tb1_month$mae,col=c("red"),type="b",ylim=y_range,xlab=xlab_text,xaxt="n")
|
|
| 606 |
ylab_text <- "MAE (C)" |
|
| 607 |
plot(1:12,tb1_month$mae,col=c("red"),type="b",ylim=y_range,xlab=xlab_text,ylab=ylab_text,xaxt="n")
|
|
| 623 | 608 |
lines(1:12,tb3_month$mae,col=c("blue"),type="b")
|
| 624 | 609 |
lines(1:12,tb4_month$mae,col=c("black"),type="b")
|
| 625 | 610 |
axis(1,at=1:12,labels=xlab_tick) |
| 626 | 611 |
title(main="Monthly average MAE") |
| 612 |
legend("topleft",legend=legend_text,
|
|
| 613 |
cex=0.9, pch=c(pch_t),col=c(col_t),lty=c(1,1,1),bty="n") |
|
| 627 | 614 |
|
| 615 |
#Second plot |
|
| 628 | 616 |
ylab_text<-"MAE (C)" |
| 629 | 617 |
xlab_text<-"Month" |
| 630 | 618 |
#y_range<-range(month_data_list$gam$mae,month_data_list$kriging$mae,month_data_list$gwr$mae) |
| 631 | 619 |
#y_range<-range(month_data_list$gam$mae) |
| 632 |
boxplot(mae~month,data=month_data_list$gam,main="GAM",ylab=ylab_text,outline=FALSE)
|
|
| 620 |
boxplot(mae~month,data=month_data_list$gam,main="Monthly MAE boxplot", xlab=xlab_text,ylab=ylab_text,outline=FALSE)
|
|
| 633 | 621 |
#boxplot(mae~month,data=month_data_list$kriging,ylim=y_range,main="Kriging",ylab=ylab_text,outline=FALSE) |
| 634 | 622 |
#boxplot(mae~month,data=month_data_list$gwr,ylim=y_range,main="GWR",ylab=ylab_text,outline=FALSE) |
| 635 | 623 |
|
Also available in: Unified diff
modifications for paper draft, single time scale method