forked from nationalparkservice/WRST-climate-futures
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathplot-creation.Rmd
136 lines (102 loc) · 5.76 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(terra);library(raster)
SiteID <- "KOVA" #*UPDATE*
base.dir = paste0("C:/Users/brobb/OneDrive - DOI/Projects/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)
write.csv(DF, here("DF.csv"))
# 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) # Excludes historical rows
df.hist <- DF %>% filter(Period == "Historical") # Includes historical rows
# df.fut = subset(DF, Period=="Future" & GCM %in% GCMs) # Excludes historical rows
df.fut <- DF %>% filter(Period == "Future" | (GCM == "Daymet" & Period == "Historical")) # Includes historical rows
delta = df.hist[,1:4]
delta[,5:8] = df.fut[,5:8] - df.hist[,5:8]
delta <- delta[-c(9:12),] # drops historical rows that subtracted to '0'
# 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) # not running for AK for now - takes very long to run and can get similar info from northern climate reports
```
``` {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)
```
```{r metrics table, include=FALSE, results=hide, eval=TRUE}
source(here::here("Code", "Plots", "climate-futures-table.R"),echo = FALSE)
```