Created
September 18, 2019 13:52
-
-
Save paleolimbot/54773c414cdb23aad8669047d2ca84cb to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(tidyverse) | |
# best case: your data is already in the right form | |
perfect_data <- tibble( | |
day_in_year = c(1, 101, 201, 1, 101, 201, 1, 101, 201), | |
depth = c(1, 1, 1, 2, 2, 2, 3, 3, 3), | |
temp = c(15, 16, 17, 14, 15, 16, 13, 14, 15) | |
) | |
ggplot(perfect_data, aes(x = day_in_year, y = depth, fill = temp)) + | |
geom_raster() + | |
scale_y_reverse() | |
# imperfect case: your data is only equally spaced in one dimension | |
# in either the time- or depth- dimension | |
imerfect_data <- tibble( | |
day_in_year = c(1, 124, 201, 1, 124, 201, 1, 124, 201), | |
depth = c(1, 1, 1, 2, 2, 2, 3, 3, 3), | |
temp = c(15, 16, 17, 14, 15, 16, 13, 14, 15) | |
) | |
# create a function that will, given a depth, estimate the temp on any given day | |
estimate_temp <- function(target_depth, target_day_in_year) { | |
data_for_depth <- imerfect_data %>% filter(depth == target_depth) | |
# this is the easiest way to do a linear interpolation | |
approx(data_for_depth$day_in_year, data_for_depth$temp, xout = target_day_in_year)$y | |
} | |
# create a tibble with the time and depth values you want | |
perfect_data <- crossing( | |
tibble(day_in_year = seq(1, 201, by = 10)), | |
# values from the equally spaced dimension must exist in `imperfect_data` | |
tibble(depth = c(1, 2, 3)) | |
) %>% | |
# group by the equally spaced dimension | |
group_by(depth) %>% | |
# estimate the temperature (our estimator only works at a single depth) | |
mutate(temp = estimate_temp(depth[1], day_in_year)) | |
ggplot(perfect_data, aes(x = day_in_year, y = depth, fill = temp)) + | |
geom_raster() + | |
scale_y_reverse() | |
# hardest case: unequal in both time and depth (probably what you have) | |
hard_data <- tibble( | |
day_in_year = c(1, 124, 201, 1, 124, 201, 1, 124, 201), | |
depth = c(0.25, 0, 1, 2.5, 2, 2.1, 3, 3, 3), | |
temp = c(15, 16, 17, 14, 15, 16, 13, 14, 15) | |
) | |
# you can always use points to get a feel for what's going on | |
ggplot(hard_data, aes(day_in_year, depth, color = temp)) + | |
geom_point(size = 3) + | |
scale_y_reverse() | |
# create a function to estimate the temperature at any depth using a smoother | |
# you will probably have to increase "k" as you get more values | |
# roughly, k is "wigglyness"...low "k" will give you a smooth surface | |
estimate_temp <- function(target_depth, target_day_in_year) { | |
smooth_fit <- mgcv::gam(temp ~ s(day_in_year, depth, k = 4), data = hard_data) | |
predict(smooth_fit, newdata = tibble(day_in_year = target_day_in_year, depth = target_depth)) | |
} | |
# create a tibble with the time and depth values you want | |
perfect_data <- crossing( | |
tibble(day_in_year = seq(1, 201, by = 10)), | |
tibble(depth = seq(0, 3, by = 0.2)) | |
) %>% | |
# estimate the temperature | |
mutate(temp = estimate_temp(depth, day_in_year)) | |
ggplot(perfect_data, aes(x = day_in_year, y = depth, fill = temp)) + | |
geom_raster() + | |
scale_y_reverse() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
(rendered version of the above:)
Created on 2019-09-18 by the reprex package (v0.2.1)