library(tidyverse)
library(lubridate)
library(scales)
library(knitr)

noaa <- src_postgres(host="localhost", dbname="noaa")

snow <- tbl(noaa, build_sql(
   "WITH wdoy_data AS (
         SELECT dte, dte - interval '120 days' as wdte,
            tmin_c, tmax_c, (tmin_c+tmax_c)/2.0 AS tavg_c,
            prcp_mm, snow_mm, snwd_mm
         FROM ghcnd_pivot
         WHERE station_name = 'FAIRBANKS INTL AP'
         AND dte > '1950-09-01')
   SELECT dte, date_part('year', wdte) AS wyear, date_part('doy', wdte) AS wdoy,
         to_char(dte, 'Mon DD') AS mmdd,
         tmin_c, tmax_c, tavg_c, prcp_mm, snow_mm, snwd_mm
   FROM wdoy_data")) %>%
   mutate(wyear=as.integer(wyear),
            wdoy=as.integer(wdoy),
            snwd_mm=as.integer(snwd_mm)) %>%
   select(dte, wyear, wdoy, mmdd,
            tmin_c, tmax_c, tavg_c, prcp_mm, snow_mm, snwd_mm) %>% collect()

write_csv(snow, "pafa_data_with_wyear_post_1950.csv")
save(snow, file="pafa_data_with_wyear_post_1950.rdata")

cum_snow <- snow %>%
   mutate(snow_na=ifelse(is.na(snow_mm),1,0),
         snow_mm=ifelse(is.na(snow_mm),0,snow_mm)) %>%
   group_by(wyear) %>%
   mutate(snow_mm_cum=cumsum(snow_mm),
         snow_na=cumsum(snow_na)) %>%
   ungroup() %>%
   mutate(snow_in_cum=round(snow_mm_cum/25.4, 1),
         snwd_in=round(snwd_mm/25.4, 0))

nov_18_snow <- cum_snow %>%
   filter(mmdd=='Nov 18') %>%
   select(wyear, snow_in_cum, snwd_in) %>%
   arrange(snow_in_cum)

decadal_avg <- nov_18_snow %>%
   mutate(decade=as.integer(wyear/10)*10) %>%
   group_by(decade) %>%
   summarize(`Snow depth`=mean(snwd_in),
            snwd_sd=sd(snwd_in),
            `Cumulative Snowfall`=mean(snow_in_cum),
            snow_cum_sd=sd(snow_in_cum))

decadal_averages <- ggplot(decadal_avg %>%
                              gather(variable, value, -decade) %>%
                              filter(variable %in% c("Cumulative Snowfall",
                                                      "Snow depth")),
                           aes(x=as.factor(decade), y=value, fill=variable)) +
            theme_bw() +
            geom_bar(stat="identity", position="dodge") +
            scale_x_discrete(name="Decade", breaks=c(1950, 1960, 1970, 1980,
                                                   1990, 2000, 2010)) +
            scale_y_continuous(name="Inches", breaks=pretty_breaks(n=10)) +
            scale_fill_discrete(name="Measurement")

print(decadal_averages)

date_x_scale <- cum_snow %>%
   filter(grepl(' (01|15)', mmdd), wyear=='1994') %>%
   select(wdoy, mmdd)

cumulative_snowfall <-
   ggplot(cum_snow %>% filter(wyear %in% c(1953, 1954, 2014, 2006, 1990),
                              wdoy>183,
                              wdoy<320),
            aes(x=wdoy, y=snow_in_cum, colour=as.factor(wyear))) +
   theme_bw() +
   geom_smooth(data=cum_snow %>% filter(wdoy>183, wdoy<320),
               aes(x=wdoy, y=snow_in_cum),
               size=0.5, colour="darkcyan",
               inherit.aes=FALSE,
               se=FALSE) +
   geom_line(position="jitter") +
   scale_x_continuous(name="",
                     breaks=date_x_scale$wdoy,
                     labels=date_x_scale$mmdd) +
   scale_y_continuous(name="Cumulative snowfall (in)",
                     breaks=pretty_breaks(n=10)) +
   scale_color_discrete(name="Winter year")

print(cumulative_snowfall)

snow_depth <-
   ggplot(cum_snow %>% filter(wyear %in% c(1953, 1954, 1962, 2014, 1990),
                              wdoy>183,
                              wdoy<320),
            aes(x=wdoy, y=snwd_in, colour=as.factor(wyear))) +
   theme_bw() +
   geom_smooth(data=cum_snow %>% filter(wdoy>183, wdoy<320),
               aes(x=wdoy, y=snwd_in),
               size=0.5, colour="darkcyan",
               inherit.aes=FALSE,
               se=FALSE) +
   geom_line(position="jitter") +
   scale_x_continuous(name="",
                     breaks=date_x_scale$wdoy,
                     labels=date_x_scale$mmdd) +
   scale_y_continuous(name="Snow Depth (in)",
                     breaks=pretty_breaks(n=10)) +
   scale_color_discrete(name="Winter year")

print(snow_depth)

library(tidyverse)
library(lubridate)
library(broom)

races_db <- src_postgres(host="localhost", dbname="races")

combined_races <- tbl(races_db, build_sql(
    "SELECT race, year, lower(name) AS name, finish_time,
        year - age AS birth_year, sex
     FROM chronotrack
     UNION
     SELECT race, year, lower(name) AS name, finish_time,
        birth_year,
        CASE WHEN age_class ~ 'M' THEN 'M' ELSE 'F' END AS sex
     FROM sportalaska
     UNION
     SELECT race, year, lower(name) AS name, finish_time,
        NULL AS birth_year, NULL AS sex
     FROM other"))

races <- tbl(races_db, build_sql(
    "SELECT race,
        lower(regexp_replace(race, '[ ’]', '_', 'g')) AS race_str,
        date_part('year', date) AS year,
        miles
     FROM races"))

racing_data <- combined_races %>%
    inner_join(races) %>%
    filter(!is.na(finish_time))

racers <- racing_data %>%
    group_by(name) %>%
    summarize(races=n(),
              birth_year=min(birth_year),
              gender_filter=ifelse(sum(ifelse(sex=='M',1,0))==
                                   sum(ifelse(sex=='F',1,0)),
                                   FALSE, TRUE),
              gender=ifelse(sum(ifelse(sex=='M',1,0))>
                            sum(ifelse(sex=='F',1,0)),
                            'M', 'F')) %>%
    ungroup() %>%
    filter(gender_filter) %>%
    select(-gender_filter)

racing_data_filled <- racing_data %>%
    inner_join(racers, by="name") %>%
    mutate(birth_year=birth_year.y) %>%
    select(name, birth_year, gender, race_str, year, miles, finish_time) %>%
    group_by(name, race_str, year) %>%
    mutate(n=n()) %>%
    filter(!is.na(birth_year), n==1) %>%
    ungroup() %>%
    collect() %>%
    mutate(fixed=ifelse(grepl('[0-9]+:[0-9]+:[0-9.]+', finish_time),
                        finish_time,
                        paste0('00:', finish_time)),
           minutes=as.numeric(seconds(hms(fixed)))/60.0,
           pace=minutes/miles,
           age=year-birth_year,
           age_class=as.integer(age/10)*10,
           group=paste0(gender, age_class),
           gender=as.factor(gender)) %>%
    filter((miles<10 & pace<16) | (miles>=10 & pace<20)) %>%
    select(-fixed, -finish_time, -n)

speeds_combined <- racing_data_filled %>%
    select(name, gender, age, age_class, group, race_str, year, pace) %>%
    spread(race_str, pace)

main_races <- c('beat_beethoven', 'chena_river_run', 'midnight_sun_run',
                'run_of_the_valkyries', 'gold_discovery_run',
                'santa_claus_half_marathon', 'golden_heart_trail_run',
                'equinox_marathon', 'hoodoo_half_marathon')

race_formula_str <-
    lapply(seq(1, length(main_races)-1),
           function(i)
               lapply(seq(i+1, length(main_races)),
                      function(j) paste(main_races[[j]], '~',
                                        main_races[[i]],
                                        '+ gender', '+ age'))) %>%
    unlist()

race_formulas <- lapply(race_formula_str, function(i) as.formula(i)) %>%
    unlist()

lm_models <- map(race_formulas, ~ lm(.x, data=speeds_combined))

models <- tibble(start_race=factor(gsub('.* ~ ([^ ]+).*',
                                        '\\1',
                                        race_formula_str),
                                   levels=main_races),
                 predicted_race=factor(gsub('([^ ]+).*',
                                            '\\1',
                                            race_formula_str),
                                       levels=main_races),
                 lm_models=lm_models) %>%
    arrange(start_race, predicted_race)

model_stats <- glance(models %>% rowwise(), lm_models)
model_coefficients <- tidy(models %>% rowwise(), lm_models)

reduced_formula_str <- model_coefficients %>%
    ungroup() %>%
    filter(p.value<0.05, term!='(Intercept)') %>%
    mutate(term=gsub('genderM', 'gender', term)) %>%
    group_by(predicted_race, start_race) %>%
    summarize(independent_vars=paste(term, collapse=" + ")) %>%
    ungroup() %>%
    transmute(reduced_formulas=paste(predicted_race, independent_vars, sep=' ~ '))

reduced_formula_str <- reduced_formula_str$reduced_formulas

reduced_race_formulas <- lapply(reduced_formula_str,
                                function(i) as.formula(i)) %>% unlist()

reduced_lm_models <- map(reduced_race_formulas, ~ lm(.x, data=speeds_combined))

n_from_lm <- function(model) {
    summary_object <- summary(model)

    summary_object$df[1] + summary_object$df[2]
}

reduced_models <- tibble(start_race=factor(gsub('.* ~ ([^ ]+).*', '\\1', reduced_formula_str),
                                           levels=main_races),
                         predicted_race=factor(gsub('([^ ]+).*', '\\1', reduced_formula_str),
                                               levels=main_races),
                         lm_models=reduced_lm_models) %>%
    arrange(start_race, predicted_race) %>%
    rowwise() %>%
    mutate(n=n_from_lm(lm_models))

reduced_model_stats <- glance(reduced_models %>% rowwise(), lm_models)
reduced_model_coefficients <- tidy(reduced_models %>% rowwise(), lm_models) %>%
    ungroup()

coefficients_and_stats <- reduced_model_stats %>%
    inner_join(reduced_model_coefficients,
               by=c("start_race", "predicted_race", "n")) %>%
    select(start_race, predicted_race, n, r.squared, term, estimate)

write_csv(coefficients_and_stats,
          "coefficients.csv")

make_scatterplot <- function(start_race, predicted_race) {
   age_limits <- speeds_combined %>%
      filter_(paste("!is.na(", start_race, ")"),
               paste("!is.na(", predicted_race, ")")) %>%
      summarize(min=min(age), max=max(age)) %>%
      unlist()

   q <- ggplot(data=speeds_combined,
               aes_string(x=start_race, y=predicted_race)) +
            # plasma works better with a grey background
            # theme_bw() +
            geom_abline(slope=1, color="darkred", alpha=0.5) +
            geom_smooth(method="lm", se=FALSE) +
            geom_point(aes(shape=gender, color=age)) +
            scale_color_viridis(option="plasma",
                              limits=age_limits) +
            scale_x_continuous(breaks=pretty_breaks(n=10)) +
            scale_y_continuous(breaks=pretty_breaks(n=6))

   svg_filename <- paste0(paste(start_race, predicted_race, sep="-"), ".svg")

   height <- 9
   width <- 16
   resize <- 0.75

   svg(svg_filename, height=height*resize, width=width*resize)
   print(q)
   dev.off()
}

lapply(seq(1, length(main_races)-1),
      function(i)
            lapply(seq(i+1, length(main_races)),
                  function(j)
                        make_scatterplot(main_races[[i]], main_races[[j]])
                  )

 library(dplyr)
 library(readr)
 library(lubridate)
 library(ggplot2)
 library(scales)
 library(grid)
 library(gtable)

 race_dates <- read_fwf("equinox_marathon_dates.rst", skip=5, n_max=54,
                        fwf_positions(c(4, 6), c(9, 19), c("number", "race_date")))

 noaa <- src_postgres(host="localhost", dbname="noaa")
 # pivot <- tbl(noaa, build_sql("SELECT * FROM ghcnd_pivot
 #                               WHERE station_name = 'UNIVERSITY EXP STN'"))
 # pivot <- tbl(noaa, build_sql("SELECT * FROM ghcnd_pivot
 #                               WHERE station_name = 'COLLEGE OBSY'"))
 pivot <- tbl(noaa, build_sql("SELECT * FROM ghcnd_pivot
                               WHERE station_name = 'FAIRBANKS INTL AP'"))

 race_day_wx <- pivot %>%
     inner_join(race_dates, by=c("dte"="race_date"), copy=TRUE) %>%
     collect() %>%
     mutate(tmin_f=round((tmin_c*9/5.0)+32, 1), tmax_f=round((tmax_c*9/5.0)+32, 1),
            prcp_in=round(prcp_mm/25.4, 2),
            snow_in=round(snow_mm/25.4, 1), snwd_in=round(snow_mm/25.4, 1),
            awnd_mph=round(awnd_mps*2.2369, 1),
            wsf2_mph=round(wsf2_mps*2.2369), 1) %>%
     select(number, race_date, tmin_f, tmax_f, prcp_in, snow_in,
            snwd_in, awnd_mph, wsf2_mph)

 week_before_race_day_wx <- pivot %>%
     mutate(year=date_part("year", dte)) %>%
     inner_join(race_dates %>%
                    mutate(year=year(race_date)),
                copy=TRUE) %>%
     collect() %>%
     mutate(tmin_f=round((tmin_c*9/5.0)+32, 1), tmax_f=round((tmax_c*9/5.0)+32, 1),
            prcp_in=round(prcp_mm/25.4, 2),
            snow_in=round(snow_mm/25.4, 1), snwd_in=round(snow_mm/25.4, 1),
            awnd_mph=round(awnd_mps*2.2369, 1), wsf2_mph=round(wsf2_mps*2.2369, 1)) %>%
     select(number, year, race_date, dte, prcp_in, snow_in) %>%
     mutate(week_before=race_date-days(7)) %>%
     filter(dte<race_date, dte>=week_before) %>%
     group_by(number, year, race_date) %>%
     summarize(pweek_prcp_in=sum(prcp_in),
               pweek_snow_in=sum(snow_in))

 all_wx <- race_day_wx %>%
     inner_join(week_before_race_day_wx) %>%
     mutate(tavg_f=(tmin_f+tmax_f)/2.0,
            snow_label=ifelse(snow_in>0, '*', NA),
            pweek_snow_label=ifelse(pweek_snow_in>0, '*', NA)) %>%
     select(number, year, race_date, tmin_f, tmax_f, tavg_f,
            prcp_in, snow_in, snwd_in, awnd_mph, wsf2_mph,
            pweek_prcp_in, pweek_snow_in,
            snow_label, pweek_snow_label);

 write_csv(all_wx, "all_wx.csv")

 madis <- src_postgres(host="localhost", dbname="madis")

 pafa_fbsa <- tbl(madis,
                  build_sql("
   WITH pafa AS (
     SELECT dt_local, temp_f, wspd_mph
     FROM observations
     WHERE station_id = 'PAFA' AND date_part('month', dt_local) = 9),
   fbsa AS (
     SELECT dt_local, temp_f, wspd_mph
     FROM observations
     WHERE station_id = 'FBSA2' AND date_part('month', dt_local) = 9)
   SELECT pafa.dt_local, pafa.temp_f AS pafa_temp_f, pafa.wspd_mph as pafa_wspd_mph,
     fbsa.temp_f AS ester_dome_temp_f, fbsa.wspd_mph as ester_dome_wspd_mph
   FROM pafa
     INNER JOIN fbsa ON
       pafa.dt_local BETWEEN fbsa.dt_local - interval '15 minutes'
         AND fbsa.dt_local + interval '15 minutes'")) %>% collect()

 write_csv(pafa_fbsa, "pafa_fbsa.csv")

 ester_dome_temps <- lm(data=pafa_fbsa,
                        ester_dome_temp_f ~ pafa_temp_f)

 summary(ester_dome_temps)
 # Model and coefficients are significant, r2 = 0.794
 # intercept = -2.69737, slope = 0.94268

 all_wx_with_ed <- all_wx %>%
   mutate(ed_min_temp_f=round(ester_dome_temps$coefficients[1]+
                              tmin_f*ester_dome_temps$coefficients[2], 1),
          ed_max_temp_f=round(ester_dome_temps$coefficients[1]+
                              tmax_f*ester_dome_temps$coefficients[2], 1))

 make_gt <- function(outside, instruments, chamber, width, heights) {
     gt1 <- ggplot_gtable(ggplot_build(outside))
     gt2 <- ggplot_gtable(ggplot_build(instruments))
     gt3 <- ggplot_gtable(ggplot_build(chamber))
     max_width <- unit.pmax(gt1$widths[2:3], gt2$widths[2:3], gt3$widths[2:3])
     gt1$widths[2:3] <- max_width
     gt2$widths[2:3] <- max_width
     gt3$widths[2:3] <- max_width
     gt <- gtable(widths = unit(c(width), "in"), heights = unit(heights, "in"))
     gt <- gtable_add_grob(gt, gt1, 1, 1)
     gt <- gtable_add_grob(gt, gt2, 2, 1)
     gt <- gtable_add_grob(gt, gt3, 3, 1)

     gt
 }

temps <- ggplot(data=all_wx_with_ed, aes(x=year, ymin=tmin_f, ymax=tmax_f, y=tavg_f)) +
   # geom_abline(intercept=32, slope=0, color="blue", alpha=0.25) +
   geom_rect(data=all_wx_with_ed %>% head(n=1),
            aes(xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=32),
            fill="darkcyan", alpha=0.25) +
   geom_abline(aes(slope=0,
                  intercept=mean(all_wx_with_ed$tmin_f)),
               color="darkorange", alpha=0.50, linetype=2) +
   geom_abline(aes(slope=0,
                  intercept=mean(all_wx_with_ed$tmax_f)),
               color="darkorange", alpha=0.50, linetype=2) +
   geom_abline(aes(slope=0,
                  intercept=mean(all_wx_with_ed$ed_min_temp_f)),
               color="darkorange", alpha=0.50, linetype=1) +
   geom_abline(aes(slope=0,
                  intercept=mean(all_wx_with_ed$ed_max_temp_f)),
               color="darkorange", alpha=0.50, linetype=1) +
   geom_linerange(aes(ymin=ed_min_temp_f, ymax=ed_max_temp_f)) +
   # geom_smooth(method="lm", se=FALSE) +
   geom_linerange(size=3, color="grey30") +
   scale_x_continuous(name="", limits=c(1963, 2015), breaks=seq(1963, 2015, 2)) +
   scale_y_continuous(name="Temperature (deg F)", breaks=pretty_breaks(n=10)) +
   theme_bw() +
   theme(plot.margin=unit(c(1, 1, 0, 0.5), 'lines')) +  # t, r, b, l
   theme(axis.text.x=element_blank(), axis.title.x=element_blank(),
         axis.ticks.x=element_blank(), panel.grid.minor.x=element_blank()) +
   ggtitle("Weather during and in the week prior to the Equinox Marathon
            Fairbanks Airport Station")

 prcp <- ggplot(data=all_wx, aes(x=year, y=prcp_in)) +
     geom_bar(stat="identity") +
     geom_text(aes(y=prcp_in+0.025, label=snow_label)) +
     scale_x_continuous(name="", limits=c(1963, 2015), breaks=seq(1963, 2015)) +
     scale_y_continuous(name="Precipitation (inches)", breaks=pretty_breaks(n=5)) +
     theme_bw() +
     theme(plot.margin=unit(c(0, 1, 0, 0.5), 'lines')) +  # t, r, b, l
     theme(axis.text.x=element_blank(), axis.title.x=element_blank(),
           axis.ticks.x=element_blank(), panel.grid.minor.x=element_blank())

 pweek_prcp <- ggplot(data=all_wx, aes(x=year, y=pweek_prcp_in)) +
     geom_bar(stat="identity") +
     geom_text(aes(y=pweek_prcp_in+0.1, label=pweek_snow_label)) +
     scale_x_continuous(name="", limits=c(1963, 2015), breaks=seq(1963, 2015)) +
     scale_y_continuous(name="Pre-week precip (inches)", breaks=pretty_breaks(n=5)) +
     theme_bw() +
     theme(plot.margin=unit(c(0, 1, 0.5, 0.5), 'lines'),
           axis.text.x=element_text(angle=45, hjust=1, vjust=1),
           panel.grid.minor.x=element_blank())

 rescale <- 0.75
 full_plot <- make_gt(temps, prcp, pweek_prcp,
                      16*rescale,
                      c(7.5*rescale, 2.5*rescale, 3.0*rescale))
 pdf("equinox_weather_grid.pdf", height=13*rescale, width=16*rescale)
 grid.newpage()
 grid.draw(full_plot)
 dev.off()

 fai_ed_temps <- ggplot(data=pafa_fbsa, aes(x=pafa_temp_f, y=ester_dome_temp_f)) +
   geom_rect(data=pafa_fbsa %>% head(n=1),
               aes(xmin=-Inf, ymin=-Inf, xmax=(32+2.69737)/0.94268, ymax=32),
               color="black", fill="darkcyan", alpha=0.25) +
   geom_point(position=position_jitter()) +
   geom_smooth(method="lm", se=FALSE) +
   scale_x_continuous(name="Fairbanks Airport Temperature (degrees F)") +
   scale_y_continuous(name="Ester Dome Temperature (degrees F)") +
   theme_bw() +
   ggtitle("Relationship between Fairbanks Airport and Ester Dome Temperatures
           September, 2008-2013")

 pdf("pafa_fbsa_sept_temps.pdf", height=10.5, width=10.5)
 print(fai_ed_temps)
 dev.off()

 fai_ed_wspds <- ggplot(data=pafa_fbsa, aes(x=pafa_wspd_mph, y=ester_dome_wspd_mph)) +
   geom_point(position=position_jitter()) +
   geom_smooth(method="lm", se=FALSE) +
   scale_x_continuous(name="Fairbanks Airport Wind Speed (MPH)") +
   scale_y_continuous(name="Ester Dome Wind (MPH)") +
   theme_bw() +
   ggtitle("Relationship between Fairbanks Airport and Ester Dome Wind Speeds
           September, 2008-2013")

 pdf("pafa_fbsa_sept_wspds.pdf", height=10.5, width=10.5)
 print(fai_ed_wspds)
 dev.off()