-
Notifications
You must be signed in to change notification settings - Fork 2
/
plot-creation.Rmd
135 lines (102 loc) · 5.44 KB
/
plot-creation.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
---
title: "Running plotting scripts"
author: "Amber Runyon"
date: "8/13/2021"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r set directories, include=FALSE, results=hide, eval=TRUE}
rm(list = ls())
library(stars);library(dplyr);library(ggplot2);library(ggthemes);library(viridis);library(here);library(ggrepel);library(rlang);library(units); library(tidyr); library(lemon);library(ggpubr);library(gridExtra);library(grid); library(gtable); library(lubridate);library(raster)
base.dir = "C:/Users/arunyon/DOI/Climate Change Scenario Planning - Data"
data.dir = paste0(base.dir)
plot.dir = "./Data/figures"
boundary.dir <- paste0(base.dir,"/Boundary_shapefiles/")
wrst <- st_read(paste0(boundary.dir, "wrst_simple.shp")) # Wrangell Mountains
wrst <- st_transform(wrst, 3338)
# insert topo
topo <- raster::stack(paste0(boundary.dir, 'HYP_HR_SR_W/HYP_HR_SR_W.tif')) # read in as stack so can see RBG layers
ext <- extent(-147, -139, 59.3, 63) # extent defined by lat/long
ak <- crop(topo, ext)
# plotRGB(ak)
ak2 <- projectRaster(ak, crs = CRS('+init=EPSG:3338')) # Alaska Albers
ak_df <- as.data.frame(ak2, xy = TRUE) # this step is important to get it to plot in ggplot
GCMs <- c("inmcm4.rcp85","ACCESS1-3.rcp45","CanESM2.rcp85")
CFs <- c("Climate Future 1", "Climate Future 2", "Climate Future 3")
cols <- c("#12045C","#ffcccc","#E10720")
CF_GCM <- data.frame(CF=CFs,GCM=GCMs,CF_col=cols)
```
```{r area map, include=FALSE,results=hide, eval=FALSE}
# source(here::here("Code", "Plots", "area-map.R"),echo = FALSE)
# holding off on script - using ArcGIS plot for now
```
```{r scatterplot, echo=FALSE, message=FALSE, warning=FALSE,eval=FALSE}
dir = paste0(data.dir,"/WRST_simple/")
df = read.csv(paste0(dir,"/Monthly_met.csv"))
tmax = cbind(df[,1:3],df[grep("tmax", colnames(df))]);names(tmax)[4:7] = gsub("\\..*","",names(tmax[,4:7]))
tmax.long = gather(tmax, season, var, -c("CF","GCM","Period")); names(tmax.long)[5]="tmax"
tmin = cbind(df[,1:3],df[grep("tmin", colnames(df))]);names(tmin)[4:7] = gsub("\\..*","",names(tmin[,4:7]))
tmin.long = gather(tmin, season, var, -c("CF","GCM","Period")); names(tmin.long)[5]="tmin"
DF = merge(tmax.long,tmin.long,by=c("CF","GCM","Period","season"));rm(tmax,tmax.long,tmin,tmin.long)
tmean = cbind(df[,1:3],df[grep("tmean", colnames(df))]);names(tmean)[4:7] = gsub("\\..*","",names(tmean[,4:7]))
tmean.long = gather(tmean, season, var, -c("CF","GCM","Period")); names(tmean.long)[5]="tmean"
DF = merge(DF, tmean.long,by=c("CF","GCM","Period","season"));rm(tmean,tmean.long)
pcp = cbind(df[,1:3],df[grep("precip", colnames(df))]);names(pcp)[4:7] = gsub("\\..*","",names(pcp[,4:7]))
pcp.long = gather(pcp, season, var, -c("CF","GCM","Period")); names(pcp.long)[5]="pcp"
DF = merge(DF, pcp.long,by=c("CF","GCM","Period","season"));rm(pcp,pcp.long)
# factors
DF$CF = factor(DF$CF, levels=CFs)
DF$GCM = factor(DF$GCM, levels=GCMs)
DF$season <- factor(DF$season, levels=c("DJF","MAM","JJA","SON"))
# create delta df
df.hist = subset(DF, Period=="Historical" & GCM %in% GCMs)
df.fut = subset(DF, Period=="Future" & GCM %in% GCMs)
delta = df.hist[,1:4]
delta[,5:8] = df.fut[,5:8] - df.hist[,5:8]
# source(here::here("Code", "Plots", "scatterplot.R"),echo = FALSE)
```
```{r seasonal plots, include=FALSE,results=hide, eval=TRUE}
source(here::here("Code", "Plots", "map_monthly_dotplots_tmax.R"),echo = FALSE)
source(here::here("Code", "Plots", "map_monthly_dotplots_tmin.R"),echo = FALSE)
source(here::here("Code", "Plots", "map_monthly_dotplots_tmean.R"),echo = FALSE)
source(here::here("Code", "Plots", "map_monthly_dotplots_pcp.R"),echo = FALSE)
```
```{r annual map-timeseries plots, include=FALSE,results=hide, eval=TRUE}
dir = paste0(data.dir,"/WRST_simple/")
source(here::here("Code", "Plots", "maps_ts_plots_tmean.R"),echo = FALSE)
source(here::here("Code", "Plots", "maps_ts_plots_pcp.R"),echo = FALSE)
source(here::here("Code", "Plots", "maps_ts_plots_water.balance.R"),echo = FALSE)
source(here::here("Code", "Plots", "maps_ts_plots_SWE.precip.R"),echo = FALSE)
source(here::here("Code", "Plots", "maps_ts_plots_max.SWE.R"),echo = FALSE)
source(here::here("Code", "Plots", "maps_ts_plots_MAMSON.SWE.R"),echo = FALSE)
source(here::here("Code", "Plots", "maps_ts_plots_soil.temp.R"),echo = FALSE)
source(here::here("Code", "Plots", "maps_ts_plots_Pr99v2.R"),echo = FALSE)
```
``` {r timeseries stack, include=FALSE,results=hide, eval=TRUE}
areas <- c("wrst_simple", 'wrst_mtns', 'copper_plateau', 'coastal_mtns', 'eastern_coast')
area_names = c("WRST", "WRST mountains", "Copper plateau", "Coastal mountains", "eastern coast")
for (a in 1:length(areas)) {
area=areas[a]
area_name = area_names[a]
dir = paste0(data.dir,"/", area, "/")
DF = read.csv(paste0(dir,"/Daily_met.csv"))
source(here::here("Code", "Plots", "daily_ts_stack.R"),echo = FALSE)
}
```
```{r daily plots, include=FALSE,results=hide, eval=TRUE}
areas <- c('wrst_simple', 'wrst_mtns', 'copper_plateau', 'coastal_mtns', 'eastern_coast')
area_names = c("WRST", "WRST mountains", "Copper plateau", "Coastal mountains", "eastern coast")
for (a in 1:length(areas)) {
area=areas[a]
area_name = area_names[a]
dir = paste0(data.dir,"/", area, "/")
runoff = read.csv(paste0(dir,"/runoff_DAY.csv"))
SWE = read.csv(paste0(dir,"/SWE_DAY.csv"))
source(here::here("Code", "Plots", "SWE-runoff-plots.R"),echo = FALSE)
}
```
```{r SWE.precip.maps, include=FALSE, results=hide, eval=TRUE}
source(here::here("Code", "Plots", "SWE.precip.maps.R"),echo = FALSE)
```