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)
png(filename=paste0(figurepath, "elevation_hubbardbrook.png"))
ggsn::north2(p9, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
# [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")
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(outputpath, "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
## Min. : 1 Min. :277345 Min. :4866025 Min. : 1.000
## 1st Qu.: 7906 1st Qu.:278415 1st Qu.:4866545 1st Qu.: 4.000
## Median :17341 Median :279265 Median :4867005 Median : 7.000
## Mean :17414 Mean :279571 Mean :4868184 Mean : 6.136
## 3rd Qu.:26795 3rd Qu.:280785 3rd Qu.:4870315 3rd Qu.: 8.000
## Max. :35487 Max. :283225 Max. :4871045 Max. :10.000
##
## year mean sd low up
## 2015 :1510 Min. :24.82 Min. :1.394 Min. :19.89 Min. :28.36
## 1997 : 0 1st Qu.:25.91 1st Qu.:2.017 1st Qu.:22.01 1st Qu.:29.98
## 1998 : 0 Median :27.12 Median :2.096 Median :23.03 Median :31.04
## 1999 : 0 Mean :26.84 Mean :2.091 Mean :22.87 Mean :30.89
## 2000 : 0 3rd Qu.:27.67 3rd Qu.:2.172 3rd Qu.:23.69 3rd Qu.:31.75
## 2001 : 0 Max. :28.71 Max. :2.511 Max. :25.34 Max. :33.39
## (Other): 0
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)
png(file=paste0(figurepath, "rain2015_total.png"))
ggsn::north2(p1, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
### Repeat for annual sd
year <- 2015
dplot <- predandfit[predandfit$year==year, ]
dim(dplot)
## [1] 1510 9
## annual dplot
summary(dplot)
## pid utmx utmy ws
## Min. : 1 Min. :277345 Min. :4866025 Min. : 1.000
## 1st Qu.: 7906 1st Qu.:278415 1st Qu.:4866545 1st Qu.: 4.000
## Median :17341 Median :279265 Median :4867005 Median : 7.000
## Mean :17414 Mean :279571 Mean :4868184 Mean : 6.136
## 3rd Qu.:26795 3rd Qu.:280785 3rd Qu.:4870315 3rd Qu.: 8.000
## Max. :35487 Max. :283225 Max. :4871045 Max. :10.000
##
## year mean sd low up
## 2015 :1510 Min. :24.82 Min. :1.394 Min. :19.89 Min. :28.36
## 1997 : 0 1st Qu.:25.91 1st Qu.:2.017 1st Qu.:22.01 1st Qu.:29.98
## 1998 : 0 Median :27.12 Median :2.096 Median :23.03 Median :31.04
## 1999 : 0 Mean :26.84 Mean :2.091 Mean :22.87 Mean :30.89
## 2000 : 0 3rd Qu.:27.67 3rd Qu.:2.172 3rd Qu.:23.69 3rd Qu.:31.75
## 2001 : 0 Max. :28.71 Max. :2.511 Max. :25.34 Max. :33.39
## (Other): 0
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)
png(filename=paste0(figurepath, "rain2015_total_sd.png"))
ggsn::north2(p, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
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)
png(filename=paste0(figurepath, "rain2001_hbes_wshed.png"))
ggsn::north2(p1, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
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)
png(filename=paste0(figurepath, "rain2011_hbes_wshed.png"))
ggsn::north2(p11, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
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)
png(filename=paste0(figurepath, "rain2001_hbes_wshed_sd.png"))
ggsn::north2(p3, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
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)
png(filename=paste0(figurepath, "rain2011_hbes_wshed_sd.png"))
ggsn::north2(p4, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
load(file=paste0(outputpath, "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)
png(filename=paste0(figurepath, "rolling_rain2010_hbes.png"))
ggsn::north2(p, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
# 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)
png(filename=paste0(figurepath, "rolling_rain2010_sd_hbes.png"))
ggsn::north2(p1, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
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)
png(filename=paste0(figurepath,"trend_rolling_rain_05-15_hbes.png"))
ggsn::north2(p, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
## 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)
png(filename=paste0(figurepath,"sd_trend_rolling_rain_05-15_hbes.png"))
ggsn::north2(p, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
### 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)
png(filename=paste0(figurepath,"allpred_rolling_rain2010.png"))
ggsn::north2(p1, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
# 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)
png(filename=paste0(figurepath,"allpred_rolling_rain2010_sd.png"))
ggsn::north2(p1, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
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)
png(filename=paste0(figurepath,"rolling_rain2010_hbes_wshed.png"))
ggsn::north2(p, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
## 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)
png(filename=paste0(figurepath,"rolling_rain2010_sd_hbes_wshed.png"))
ggsn::north2(p1, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
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)
png(filename=paste0(figurepath,"trend_wshed05-15.png"))
ggsn::north2(p, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
## 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)
png(filename=paste0(figurepath,"trend_wshed05-15_sd.png"))
ggsn::north2(p, x=0.75, y=0.25, symbol=12)
dev.off()
## png
## 2
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.9083