R/get_depths.R

In fawda123/CTDplot: Plot CTD water quality data along an estuarine axis

#' Match CTD points with closest hydro soundings
#'
#' Match CTD points with closest hydro soundings
#'
#' @param ctd_pts data.frame of CTD station locations, see details
#' @param dep_pts data.frame of depth sounding points, see details
#' @param expand numeric for expanding sample size of CTD points
#' @param plot logical to return plot of matched points
#' @param show_bath logical to show bathymetric depth points on the plot
#' @param show_stat logical to show stations on the plot
#' @param zoom numeric for map zoom, passed to \code{\link[ggmap]{get_stamenmap}}
#'
#' @details
#' CTD station locations are linearly expanded and then matched to the nearests points in the depth sounding layer.  This is done to create a smoother depth gradient in the contour plot.  The \code{ctd_pts} data.frame must have three columns for longitude (\code{Long}), latitude (\code{Latitude}), and station name (\code{Station}).  The \code{dep_pts} data.frame should have three columns for longitude (\code{Long}), latitude (\code{Latitude}), and depth (\code{Depth}, non-negative).  The \code{plot} argument can be used to verify the match of the interpolated points with the depth soundings.
#'
#' All coordinates are assumed to be geographic decimal degrees using the WGS 1984 projection, negative longitude is west of the Prime Meridian.
#'
#' @return
#' The expanded data.frame from \code{ctd_pts} including depth (\code{Depth}) and cumulative distance (\code{Dist}, in km) columns.  If \code{plot = TRUE}, a ggmap plot is returned showing the matched points.
#'
#' @export
#'
#' @import dplyr geosphere ggmap ggplot2 ggrepel
#'
#' @examples
#' \dontrun{
#' dt <- as.Date('2014-04-21')
#' toplo <- ctd[ctd$Date %in% dt, ]
#' get_depths(toplo, PB_dep_pts)
#' get_depths(toplo, PB_dep_pts, plot = TRUE)
#' }
get_depths <- function(ctd_pts, dep_pts, expand = 200, plot = FALSE, show_bath = TRUE, show_stat = TRUE, zoom = 12){

  # get unique locations, max depth at each location
  ctd_pts <- select(ctd_pts, Station, Long, Lat, Depth) %>%
    group_by(Station) %>%
    filter(Depth == max(Depth)) %>%
    ungroup

  # expand value has to be changed to account for adding stations to the interpolated list
  expand <- expand - nrow(ctd_pts) + 2

  # expand ctd points by linear interp
  xint <- stats::approx(x = ctd_pts$Long, n = expand)$y
  yint <- stats::approx(x = ctd_pts$Lat, n = expand)$y
  locs <- data.frame(
  	Long = xint,
  	Lat = yint
  	)

  # get distance of new points from all depth points
  # find the depth points with the minimum distance to each point
  clo_dep <- distm(rbind(locs, dep_pts[, c('Long', 'Lat')])) %>%
    .[1:nrow(locs), ((1 + nrow(locs)):ncol(.))] %>%
    apply(1, function(x) which.min(x)[1])
  locs <- data.frame(locs, dep_pts[clo_dep, ])

  # estimate cumulative distance between points from geosphere (km)
  # extract top part of lower triangle
  dist <- distm(locs[, c('Long', 'Lat')])/1000
  d <- row(dist) - col(dist)
  dist <- split(dist, d)$`1`
  locs$Dist <- c(0, cumsum(dist))

  # get distance between original survey stations
  distorig <- distm(ctd_pts[, c('Long', 'Lat')])/1000
  d <- row(distorig) - col(distorig)
  distorig <- split(distorig, d)$`1`
  ctd_pts$Dist <- c(0, cumsum(distorig))

  # for join segments on plot
  mtchs <- select(locs, Long, Lat, Long.1, Lat.1)

  # format output
  # add scale, smooth depth
  out <- mutate(locs,
    Station = NA
    ) %>%
    select(Station, Long, Lat, Depth, Dist) %>%
    .[-c(1, nrow(.)), ] %>%
    rbind(., ctd_pts) %>%
    arrange(Dist)

  # return plot if T
  if(plot){

    # extent
    ext <- ggmap::make_bbox(out$Long, out$Lat)
    map <- ggmap::get_stamenmap(ext, zoom = zoom, maptype = "toner-lite")

    # base map
    pbase <- ggmap::ggmap(map) +
      theme_bw() +
      theme(
        axis.title.x = element_blank(),
        axis.title.y = element_blank()
      )

    # add depth points to interped
    if(show_bath){

      p <- pbase +
        geom_segment(data = mtchs, aes(x = Long, y = Lat, xend = Long.1, yend = Lat.1)) +
        geom_point(data = dep_pts, aes(x = Long, y = Lat), alpha = 0.8) +
        geom_point(data = out, aes(x = Long, y = Lat, fill = Depth), pch = 21, size = 3, colour = 'black', alpha = 0.8) +
        scale_fill_distiller(palette = 'Spectral')

    # otherwise just interped
    } else {

      p <- pbase +
        geom_point(data = out, aes(x = Long, y = Lat, fill = Depth), pch = 21, size = 3, colour = 'black', alpha = 0.8) +
        scale_fill_distiller(palette = 'Spectral')

    }

    # add stations
    if(show_stat){

      p <- p +
        geom_label_repel(data = ctd_pts, aes(x = Long, y = Lat, label = Station),
          point.padding = grid::unit(0, "lines"),
          box.padding = unit(3, "lines"),
          alpha = 0.8,
          force = 2
        )

    }

    return(p)

  }

  return(out)

}