Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
downloads temperature data for 4 cities in 2015 then calculates and makes plots of GP and GDD and permutations
# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment