forked from nationalparkservice/WRST-climate-futures
-
Notifications
You must be signed in to change notification settings - Fork 0
/
plot-creation.Rmd
127 lines (95 loc) · 5.1 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
---
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)
SiteID <- "CAKR"
base.dir = paste0("C:/Users/arunyon/OneDrive - DOI/AKR-CFs/",SiteID) #*UPDATE*
data.dir = paste0(base.dir,"/Data")
plot.dir = paste0(base.dir,"/Figures")
# plot.dir = "./Data/figures"
nps_boundary <- st_read('./Data/nps_boundary/nps_boundary.shp')
park <- filter(nps_boundary, UNIT_CODE == SiteID)
shp <- st_transform(park, 3338)
park_ext <- st_bbox(filter(nps_boundary, UNIT_CODE == SiteID))
# #Unique
# huc <- st_read("C:/Users/arunyon/OneDrive - DOI/Documents/GIS/AK_HUC8/Kachemak_Tuxedni.shp")
# shp <- st_transform(huc, 3338)
# insert topo
topo <- stack('./Data/nps_boundary/HYP_HR_SR_W/HYP_HR_SR_W.tif') # read in as stack so can see RBG layers
# shp <- st_transform(huc, st_crs(topo))
ext <- extent(park_ext$xmin, park_ext$xmax, park_ext$ymin, park_ext$ymax) # 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("MRI-CGCM3.rcp45","CCSM4.rcp85")
CFs <- c("Climate Future 1", "Climate Future 2")
cols <- c("#6EB2D4","#CA0020")
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)
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)
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}
DF = read.csv(paste0(data.dir,"/Daily_met.csv"))
source(here::here("Code", "Plots", "daily_ts_stack.R"),echo = FALSE)
```
```{r daily plots, include=FALSE,results=hide, eval=TRUE}
runoff = read.csv(paste0(data.dir,"/runoff_DAY.csv"))
SWE = read.csv(paste0(data.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)
```