upload codebase

Author: joshhjacobson

.gitignore

@@ -0,0 +1,4 @@
+*.DS_Store
+*.nc
+*.RData
+data

README.md

@@ -1,2 +1,6 @@
-# FTE
-The Fraction of Threshold Exceedance (FTE) metric for verification of spatial structure in ensembles of forecast fields
+# Fraction of Threshold Exceedance
+
+This project is licensed under the MIT License - see the [LICENSE.md](LICENSE.md) file for details.
+
+## Suggested citation
+

beta_params_range.R

@@ -0,0 +1,44 @@
+
+## Demonstrate dependence on domain size by analyzing range of ratios for a1 = 1,3 and tau fixed at 1
+
+library(tidyverse)
+library(reshape2)
+
+## Beta param plot for a1 = 2
+cont_fit_tab <- read.table('./data/cont_fit_tab.RData')
+
+df <- cont_fit_tab %>% 
+  filter((tau==1 & s1==1) |  (tau==1 & s1==3)) %>% 
+  dplyr::select(-tau) %>% 
+  melt(id.vars=c("s1", "ratio"))
+
+range_labs <- c(
+  `1` = "a0 = 1",
+  `3` = "a0 = 3"
+)
+
+png("beta_params_range.png", units="in", height=3.3, width=6.4, res=200, pointsize=9)
+
+ggplot(data=df, aes(x=ratio, y=value, color=variable, linetype=variable)) +
+  facet_grid(~s1, labeller = as_labeller(range_labs)) +
+  geom_hline(yintercept=1, linetype=3, size=0.3) +
+  geom_line(size=0.3) + 
+  geom_point(size=0.8) +
+  scale_colour_manual(values=c(a="darkred", b="skyblue3")) +
+  scale_linetype_manual(breaks=c("a","b"), values=c(2, 1)) +
+  scale_y_continuous(breaks = c(0.75, 1.0, 1.25)) +
+  labs(x="Ratio", y="Parameter") +
+  theme_bw() +
+  theme(legend.title = element_blank(),
+        strip.background = element_blank(),
+        text = element_text(color="black"),
+        strip.text= element_text(size=12),
+        axis.text = element_text(size=9, color="black"),
+        legend.text = element_text(size=10),
+        panel.grid.minor = element_blank(),
+        panel.grid.major = element_blank(),
+        aspect.ratio = 1/1,
+        plot.margin = unit(c(0,0,0,0.1), "cm"))
+
+dev.off()
+

beta_params_tau.R

@@ -0,0 +1,46 @@
+
+## Demonstrate dependence on tau by analyzing range of taus at low, even, and high ratios (a1 = 2)
+
+library(tidyverse)
+library(reshape2)
+
+cont_fit_tab <- read.table('./data/cont_fit_tab.RData')
+
+df <- cont_fit_tab %>% 
+  filter(s1==2, ratio==0.5 | ratio==0.9 | ratio==1.0 | ratio == 1.1 | ratio==1.5) %>% 
+  select(-s1) %>% 
+  melt(id.vars=c("tau", "ratio"))
+
+ratio_labs <- c(
+  `0.5` = "ratio = 0.5",
+  `0.9` = "ratio = 0.9",
+  `1` = "ratio = 1.0",
+  `1.1` = "ratio = 1.1",
+  `1.5` = "ratio = 1.5"
+)
+
+png("beta_params_tau.png", units="in", height=2.2, width=8.4, res=200, pointsize=10)
+
+ggplot(data=df, aes(x=tau, y=value, color=variable, linetype=variable)) +
+  facet_grid(~ratio, labeller = as_labeller(ratio_labs)) +
+  geom_hline(yintercept=1, linetype=3, size=0.2) +
+  geom_line(size=0.25) + 
+  geom_point(size=0.6) +
+  scale_colour_manual(values=c(a="darkred", b="skyblue3")) +
+  scale_linetype_manual(breaks=c("a","b"), values=c(2, 1)) +
+  scale_y_continuous(breaks = c(0.5, 1.0, 1.5)) +
+  labs(x="Threshold", y="Parameter") +
+  theme_bw() +
+  theme(legend.title = element_blank(),
+        strip.background = element_blank(),
+        text = element_text(color="black"),
+        strip.text= element_text(size=12),
+        axis.text = element_text(size=9, color="black"),
+        legend.text = element_text(size=10),
+        panel.grid.minor = element_blank(),
+        panel.grid.major = element_blank(),
+        aspect.ratio = 1/1,
+        plot.margin = unit(c(0,0,0,0.1), "cm"))
+
+dev.off()
+

fte_downscaled.R

@@ -0,0 +1,113 @@
+
+## Produce FTE histograms for downscaled GEFS data
+
+library(lubridate)
+library(dplyr)
+library(tidyr)
+library(reshape2)
+library(ggplot2)
+library(fitdistrplus)
+library(latex2exp)
+
+
+# Threshold exceedence ranking --------------------------------------------
+
+## import fields_dat dim(lat x lon x time x mem)
+load("./data/gefs_downscaled_fields.RData")
+
+## compute the mean threshold exceedence of fields_df at array of thresholds
+## and return analysis ranks
+exceed_ranks <- function(dat_arr, tau){
+  # dat_arr: 3d array (lon x lat x member)
+  # tau: vector of thresholds
+  ranks <- array(dim = length(tau))
+  for (i in 1:length(tau)) {
+    m <- apply(dat_arr, 3, function(field) mean(as.vector(field) > tau[i]))
+    if(length(unique(m)) != 1) { 
+      ranks[i] <- rank(m, ties.method = "random")[1]
+    } else {
+      # exclude exact ties
+      ranks[i] <- NA
+    }
+  }
+  return(ranks)
+}
+
+disagg_rank <- function(r) {
+  return(runif(1, r-1/24, r+1/24))
+}
+
+## iterate over time, compute ranks at different thresholds
+set.seed(10)
+tau <- c(5,10,20)
+ranks_df <- data.frame(t(apply(field_dat, 3, exceed_ranks, tau=tau))) # (day x tau)
+
+
+# FTE Histograms ----------------------------------------------------------
+
+## stratify days by month and build histograms for each month-threshold pair
+names(ranks_df) <- paste('tau', tau, sep='')
+
+m <- c(1, 4, 7, 10) # Jan, Apr, Jul, Oct
+dates <- seq.Date(as.Date('2002-01-02'), as.Date('2015-12-30'), by='day')
+date_idx <- (month(dates) %in% m)
+
+## fit beta parameters to density fte hists
+set.seed(20)
+down_fit_tab <- ranks_df %>%
+  mutate(month = month(dates[date_idx])) %>%
+  melt(id.vars='month', variable.name='tau', value.name='rank') %>%
+  mutate(rank = sapply((rank-0.5)/12, disagg_rank)) %>%
+  group_by(tau, month) %>%
+  drop_na() %>%
+  summarise(params=paste(fitdist(rank,'beta')$estimate, collapse=" ")) %>%
+  separate(params, c('a', 'b'), sep=" ") %>%
+  mutate(a=round(as.numeric(a), 3), b=round(as.numeric(b),3)) %>%
+  unite(params, a:b, sep = ", ")
+
+## build dataframe with months and tau as factors to facet over
+month_labs <- rep("", 12)
+month_labs[c(1,4,7,10)] <- c("January", "April", "July", "October")
+tau_labs <- c(tau5=TeX(paste("$\\tau = $", "5 mm")), 
+              tau10=TeX(paste("$\\tau = $", "10 mm")), 
+              tau20=TeX(paste("$\\tau = $", "20 mm")))
+
+ranks_df <- ranks_df %>%
+  mutate(month = month(dates[date_idx])) %>%
+  melt(id.vars='month', variable.name='tau', value.name='rank') %>%
+  mutate(month = as.factor(month_labs[month])) %>%
+  mutate(rank = (rank-0.5)/12)
+
+ranks_df$month <- factor(ranks_df$month, 
+                     levels = c("January", "April", "July", "October"),
+                     labels = c("bold(January)", "bold(April)", "bold(July)", "bold(October)"))
+levels(ranks_df$tau) <- tau_labs
+
+
+png("fte_downscaled.png", units="in", height=6.2, width=8, res=200, pointsize=10)
+
+ranks_df %>%
+  ggplot(aes(x=rank)) +
+  geom_hline(yintercept=1, linetype=3, size=0.3, color="grey") +
+  geom_histogram(aes(y=..density..), bins=12, fill="black", color="white") +
+  scale_y_continuous(breaks=seq(0,2)) +
+  facet_grid(rows=vars(tau), cols=vars(month), labeller=label_parsed) +
+  annotate("text", x=0.48, y=2.5, size=4, label=down_fit_tab$params) +
+  labs(y="", x="") +
+  theme_bw() +
+  theme(legend.title = element_blank(),
+        strip.background = element_blank(),
+        text = element_text(color="black"),
+        strip.text = element_text(size=12),
+        axis.text = element_text(size=9, color="black"),
+        axis.text.x = element_blank(),
+        axis.ticks.x = element_blank(),
+        panel.grid.minor = element_blank(),
+        panel.grid.major.x = element_blank(),
+        panel.grid.major = element_blank(),
+        aspect.ratio = 1/1,
+        plot.margin = unit(c(0,0,0,0), "cm"))
+
+dev.off()
+
+

gefs_downscaled.R

@@ -0,0 +1,85 @@
+
+## Produce FTE histograms for downscaled GEFS data
+
+library(ncdf4)
+library(lubridate)
+library(dplyr)
+library(reshape2)
+library(ggplot2)
+
+
+# Analyses ----------------------------------------------------------------
+
+## NOTE: dates range from 2002010200 to 2015123000
+
+## grab analyses and time data
+ncin <- nc_open('./data/gefs/original/refcstv2_precip_ccpav3_subset_066_to_072.nc')
+apcp_anal <- ncvar_get(ncin, 'apcp_anal')
+
+## subset data to same region
+lons_anal <- ncvar_get(ncin, 'lons_anal')
+lats_anal <- ncvar_get(ncin, 'lats_anal')
+lon_idx <- which(lons_anal[,1] >= -91 &  lons_anal[,1] <= -81)
+lat_idx <- which(lats_anal[1,] >= 30 &  lats_anal[1,] <= 40)
+apcp_anal <- apcp_anal[lon_idx, lat_idx, ]
+rm(lons_anal, lats_anal, lon_idx, lat_idx)
+
+## create an array to store all each analysis-ensemble pair dim(lon, lat, time, obs/ens_n)
+field_dat <- array(dim = c(dim(apcp_anal), 12))
+field_dat[,,,1] <- apcp_anal
+
+nc_close(ncin)
+rm(ncin, apcp_anal)
+
+
+# Ensembles (downscaled) --------------------------------------------------
+
+## NOTES: Each array has space for 31 days, but 
+#   - January 2002 starts from 01/02 (start timeframe one day earlier, and trim)
+#   - April always has 30 days
+
+m <- c(1, 4, 7, 10) # Jan, Apr, Jul, Oct
+dates <- seq.Date(as.Date('2002-01-02'), as.Date('2015-12-30'), by='day')
+
+files <- list.files(path='./data', pattern='*.nc', full.names=TRUE)
+lapply(seq_along(files), function(i) {
+  ncin <- nc_open(files[i])
+  fcsts <- ncvar_get(ncin, 'downscaled')
+  
+  lons_fcst <- ncvar_get(ncin, 'Longitude')
+  lats_fcst <- ncvar_get(ncin, 'Latitude')
+  lon_idx <- which(lons_fcst[,1] >= -91 &  lons_fcst[,1] <= -81)
+  lat_idx <- which(lats_fcst[1,] >= 30 &  lats_fcst[1,] <= 40)
+  
+  # april
+  if(m[i]==4){
+    fcsts <- fcsts[,,,1:30,]
+  }
+  
+  d <- dim(fcsts)
+  dim(fcsts) <- c(d[1], d[2], d[3], d[4]*d[5])
+  
+  # january
+  if(m[i]==1){
+    # nil <- seq(31, dim(fcsts)[4], by=31)
+    fcsts <- fcsts[,,, -31]
+  }
+  
+  date_idx <- which(month(dates)==m[i])
+  for(j in 1:11) {
+    jj <- j+1
+    # print(dim(field_dat[,, date_idx, jj]))
+    # print(dim(fcsts[,,j,]))
+    field_dat[,, date_idx, jj] <<- fcsts[lon_idx, lat_idx, j,]
+  }
+  
+  nc_close(ncin)
+  rm(ncin)
+})
+rm(files)
+
+date_idx <- (month(dates) %in% m)
+field_dat <- field_dat[,, date_idx,] 
+save(field_dat, file="gefs_downscaled_fields.RData")
+
+

gefs_upscaled.R

@@ -0,0 +1,36 @@
+
+## Produce FTE histograms for upscaled GEFS reanalysis data
+
+library(ncdf4)
+library(lubridate)
+library(dplyr)
+library(reshape2)
+library(ggplot2)
+
+
+# Fcst + Analysis ---------------------------------------------------------
+
+## NOTE: dates range from 2002010200 to 2015123000
+
+## grab field and time data
+ncin <- nc_open('./data/refcstv2_precip_ccpav3_subset_066_to_072.nc')
+fcst_ens <- ncvar_get(ncin, 'apcp_fcst_ens')
+anal_upsc <- ncvar_get(ncin, 'apcp_anal_upsc')
+init_anal <- ncvar_get(ncin, 'yyyymmddhh_init')
+
+## get indices of data in upscaled region
+lons_fcst <- ncvar_get(ncin, 'lons_fcst')
+lats_fcst <- ncvar_get(ncin, 'lats_fcst')
+lon_idx <- which(lons_fcst[,1] >= -91 &  lons_fcst[,1] <= -81)
+lat_idx <- which(lats_fcst[1,] >= 30 &  lats_fcst[1,] <= 40)
+
+## create an array to store each analysis-ensemble pair dim(lon, lat, time, obs/ens_n)
+field_dat <- array(dim = c(length(lon_idx), length(lat_idx), length(init_anal), 12))
+
+## populate and trim to same regions
+field_dat[,,,1] <- anal_upsc[lon_idx, lat_idx, ]
+field_dat[,,,2:12] <- fcst_ens[lon_idx, lat_idx, , ]
+
+nc_close(ncin)
+rm(ncin, anal_upsc, fcst_ens, init_anal, lats_fcst, lons_fcst, lat_idx, lon_idx)
+