Skip to content

Instantly share code, notes, and snippets.

@micahwoods
Created June 25, 2016 12:01
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save micahwoods/f2b4791bcfbae1f5ebe8bd66264adf77 to your computer and use it in GitHub Desktop.
Save micahwoods/f2b4791bcfbae1f5ebe8bd66264adf77 to your computer and use it in GitHub Desktop.
downloads Holly Springs data and plots it, specifically looking at air temperature, soil 10 cm temp, and soil 10 cm H2O
# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment