downloads temperature data for 4 cities in 2015 then calculates and makes plots of GP and GDD and permutations

gp_vs_gdd.R

# compare C3 growth potential (GP) and growing degree days (GDD)
# at a few locations, start with Sydney, Minneapolis, Tokyo, London
# these data and charts are shown at http://www.blog.asianturfgrass.com/2016/01/g.html

# load libraries
library("dplyr")
library("lubridate")
library("ggplot2")
library("cowplot")
library("RColorBrewer")

# functions for conversion F to C and for calculating GP
f2c <- function(x) {
  tempC <- (x - 32) * (5/9)
  return(tempC)
}

# function to calculate gp C3
c3gp <- function(x) {
  GP <- 2.71828 ^ (-0.5 * ((x - 20) / 5.5) ^ 2)
  return(GP)
}

calcValues <- function(x) {
  tempC <- f2c(x)
  gp <- c3gp(tempC)
  gdd0 <- ifelse(tempC <= 0, 0, tempC)
  gdd10 <- ifelse(tempC <= 10, 0, tempC - 10)
  gpSum <- cumsum(gp)
  gddSum0 <- cumsum(gdd0)
  gddSum10 <- cumsum(gdd10)
  return(list(tempC = tempC, gp = gp, 
              gdd0 = gdd0, gdd10 = gdd10,
              gpSum = gpSum, gddSum0 = gddSum0, gddSum10 = gddSum10))
}

# get data for 2015, for simplicity get weather underground airport data for SYD, MSP, LHR, HND
syd <- read.csv("http://www.wunderground.com/history/airport/YSSY/2015/1/1/CustomHistory.html?dayend=31&monthend=12&yearend=2015&req_city=NA&req_state=NA&req_statename=NA%22&format=1",
                header = TRUE)

msp <- read.csv("http://www.wunderground.com/history/airport/KMSP/2015/1/1/CustomHistory.html?dayend=31&monthend=12&yearend=2015&req_city=NA&req_state=NA&req_statename=NA%22&format=1",
                header = TRUE)
  
lhr <- read.csv("http://www.wunderground.com/history/airport/EGLL/2015/1/1/CustomHistory.html?dayend=31&monthend=12&yearend=2015&req_city=NA&req_state=NA&req_statename=NA%22&format=1",
                header = TRUE)
  
hnd <- read.csv("http://www.wunderground.com/history/airport/RJTT/2015/1/1/CustomHistory.html?dayend=31&monthend=12&yearend=2015&req_city=NA&req_state=NA&req_statename=NA%22&format=1",
                header = TRUE)

syd$date <- ymd(as.character(syd$AEDT))
msp$date <- ymd(as.character(msp$CST))
lhr$date <- ymd(as.character(lhr$GMT))
hnd$date <- ymd(as.character(hnd$JST))

syd$city <- "Sydney"
msp$city <- "Minneapolis"
lhr$city <- "London"
hnd$city <- "Tokyo"

syd <- select(syd, date, Mean.TemperatureF, city)
msp <- select(msp, date, Mean.TemperatureF, city)
lhr <- select(lhr, date, Mean.TemperatureF, city)
hnd <- select(hnd, date, Mean.TemperatureF, city)

sydOut <- as.data.frame(calcValues(syd$Mean.TemperatureF))
mspOut <- as.data.frame(calcValues(msp$Mean.TemperatureF))
lhrOut <- as.data.frame(calcValues(lhr$Mean.TemperatureF))
hndOut <- as.data.frame(calcValues(hnd$Mean.TemperatureF))

syd <- cbind.data.frame(syd, sydOut)
msp <- cbind.data.frame(msp, mspOut)
lhr <- cbind.data.frame(lhr, lhrOut)
hnd <- cbind.data.frame(hnd, hndOut)

all <- rbind.data.frame(syd, msp, lhr, hnd)

# plot temperature
p <- ggplot(data = all, aes(x = date, y = tempC))
t <- p + geom_point(aes(colour = city), alpha = 0.2) + 
  geom_smooth(se = FALSE, aes(colour = city), span = 0.25,
              size = 0.6) +
  scale_color_brewer(palette = "Set1") +
  background_grid(major = "xy") +
  theme(legend.title=element_blank(),
        legend.position="top",
        plot.margin=unit(c(0,0.5,0,0), "cm")) +
  labs(x = "",
       y = "Mean daily temperature (°C)")
save_plot("temperature.svg", t,
          base_aspect_ratio = 1.3)

# plot GP
p <- ggplot(data = all, aes(x = date, y = gp))
t <- p + geom_point(aes(colour = city), alpha = 0.2) + 
  geom_smooth(se = FALSE, aes(colour = city), span = 0.25,
              size = 0.6) +
  scale_color_brewer(palette = "Set1") +
  background_grid(major = "xy") +
  theme(legend.title=element_blank(),
        legend.position="top",
        plot.margin=unit(c(0,0.5,0,0), "cm")) +
  labs(x = "",
       y = bquote("Temperature-based"~C[3]~"growth potential (GP)"))
save_plot("gp.svg", t,
          base_aspect_ratio = 1.3)

# plot GDD zero
p <- ggplot(data = all, aes(x = date, y = gdd0))
t <- p + geom_point(aes(colour = city), alpha = 0.2) + 
  geom_smooth(se = FALSE, aes(colour = city), span = 0.25,
              size = 0.6) +
  scale_color_brewer(palette = "Set1") +
  background_grid(major = "xy") +
  theme(legend.title=element_blank(),
        legend.position="top",
        plot.margin=unit(c(0,0.5,0,0), "cm")) +
  labs(x = "",
       y = bquote("Growing degree days ("*GDD[0]*")"))
save_plot("gdd0.svg", t,
          base_aspect_ratio = 1.3)

# plot GDD 10
p <- ggplot(data = all, aes(x = date, y = gdd10))
t <- p + geom_point(aes(colour = city), alpha = 0.2) + 
  geom_smooth(se = FALSE, aes(colour = city), span = 0.25,
              size = 0.6) +
  scale_color_brewer(palette = "Set1") +
  background_grid(major = "xy") +
  theme(legend.title=element_blank(),
        legend.position="top",
        plot.margin=unit(c(0,0.5,0,0), "cm")) +
  labs(x = "",
       y = bquote("Growing degree days ("*GDD[10]*")"))
save_plot("gdd10.svg", t,
          base_aspect_ratio = 1.3)

# look at just two cities to avoid too much overlap
twoCity <- filter(all, city == "London" |
                    city == "Minneapolis")

# plot Temp vs GP 
p <- ggplot(data = twoCity, aes(x = tempC, y = gp))
t <- p + geom_jitter(size = 1,
                aes(colour = city), alpha = 0.3, height = 0.025, width = 0.5) + 
  scale_color_brewer(palette = "Set1") +
  background_grid(major = "xy") +
  theme(legend.title=element_blank(),
        legend.position="top",
        plot.margin=unit(c(0,0.5,0,0), "cm")) +
  labs(x = "Mean daily temperature (°C)",
       y = bquote("Temperature-based"~C[3]~"growth potential (GP)"))
save_plot("tVsGp.svg", t,
          base_aspect_ratio = 1.3)

# plot temp vs GDD 
p <- ggplot(data = twoCity, aes(x = tempC, y = gdd0))
t <- p + geom_jitter(size = 1, aes(colour = city), alpha = 0.3) + 
  scale_color_brewer(palette = "Set1") +
  background_grid(major = "xy") +
  theme(legend.title=element_blank(),
        legend.position="top",
        plot.margin=unit(c(0,0.5,0,0), "cm")) +
  labs(x = "Mean daily temperature (°C)",
       y = bquote("Growing degree days ("*GDD[0]*")"))
save_plot("tVsGdd0.svg", t,
          base_aspect_ratio = 1.3)

# plot cumSum of GP
p <- ggplot(data = all, aes(x = date, y = gpSum))
t <- p + geom_step(aes(colour = city)) +
  scale_color_brewer(palette = "Set1") +
  background_grid(major = "xy") +
  theme(legend.title=element_blank(),
        legend.position="top",
        plot.margin=unit(c(0,0.5,0,0), "cm")) +
  labs(x = "",
       y = "Cumulative sum of GP")
save_plot("cumsumGP.svg", t,
          base_aspect_ratio = 1.3)
  
# plot cumSum of gdd0
p <- ggplot(data = all, aes(x = date, y = gddSum0))
t <- p + geom_step(aes(colour = city)) +
  scale_color_brewer(palette = "Set1") +
  background_grid(major = "xy") +
  theme(legend.title=element_blank(),
        legend.position="top",
        plot.margin=unit(c(0,0.5,0,0), "cm")) +
  labs(x = "",
       y = bquote("Cumulative sum of"~GDD[0]))
save_plot("cumsumGDD0.svg", t,
          base_aspect_ratio = 1.3)

# plot cumSum of gdd10
p <- ggplot(data = all, aes(x = date, y = gddSum10))
t <- p + geom_step(aes(colour = city)) +
  scale_color_brewer(palette = "Set1") +
  background_grid(major = "xy") +
  theme(legend.title=element_blank(),
        legend.position="top",
        plot.margin=unit(c(0,0.5,0,0), "cm")) +
  labs(x = "",
       y = bquote("Cumulative sum of"~GDD[10]))
save_plot("cumsumGDD10.svg", t,
          base_aspect_ratio = 1.3)

p <- ggplot(data = all, aes(x = gpSum, y = gddSum0))
t <- p + geom_step(aes(colour = city)) +
  scale_color_brewer(palette = "Set1") +
  background_grid(major = "xy") +
  theme(legend.title=element_blank(),
        legend.position="top",
        plot.margin=unit(c(0,0.5,0,0), "cm")) +
  labs(x = "Cumulative sum of GP",
       y = bquote("Cumulative sum of"~GDD[0]))
save_plot("gpVsGDD0.svg", t,
          base_aspect_ratio = 1.3)

p <- ggplot(data = all, aes(x = gpSum, y = gddSum10))
t <- p + geom_step(aes(colour = city)) +
  scale_color_brewer(palette = "Set1") +
  background_grid(major = "xy") +
  theme(legend.title=element_blank(),
        legend.position="top",
        plot.margin=unit(c(0,0.5,0,0), "cm")) +
  labs(x = "Cumulative sum of GP",
       y = bquote("Cumulative sum of"~GDD[10]))
save_plot("gpVsGDD10.svg", t,
          base_aspect_ratio = 1.3)