Bug fixes for variable calculations

nationalparkservice · Jan 19, 2023 · 59b82bb · 59b82bb
1 parent c429b75
commit 59b82bb
Show file tree

Hide file tree

Showing 7 changed files with 15 additions and 15 deletions.
diff --git a/Code/Metric-development/vic/MAMSON.SWE.R b/Code/Metric-development/vic/MAMSON.SWE.R
@@ -97,7 +97,7 @@ for (G in 1:length(GCMs)){
 
 
   fut <- aggregate(fut_var_stars, by = by_t, FUN = function(x) mean(x)) # Doesn't work in lat/long. Must be projected. Removes units from tmax. Also aggregates to a lower resolution.
-  fut <- fut[,1:31,,]
+  fut <- fut[,1:length(future.period),,]
   fut1 <- split(fut, "time")
 
   df<-data.frame(year=future.period,mean=NA)
@@ -109,7 +109,7 @@ for (G in 1:length(GCMs)){
   DF<-rbind(DF,df)
 
 
-  mean_hist <- st_apply(hist, c("x", "y"), mean) # find mean
+  # mean_hist <- st_apply(hist, c("x", "y"), mean) # find mean
   mean_fut <- st_apply(fut, c("x", "y"), mean)
   delta <- mean_fut - mean_grid
   saveRDS(delta, file = paste(data.dir,paste(var,gcm,rcp,sep="_"),sep="/"))

diff --git a/Code/Metric-development/vic/SWE.precip.R b/Code/Metric-development/vic/SWE.precip.R
@@ -48,8 +48,8 @@ DF = data.frame()
   s <- aggregate(var_swe, by = by_t, FUN = function(x) max(x)) #Don't need to divide by #yrs b/c by year
   s1 <- split(s, "time")
 
-  ratio1 <- s1[1:37,,]/p1[1:37,,]
-  ratio  <- s[,1:37,,]/p[,1:37,,]
+  ratio1 <- s1[1:length(daymet.period),,]/p1[1:length(daymet.period),,]
+  ratio  <- s[,1:length(daymet.period),,]/p[,1:length(daymet.period),,]
 
 
   df<-data.frame(year=daymet.period,mean=NA)
@@ -117,8 +117,8 @@ for (G in 1:length(GCMs)){
   s <- aggregate(var_swe, by = by_t, FUN = function(x) max(x)) #Don't need to divide by #yrs b/c by year
   s1 <- split(s, "time")
 
-  ratio1 <- s1[1:31,,]/p1[1:31,,]
-  ratio  <- s[,1:31,,]/p[,1:31,,]
+  ratio1 <- s1[1:length(future.period),,]/p1[1:length(future.period),,]
+  ratio  <- s[,1:length(future.period),,]/p[,1:length(future.period),,]
   mean_fut<-st_apply(ratio, c("x", "y"), mean)
 
   df<-data.frame(year=future.period,mean=NA)

diff --git a/Code/Metric-development/vic/max.SWE.R b/Code/Metric-development/vic/max.SWE.R
@@ -82,7 +82,7 @@ for (G in 1:length(GCMs)){
 
   by_t = "1 year"
   hist <- aggregate(hist_var_stars, by = by_t, FUN = function(x) max(x)) #Don't need to divide by #yrs b/c by year
-  hist <- hist[,2:51,,]
+  hist <- hist[,1:length(historical.period),,]
   # hist1 <- split(hist, "time")
   # 
   # 
@@ -96,7 +96,7 @@ for (G in 1:length(GCMs)){
 
 
   fut <- aggregate(fut_var_stars, by = by_t, FUN = function(x) max(x)) # Doesn't work in lat/long. Must be projected. Removes units from tmax. Also aggregates to a lower resolution.
-  fut <- fut[,2:32,,]
+  fut <- fut[,1:length(future.period),,]
   fut1 <- split(fut, "time")
 
   df<-data.frame(year=future.period,mean=NA)

diff --git a/Code/Metric-development/vic/soil.temp.R b/Code/Metric-development/vic/soil.temp.R
@@ -81,7 +81,7 @@ for (G in 1:length(GCMs)){
 
   by_t = "1 year"
   hist <- aggregate(hist_var_stars, by = by_t, FUN = function(x) mean(x)) #Don't need to divide by #yrs b/c by year
-  hist <- hist[,2:51,,]
+  hist <- hist[,1:length(historical.period),,]
   # hist1 <- split(hist, "time")
   # 
   # 
@@ -95,7 +95,7 @@ for (G in 1:length(GCMs)){
 
 
   fut <- aggregate(fut_var_stars, by = by_t, FUN = function(x) mean(x)) # Doesn't work in lat/long. Must be projected. Removes units from tmax. Also aggregates to a lower resolution.
-  fut <- fut[,2:32,,]
+  fut <- fut[,1:length(future.period),,]
   fut1 <- split(fut, "time")
 
   df<-data.frame(year=future.period,mean=NA)

diff --git a/Code/Metric-development/vic/water.balance.R b/Code/Metric-development/vic/water.balance.R
@@ -83,7 +83,7 @@ for (G in 1:length(GCMs)){
 
   by_t = "1 year"
   hist <- aggregate(hist_var_stars, by = by_t, FUN = function(x) sum(x)) #Don't need to divide by #yrs b/c by year
-  hist <- hist[,2:51,,]
+  hist <- hist[,1:length(historical.period),,]
   # hist1 <- split(hist, "time")
   # 
   # 
@@ -97,7 +97,7 @@ for (G in 1:length(GCMs)){
 
 
   fut <- aggregate(fut_var_stars, by = by_t, FUN = function(x) sum(x)) # Doesn't work in lat/long. Must be projected. Removes units from tmax. Also aggregates to a lower resolution.
-  fut <- fut[,2:32,,]
+  fut <- fut[,1:length(future.period),,]
   fut1 <- split(fut, "time")
 
   df<-data.frame(year=future.period,mean=NA)

diff --git a/Code/example-scripts/metric_test_code.R b/Code/example-scripts/metric_test_code.R
@@ -8,7 +8,7 @@ test <- readRDS("C:/Users/achildress/Documents/wrst_temp/Data/WRST_simple/inmcm4
 
 
 ggplot() + # Resolution is course
-  geom_stars(data = hist1[7,], alpha = 0.8) + 
+  geom_stars(data = mean_grid, alpha = 0.8) + 
   # facet_wrap("time") +
   # scale_fill_viridis() + 
   #coord_equal() + 

diff --git a/Variable_calculations.Rmd b/Variable_calculations.Rmd
@@ -89,11 +89,11 @@ source(here::here("Code", "Metric-development", "vic","soil.temp.R"),echo = FALS
 ```{r Create daily variables, results=hide, eval=TRUE}
 gc()
 memory.limit(60000)
-source(here::here("Code", "Metric-development", "met","met-daily-stars_by_year.R"),echo = FALSE) 
+source(here::here("Code", "Metric-development", "met","met-daily-stars_by_year.R"),echo = FALSE) #Takes ~1 hrs to run
 source(here::here("Code", "Metric-development", "vic","SWE-runoff-Pr99.R"),echo = FALSE)
 source(here::here("Code", "Metric-development", "met","Pr99v2.R"),echo = FALSE)
 
-source(here::here("Code", "Metric-development","daily-all-geographies.R"),echo = FALSE,eval=FALSE)
+# source(here::here("Code", "Metric-development","daily-all-geographies.R"),echo = FALSE,eval=FALSE)
 
 ```