Depth-time visualization of lake profile data

General
,
1

Greetings R Community.

I have a profiles of oxygen from an ice covered lake from fall, winter, and spring. I would like to make a 2-D contour plot with interpolated values between those that I have measured. I followed the excellent tutorial/blog posted by @paleolimbot, but my data are imperfect. First, since there is ice cover growth, the profiles do not always start at the surface (0m), so interpolating between depths needs to be done in a staggered way to account for the loss of samples in the upper water column.

library(tidyverse)
library(lubridate)
library(readxl)

truelove<-read_excel("truelove.xlsx")

print(truelove)

A tibble: 15 x 3

date depth oxygen_mL

1 1961-09-10 00:00:00 0 9.75
2 1961-09-10 00:00:00 2 9.69
3 1961-09-10 00:00:00 4 9.63
4 1961-09-10 00:00:00 6 9.56
5 1961-10-25 00:00:00 1 10.8
6 1961-10-25 00:00:00 2 10.5
7 1961-10-25 00:00:00 4 9.67
8 1961-10-25 00:00:00 6 6.62
9 1961-11-21 00:00:00 1 10.3
10 1961-11-21 00:00:00 2 10.0
11 1961-11-21 00:00:00 4 8.36
12 1961-11-21 00:00:00 6 6.39
13 1962-01-10 00:00:00 2 9.58
14 1962-01-10 00:00:00 4 8.31
15 1962-01-10 00:00:00 6 5.09

estimate_oxy_by_date <- function(target_date, target_depth) {
data_for_date <- truelove %>%
filter(date == target_date) %>%
arrange(depth)
approx(data_for_date$depth, data_for_date$oxygen_mL, xout = target_depth)$y
}

estimate_oxy_by_date(ymd("1962-05-07"), c(0, 2, 4, 5, 5.5, 6)) here the interpolation with depth works at each date samples in a profile were taken.

oxy_interp_depth<-crossing(tibble(date = unique(truelove$date)),
tibble(depth = seq(0, 6, length.out = 24))) %>%
group_by(date) %>%
mutate(oxygen_mL = estimate_oxy_by_date(date[1], depth)) #this also worked and I was able to plot interpolated profiles of oxygen

Now trying to interpolate between dates sampled

estimate_oxy_by_depth <- function(target_depth, target_date) {
data_for_depth <- oxy_interp_depth %>%
filter(depth == target_depth) %>%
arrange(date)
approx(data_for_depth$date, data_for_depth$oxygen_mL, xout = target_date)$y
}

estimate_oxy_by_date(target_depth = 4, target_date = seq(ymd("1961-10-25"), ymd("1962-11-21"), by = 1))

oxy_raster <- crossing(tibble(date = seq(ymd("1961-09-10"), ymd("1962-05-07"), by=2)),
tibble(depth = unique(oxy_interp_depth$depth))) %>%
group_by(depth) %>%
mutate(oxygen_mL = estimate_oxy_by_depth(depth[1], date))

This crossing() gave me the error:

Error in approx(data_for_depth$date, data_for_depth$oxygen_mL, xout = target_date) :
need at least two non-NA values to interpolate
Called from: approx(data_for_depth$date, data_for_depth$oxygen_mL, xout = target_date)First time posting in an R help, so please let me know what other information is helpful for you all.

I think this has to do with the fact that the profiles start at different depths through time. Just an idea though. Any help is greatly appreciated!

Thanks!

2

Thanks for the post and glad you enjoyed the read!

It's difficult to help with this extensive of an example - you will have more luck here and elsewhere if you create a smaller version of your problem and render it with the reprex package. The main feature of this is that all I have to do to help you is copy your code and run it (way faster for me). In your case, it would look like this:


library(tidyverse)
library(lubridate)
#> 
#> Attaching package: 'lubridate'
#> The following object is masked from 'package:base':
#> 
#>     date

truelove <- tibble(
    date = as.Date(c("1961-09-10", "1961-09-10", "1961-09-10", "1961-09-10",
             "1961-10-25", "1961-10-25", "1961-10-25", "1961-10-25",
             "1961-11-21", "1961-11-21", "1961-11-21", "1961-11-21", "1962-01-10",
             "1962-01-10", "1962-01-10")),
   depth = c(0, 2, 4, 6, 1, 2, 4, 6, 1, 2, 4, 6, 2, 4, 6),
   oxygen_mL = c(9.75, 9.69, 9.63, 9.56, 10.8, 10.5, 9.67, 6.62, 10.3, 10,
                 8.36, 6.39, 9.58, 8.31, 5.09)
)

estimate_oxy_by_date <- function(target_date, target_depth) {
  data_for_date <- truelove %>%
    filter(date == target_date) %>%
    arrange(depth)
  
  approx(data_for_date$depth, data_for_date$oxygen_mL, xout = target_depth)$y
}

oxy_interp_depth <- crossing(
  tibble(date = unique(truelove$date)),
  tibble(depth = seq(0, 6, length.out = 24))
) %>%
  group_by(date) %>%
  mutate(oxygen_mL = estimate_oxy_by_date(date[1], depth))

estimate_oxy_by_depth <- function(target_depth, target_date) {
  data_for_depth <- oxy_interp_depth %>%
    filter(depth == target_depth) %>%
    arrange(date)
  
  approx(data_for_depth$date, data_for_depth$oxygen_mL, xout = target_date)$y
}

oxy_raster <- crossing(tibble(date = seq(ymd("1961-09-10"), ymd("1962-05-07"), by = 2)),
                       tibble(depth = unique(oxy_interp_depth$depth))) %>%
  group_by(depth) %>%
  mutate(oxygen_mL = estimate_oxy_by_depth(depth[1], date))
#> Error in approx(data_for_depth$date, data_for_depth$oxygen_mL, xout = target_date): need at least two non-NA values to interpolate

Created on 2020-05-18 by the reprex package (v0.3.0)

As you suspected, the error is because there's only one measurement at zero depth, so approx() really is trying to interpolate where there isn't enough data to do so. The solution is to limit the depth in oxy_raster to only depth where this makes sense (greater than one). I'd insert filter(depth > 1) to make this work:

library(tidyverse)
library(lubridate)
#> 
#> Attaching package: 'lubridate'
#> The following object is masked from 'package:base':
#> 
#>     date

truelove <- tibble(
    date = as.Date(c("1961-09-10", "1961-09-10", "1961-09-10", "1961-09-10",
             "1961-10-25", "1961-10-25", "1961-10-25", "1961-10-25",
             "1961-11-21", "1961-11-21", "1961-11-21", "1961-11-21", "1962-01-10",
             "1962-01-10", "1962-01-10")),
   depth = c(0, 2, 4, 6, 1, 2, 4, 6, 1, 2, 4, 6, 2, 4, 6),
   oxygen_mL = c(9.75, 9.69, 9.63, 9.56, 10.8, 10.5, 9.67, 6.62, 10.3, 10,
                 8.36, 6.39, 9.58, 8.31, 5.09)
)

estimate_oxy_by_date <- function(target_date, target_depth) {
  data_for_date <- truelove %>%
    filter(date == target_date) %>%
    arrange(depth)
  
  approx(data_for_date$depth, data_for_date$oxygen_mL, xout = target_depth)$y
}

oxy_interp_depth <- crossing(
  tibble(date = unique(truelove$date)),
  tibble(depth = seq(0, 6, length.out = 24))
) %>%
  group_by(date) %>%
  mutate(oxygen_mL = estimate_oxy_by_date(date[1], depth))

estimate_oxy_by_depth <- function(target_depth, target_date) {
  data_for_depth <- oxy_interp_depth %>%
    filter(depth == target_depth) %>%
    arrange(date)
  
  approx(data_for_depth$date, data_for_depth$oxygen_mL, xout = target_date)$y
}

oxy_raster <- crossing(tibble(date = seq(ymd("1961-09-10"), ymd("1962-05-07"), by = 2)),
                       tibble(depth = unique(oxy_interp_depth$depth))) %>%
  filter(depth > 1) %>% 
  group_by(depth) %>%
  mutate(oxygen_mL = estimate_oxy_by_depth(depth[1], date))

oxy_raster
#> # A tibble: 2,400 x 3
#> # Groups:   depth [20]
#>    date       depth oxygen_mL
#>    <date>     <dbl>     <dbl>
#>  1 1961-09-10  1.04      9.72
#>  2 1961-09-10  1.30      9.71
#>  3 1961-09-10  1.57      9.70
#>  4 1961-09-10  1.83      9.70
#>  5 1961-09-10  2.09      9.69
#>  6 1961-09-10  2.35      9.68
#>  7 1961-09-10  2.61      9.67
#>  8 1961-09-10  2.87      9.66
#>  9 1961-09-10  3.13      9.66
#> 10 1961-09-10  3.39      9.65
#> # … with 2,390 more rows

Created on 2020-05-18 by the reprex package (v0.3.0)

3

Thanks, @paleolimbot for the help. And a big thanks for introducing me to the reprex package. This will save time when posting to the internet help boards and also fixing code with colleagues. Your code worked well, but I could not figure out why mine (which was identical) did not work. Then I realized that it was due to the read_excel(). My dates were not coming in as dates. I am a recovering PC user and so this mistake should be a warning to all other recovering PC users and beginning R users, while handy, read_excel() is not ideal.

4

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.