downloads Holly Springs data and plots it, specifically looking at air temperature, soil 10 cm temp, and soil 10 cm H2O

20160625_holly_springs_hourly.R

# get hourly Holly Springs and plot it

library("ggplot2")
library("lubridate")
library("dplyr")
library("cowplot")
library("reshape2")

hollyHourly2015 <- read.table("http://www1.ncdc.noaa.gov/pub/data/uscrn/products/hourly02/2015/CRNH0203-2015-MS_Holly_Springs_4_N.txt",
                            header = FALSE)

# right, got 8760 hours now
# get columns of interest

hollyH <- hollyHourly2015[, c(4, 5, 10, 29, 30, 34, 35)]

colnames(hollyH) <- c("date", "time", "air",
                    "water5", "water10", "soil5", "soil10")

hollyH$date <- ymd(hollyH$date)

startTime <- ymd_hm(201412311900)

timeVector <- startTime + dhours(0:(length(hollyH$time) - 1))

hollyH$dateTime <- timeVector

hollyH$hour <- hour(hollyH$dateTime)

hollyH$month <- month(hollyH$dateTime)

hollyH <- filter(hollyH, air > -9000)

hollyH <- filter(hollyH, soil10 > -9000)

hollyH <- filter(hollyH, soil5 > -9000)

forPlotholly <- melt(hollyH, id.vars = "dateTime")

forPlotholly <- filter(forPlotholly, variable == "air" |
                       variable == "soil10")

forPlotholly$variable <- ifelse(forPlotholly$variable == "air", "Air", "Soil 10 cm")

# look at air temp and soil temp at 10 cm by hour

p <- ggplot(data = forPlotholly, aes(x = dateTime, y = value))
holly1 <- p + geom_point(aes(shape = variable, alpha = variable,
                             colour = variable)) +
  scale_shape_manual(values=c(1, 17)) +   
  scale_alpha_manual(values = c(0.2, 0.4)) +
  background_grid(major = "xy") +
  scale_x_datetime(date_labels = "%m", date_breaks = "1 month") +
  theme(legend.title = element_blank()) +
  labs(x = "Month",
       y = "Average hourly temperature (°C)",
       title = "Holly Springs, Mississippi, 2015")

holly1

save_plot("Desktop/holly_year.png", holly1, base_height = 6, base_aspect_ratio = 1.4)

p <- ggplot(data = hollyH, aes(x = dateTime, y = water10 * 100))
waterHol <- p + geom_point(shape = 1, colour = "blue", alpha = 0.3) +
  background_grid(major = "xy") +
  scale_y_continuous(limits = c(0, 50), breaks = seq(0, 50, 10)) +
  labs(x = "Date",
       y = "Average hourly Soil water\nat 10 cm depth (% by volume)",
       title = "Holly Springs, Mississippi, 2015")

waterHol

save_plot("Desktop/holly_water_year.png", waterHol, base_height = 6, base_aspect_ratio = 1.4)

# look at hot soil
hotSoil <- filter(hollyH, soil10 >= 23 & month >= 6 & month <= 9)

hotSoil$delta10 <- hotSoil$air - hotSoil$soil10

median(hotSoil$water10)

hotSoil$soilWaterLevel <- ifelse(hotSoil$water10 < .27, "Soil water < 27%", "Soil water >= 27%")

# do boxplot by time of day -- facet it by soil moisture
p <- ggplot(data = hotSoil, aes(x = as.factor(hour), y = delta10))
facetDeltaholly <- p + geom_boxplot(outlier.shape = NA) + 
  background_grid(major = "xy") +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_jitter(alpha = 0.1, shape = 1) +
  facet_wrap(~soilWaterLevel, ncol = 1) +
#  scale_y_continuous(limits = c(-8, 6)) +
  labs(x = "For the hour ending at this time",
       y = "Average air temperature - 10 cm average soil temperature, °C",
       title = "Holly Springs, MS, USA, 2015\nall hours from June 1 to Sept 30 when 10 cm soil temperature >= 23°C")

facetDeltaholly

save_plot("Desktop/figure3.png", 
          facetDeltaholly, base_height = 6, base_aspect_ratio = 1.4)

# try to ascertain something by taking diff of mean Delta air-soil of high
# and low water, done by hour of the day

eto <- hotSoil %>%
  group_by(soilWaterLevel, hour) %>%
  summarise(meanDelta = mean(delta10))

eto2 <- hotSoil %>%
  group_by(soilWaterLevel) %>%
  summarise(meanSoil = mean(soil10),
            meanAir = mean(air))

dry <- filter(eto, soilWaterLevel == "Soil water < 27%")
wet <- filter(eto, soilWaterLevel == "Soil water >= 27%")

colnames(dry) <- c("level", "hour", "deltaDry")
colnames(wet) <- c("level", "hour", "deltaWet")

forAno <- cbind.data.frame(dry, wet$deltaWet)

colnames(forAno) <- c("level", "hour", "deltaDry", "deltaWet")

p <- ggplot(data = eto, aes(x = hour, y = meanDelta))
meanDiffHolly <- p + geom_point(aes(colour = soilWaterLevel, shape = soilWaterLevel)) +
  scale_x_continuous(breaks = 0:23) +
  scale_y_continuous(breaks = -5:3) +
  background_grid(major = "xy") +
  theme(legend.title = element_blank()) +
  labs(title = "Holly Springs, MS, USA, 2015\nall hours from June 1 to Sep 30 when 10 cm soil temperature >= 23°C",
       x = "Hour ending at",
       y = "Avg difference of air temperature - 10 cm soil temperature, °C")

meanDiffHolly

save_plot("Desktop/figure4.png", 
          meanDiffHolly, base_height = 6, base_aspect_ratio = 1.4)