## Note the start time
start.time<-proc.time()[3]
# figurepath <- "YourPath" # Folder to save or read the figures from
# filepath <- "YourFilePath" # Folder to save or read the tables from
longrun <- FALSE # Should we run the model to check the results.
# If FALSE we use the results, e.g. tables, figures and model fits from previous runs
# mappath <- "mapfiles"
# dpath <- "datafiles/"
# modpath <- "precip_modelfits/"
# load(file=paste0(dpath, "C8rainfall.RData"))library(bmstdr)
library(akima)
library(ggplot2)
library(rgdal)
library(spdep)
library(broom)
library(tidyr)
library(plyr)
library(RColorBrewer)
library(ggpubr)
library(spdep)
library(GGally)
library(geoR)‘RandomFieldsUtils’ will use OMP ‘RandomFields’ will use OMP
library(spTimer)
library(sp)
library(doBy)
library(grid)
library(rgeos)1 Code to reproduce Table 8.3
head(ds3)
## id site utmx utmy longitude latitude elevation aspect slope facing wkno wknoofY startdate
## 1 1 RG1 281359.8 4870165 -71.72541 43.95233 521 145.8855 12.66026 South-facing 1 2 1997-01-06
## 2 1 RG1 281359.8 4870165 -71.72541 43.95233 521 145.8855 12.66026 South-facing 2 3 1997-01-13
## 3 1 RG1 281359.8 4870165 -71.72541 43.95233 521 145.8855 12.66026 South-facing 3 4 1997-01-20
## 4 1 RG1 281359.8 4870165 -71.72541 43.95233 521 145.8855 12.66026 South-facing 4 5 1997-01-27
## 5 1 RG1 281359.8 4870165 -71.72541 43.95233 521 145.8855 12.66026 South-facing 5 6 1997-02-03
## 6 1 RG1 281359.8 4870165 -71.72541 43.95233 521 145.8855 12.66026 South-facing 6 7 1997-02-10
## enddate rainfall year month startday index sin1 cos1 sin2 cos2 xutmx xutmy xelevation
## 1 1997-01-12 26.4 1997 January 6 1 0.0344 0.9994 0.0688 0.9976 1.019005 0.8558666 -1.19204
## 2 1997-01-19 23.4 1997 January 13 2 0.0516 0.9987 0.1031 0.9947 1.019005 0.8558666 -1.19204
## 3 1997-01-26 28.9 1997 January 20 3 0.0688 0.9976 0.1373 0.9905 1.019005 0.8558666 -1.19204
## 4 1997-02-02 38.9 1997 January 27 4 0.0860 0.9963 0.1713 0.9852 1.019005 0.8558666 -1.19204
## 5 1997-02-09 29.3 1997 February 3 5 0.1031 0.9947 0.2051 0.9787 1.019005 0.8558666 -1.19204
## 6 1997-02-16 21.8 1997 February 10 6 0.1202 0.9927 0.2387 0.9711 1.019005 0.8558666 -1.19204
## xaspect xslope xsin1 xcos1 xsin2 xcos2 ts ws
## 1 0.01548894 -0.3160637 -1.724396 1.146756 -2.044907 1.212785 1 1
## 2 0.01548894 -0.3160637 -1.649013 1.140697 -1.932230 1.205267 1 1
## 3 0.01548894 -0.3160637 -1.573629 1.131176 -1.819882 1.194379 1 1
## 4 0.01548894 -0.3160637 -1.498245 1.119923 -1.708190 1.180639 1 1
## 5 0.01548894 -0.3160637 -1.423300 1.106075 -1.597156 1.163789 1 1
## 6 0.01548894 -0.3160637 -1.348355 1.088763 -1.486779 1.144087 1 1
dim(ds3)
## [1] 23784 34
site <- unique(ds3[, c("id", "site", "utmx", "utmy", "longitude", "latitude", "elevation", "ws")])
sum_annual_av <- summaryBy(rainfall ~ id+year, data = ds3, FUN = function(x) { c(m = sum(x))})
#sum_annual_av$rainfall.m <- sum_annual_av$rainfall.m
head(sum_annual_av )
## id year rainfall.m
## 1 1 1997 1243.6
## 2 1 1998 1418.1
## 3 1 1999 1371.4
## 4 1 2000 1407.1
## 5 1 2001 993.3
## 6 1 2002 1327.9
colnames(sum_annual_av)[3] <- "precipitation"
rgstats <- summaryBy(precipitation ~ id, data = sum_annual_av, FUN = function(x) { c(m = mean(x), c(s=sd(x)))})
# rgstats
# head(rgstats)
dim(rgstats)
## [1] 24 3
week_av <- summaryBy(rainfall ~ id, data = ds3, FUN = function(x) { c(m = mean(x), s=sd(x))})
v <- cbind.data.frame(Gauge=site$site, Watershed=site$ws, mean.week=round(week_av[,2], 2),
sd.week= round(week_av[,3], 2), mean.annual=round(rgstats[,2], 2),
sd.annual= round(rgstats[,3], 2))
rownames(v) <- site$id
# v
table8.3 <- v
table8.3
## Gauge Watershed mean.week sd.week mean.annual sd.annual
## 1 RG1 1 27.80 25.72 1450.02 176.36
## 2 RG2 1 28.43 26.30 1483.07 186.28
## 3 RG3 1 26.91 24.67 1403.48 165.08
## 4 RG4 3 27.48 25.12 1433.26 166.40
## 5 RG5 3 27.68 24.98 1443.84 164.10
## 6 RG6 4 28.39 25.76 1480.65 177.31
## 7 RG7 4 28.19 25.71 1470.43 171.42
## 8 RG8 4 28.18 25.85 1469.97 173.39
## 9 RG9 6 28.58 26.11 1490.80 187.80
## 10 RG10 6 29.23 26.57 1524.58 184.56
## 11 RG11 6 28.25 26.14 1473.36 174.79
## 12 RG12 7 29.14 26.48 1519.93 192.31
## 13 RG13 7 29.33 26.83 1529.91 193.39
## 14 RG14 7 29.72 27.38 1550.19 208.74
## 15 RG15 7 29.37 26.96 1531.77 197.50
## 16 RG16 7 29.59 26.69 1543.55 192.02
## 17 RG17 8 29.10 26.52 1517.77 193.40
## 19 RG19 8 29.29 27.45 1527.55 203.24
## 20 RG20 8 29.28 27.17 1527.14 194.59
## 21 RG21 9 29.81 27.13 1554.59 195.25
## 22 RG22 10 24.63 23.07 1284.78 161.22
## 23 RG23 9 29.61 27.62 1544.50 189.77
## 24 RG24 9 30.83 28.27 1608.17 207.25
## 25 RG25 9 28.87 26.33 1505.89 186.65
dput(table8.3, file=paste0(filepath, "table8.3.txt"))
wsheds <- readOGR(dsn=path.expand(mappath), layer = "wsheds_polygon")
## OGR data source with driver: ESRI Shapefile
## Source: "/home/sks/Dropbox/sks/bookproject/bmstdr_Sahu/chapter_code_files/mapfiles", layer: "wsheds_polygon"
## with 9 features
## It has 10 fields
## Integer64 fields read as strings: WSHEDS_ WSHEDS_ID WSHEDS1_ WSHEDS1_ID SMLBND_ SMLBND_ID
udat <- tidy(wsheds)
## Regions defined for each Polygons
centers <- gCentroid(wsheds, byid=TRUE)
centers <- as.data.frame(centers)
centers$ws <- paste0("W", c(3, 2, 1, 4, 6, 5, 8, 7, 9))
xo <- seq(from=min(ds3$utmx)-10, to = max(ds3$utmx)+10, length=200)
yo <- seq(from=min(ds3$utmy)-10, to = max(ds3$utmy)+ 10, length=200)
a <- c(277325, 281995)
b <- c( 4866015, 4871045)
xminn <- a[1]
xmaxx <- a[2]
yminn <- b[1]
ymaxx <- b[2] 2 Reproduce the EDA figures
summary(ds3$rainfall)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 8.80 22.40 28.65 41.20 197.30
count(ds3$rainfall == 0) ## 1247
## x freq
## 1 FALSE 22537
## 2 TRUE 1247
100*count(ds3$rainfall == 0)/nrow(ds3) ## 5.24% 0 rain fall
## x freq
## 1 0.000000000 94.756979
## 2 0.004204507 5.243021
100*count(ds3$rainfall < 0.9)/nrow(ds3) ## 6.70% less than 0.9
## x freq
## 1 0.000000000 93.30222
## 2 0.004204507 6.69778
## Rain histogram
p1 <- ggplot(ds3, aes(x = rainfall))+
geom_histogram(bins = 100, color = "blue", fill = "white") +
geom_vline(aes(xintercept = mean(rainfall)), linetype = "dashed", size = 0.6,
colour = "red")+
labs(x = "Precipitation (mm)", y = "Frequency")+
scale_x_continuous(breaks = seq(0,200,20))+
scale_y_continuous(breaks = seq(0,1600,200))+
theme(legend.position = "right", axis.text.x = element_text(size=12),
axis.title.y = element_text(size = 16), axis.title.x = element_text(size=16),
axis.text.y = element_text(size = 12))
p1
ggsave(filename=paste0(figurepath, "rain_hist_orig.png"))
## Saving 7 x 5 in image
ds3$sq_rain <- sqrt(ds3$rainfall)
ds3$log_rain <- log(ds3$rainfall)
# Frequency distribution (square-root scale)
p2 <- ggplot(ds3, aes(x = sq_rain)) +
geom_histogram(bins = 100, color = "blue", fill = "white") +
geom_vline(aes(xintercept = mean(sq_rain)), linetype = "dashed", size = 0.6,
colour = "red")+
labs(x = "Square-root precipitation (mm)", y = "Frequency")+
scale_x_continuous(breaks = seq(0,18,2))+
scale_y_continuous(breaks = seq(0,1400,200))+
theme(legend.position = "right", axis.text.x = element_text(size=12),
axis.title.y = element_text(size = 16), axis.title.x = element_text(size=16),
axis.text.y = element_text(size = 12))
p2
ggsave(filename=paste0(figurepath, "rain_hist_sqrt.png"))
## Saving 7 x 5 in image
## Mean vs variance on original scale
sum_original <- summaryBy(rainfall ~ id + year, data = ds3,
FUN = function(x) { c(m = mean(x), s = sd(x))})
sum_original$rainfall.v <- (sum_original$rainfall.s)^2
u <- lm(data=sum_original, formula=rainfall.v~rainfall.m)
coef(u)
## (Intercept) rainfall.m
## -650.01764 46.83705
p3 <- ggplot(sum_original, aes(x=rainfall.m, y=rainfall.v, shape=id, color=year)) +
geom_point(size=2) + scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,11,12,
13,14,15,16,17,18,19,20,21,22,23,24))+
labs(#title="Mean versus variance by site and year (original scale)",#
x = "Mean original (mm)", y = "Variance")+
geom_abline(intercept = -650, slope = 46.8, color="black", linetype="solid", size=1)+
theme(legend.position = NULL, legend.key.size = unit(0.4, "cm"),
axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 12))
p3
ggsave(filename=paste0(figurepath, "mean_var_orig.png"))
## Saving 7 x 5 in image
head(sum_original)
## id year rainfall.m rainfall.s rainfall.v
## 1 1 1997 23.91538 19.43387 377.6754
## 2 1 1998 27.27115 26.60934 708.0570
## 3 1 1999 26.37308 28.79311 829.0432
## 4 1 2000 27.05962 24.62686 606.4821
## 5 1 2001 18.74151 19.77416 391.0175
## 6 1 2002 25.53654 17.67785 312.5063
# on Log scale
sum_log <- summaryBy(log_rain ~ id + year , data = ds3,
FUN = function(x) { c(m = mean(x), s = sd(x))})
sum_log$log_rain.v <- (sum_log$log_rain.s)^2
p4 <- ggplot(sum_log, aes(x=log_rain.m, y=log_rain.v, shape=id, color=year))+
geom_point(size=2) + scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,11,12,
13,14,15,16,17,18,19,20,21,22,23,24))+
labs(#title="Mean versus variance by site and year (log scale)"#,
x = "Mean log precipitation (mm)", y = "Variance")+
theme(legend.position = NULL, legend.key.size = unit(0.4, "cm"),
axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 12))
p4
## Warning: Removed 435 rows containing missing values (geom_point).
ggsave(filename=paste0(figurepath, "mean_var_log.png"))
## Saving 7 x 5 in image
## Warning: Removed 435 rows containing missing values (geom_point).
# On Square-root scale
sum_sqrt <- summaryBy(sq_rain ~ id + year , data = ds3,
FUN = function(x) { c(m = mean(x), s = sd(x))})
sum_sqrt$sq_rain.v <- (sum_sqrt$sq_rain.s)^2
head(sum_sqrt)
## id year sq_rain.m sq_rain.s sq_rain.v
## 1 1 1997 4.322541 2.309457 5.333593
## 2 1 1998 4.602614 2.491275 6.206452
## 3 1 1999 4.466700 2.558824 6.547580
## 4 1 2000 4.699064 2.253004 5.076026
## 5 1 2001 3.799500 2.094780 4.388102
## 6 1 2002 4.706437 1.858059 3.452382
u <- lm(data=sum_sqrt, formula=sq_rain.v~sq_rain.m)
# coef(u)
p5 <- ggplot(sum_sqrt, aes(x=sq_rain.m, y=sq_rain.v, shape=id, color=year))+
geom_point(size=2) + scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,11,12,
13,14,15,16,17,18,19,20,21,22,23,24))+
geom_abline(intercept = 0.30, slope = 1.25, color="black", linetype="solid", size=1)+
labs(#title="Mean versus variance by site and year (square-root scale)"#,
x = "Mean square-root", y = "Variance")+
theme(legend.position = NULL, legend.key.size = unit(0.4, "cm"),
axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
axis.text.x = element_text(size = 12), axis.text.y = element_text(size = 12))
p5
ggsave(file=paste0(figurepath, "mean_var_sqrt.png"))
## Saving 7 x 5 in image
ggarrange(p3, p5, common.legend = TRUE, legend = "bottom", nrow = 1, ncol = 2)
ggsave(filename=paste0(figurepath, "rain_mean_vars.png"), width=6, height=5, dpi=300, device = "png")
## rain v elevation
sum_sqrt_id <- summaryBy(sq_rain ~ id + elevation + longitude + latitude, data = ds3,
FUN = function(x) { c(m = mean(x), s = sd(x))})
p8 <- ggplot(sum_sqrt_id, aes(x=elevation, y=sq_rain.m, color=id))+
geom_point(shape=16, size=2)+
stat_smooth(method = "lm", col = "red")+
labs(x = "Elevation (m)", y = "Average square-root precipitation (mm)")+
geom_text(aes(label=id), hjust = -0.7, show.legend = F)+
theme(legend.position = "none", axis.text.x = element_text(size=12),
axis.title.y = element_text(size = 16), axis.title.x = element_text(size=16),
axis.text.y = element_text(size = 12))
p8
## `geom_smooth()` using formula 'y ~ x'
ggsave(filename=paste0(figurepath, "rain_v_elevation.png"))
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
## Drawing maps of elevation
u <- pgrid
dimnames(u)[[2]] <- c("utmx", "utmy", "ws", "ele", "asp", "slp")
wiw <- u
a <- table(pgrid$Watershednum)
wnum <- 0.04 * a
wnum
##
## 1 2 3 4 5 6 7 8 9
## 56.48 71.80 160.96 155.20 83.40 59.68 306.40 244.32 281.24
utm <- data.frame(x=wiw$utmx,y=wiw$utmy)
coordinates(utm) <- ~x+y
class(utm)
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"
proj4string(utm) <- CRS("+proj=utm +zone=19 +datum=WGS84 +units=m +ellps=WGS84")
lonlat <- spTransform(utm,CRS("+proj=longlat +datum=WGS84"))
wiw$long <- lonlat$x
wiw$lat <- lonlat$y
wshedlonglat <- spTransform(wsheds,CRS("+proj=longlat +datum=WGS84"))
class(wshedlonglat)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
p9 <- ggplot(data=wiw, aes(x=utmx, y=utmy)) +
geom_tile(data=wiw, aes(x=utmx, y=utmy, fill=ele, group=ws)) +
geom_polygon(data = wsheds, aes(x = long, y = lat, group = group), colour = "black", fill = NA) +
guides(fill=guide_colorbar(title="Elevation"))+
scale_fill_gradientn(colours=com,na.value="black",limits=c(358, 844))+
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
geom_text(data=site, aes(x=utmx, y=utmy-275, label=site), size=2.5, nudge_x=100, nudge_y=100) +
geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
# geom_text(data=wlabel, aes(x, y, label=wsl), size=3) +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank()) +
theme(legend.title = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx )
## Regions defined for each Polygons
ggsn::north2(p9, x=0.75, y=0.25, symbol=12) 
ggsave(file=paste0(figurepath, "elevation_hubbardbrook.png"))
## Saving 7 x 5 in image
# [Figure 13(b)]: Longitude
p81 <- ggplot(sum_sqrt_id, aes(x=longitude, y=sq_rain.m, color=id))+
geom_point(shape=16, size=2)+
stat_smooth(method = "lm", col = "red")+
labs(x = "Longitude", y = "Average square-root precipitation (mm)")+
geom_text(aes(label=id), hjust = -0.7, show.legend = F)+
theme(legend.position = "none", axis.text.x = element_text(size=12),
axis.title.y = element_text(size = 16), axis.title.x = element_text(size=16),
axis.text.y = element_text(size = 12))
p81
## `geom_smooth()` using formula 'y ~ x'
ggsave(filename=paste0(figurepath, "rain_v_longitude.png"))
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
# [Figure 13(c)]: Latitude
p91 <- ggplot(sum_sqrt_id, aes(x=latitude, y=sq_rain.m, color=id))+
geom_point(shape=16, size=2)+
stat_smooth(method = "lm", col = "red")+
labs(x = "Latitude", y = "Average square-root precipitation (mm)")+
geom_text(aes(label=id), hjust = -0.7, show.legend = F)+
theme(legend.position = "none", axis.text.x = element_text(size=12),
axis.title.y = element_text(size = 16), axis.title.x = element_text(size=16),
axis.text.y = element_text(size = 12))
p91
## `geom_smooth()` using formula 'y ~ x'
ggsave(file=paste0(figurepath, "rain_v_latitude.png"))
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
ggarrange(p81, p91, common.legend = TRUE, legend = "none", nrow = 1, ncol = 2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
ggsave(filename=paste0(figurepath, "ann_and_rolling.png"), width=6, height=3, dpi=300, device = "png")
# ggsave(paste0(figurepath, "rain_v_latitude.png"), width=4, height=3.44, dpi=300)
## Temporal pattern
# 1997 -- 2006
p9 <- ggplot(ds3[ds3$year %in% c("1997","1998","1999","2000","2001","2002","2003","2004",
"2005","2006"),], aes(x=month, y=sq_rain, color=month)) +
geom_boxplot() + facet_wrap(~year, ncol=5) +
stat_summary(fun=mean, geom = "line", aes(group=1))+
stat_summary(fun=mean, geom="point")+
labs(#title="Boxplot of square-root precipitation against months by year",
x = "Month", y = "Square-root precipitation (mm)")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 90, hjust = 1, size=11),
axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
axis.text.y = element_text(size = 12))
ggsave(filename=paste0(figurepath, "9a.png"))
## Saving 7 x 5 in image
# 2007 -- 2015
p10 <- ggplot(ds3[ds3$year %in% c("2007","2008","2009","2010","2011","2012","2013","2014","2015"),],
aes(x=month, y=sq_rain, color=month)) + geom_boxplot() + facet_wrap(~year, ncol=5)+
stat_summary(fun=mean, geom = "line", aes(group=1))+
stat_summary(fun=mean, geom="point")+
labs(#title="Boxplot of square-root rain against months by year",
x = "Month", y = "")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 90, hjust = 1, size=11),
axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
axis.text.y = element_text(size = 12))
p10
ggsave(filename=paste0(figurepath, "9b.png"))
## Saving 7 x 5 in image
class(ds3$month)
## [1] "factor"
ds3$nmonth <- ds3$month
levels(ds3$nmonth) <- as.character(1:12)
table(ds3$nmonth)
##
## 1 2 3 4 5 6 7 8 9 10 11 12
## 1992 1824 2064 1944 1992 1992 1992 2016 1968 1992 1968 2040
p11 <- ggplot(ds3, aes(x=nmonth, y=sq_rain, color=nmonth)) +
geom_boxplot() + facet_wrap(~year, ncol=4) +
stat_summary(fun=mean, geom = "line", aes(group=1))+
stat_summary(fun=mean, geom="point")+
labs(x = "Month", y = "Square-root precipitation (mm)")+
theme(legend.position = "none") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=11),
axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
axis.text.y = element_text(size = 12))
p11
ggsave(filename=paste0(figurepath, "9ab.png"))
## Saving 7 x 5 in image
## Annual rain falls
# load(file="ds3.RData")
sumsite <- summaryBy(rainfall ~ id +year, data = ds3,
FUN = function(x) { c(s = mean(x))})
dim(sumsite) ## 19 years time 24 sites =456
## [1] 456 3
summary(sumsite$rainfall.s)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 17.10 26.43 28.23 28.66 31.08 37.62
# summary(sumsite$rainfall.s)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 17.10 26.43 28.23 28.66 31.08 37.62
## Work with mean
sumsite$year <- as.numeric(as.character(sumsite$year))
summary(sumsite)
## id year rainfall.s
## 1 : 19 Min. :1997 Min. :17.10
## 2 : 19 1st Qu.:2001 1st Qu.:26.43
## 3 : 19 Median :2006 Median :28.23
## 4 : 19 Mean :2006 Mean :28.66
## 5 : 19 3rd Qu.:2011 3rd Qu.:31.08
## 6 : 19 Max. :2015 Max. :37.62
## (Other):342
pannual <- ggplot(data=sumsite, aes(y=rainfall.s*52, x=year, color=id)) +
labs(x = "Year", y = "Average annual precipitation (mm)")+
geom_line()
pannual
ggsave(filename=paste0(figurepath, "average_annual.png"))
## Saving 7 x 5 in image
## Provides 3 year rolling averages
three_year_rolling_av <- function(x) {
n <- length(x)
y <- rep(NA, n-2)
for (k in 3:n) {
y[k-2] <- mean(x[(k-2):k])
}
y
}
a <- summaryBy(rainfall.s ~ id, data = sumsite, FUN = three_year_rolling_av)
# b <- matrix(1:12, 3, 4)
b <- as.vector(t(a[,-1]))
rolling <- data.frame(id=rep(a$id, each=17), year=rep(1999:2015, 24), rainfall=b)
###
head(rolling)
## id year rainfall
## 1 1 1999 25.85321
## 2 1 2000 26.90128
## 3 1 2001 24.05807
## 4 1 2002 23.77922
## 5 1 2003 25.20166
## 6 1 2004 27.41538
dim(rolling) # 408 = 17 * 24
## [1] 408 3
table(rolling$id)
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 19 20 21 22 23 24 25
## 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17
# site <- unique(ds3[, c("id", "utmx", "utmy", "ws")])
head(site)
## id site utmx utmy longitude latitude elevation ws
## 1 1 RG1 281359.8 4870165 -71.72541 43.95233 521 1
## 9911 2 RG2 281196.7 4870489 -71.72737 43.95493 626 1
## 16848 3 RG3 281109.7 4870900 -71.72866 43.95846 697 1
## 17839 4 RG4 281691.7 4870494 -71.72088 43.95496 606 3
## 18830 5 RG5 281907.4 4870948 -71.71883 43.95953 651 3
## 19821 6 RG6 280574.1 4870760 -71.73612 43.95690 724 4
rollingsites <- merge(rolling, site)
head(rollingsites)
## id year rainfall site utmx utmy longitude latitude elevation ws
## 1 1 1999 25.85321 RG1 281359.8 4870165 -71.72541 43.95233 521 1
## 2 1 2000 26.90128 RG1 281359.8 4870165 -71.72541 43.95233 521 1
## 3 1 2001 24.05807 RG1 281359.8 4870165 -71.72541 43.95233 521 1
## 4 1 2002 23.77922 RG1 281359.8 4870165 -71.72541 43.95233 521 1
## 5 1 2003 25.20166 RG1 281359.8 4870165 -71.72541 43.95233 521 1
## 6 1 2004 27.41538 RG1 281359.8 4870165 -71.72541 43.95233 521 1
wsaverage <- summaryBy(rainfall ~ ws+year, data = rollingsites,
FUN = function(x) { c(s = mean(x))})
head(wsaverage)
## ws year rainfall.s
## 1 1 1999 25.70962
## 2 1 2000 26.87543
## 3 1 2001 24.13895
## 4 1 2002 23.95134
## 5 1 2003 25.24365
## 6 1 2004 27.26474
wsy <- wsaverage[wsaverage$year==2015, ]
wsy$rainfall.s <- wsy$rainfall.s * 52
wsy
## ws year rainfall.s
## 17 1 2015 1385.856
## 34 3 2015 1373.567
## 51 4 2015 1396.778
## 68 6 2015 1408.589
## 85 7 2015 1441.453
## 102 8 2015 1443.244
## 119 9 2015 1480.100
## 136 10 2015 1206.433
a <- rolling[rolling$year==2010,]
a <- a[order(a$rainfall), ]
twosites <- a[c(1, 2, 23, 24), ]
twosites
## id year rainfall
## 352 22 2010 27.14167
## 46 3 2010 28.90897
## 233 14 2010 32.55833
## 386 24 2010 33.34808
proll <- ggplot(data=rolling, aes(x=year, y=rainfall*52, color=id)) +
geom_line() +
geom_text(data=twosites, nudge_y = -20, aes(label=id, x=year, y=rainfall*52, color=id, size=2))+
labs(x = "Year", y = "Three rolling average precipitation (mm)")
proll
ggsave(filename=paste0(figurepath, "rolling_average_annual.png"))
## Saving 7 x 5 in image
ggarrange(pannual, proll, common.legend = TRUE, legend = "none", nrow = 1, ncol = 2)
ggsave(filename=paste0(figurepath, "ann_and_rolling.png"), width=6, height=3, dpi=300, device = "png")3 Reproduce Figure 8.2 Phi choice
sn <- length(unique(ds3$id))
tn <- nrow(ds3)/sn
sn
tn
f2 <- rainfall ~ xutmx + xelevation + xsin1 + xcos1
valids <- c(21, 22)
# Find the validation rows
vs <- getvalidrows(sn=sn, tn=tn, valids=c(21, 22), allt=T)
# 24 sites
## want to validate at RG22 and RG23 which are site numbers
## 21 and 22.
if (longrun) {
nItr <- 6000 # number of MCMC samples for each model
nBurn <- 1000 # number of burn-in from the MCMC samples
## Select phis
phis <- seq(from=0.01, to=1.5, by=0.03)
k <- length(phis)
u2 <- matrix(NA, nrow=k, ncol=4)
vs <- getvalidrows(sn=sn, tn=tn, valids=c(21), allt=T)
for (i in 1:k) {
cat("i=", i, "to go to", k)
b3 <- Bsptime(data=ds3, formula=f2, package="spTimer", scale.transform ="NONE",
time.data=time.data, prior.phi = "FIXED", prior.phi.param =phis[i],
truncation.para = list(at = 0,lambda =2), validrows=vs,
model="truncatedGP", coordtype="utm", coords=3:4, n.report= 1,
N=nItr, burn.in = nBurn, mchoice = FALSE)
u2[i, ] <- unlist(b3$stats)
}
phiselect <- cbind.data.frame(phis=phis, u2)
colnames(phiselect) <- c("phis", "rmse", "mae", "crps", "cvg")
phiselect
plot(phiselect$phis, phiselect$rmse)
k <- nrow(phiselect)
aphi <- data.frame(phis=rep(phiselect$phis, 3),
Criteria=c(rep("rmse", k), rep("mae", k), rep("crps", k)),
vcrit=c(phiselect$rmse, phiselect$mae, phiselect$crps))
aphi$Criteria <- as.factor(aphi$Criteria)
aphi$Criteria <- factor(aphi$Criteria, levels=rev(levels(aphi$Criteria)))
save(phiselect, aphi, file=paste0(modpath, "precip_example_phiselect.RData"))
} else {
load(file=paste0(modpath, "precip_example_phiselect.RData"))
}
p <- ggplot() +
geom_point(data=aphi, aes(x=phis, y=vcrit, group=Criteria, color=Criteria, shape=Criteria)) +
labs(x =expression(phi), y = "criteria value")+
theme(legend.position=c(0.2, 0.85))
p
ggsave(paste0(figurepath, "rainfall_phichoice.png"))4 Table 8.4 Parameter estimates
if (longrun) {
nItr <- 6000 # number of MCMC samples for each model
nBurn <- 1000 # number of burn-in from the MCMC samples
f2 <- rainfall ~ xutmx + xelevation + xsin1 + xcos1
### Parameter estimates
M1 <- Bsptime(data=ds3, formula=f2, package="spTimer", scale.transform ="NONE",
time.data=time.data, prior.phi = "FIXED", prior.phi.param =1,
truncation.para = list(at = 0,lambda =2),
model="truncatedGP",coordtype="utm", coords=3:4, n.report= 2,
N=nItr, burn.in = nBurn, mchoice = T)
table8.4 <- round(M1$params[1:7, ], 4)
table8.4
dput(table8.4, file=paste0(filepath, "table8.4.txt"))
# gof penalty pmcc
# 1387.94 9449.32 10837.26
} else {
table8.4 <- dget(file=paste0(filepath, "table8.4.txt"))
}
table8.4
## mean sd 2.5% 97.5%
## (Intercept) 10.6824 0.0101 10.6625 10.7019
## xutmx -0.0466 0.0089 -0.0645 -0.0294
## xelevation 0.0049 0.0060 -0.0066 0.0167
## xsin1 0.3090 0.0353 0.2391 0.3784
## xcos1 0.2093 0.0354 0.1398 0.2785
## sig2eps 0.0051 0.0001 0.0050 0.0052
## sig2eta 0.5130 0.0047 0.5037 0.52255 Table 8.5 (validation)
if (longrun) {
f2 <- rainfall ~ xutmx + xelevation + xsin1 + xcos1
nItr <- 6000 # number of MCMC samples for each model
nBurn <- 1000 # number of burn-in from the MCMC samples
u1 <- matrix(NA, nrow=24, ncol=4)
for (i in 1:24) {
vs <- getvalidrows(sn=sn, tn=tn, valids=c(i), allt=T)
cat("Out of 24, i =", i, "\n")
b <- Bsptime(data=ds3, formula=f2, package="spTimer", scale.transform ="NONE",
time.data=time.data, validrows=vs,
prior.phi = "FIXED", prior.phi.param = 1.0,
truncation.para = list(at = 0,lambda =2),
model="truncatedGP", coordtype="utm", coords=3:4, n.report= 2,
N=nItr, burn.in = nBurn, mchoice = F)
u1[i, ] <- unlist(b$stats)
}
head(ds3)
sum_annual_av <- summaryBy(rainfall ~ id+year, data = ds3, FUN = function(x) { c(m = sum(x))})
head(sum_annual_av )
#sum_annual_av$rainfall.m <- sum_annual_av$rainfall.m
head(sum_annual_av )
colnames(sum_annual_av)[3] <- "precipitation"
rgstats <- summaryBy(precipitation ~ id, data = sum_annual_av, FUN = function(x) { c(m = mean(x), c(s=sd(x)))})
rgstats
dim(rgstats)
site <- unique(ds3[, c("id", "site")])
v <- cbind.data.frame(Gauge=site$site, rmse=round(u1[,1], 2),
mae= round(u1[,2], 2), crps= round(u1[,3], 2), cvg = round(u1[,4], 2))
anntable <- v
rownames(anntable) <- site$id
dput(anntable, file=paste0(filepath, "table8.5.txt"))
} else
table8.5 <- dget(file=paste0(filepath, "table8.5.txt")
table8.5
load(file=paste0(modpath, "pred_original_frames.RData"))
year <- 2015
dplot <- predandfit[predandfit$year==year, ]
# dplot <- obigdat[obigdat$year==year, ] ## predfit better
dim(dplot)
## [1] 1510 9
## annual dplot
summary(dplot)
## pid utmx utmy ws year mean
## Min. : 1 Min. :277345 Min. :4866025 Min. : 1.000 2015 :1510 Min. :24.82
## 1st Qu.: 7906 1st Qu.:278415 1st Qu.:4866545 1st Qu.: 4.000 1997 : 0 1st Qu.:25.91
## Median :17341 Median :279265 Median :4867005 Median : 7.000 1998 : 0 Median :27.12
## Mean :17414 Mean :279571 Mean :4868184 Mean : 6.136 1999 : 0 Mean :26.84
## 3rd Qu.:26795 3rd Qu.:280785 3rd Qu.:4870315 3rd Qu.: 8.000 2000 : 0 3rd Qu.:27.67
## Max. :35487 Max. :283225 Max. :4871045 Max. :10.000 2001 : 0 Max. :28.71
## (Other): 0
## sd low up
## Min. :1.394 Min. :19.89 Min. :28.36
## 1st Qu.:2.017 1st Qu.:22.01 1st Qu.:29.98
## Median :2.096 Median :23.03 Median :31.04
## Mean :2.091 Mean :22.87 Mean :30.89
## 3rd Qu.:2.172 3rd Qu.:23.69 3rd Qu.:31.75
## Max. :2.511 Max. :25.34 Max. :33.39
##
dplot$mean <- dplot$mean*52
surf <- interp(dplot$utmx, dplot$utmy, dplot$mean, xo=xo, yo=yo)
names(surf)
## [1] "x" "y" "z"
surf <- list(x=surf$x, y=surf$y, z=surf$z)
# pnts <- data.frame(long=surf$x, lat=surf$y)
pnts <- expand.grid(surf$x, surf$y)
colnames(pnts) <- c("long", "lat")
pts1 <- pnts
coordinates(pnts)<- ~long + lat
pnts <- SpatialPoints(coords = pnts)
proj4string(pnts) <- proj4string(wsheds)
## Warning in proj4string(wsheds): CRS object has comment, which is lost in output
proj4string(wsheds)
## Warning in proj4string(wsheds): CRS object has comment, which is lost in output
## [1] "+proj=utm +zone=19 +datum=NAD27 +units=m +no_defs"
proj4string(pnts)
## Warning in proj4string(pnts): CRS object has comment, which is lost in output
## [1] "+proj=utm +zone=19 +datum=NAD27 +units=m +no_defs"
identicalCRS(wsheds, pnts)
## [1] TRUE
a1 <- sp::over(pnts, wsheds, fn = NULL) #, returnList = TRUE)
# summary(a1)
v <- pts1[!is.na(a1), ]
z <- as.vector(surf$z)
length(z)
## [1] 40000
summary(z)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1295 1351 1393 1390 1429 1477 18159
z[is.na(a1)] <- NA
z <- matrix(z, ncol=200, nrow=200)
interp1 <- data.frame(long = surf$x, z)
names(interp1)[1:length(surf$y)+1] <- surf$y
interp1 <- gather(interp1, key = lat, value =precipitation,-long,convert = TRUE)
zr <- range(interp1$precipitation, na.rm=T)
# a <- range(site$utmx)
# b <- range(site$utmy)
#library(rgeos)
a <- range(interp1$long)
b <- range(interp1$lat)
#library(RColorBrewer)
pcolors <- brewer.pal(8,"Greens")
# pcolors <- brewer.pal(9,"YlOrRd")
# display.brewer.all()
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
p1 <- ggplot(data=interp1, aes(x=long, y=lat)) +
geom_raster(aes(x=long, y=lat, fill=precipitation)) +
geom_polygon(data = wsheds, aes(x = long, y = lat, group = group), colour = "black", fill = NA) +
guides(fill=guide_colorbar(title="Precipitation"))+
scale_fill_gradientn(colours=pcolors, na.value="gray95", limits=zr) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("Annual total precipitation in ", year), x="", y = "", size=2.5)
## Regions defined for each Polygons
ggsn::north2(p1, x=0.75, y=0.25, symbol=12) 
ggsave(file=paste0(figurepath, "rain2015_total.png"))
## Saving 7 x 5 in image
### Repeat for annual sd
year <- 2015
dplot <- predandfit[predandfit$year==year, ]
dim(dplot)
## [1] 1510 9
## annual dplot
summary(dplot)
## pid utmx utmy ws year mean
## Min. : 1 Min. :277345 Min. :4866025 Min. : 1.000 2015 :1510 Min. :24.82
## 1st Qu.: 7906 1st Qu.:278415 1st Qu.:4866545 1st Qu.: 4.000 1997 : 0 1st Qu.:25.91
## Median :17341 Median :279265 Median :4867005 Median : 7.000 1998 : 0 Median :27.12
## Mean :17414 Mean :279571 Mean :4868184 Mean : 6.136 1999 : 0 Mean :26.84
## 3rd Qu.:26795 3rd Qu.:280785 3rd Qu.:4870315 3rd Qu.: 8.000 2000 : 0 3rd Qu.:27.67
## Max. :35487 Max. :283225 Max. :4871045 Max. :10.000 2001 : 0 Max. :28.71
## (Other): 0
## sd low up
## Min. :1.394 Min. :19.89 Min. :28.36
## 1st Qu.:2.017 1st Qu.:22.01 1st Qu.:29.98
## Median :2.096 Median :23.03 Median :31.04
## Mean :2.091 Mean :22.87 Mean :30.89
## 3rd Qu.:2.172 3rd Qu.:23.69 3rd Qu.:31.75
## Max. :2.511 Max. :25.34 Max. :33.39
##
dplot$sd <- dplot$sd*52
xo <- seq(from=min(dplot$utmx)-10, to = max(dplot$utmx)+10, length=200)
yo <- seq(from=min(dplot$utmy)-10, to = max(dplot$utmy)+ 10, length=200)
surf <- interp(dplot$utmx, dplot$utmy, dplot$sd, xo=xo, yo=yo)
names(surf)
## [1] "x" "y" "z"
surf <- list(x=surf$x, y=surf$y, z=surf$z)
# pnts <- data.frame(long=surf$x, lat=surf$y)
pnts <- expand.grid(surf$x, surf$y)
colnames(pnts) <- c("long", "lat")
pts1 <- pnts
coordinates(pnts)<- ~long + lat
pnts <- SpatialPoints(coords = pnts)
proj4string(pnts) <- proj4string(wsheds)
## Warning in proj4string(wsheds): CRS object has comment, which is lost in output
proj4string(wsheds)
## Warning in proj4string(wsheds): CRS object has comment, which is lost in output
## [1] "+proj=utm +zone=19 +datum=NAD27 +units=m +no_defs"
proj4string(pnts)
## Warning in proj4string(pnts): CRS object has comment, which is lost in output
## [1] "+proj=utm +zone=19 +datum=NAD27 +units=m +no_defs"
identicalCRS(wsheds, pnts)
## [1] TRUE
a1 <- sp::over(pnts, wsheds, fn = NULL) #, returnList = TRUE)
# summary(a1)
v <- pts1[!is.na(a1), ]
z <- as.vector(surf$z)
length(z)
## [1] 40000
summary(z)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 79.29 95.12 100.40 100.81 106.38 128.25 20299
z[is.na(a1)] <- NA
z <- matrix(z, ncol=200, nrow=200)
interp1 <- data.frame(long = surf$x, z)
names(interp1)[1:length(surf$y)+1] <- surf$y
interp1 <- gather(interp1, key = lat, value =sd,-long,convert = TRUE)
zr <- range(interp1$sd, na.rm=T)
# a <- range(site$utmx)
# b <- range(site$utmy)
a <- range(interp1$long)
b <- range(interp1$lat)
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
pcolors <- brewer.pal(9,"YlOrRd")
p <- ggplot(data=interp1, aes(x=long, y=lat)) +
geom_raster(aes(x=long, y=lat, fill=sd)) +
geom_polygon(data = wsheds, aes(x = long, y = lat, group = group), colour = "black", fill = NA) +
guides(fill=guide_colorbar(title="sd"))+
scale_fill_gradientn(colours=pcolors, na.value="gray95", limits=zr) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("sd of annual total precipitation in ", year), x="", y = "", size=2.5)
## Regions defined for each Polygons
ggsn::north2(p, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath, "rain2015_total_sd.png"))
## Saving 7 x 5 in image
head(owsdat)
## ws year mean sd low up
## 1 1 1997 25.62869 0.3317683 24.96394 26.28748
## 2 1 1998 25.90971 0.3317794 25.28463 26.57166
## 3 1 1999 25.75004 0.3728270 25.03136 26.45557
## 4 1 2000 25.77275 0.3558085 25.14459 26.48592
## 5 1 2001 25.25588 0.3270117 24.57637 25.87958
## 6 1 2002 25.71273 0.3752823 24.96885 26.39994
# 2001 and 2011
d0111 <- rbind(owsdat[owsdat$year==2001, ], owsdat[owsdat$year==2011, ])
zr <- range(d0111$mean) * 52
year <- 2001
dplot <- owsdat[owsdat$year==year, ]
k <- c(3, 2, 1, 4, 6, 5, 8, 7, 9)
udat$rain <- NA
for(i in 1:9){
udat$rain[which(udat$id==(i-1))] <- dplot$mean[k[i]] * 52
}
# zr <- range(dplot$mean) * 52
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
p1 <- ggplot() +
geom_polygon(data = udat, aes(x = long, y = lat, group = group, fill=rain), colour = "black") +
guides(fill=guide_colorbar(title="Precipitation"))+
scale_fill_gradientn(colours=com, na.value="gray95", limits=zr) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("Annual total precipitation in ", year), x="", y = "", size=2.5)
ggsn::north2(p1, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath, "rain2001_hbes_wshed.png"))
## Saving 7 x 5 in image
year <- 2011
dplot <- owsdat[owsdat$year==year, ]
k <- c(3, 2, 1, 4, 6, 5, 8, 7, 9)
udat$rain <- NA
for(i in 1:9){
udat$rain[which(udat$id==(i-1))] <- dplot$mean[k[i]] * 52
}
#zr <- range(dplot$mean) * 52
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
p11 <- ggplot() +
geom_polygon(data = udat, aes(x = long, y = lat, group = group, fill=rain), colour = "black") +
guides(fill=guide_colorbar(title="Precipitation"))+
scale_fill_gradientn(colours=com, na.value="gray95", limits=zr) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
ylab("")+xlab("") + theme(axis.text = element_blank(), axis.ticks = element_blank()) +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("Annual total precipitation in ", year), x="", y = "", size=2.5)
ggsn::north2(p11, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath, "rain2011_hbes_wshed.png"))
## Saving 7 x 5 in image
d0111 <- rbind(owsdat[owsdat$year==2001, ], owsdat[owsdat$year==2011, ])
zr <- range(d0111$sd) * 52
year <- 2001
dplot <- owsdat[owsdat$year==year, ]
k <- c(3, 2, 1, 4, 6, 5, 8, 7, 9)
udat$sd <- NA
for(i in 1:9){
udat$sd[which(udat$id==(i-1))] <- dplot$sd[k[i]] * 52
}
# zr <- range(dplot$mean) * 52
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
p3 <- ggplot() +
geom_polygon(data = udat, aes(x = long, y = lat, group = group, fill=sd), colour = "black") +
guides(fill=guide_colorbar(title="sd"))+
scale_fill_gradientn(colours=com, na.value="gray95", limits=zr) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
ylab("")+xlab("") +theme(axis.text = element_blank(), axis.ticks = element_blank()) +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("sd of annual total precipitation in ", year), x="", y = "", size=2.5)
ggsn::north2(p3, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath, "rain2001_hbes_wshed_sd.png"))
## Saving 7 x 5 in image
year <- 2011
dplot <- owsdat[owsdat$year==year, ]
k <- c(3, 2, 1, 4, 6, 5, 8, 7, 9)
udat$sd <- NA
for(i in 1:9){
udat$sd[which(udat$id==(i-1))] <- dplot$sd[k[i]] * 52
}
#zr <- range(dplot$mean) * 52
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
p4 <- ggplot() +
geom_polygon(data = udat, aes(x = long, y = lat, group = group, fill=sd), colour = "black") +
guides(fill=guide_colorbar(title="sd"))+
scale_fill_gradientn(colours=com, na.value="gray95", limits=zr) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("sd of annual total precipitation in ", year), x="", y = "", size=2.5)
ggsn::north2(p4, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath, "rain2011_hbes_wshed_sd.png"))
## Saving 7 x 5 in image
load(file=paste0(modpath, "pred_rolling_frames.RData"))
year <- 2010
# dplot <- allpred[allpred$year==year, ]
dplot <- bigdat[bigdat$year==year, ]
dim(dplot)
## [1] 1486 9
head(dplot)
## pid utmx utmy ws year mean sd low up
## 16298 22717 281225 4869955 1 2010 25.56022 1.244592 23.17574 27.88903
## 16400 22892 281225 4869975 1 2010 25.47505 1.169797 23.49283 27.55400
## 16757 23558 281195 4870045 1 2010 25.65757 1.245526 23.43243 28.07581
## 16910 23657 281175 4870055 1 2010 25.71080 1.327097 23.31293 28.63666
## 17114 23963 281215 4870085 1 2010 25.61617 1.251508 23.22060 28.00790
## 17726 24848 281135 4870165 1 2010 25.72020 1.200624 23.47980 28.14892
## annual dplot
dplot$mean <- dplot$mean*52
dim(dplot)
## [1] 1486 9
xo <- seq(from=min(dplot$utmx)-10, to = max(dplot$utmx)+10, length=200)
yo <- seq(from=min(dplot$utmy)-10, to = max(dplot$utmy)+ 10, length=200)
surf <- interp(dplot$utmx, dplot$utmy, dplot$mean, xo=xo, yo=yo)
names(surf)
## [1] "x" "y" "z"
surf <- list(x=surf$x, y=surf$y, z=surf$z)
# pnts <- data.frame(long=surf$x, lat=surf$y)
pnts <- expand.grid(surf$x, surf$y)
colnames(pnts) <- c("long", "lat")
pts1 <- pnts
coordinates(pnts)<- ~long + lat
pnts <- SpatialPoints(coords = pnts)
proj4string(pnts) <- proj4string(wsheds)
## Warning in proj4string(wsheds): CRS object has comment, which is lost in output
proj4string(wsheds)
## Warning in proj4string(wsheds): CRS object has comment, which is lost in output
## [1] "+proj=utm +zone=19 +datum=NAD27 +units=m +no_defs"
proj4string(pnts)
## Warning in proj4string(pnts): CRS object has comment, which is lost in output
## [1] "+proj=utm +zone=19 +datum=NAD27 +units=m +no_defs"
identicalCRS(wsheds, pnts)
## [1] TRUE
a1 <- sp::over(pnts, wsheds, fn = NULL) #, returnList = TRUE)
# summary(a1)
v <- pts1[!is.na(a1), ]
z <- as.vector(surf$z)
length(z)
## [1] 40000
summary(z)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1291 1355 1384 1384 1412 1475 21741
z[is.na(a1)] <- NA
z <- matrix(z, ncol=200, nrow=200)
interp1 <- data.frame(long = surf$x, z)
names(interp1)[1:length(surf$y)+1] <- surf$y
interp1 <- gather(interp1, key = lat, value =precipitation,-long,convert = TRUE)
zr <- range(interp1$precipitation, na.rm=T)
# a <- range(site$utmx)
# b <- range(site$utmy)
a <- range(interp1$long)
b <- range(interp1$lat)
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
p <- ggplot(data=interp1, aes(x=long, y=lat)) +
geom_raster(aes(x=long, y=lat, fill=precipitation)) +
geom_polygon(data = wsheds, aes(x = long, y = lat, group = group), colour = "black", fill = NA) +
guides(fill=guide_colorbar(title="Precipitation"))+
scale_fill_gradientn(colours=com, na.value="gray95", limits=zr) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("Three year rolling average precipitation in ", year), x="", y = "", size=2.5)
## Regions defined for each Polygons
ggsn::north2(p, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath, "rolling_rain2010_hbes.png"))
## Saving 7 x 5 in image
# wsy <- wsaverage[wsaverage$year==2008, ]
# wsy$rainfall.s <- wsy$rainfall.s * 52
# wsy
### Repeat for annual sd
year <- 2010
# dplot <- allpred[allpred$year==year, ]
dplot <- bigdat[bigdat$year==year, ]
dim(dplot)
## [1] 1486 9
## annual dplot
dplot$sd <- dplot$sd*52
surf <- interp(dplot$utmx, dplot$utmy, dplot$sd, xo=xo, yo=yo)
names(surf)
## [1] "x" "y" "z"
surf <- list(x=surf$x, y=surf$y, z=surf$z)
z <- as.vector(surf$z)
length(z)
## [1] 40000
summary(z)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 55.33 62.72 64.20 64.08 65.40 71.64 21741
z[is.na(a1)] <- NA
z <- matrix(z, ncol=200, nrow=200)
interp1 <- data.frame(long = surf$x, z)
names(interp1)[1:length(surf$y)+1] <- surf$y
interp1 <- gather(interp1, key = lat, value =sd,-long,convert = TRUE)
zr <- range(interp1$sd, na.rm=T)
pcolors <- brewer.pal(9,"YlOrRd")
p1 <- ggplot(data=interp1, aes(x=long, y=lat)) +
geom_raster(aes(x=long, y=lat, fill=sd)) +
geom_polygon(data = wsheds, aes(x = long, y = lat, group = group), colour = "black", fill = NA) +
guides(fill=guide_colorbar(title="sd"))+
scale_fill_gradientn(colours=pcolors, na.value="gray95", limits=zr) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("sd of rolling three year rolling average precipitation in ", year), x="", y = "", size=2.5)
## Regions defined for each Polygons
ggsn::north2(p1, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath, "rolling_rain2010_sd_hbes.png"))
## Saving 7 x 5 in image
head(wsdat)
## ws year mean sd low up
## 1 1 1999 25.77634 0.2768485 25.26947 26.28308
## 2 1 2000 25.83333 0.2775178 25.25166 26.37674
## 3 1 2001 25.62401 0.2730585 25.07344 26.19294
## 4 1 2002 25.61009 0.2668279 25.10405 26.14883
## 5 1 2003 25.66334 0.2705845 25.16294 26.18858
## 6 1 2004 25.83028 0.2879703 25.28787 26.37979
dim(wsdat)
## [1] 153 6
## Plotting trend
head(trenddat)
## pid utmx utmy ws year2 year1 mean sd low up
## 1 16 278975 4866035 9 2015 2010 0.8475066 6.438216 -11.00855 13.10672
## 2 29 278945 4866045 9 2015 2010 0.6106068 6.360907 -11.55210 12.76845
## 3 57 278965 4866055 9 2015 2010 -0.8388933 6.789405 -11.90880 14.44912
## 4 127 278035 4866075 7 2015 2010 0.4302273 6.972584 -12.00067 13.75223
## 5 159 278945 4866075 9 2015 2010 -0.4430621 5.866910 -11.01959 11.47181
## 6 182 277845 4866085 7 2015 2010 0.5254033 5.673035 -10.33862 10.75350
dim(trenddat)
## [1] 1486 10
dplot <- trenddat
dim(dplot)
## [1] 1486 10
head(dplot)
## pid utmx utmy ws year2 year1 mean sd low up
## 1 16 278975 4866035 9 2015 2010 0.8475066 6.438216 -11.00855 13.10672
## 2 29 278945 4866045 9 2015 2010 0.6106068 6.360907 -11.55210 12.76845
## 3 57 278965 4866055 9 2015 2010 -0.8388933 6.789405 -11.90880 14.44912
## 4 127 278035 4866075 7 2015 2010 0.4302273 6.972584 -12.00067 13.75223
## 5 159 278945 4866075 9 2015 2010 -0.4430621 5.866910 -11.01959 11.47181
## 6 182 277845 4866085 7 2015 2010 0.5254033 5.673035 -10.33862 10.75350
xo <- seq(from=min(dplot$utmx)-10, to = max(dplot$utmx)+10, length=200)
yo <- seq(from=min(dplot$utmy)-10, to = max(dplot$utmy)+ 10, length=200)
surf <- interp(dplot$utmx, dplot$utmy, dplot$mean, xo=xo, yo=yo)
names(surf)
## [1] "x" "y" "z"
surf <- list(x=surf$x, y=surf$y, z=surf$z)
# pnts <- data.frame(long=surf$x, lat=surf$y)
pnts <- expand.grid(surf$x, surf$y)
colnames(pnts) <- c("long", "lat")
pts1 <- pnts
coordinates(pnts)<- ~long + lat
pnts <- SpatialPoints(coords = pnts)
proj4string(pnts) <- proj4string(wsheds)
## Warning in proj4string(wsheds): CRS object has comment, which is lost in output
proj4string(wsheds)
## Warning in proj4string(wsheds): CRS object has comment, which is lost in output
## [1] "+proj=utm +zone=19 +datum=NAD27 +units=m +no_defs"
proj4string(pnts)
## Warning in proj4string(pnts): CRS object has comment, which is lost in output
## [1] "+proj=utm +zone=19 +datum=NAD27 +units=m +no_defs"
identicalCRS(wsheds, pnts)
## [1] TRUE
a1 <- sp::over(pnts, wsheds, fn = NULL) #, returnList = TRUE)
# summary(a1)
v <- pts1[!is.na(a1), ]
z <- as.vector(surf$z)
length(z)
## [1] 40000
summary(z)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -1.166 0.040 0.250 0.221 0.417 1.522 21741
z[is.na(a1)] <- NA
z <- matrix(z, ncol=200, nrow=200)
interp1 <- data.frame(long = surf$x, z)
names(interp1)[1:length(surf$y)+1] <- surf$y
interp1 <- gather(interp1, key = lat, value =trend,-long,convert = TRUE)
zr <- range(interp1$trend, na.rm=T)
a <- range(interp1$long)
b <- range(interp1$lat)
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
pcolors <- brewer.pal(8,"Greens")
p <- ggplot(data=interp1, aes(x=long, y=lat)) +
geom_raster(aes(x=long, y=lat, fill=trend)) +
geom_polygon(data = wsheds, aes(x = long, y = lat, group = group), colour = "black", fill = NA) +
guides(fill=guide_colorbar(title="Trend"))+
scale_fill_gradientn(colours=pcolors, na.value="gray95", limits=zr) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("Percentage trend between 2005 and 2015"), x="", y = "", size=2.5)
## Regions defined for each Polygons
ggsn::north2(p, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath,"trend_rolling_rain_05-15_hbes.png"))
## Saving 7 x 5 in image
## Repeated the above for sd
head(trenddat)
## pid utmx utmy ws year2 year1 mean sd low up
## 1 16 278975 4866035 9 2015 2010 0.8475066 6.438216 -11.00855 13.10672
## 2 29 278945 4866045 9 2015 2010 0.6106068 6.360907 -11.55210 12.76845
## 3 57 278965 4866055 9 2015 2010 -0.8388933 6.789405 -11.90880 14.44912
## 4 127 278035 4866075 7 2015 2010 0.4302273 6.972584 -12.00067 13.75223
## 5 159 278945 4866075 9 2015 2010 -0.4430621 5.866910 -11.01959 11.47181
## 6 182 277845 4866085 7 2015 2010 0.5254033 5.673035 -10.33862 10.75350
dim(trenddat)
## [1] 1486 10
dplot <- trenddat
dim(dplot)
## [1] 1486 10
surf <- interp(dplot$utmx, dplot$utmy, dplot$sd, xo=xo, yo=yo)
names(surf)
## [1] "x" "y" "z"
surf <- list(x=surf$x, y=surf$y, z=surf$z)
# pnts <- data.frame(long=surf$x, lat=surf$y)
pnts <- expand.grid(surf$x, surf$y)
colnames(pnts) <- c("long", "lat")
pts1 <- pnts
coordinates(pnts)<- ~long + lat
pnts <- SpatialPoints(coords = pnts)
proj4string(pnts) <- proj4string(wsheds)
## Warning in proj4string(wsheds): CRS object has comment, which is lost in output
proj4string(wsheds)
## Warning in proj4string(wsheds): CRS object has comment, which is lost in output
## [1] "+proj=utm +zone=19 +datum=NAD27 +units=m +no_defs"
proj4string(pnts)
## Warning in proj4string(pnts): CRS object has comment, which is lost in output
## [1] "+proj=utm +zone=19 +datum=NAD27 +units=m +no_defs"
identicalCRS(wsheds, pnts)
## [1] TRUE
a1 <- sp::over(pnts, wsheds, fn = NULL) #, returnList = TRUE)
# summary(a1)
v <- pts1[!is.na(a1), ]
z <- as.vector(surf$z)
length(z)
## [1] 40000
summary(z)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 5.360 6.302 6.536 6.495 6.694 7.401 21741
z[is.na(a1)] <- NA
z <- matrix(z, ncol=200, nrow=200)
interp1 <- data.frame(long = surf$x, z)
names(interp1)[1:length(surf$y)+1] <- surf$y
interp1 <- gather(interp1, key = lat, value =sd,-long,convert = TRUE)
zr <- range(interp1$sd, na.rm=T)
# a <- range(site$utmx)
# b <- range(site$utmy)
a <- range(interp1$long)
b <- range(interp1$lat)
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
p <- ggplot(data=interp1, aes(x=long, y=lat)) +
geom_raster(aes(x=long, y=lat, fill=sd)) +
geom_polygon(data = wsheds, aes(x = long, y = lat, group = group), colour = "black", fill = NA) +
guides(fill=guide_colorbar(title="sd"))+
scale_fill_gradientn(colours=com, na.value="gray95", limits=zr) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("sd of percentage trend between 2005 and 2015"), x="", y = "", size=2.5)
## Regions defined for each Polygons
ggsn::north2(p, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath,"sd_trend_rolling_rain_05-15_hbes.png"))
## Saving 7 x 5 in image
### All pred is including data from the fitting sites
year <- 2010
dplot <- allpred[allpred$year==year, ]
# dplot <- bigdat[bigdat$year==year, ]
dim(dplot)
## [1] 1510 9
## annual dplot
dplot$mean <- dplot$mean * 52
surf <- interp(dplot$utmx, dplot$utmy, dplot$mean, xo=xo, yo=yo)
names(surf)
## [1] "x" "y" "z"
surf <- list(x=surf$x, y=surf$y, z=surf$z)
z <- as.vector(surf$z)
length(z)
## [1] 40000
summary(z)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1293 1408 1467 1489 1581 1713 17365
z[is.na(a1)] <- NA
z <- matrix(z, ncol=200, nrow=200)
pcolors <- brewer.pal(8,"Greens")
interp1 <- data.frame(long = surf$x, z)
names(interp1)[1:length(surf$y)+1] <- surf$y
interp1 <- gather(interp1, key = lat, value =precipitation,-long,convert = TRUE)
zr <- range(interp1$precipitation, na.rm=T)
p1 <- ggplot(data=interp1, aes(x=long, y=lat)) +
geom_raster(aes(x=long, y=lat, fill=precipitation)) +
geom_polygon(data = wsheds, aes(x = long, y = lat, group = group), colour = "black", fill = NA) +
guides(fill=guide_colorbar(title="Precipitation"))+
scale_fill_gradientn(colours=pcolors, na.value="gray95", limits=zr) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("Rolling three year average precipitation in ", year), x="", y = "", size=2.5)
## Regions defined for each Polygons
ggsn::north2(p1, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath,"allpred_rolling_rain2010.png"))
## Saving 7 x 5 in image
# Repeat for sd
year <- 2010
dplot <- allpred[allpred$year==year, ]
# dplot <- bigdat[bigdat$year==year, ]
dim(dplot)
## [1] 1510 9
## annual dplot
dplot$sd <- dplot$sd * 52
surf <- interp(dplot$utmx, dplot$utmy, dplot$sd, xo=xo, yo=yo)
names(surf)
## [1] "x" "y" "z"
surf <- list(x=surf$x, y=surf$y, z=surf$z)
z <- as.vector(surf$z)
length(z)
## [1] 40000
summary(z)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 51.75 58.83 60.77 61.14 63.41 71.64 17365
z[is.na(a1)] <- NA
z <- matrix(z, ncol=200, nrow=200)
interp1 <- data.frame(long = surf$x, z)
names(interp1)[1:length(surf$y)+1] <- surf$y
interp1 <- gather(interp1, key = lat, value =sd,-long,convert = TRUE)
zr <- range(interp1$sd, na.rm=T)
p1 <- ggplot(data=interp1, aes(x=long, y=lat)) +
geom_raster(aes(x=long, y=lat, fill=sd)) +
geom_polygon(data = wsheds, aes(x = long, y = lat, group = group), colour = "black", fill = NA) +
guides(fill=guide_colorbar(title="sd"))+
scale_fill_gradientn(colours=com, na.value="gray95", limits=zr) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("sd of rolling three year average precipitation in ", year), x="", y = "", size=2.5)
## Regions defined for each Polygons
ggsn::north2(p1, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath,"allpred_rolling_rain2010_sd.png"))
## Saving 7 x 5 in image
year <- 2010
dplot <- wsdat[wsdat$year==year, ]
k <- c(3, 2, 1, 4, 6, 5, 8, 7, 9)
udat$rain <- NA
for(i in 1:9){
udat$rain[which(udat$id==(i-1))] <- dplot$mean[k[i]] * 52
}
zr <- range(dplot$mean) * 52
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
pcolors <- brewer.pal(8,"Greens")
p <- ggplot() +
geom_polygon(data = udat, aes(x = long, y = lat, group = group, fill=rain), colour = "black") +
guides(fill=guide_colorbar(title="Precipitation"))+
scale_fill_gradientn(colours=pcolors, na.value="gray95", limits=zr) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("Rolling three year average precipitation in ", year), x="", y = "", size=2.5)
ggsn::north2(p, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath,"rolling_rain2010_hbes_wshed.png"))
## Saving 7 x 5 in image
## Repeat for sd
year <- 2010
dplot <- wsdat[wsdat$year==year, ]
head(dplot)
## ws year mean sd low up
## 12 1 2010 25.95719 0.2699527 25.37861 26.45994
## 29 2 2010 25.49807 0.2862042 24.90818 26.00943
## 46 3 2010 25.28151 0.2853951 24.74680 25.80540
## 63 4 2010 25.91094 0.2482633 25.41130 26.37017
## 80 5 2010 26.07174 0.2630972 25.48658 26.52135
## 97 6 2010 26.52863 0.2393418 26.08768 26.93052
k <- c(3, 2, 1, 4, 6, 5, 8, 7, 9)
udat$sd <- NA
for(i in 1:9){
udat$sd[which(udat$id==(i-1))] <- dplot$sd[k[i]] * 52
}
zr <- range(dplot$sd) * 52
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
p1 <- ggplot() +
geom_polygon(data = udat, aes(x = long, y = lat, group = group, fill=sd), colour = "black") +
guides(fill=guide_colorbar(title="sd"))+
scale_fill_gradientn(colours=com, na.value="gray95", limits=zr) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= paste0("sd of rolling annual total precipitation in ", year), x="", y = "", size=2.5)
ggsn::north2(p1, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath,"rolling_rain2010_sd_hbes_wshed.png"))
## Saving 7 x 5 in image
dplot <- wstrenddat
k <- c(3, 2, 1, 4, 6, 5, 8, 7, 9)
udat$rain <- NA
for(i in 1:9){
udat$rain[which(udat$id==(i-1))] <- dplot$mean[k[i]]
}
zr <- range(dplot$mean)
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
pcolors <- rev(brewer.pal(9,"YlOrRd"))
p <- ggplot() +
geom_polygon(data = udat, aes(x = long, y = lat, group = group, fill=rain), colour = "black") +
guides(fill=guide_colorbar(title="Trend"))+
scale_fill_gradientn(colours=pcolors, na.value="gray95", limits=zr) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= "Percentage trend between 2005 and 2015", x="", y = "", size=2.5)
ggsn::north2(p, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath,"trend_wshed05-15.png"))
## Saving 7 x 5 in image
## Repeat for sd
dplot <- wstrenddat
k <- c(3, 2, 1, 4, 6, 5, 8, 7, 9)
udat$sd <- NA
for(i in 1:9){
udat$sd[which(udat$id==(i-1))] <- dplot$sd[k[i]]
}
zr <- range(dplot$sd)
com <- rev(c("firebrick4","firebrick2","green","dodgerblue2","dodgerblue4")) #colour palette
p <- ggplot() +
geom_polygon(data = udat, aes(x = long, y = lat, group = group, fill=sd), colour = "black") +
guides(fill=guide_colorbar(title="sd"))+
scale_fill_gradientn(colours=com, na.value="gray95", limits=zr) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
ylab("")+xlab("") +
theme_bw()+theme(text=element_text(family="Times")) +
#geom_point(data=site, aes(x=utmx, y=utmy-275), size=1) +
geom_text(data=centers, aes(x, y, label=ws), size=3.5) +
ggsn::scalebar(dist =1, location = "bottomright", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T,
x.min=xminn, x.max=xmaxx, y.min=yminn, y.max=ymaxx ) +
labs(title= "sd of percentage trend between 2005 and 2015", x="", y = "", size=2.5)
ggsn::north2(p, x=0.75, y=0.25, symbol=12) 
ggsave(filename=paste0(figurepath,"trend_wshed05-15_sd.png"))
## Saving 7 x 5 in image
wsdat$ws <- as.factor(wsdat$ws)
p <- ggplot(data=wsdat, aes(x=year, y=mean*52, color=ws)) +
geom_line() +
labs(x = "Year", y = "Three rolling average precipitation (mm)")
p
# All done
end.time <- proc.time()[3]
comp.time<- (end.time-start.time)/60
# comp.time<-fancy.time(comp.time)
print(comp.time)
## elapsed
## 0.6817833