R/footprint.R

In grahamstewart12/fluxtools:

## =============================================================================
#' Calculate Fetch Distances for a Dataset
#'
#' @description
#'
#' @param data
#' @param height
#' @param displacement
#' @param roughness_length
#' @param percent
#' @param method
#'
#' @return
#' @export
#'
#' @examples
calculate_fetch <- function(data, ws = "ws", ustar = "ustar", zeta = "zeta",
                            l = "l", height, displacement, roughness_length,
                            percent, method = c("KM01", "H00")) {
  cat(paste0("\n", "Calculating fetch lengths", "\n"))
  # Simplify variables names
  z <- height
  zd <- displacement
  zo <- roughness_length
  x <- rep(NA, nrow(data))
  p <- sort(percent)
  method <- match.arg(method)
  pbar <- dplyr::progress_estimated(length(x))
  if (method == "KM01") {
    vars <- c(ws, ustar, zeta)
    if (!all(vars %in% names(data))) {
      stop("Missing ", paste0(
        vars[!(vars %in% names(data))], collapse = ", "
      ), ".")
    }
    out <- data[, 0]
    names <- dplyr::if_else(
      !suppressWarnings(is.na(as.numeric(p))),
      paste0("x_", p, "perc"),
      paste0("x_", p)
    )
    out[, names] <- NA
    for (i in seq_along(x)) {
      # Set model inputs - Kormann & Meixner 2001
      ws_ <- data[i, ws]
      ustar_ <- data[i, ustar]
      zeta_ <- data[i, zeta]
      # Execute model function
      out[i, ] <- kormann_ftp_dist(ws_, ustar_, zeta_, z, zd, p)
      pbar$tick()$print()
    }
    # Set varnames and units
    varnames(out) <- paste0(
      "calculate_fetch(percent = ", names(out), ", method = ", method, ")"
    )
    openeddy::units(out) <- rep("m", ncol(out))
  } else if (method == "H00") {
    if (!l %in% names(data)) stop("Missing l.")
    out <- data[, 0]
    names <- dplyr::if_else(
      !suppressWarnings(is.na(as.numeric(p))),
      paste0("x_", p, "perc"),
      paste0("x_", p)
    )
    names <- c("xoffset", names)
    out[, names] <- NA
    for (i in seq_along(x)) {
      # Set model inputs - Hsieh 2000
      l_ <- data[i, l]
      # Execute model function
      out[i, ] <- hsieh_ftp_dist(z, zd, zo, l_, p)
      pbar$tick()$print()
    }
    # Set varnames and units
    varnames(out) <- paste0(
      "calculate_fetch(percent = ", names(out), ", method = ", method, ")"
    )
    openeddy::units(out) <- rep("m", ncol(out))
  }
  # Return vector if only one p was requested
  if (ncol(out) < 2) out <- out[, 1]
  return(out)
}

## =============================================================================
#' Calculate two-dimensional footprints
#'
#' Create a two-dimensional matrix of footprint probabilities.
#'
#' @param timestamp A vector of class POSIXct
#' @param wind_dir A numeric vector representing wind direction in degrees from
#'   North.
#' @param ustar A numeric vector representing friction velocity in units
#'   ms+1s-1.
#' @param L A numeric vector representing Monin-Obukhov length in units m.
#' @param v_sd A numeric vector representing the standard deviation of the
#'   crosswind component.
#' @param wind_speed A numeric vector representing wind speed in units ms+1s-1.
#' @param fetch An integer specifying the maximum length in each direction to
#'   calculate footprint probabilities. Defaults to 100 * (z - zd).
#' @param grid An integer specifying the number of grid cells to be used in
#'   calculating footprint probabilities. Defaults to the fetch length (i.e. a
#'   1-m resolution).
#' @param model A character string naming the model to be used in footprint
#'   calculations. Can be "KM01" for the Kormann & Meixner (2001) model or "H00"
#'   for the Hsieh et al. (2000) model.
#' @param tower_height The tower height (z) in meters.
#' @param displacement The zero-plane displacement height (zd) in meters.
#' @param roughness_length The aerodynamic roughness length (z0) in meters.
#' @param percent_out
#' @param raster_out
#' @param coords
#' @param mask_area
#' @param normalize
#' @param path
#'
#' @return
#' @export
#'
#' @examples
fp_calculate <- function(timestamp, wind_dir, ustar, L, v_sd, wind_speed,
                         weights = NULL, model = c("KM01", "H00"),
                         tower_height, displacement, roughness_length,
                         fetch = (tower_height - displacement) * 100,
                         grid = 2 * fetch, percent_out = FALSE,
                         raster_out = FALSE, coords = NULL, mask_area = NULL,
                         normalize = TRUE, path = NULL) {
  #browser()
  model <- match.arg(model)
  # Simplify variable names
  wd <- wind_dir
  ws <- wind_speed
  z <- tower_height
  zd <- displacement
  zo <- roughness_length

  # Check inputs required for either model
  if (!is.null(weights)) check_input(weights, "numeric_na")
  if (any(missing(wd), missing(ustar), missing(L), missing(v_sd))) {
    stop("All models require 'wd', 'ustar', 'L', and 'v_sd'.", call. = FALSE)
  }
  # Check inputs based on selected model
  if (model == "KM01") {
    # Requires z, zd, zo, L, wd, ustar, v_sd, ws
    if (missing(ws)) stop("KM01 model requires 'ws'.", call. = FALSE)
    lens <- c(length(ws), length(wd), length(ustar), length(L), length(v_sd))
  } else if (model == "H00") {
    # Requires z, zd, zo, L, wd, ustar, v_sd
    lens <- c(length(wd), length(L), length(v_sd), length(ustar))
  }
  # Check input lengths
  if (min(lens) == max(lens)) {
    len <- lens[1]
  } else stop("All variables must be of same length.", call. = FALSE)
  # Check that coordinates are given if raster output is selected
  if (raster_out & is.null(coords)) {
    stop("Coordinates must be provided for raster output.", call. = FALSE)
  }

  # Create general spatial grid
  dm <- fp_grid(d = 1, n = grid)
  # Create empty list to hold footprints
  out <- vector("list", length = len)

  # Begin loop for calculating footprint at each input record
  pbar <- dplyr::progress_estimated(len)
  for (i in 1:len) {
    # Create temporary object for holding footprint
    fp <- list()
    # Set output name as the timestamp
    name <- format(timestamp[i], "%Y-%m-%dT%H%M%S")

    # Execute footprint engine
    fp <- fp_engine(
      i, wd, ustar, L, v_sd, ws, z, zd, zo, weights, model, grid, dm
    )

    # Check for invalid or null output
    if (fp_is_empty(fp)) {
      # If no path, add null footprint to list
      if (is.null(path)) {
        out[[i]] <- fp
        names(out)[[i]] <- name
      }
      pbar$tick()$print()
      next
    }

    # Prepare for output differently depending on output type
    if (raster_out) {
      # Convert matrices to raster
      fp <- fp_rasterize(fp, coords, mask_area)
      phi_sum <- raster::cellStats(fp, sum)
      # Multiply by weights if indicated
      if (!is.null(weights)) {
        fp <- raster::calc(fp, function(x) x * weights[i])
        phi_sum <- raster::cellStats(fp, sum)
      }
      # Normalize *after* clipping
      if (normalize) fp <- raster::calc(fp, function(x) x / phi_sum)
      # Include phi_sum as an attribute
      attr(fp, "phi_sum") <- phi_sum
      # Set raster name as the timestamp
      names(fp) <- format(as.numeric(timestamp[i]), scientific = FALSE)
    } else {
      # (Output footprint as a list of matrices)
      phi_sum <- sum(sum(fp$z))
      #phi_sum <- sum(sum(phi))
      # Multiply by weights if indicated
      if (!is.null(weights)) {
        fp$z <- fp$z * weights[i]
        phi_sum <- sum(sum(fp$z))
      }
      # Normalize integrated footprint to 1
      if (normalize) fp$z <- fp$z / phi_sum
      # Include phi_sum in output list
      fp$phi_sum <- phi_sum
    }

    # Write footprint to file if path is given
    if (!is.null(path)) {
      full_path <- paste0(path, "/", name, "_fp", model)
      # Different handling depending on footprint format
      if (raster_out) {
        raster::writeRaster(fp, full_path, format = "raster")
      } else {
        fp <- fp_tidy(fp)
        write.csv(fp, full_path)
      }
    } else {
      # If no path, add footprint to list
      out[[i]] <- fp
      names(out)[[i]] <- name
    }
    pbar$tick()$print()
  }
  # Return footprint list if not writing to file
  if (is.null(path)) out
}

## =============================================================================
#' Title
#'
#' @param x
#' @param weights
#' @param output
#'
#' @return
#' @export
#'
#' @examples
fp_average <- function(x, weights = NULL, output = c("matrix", "raster")) {
  #browser()
  class <- class(x)
  if (class == "RasterLayer") {}
  output <- match.arg(output)
  # Get grid size from input object
  grid <- length(x[[1]][[1]][, 1])
  len <- length(x)
  if (len >= 1) {
    if (!is.null(weights)) {
      avg_ftp <- list(
        z = matrix(0, nrow = grid, ncol = grid),
        fpw = matrix(0, nrow = grid, ncol = grid)
      )
    } else avg_ftp <- list(z = matrix(0, nrow = grid, ncol = grid))
    #finput <- t(mapply(list, Ws, Wd, ustar, zol, sigma_v))
    count = 0
    pbar <- dplyr::progress_estimated(len)
    for (i in 1:len) {
      fp <- x[[i]]
      if (any(!is.na(fp$z))) {
        if (!is.null(weights)) {
          # If using weights, only include record if both fp and weight exist
          if (!is.na(weights[i])) {
            avg_ftp$z <- avg_ftp$z + fp$z
            avg_ftp$fpw <- avg_ftp$fpw + fp$z * weights[i]
            count <- count + 1
          }
        } else if (is.null(weights)) {
          # Only need fp if not using weights
          avg_ftp$z <- avg_ftp$z + fp$z
          count <- count + 1
        }
      }
      pbar$tick()$print()
    }
  }
  avg_ftp$z <- avg_ftp$z / count
  if (!is.null(weights)) {
    # Similar weighting routine to Budishchev et al. (2014)
    wsum <- sum(avg_ftp$fpw, na.rm = TRUE) # equals sum of weights
    avg_ftp$fpw <- avg_ftp$fpw / wsum
    #avg_ftp$fpw <- avg_ftp$fpw / count # this preserves the weights units
  }
  out <- list(
    x = fp$x, # distance east, from last footprint used in avg
    y = fp$y, # distance north, from last footprint used in avg
    z = if (is.null(weights)) avg_ftp$z else avg_ftp$fpw # probability
  )
  if (output == "raster") {
    out <- raster::raster(out)
  }
  out
}

## =============================================================================
#' Calculate footprint topology
#'
#' Calculates the proportion of footprint coverage within a given area of
#' interest.
#'
#' @param timestamp
#' @param wind_dir
#' @param ustar
#' @param L
#' @param v_sd
#' @param wind_speed
#' @param model
#' @param tower_height
#' @param displacement
#' @param roughness_length
#' @param fetch
#' @param grid
#' @param coords
#' @param mask_area
#'
#' @return
#' @export
#'
#' @examples
fp_topology <- function(timestamp, wind_dir, ustar, L, v_sd, wind_speed,
                        model = c("KM01", "H00"), tower_height, displacement,
                        roughness_length,
                        fetch = (tower_height - displacement) * 100,
                        grid = 2 * fetch, coords, mask_area) {
  model <- match.arg(model)
  # Simplify variable names
  wd <- wind_dir
  ws <- wind_speed
  z <- tower_height
  zd <- displacement
  zo <- roughness_length

  # Check inputs required for either model
  if (any(missing(wd), missing(ustar), missing(L), missing(v_sd))) {
    stop("All models require 'wd', 'ustar', 'L', and 'v_sd'.", call. = FALSE)
  }
  # Check inputs based on selected model
  if (model == "KM01") {
    # Requires z, zd, zo, L, wd, ustar, v_sd, ws
    if (missing(ws)) stop("KM01 model requires 'ws'.", call. = FALSE)
    lens <- c(length(ws), length(wd), length(ustar), length(L), length(v_sd))
  } else if (model == "H00") {
    # Requires z, zd, zo, L, wd, ustar, v_sd
    lens <- c(length(wd), length(L), length(v_sd), length(ustar))
  }
  # Check input lengths
  if (min(lens) == max(lens)) {
    len <- lens[1]
  } else stop("All variables must be of same length.", call. = FALSE)
  # Check that coordinates and a site area object are given
  if (missing(coords)) {
    stop("Coordinates 'coords' must be provided.", call. = FALSE)
  }
  if (missing(mask_area)) {
    stop("Area of interest 'mask_area' must be provided.", call. = FALSE)
  }

  # Create general spatial grid
  dm <- fp_grid(d = 1, n = grid)
  # Create empty vector to hold phi values
  out <- vector("double", length = len)

  # Begin loop for calculating footprint at each input record
  pbar <- dplyr::progress_estimated(len)
  for (i in 1:len) {
    # Create temporary object for holding footprint
    fp <- list()

    # Execute footprint engine
    fp <- fp_engine(
      i, wd, ustar, L, v_sd, ws, z, zd, zo, weights = NULL, model, grid, dm
    )

    # Check for invalid or null output
    if (fp_is_empty(fp)) {
      out[i] <- NA
      pbar$tick()$print()
      next
    }

    # Convert matrices to raster
    fp <- fp_rasterize(fp, coords, mask_area)
    # Calculate phi value
    phi <- raster::cellStats(fp, sum)
    out[i] <- phi

    pbar$tick()$print()
  }
  # Return output vector
  out
}

## =============================================================================
#' Fetch Distance Using Hsieh (2000) Model
#'
#' @description
#'
#' @param z
#' @param zd
#' @param zo
#' @param L
#' @param p
#'
#' @return
#' @export
#'
#' @examples
hsieh_ftp_dist <- function(z, zd, zo, L, p) {
  # Recommended to gap-fill L prior to using this model
  if (!is.na(L)) {
    # Define and calculate footprint parameters
    k <- 0.4
    di <- 5
    zm <- z - zd
    zu <- zm * (log(zm / zo) - 1 + (zo / zm))
    zL <- zu / L
    if (abs(zL) < 0.04) {
      D <- 0.97
      P <- 1
    } else if (zL < 0) {
      D <- 0.28
      P <- 0.59
    } else if (zL > 0) {
      D <- 2.44
      P <- 1.33
    }
    DP <- D * zu^P * abs(L)^(1 - P)
    # Calculate x if flux percents are given
    all_p <- p
    p <- as.numeric(all_p[!suppressWarnings(is.na(as.numeric(all_p)))])
    if (length(p) > 0) {
      # Define footprint model
      f_int <- function(x) exp(-(D * zu^P * abs(L)^(1 - P) / (k^2 * x)))
      #f <- function(x) (DP / (k^2 * x^2)) * (exp(-DP / (k^2 * x)))
      p <- c(1, sort(p))
      fx <- rep(NA, length(p))
      x <- 1
      while (x < 10000) {
        #fp <- try(integrate(f, 0, x)$value, silent = TRUE)
        #if (class(fp) == "try-error") {
        # Certain (uncommon) meteorological conditions result in a non-finite
        # integral of the footprint function: set these records to NA
        #  fp <- NA
        #  break
        #}
        fp <- f_int(x)
        fx[which(((fp > (p / 100) & is.na(fx))) == TRUE)[1]] <- x
        if (!anyNA(fx)) break # finished if all fetch lengths are found
        x <- x + 1
      }
    }
    # Calculate x if flux peak is requested
    if ("peak" %in% all_p) {
      peak <- DP / (2 * k^2)
      fx <- c(fx, peak)
    }
  } else {
    # Cannot get fetch lengths if meteorological parameters are absent
    fx <- rep(NA, length(p) + 1)
  }
  return(fx)
}

## =============================================================================
#' Fetch Distance Using Kormann & Meixner (2001) Model
#'
#' @description
#'
#' @param wind_speed
#' @param ustar
#' @param zeta
#' @param z
#' @param zd
#' @param p
#'
#' @return
#' @export
#'
#' @examples
kormann_ftp_dist <- function(wind_speed, ustar, zeta, z, zd, p) {
  # Simplify variable names
  ws <- wind_speed
  if (!is.na(ws) & !is.na(ustar) & !is.na(zeta)) {
    # Calculate model parameters
    zm <- z - zd  # measurement height
    k <- 0.4  # Von Karman constant
    if (zeta > 0) {
      # similarity relations
      phi_m <- 1 + 5 * zeta
      phi_c <- phi_m
      # exponent of diffusivity power law
      n <- 1 / phi_m
    } else if (zeta <= 0) {
      # similarity relations
      phi_m <- (1 - 16 * zeta)^-0.25
      phi_c <- (1 - 16 * zeta)^-0.5
      eta <- (1 - 16 * zeta)^0.25
      # exponent of the diffusivity power law
      n <- (1 - 24 * zeta) / (1 - 16 * zeta)
    }
    # Proportionality constant of the diffusivity power law (Eqs. 11 & 32)
    key <- k * ustar * zm / (phi_c * zm^n)
    # Exponent of the wind speed power law
    m <- ustar * phi_m / (k * ws)
    # Proportionality constant of the wind speed power law (Eqs. 11 & 31)
    U <- ws / (zm^m)
    # Intermediate parameters
    r <- 2 + m - n  # shape factor
    mu <- (1 + m) / r  # constant
    xi <- U * zm^r / (r^2 * key)  # flux length scale
    # Calculate x if flux percents are given
    all_p <- p
    p <- as.numeric(all_p[!suppressWarnings(is.na(as.numeric(all_p)))])
    if (length(p) > 0) {
      # Define footprint model as a function of x (distance from tower)
      f <- function(x) xi^mu * exp(-xi / x) / (x^(1 + mu) * gamma(mu))
      p <- sort(p)
      fx <- rep(NA, length(p))
      x <- 1
      while (x < 10000) {
        fp <- try(integrate(f, 0, x)$value, silent = TRUE)
        if (class(fp) == "try-error") {
          # Certain (uncommon) meteorological conditions result in a non-finite
          # integral of the footprint function: set these records to NA
          fp <- NA
          break
        }
        fx[which(((fp > (p / 100) & is.na(fx))) == TRUE)[1]] <- x
        if (!anyNA(fx)) break # finished if all fetch lengths are found
        x <- x + 1
      }
    }
    if ("peak" %in% all_p) {
      peak <- xi / (1 + mu)
      fx <- c(fx, peak)
    }
  } else {
    # Cannot get fetch lengths if meteorological parameters are absent
    fx <- rep(NA, length(p))
  }
  return(fx)
}

fp_engine <- function(i, wd, ustar, L, v_sd, ws, z, zd, zo, weights, model,
                      grid, dm) {
  # Check required variables
  if (model == "KM01") {
    vars <- c(ws[i], wd[i], ustar[i], L[i], v_sd[i])
  } else if (model == "H00") {
    vars <- c(wd[i], ustar[i], L[i], v_sd[i])
  }
  if (!is.null(weights)) vars <- c(vars, weights[i])

  # All required variables must be present - if not return empty matrix
  if (anyNA(vars)) {
    phi <- matrix(NA, nrow = grid, ncol = grid)
    fp <- list(x = dm$x, y = dm$y, z = phi)
    return(fp)
  }

  # Construct matrix to hold calculated values at each grid point
  phi <- matrix(0, nrow = grid + 1, ncol = grid + 1)

  # Wind direction and mathematical definitions of angles
  theta <- ((360 - wd[i]) %% 360) * (pi / 180)
  x_ <- dm$x * cos(theta) - dm$y * sin(theta)
  y_ <- dm$x * sin(theta) + dm$y * cos(theta)

  # Subset downwind area of grid for calculating contributions
  dw <- x_ > 0

  # Execute chosen model
  if (model == "KM01") {
    # Kormann model not valid if abs(zeta) > 3
    # Handle invalid cases
    if (!dplyr::between((z - zd) / L[i], -3, 3)) {
      fpp <- matrix(NA, nrow = grid, ncol = grid)
      fp <- list(fpp = fpp, fpx = fpx, fpy = fpy)
      return(fp)
    }

    # Evaluate model across grid points
    phi[dw] <- kormann_model(
      x_, y_, ws[i], ustar[i], L[i], v_sd[i], z, zd, zo, grid
    )

  } else if (model == "H00") {
    # Hsieh model
    phi[dw] <- hsieh_model(
      x_, y_, ustar[i], L[i], v_sd[i], z, zd, zo, grid
    )
  }

  # Output list of matrices holding all footprint data
  list(x = dm$x, y = dm$y, z = phi)
}

fp_grid <- function(extent, d, n) {
  if (missing(extent) & missing(d) & missing(n)) {
    extent <- c(-500, 500, -500, 500)
    dx <- 1
    dy <- 1
    nx <- 500
    ny <- 500
  } else if (!missing(extent)) {
    if (!missing(d)) {
      dx <- d
      dy <- d
      nx <- (extent[2] - extent[1]) / dx
      ny <- (extent[4] - extent[3]) / dy
    } else {
      if (missing(n)) {
        nx <- 500
        ny <- 500
      }
      dx <- (extent[2] - extent[1]) / nx
      dy <- (extent[4] - extent[3]) / ny
    }
  } else if (!missing(d) & !missing(n)) {
    dx <- d
    dy <- d
    nx <- n
    ny <- n
    extent <- c(-nx * dx / 2, nx * dx / 2, -ny * dy / 2, ny * dy / 2)
  } else if (!missing(d)) {
    extent <- c(-500, 500, -500, 500)
    dx <- d
    dy <- d
    nx <- (extent[2] - extent[1]) / dx
    ny <- (extent[4] - extent[3]) / dy
  } else if (!missing(n)) {
    nx <- n
    ny <- n
    extent <- c(-500, 500, -500, 500)
    dx <- (extent[2] - extent[1]) / nx
    dy <- (extent[4] - extent[3]) / ny
  }

  # Put extent into more convenient vars
  xmin <- extent[1]
  xmax <- extent[2]
  ymin <- extent[3]
  ymax <- extent[4]
  #extent <- raster::extent(extent)

  # Define physical domain in cartesian coordinates
  x <- seq(xmin, xmax, length.out = nx + 1)
  y <- -seq(ymin, ymax, length.out = ny + 1)

  x_2d <- matrix(rep(x, each = ny + 1), nrow = ny + 1)
  y_2d <- matrix(rep(y, nx + 1), nrow = ny + 1)
  f_2d <- matrix(0, nrow = ny + 1, ncol = nx + 1)

  list(x = x_2d, y = y_2d, f = f_2d)
}

fp_rotate <- function(a) {
  r <- a * pi / 180 # degrees to radians
  matrix(c(cos(r), sin(r), -sin(r), cos(r)), nrow = 2, ncol = 2)
}

fp_save <- function(x, path) {

  write.table(x, path, row.names = FALSE, col.names = FALSE)
  MASS::write.matrix(x, path, sep = "\t")
}

## =============================================================================
kormann_model <- function(x, y, ws, ustar, l, v_sd, z, zd, zo, grid) {
  # Not meant to be used by itself - use fp_calculate instead
  # Construct helper matrices
  h0 <- matrix(0, nrow = grid + 1, ncol = grid + 1)
  h1 <- matrix(0, nrow = grid + 1, ncol = grid + 1)
  h2 <- matrix(0, nrow = grid + 1, ncol = grid + 1)
  h3 <- matrix(0, nrow = grid + 1, ncol = grid + 1)

  ## KORMANN MODEL - via Forbrich et al. (2011)
  # Set model constants
  zm <- z - zd  # aerodynamic measurement height
  k <- 0.4  # Von Karman constant
  zmol <- zm / l
  # limited to -3 to +3 zeta
  # Dimensionless gradient functions of wind and temp profiles (Dyer 1974)
  if (l > 0) {
    phi <- 1 + 5 * zmol
    phim <- 1 + 5 * zmol
    phic <- 1 + 5 * zmol
    m <- (1 + 5 * (zmol)) / (log(zm / zo) + 5 * zmol)
    n <- 1 / (1 + 5 * zmol)
  } else if (l < 0) {
    zeta <- (1 - 16 * (zmol))^(-1 / 4)
    phi <- 1 / (zeta^2)
    phim <- (1 - 16 * (zmol))^(-1 / 4)
    phic <- (1 - 16 * (zmol))^(-1 / 2)
    m <- ustar / k / zeta / ws
    n <- (1 - 24 * zmol) / (1 - 16 * zmol)
  }
  # Coefficients of diffusivity power law (Eqs. A.4-A.5)
  # Many different ways to calculate U
  #U <- (ustar * log((zm / zo) + (phim))) / (k * zm^m)
  #U <- (ustar * phim) / (k * m * zm^m) # neftel
  #kappa <- (k * ustar * zm) / (phic * zm^n)
  kappa <- (k * ustar * zm) / (phi * zm^n)
  U <- ws / zm^n
  # Intermediate parameters (Eqs. A.6-A.7; p = m)
  r <- 2 + m - n
  mu <- (1 + m) / r
  # Flux length scale (Eq. A.8)
  xi <- (U * zm^r) / (r^2 * kappa)
  # Neftel parameters
  A <- 1 + mu
  B <- (U * (zm^r)) / (r^2 * kappa)
  C <- B^mu / gamma(mu)
  D <- v_sd / U * gamma(1 / r) / gamma(mu) * (U / (r^2 * kappa))^(m / r)
  E <- (r - m) / r

  # Subset downwind area of grid for calculating contributions
  dw <- x > 0
  # Splitting Eq 2.1 into a few helpers as per OzFlux
  h0[dw] <- x[dw]^E
  h1[dw] <- (sqrt(2 * pi) * D * (h0[dw]))
  h2[dw] <- exp(-(y[dw]^2) / (2 * (D * (h0[dw]))^2))
  h3[dw] <- C * exp(-B / x[dw]) * (x[dw]^(-A))
  # Combine helpers - two-dimensional footprint function (Eq. 2.1)
  h2[dw] * h3[dw] / h1[dw]
}

hsieh_model <- function(x, y, ustar, l, v_sd, z, zd, zo, grid) {
  # Not meant to be used by itself - use fp_calculate instead
  # Construct helper matrices
  f <- matrix(0, nrow = grid + 1, ncol = grid + 1)
  sig <- matrix(0, nrow = grid + 1, ncol = grid + 1)
  dy <- matrix(0, nrow = grid + 1, ncol = grid + 1)

  ## HSIEH MODEL - via Hsieh et al. (2000)
  # Set model constants
  zm <- z - zd  # aerodynamic measurement height
  k <- 0.4  # Von Karman constant
  # define and calculate footprint parameters
  zu <- zm * (log(zm / zo) - 1 + (zo / zm))
  zl <- zu / l
  if (abs(zl) < 0.04) {
    d <- 0.97
    p <- 1
  } else if (zl < 0) {
    d <- 0.28
    p <- 0.59
  } else if (zl > 0) {
    d <- 2.44
    p <- 1.33
  }
  dp <- d * zu^p * abs(l)^(1 - p)
  # Subset downwind area of grid for calculating contributions
  dw <- x > 0
  # calculate cross-wind integrated footprint
  f[dw] <- (dp / (k^2 * x[dw]^2)) * (exp(-dp / (k^2 * x[dw])))
  # define and calculate 2D extension parameters
  a1 <- 0.3
  p1 <- 0.86
  sig[dw] <- a1 * zo * (v_sd / ustar) * (x[dw] / zo) ^ p1
  dy[dw] <- 1 / (sqrt(2 * pi) * sig[dw]) * exp(-0.5 * (y[dw] / sig[dw]) ^ 2)
  # calculate 2D footprint adjusted by lateral diffusion
  dy[dw] * f[dw]
}

fp_tidy <- function(x, weights = FALSE) {
  # Define function for rotating grid
  rotate <- function(x) t(apply(x, 2, rev))
  fpx <- x$x
  fpy <- x$y
  if (weights == FALSE) fpz <- x$z else fpz <- x$fpw
  df <- data.frame(
    x = numeric(length(fpx)),
    y = numeric(length(fpy)),
    z = numeric(length(fpz))
  )
  df$z <- as.vector(fpz)
  df$x <- as.vector(fpx)
  df$y <- as.vector(fpy)
  t <- df[order(-df$z), ]
  m <- xtabs(z ~ x + y, data = t)
  f <- x
  f$z <- rotate(rotate(rotate(unclass(m))))
  f <- as.data.frame(as.table(f$z))
  f$y <- as.numeric(as.character(f$y))
  f$x <- as.numeric(as.character(f$x))
  names(f)[3] <- "z"
  f
}

fp_percent <- function(x, p, return = c("raster", "polygon")) {
  return <- match.arg(return)
  # Sort raster values from high to low
  rs <- sort(as.vector(x), decreasing = TRUE)
  # Calculate cumulative sums of footprint weights
  cs <- cumsum(rs)
  # Get the minimum value of cumulative sum less than p
  rp <- min(rs[cs < p])
  out <- x
  # Set all other cells to NA
  out[out < rp] <- NA
  # Return a polygon if requested
  if (return == "polygon") out <- fp_polygonize(out)
  out
}

fp_is_empty <- function(fp) {
  #browser()
  class <- class(fp)
  # If fp is a spatial object, it is not empty
  if ("SpatialLinesDataFrame" %in% class) {
    FALSE
  } else if ("SpatialPolygons" %in% class) {
    FALSE
  } else if ("RasterLayer" %in% class) {
    FALSE
  } else if (length(fp) == 1) {
    # If fp has length one, NA indicates it is empty
    is.na(fp)
    # If fp is a matrix, check contents
  } else all(!is.finite(fp$z))
}

#' Title
#'
#' @param x
#' @param coords
#' @param mask_area
#'
#' @return
#' @export
#'
#' @examples
fp_rasterize <- function(x, coords, mask_area) {
  # Add geographical coordinates to grid
  x$y <- x$y + coords[1] # lat
  x$x <- x$x + coords[2] # lon
  # Convert list of matrices to XYZ data frame
  out <- fp_tidy(x)
  # Create raster object
  out <- raster::rasterFromXYZ(
    out, crs = sp::CRS("+proj=utm +zone=18 ellps=WGS84")
  )
  # Mask by an extent if indicated
  if (!missing(mask_area)) {
    out <- suppressWarnings(raster::mask(out, mask_area))
  }
  out
}

fp_polygonize <- function(x) {
  class <- class(x)
  #browser()
  if (fp_is_empty(x)) return(NA)
  if (class == "RasterLayer") {
    # Method for RasterLayer class footprint
    # Set all raster values to 0 - they don't matter anymore
    homog <- raster::calc(x, function(x) x * 0)
    # Return polygon of the total area
    raster::rasterToPolygons(homog, dissolve = TRUE)
  } else if (class == "SpatialLinesDataFrame") {
    # Method for SpatialLinesDataFrame class footprint
    # Tidy ftp for creating polygon object
    fort <- broom::tidy(x)
    # Create spatial polygons object
    poly <- sp::Polygon(fort[, 1:2])
    polys <- sp::Polygons(list(sp::Polygon(poly)), "poly")
    sp::SpatialPolygons(
      list(polys), proj4string = sp::CRS("+proj=utm +zone=18 +ellps=GRS80")
    )
  }
}

write_footprint <- function(x, path, matrix = "fpp", save_log_vals = FALSE,
                            tag = "fp") {
  len <- length(x)
  if (purrr::vec_depth(x) == 2) {
    out <- x[[matrix]]
    path_out <- paste0(path, paste0(tag, ".csv"))
    write.table(out, path_out, sep = ",", row.names = FALSE, col.names = FALSE)
  } else if (purrr::vec_depth(x) == 3) {
    pbar <- dplyr::progress_estimated(len)
    for (i in 1:len) {
      out <- x[[i]][matrix]
      if (save_log_vals) {
        out <- log(out)
      }
      name_out <- paste0(tag, "_", names(x)[[i]], ".csv")
      path_out <- paste0(path, name_out)
      write.table(
        out, path_out, sep = ",", row.names = FALSE, col.names = FALSE
      )
      pbar$tick()$print()
    }
  }
}

read_footprint <- function(path) {
  files <- list.files(path)
  len <- length(files)
  out <- vector("list", len)
  pbar <- dplyr::progress_estimated(len)
  for (i in seq_along(files)) {
    out[[i]] <- read.table(paste0(path, files[i]), sep = ",")
    names(out)[[i]] <- stringr::str_replace(files[i], "fp_", "")
    names(out)[[i]] <- stringr::str_replace(files[i], "\\.csv", "")
    pbar$tick()$print()
  }
  out
}

spread_footprints <- function(x) {
  len <- length(x[[1]]$fpp[1, ])
  for (i in seq_along(x)) {
    fp <- as.data.frame(x[[i]]$fpp) %>%
      tidyr::gather(var, val) %>%
      dplyr::mutate(
        x = rep(seq(1:len), times = len, each = 1),
        y = rep(seq(1:len), times = 1, each = len)
      ) %>%
      tidyr::unite("xy", x, y) %>%
      select(-var) %>%
      tidyr::spread(xy, val) %>%
      dplyr::mutate(timestamp = names(x)[[i]]) %>%
      dplyr::select(timestamp, dplyr::everything())
    out <- dplyr::bind_rows(out, fp)
  }
}

project_footprint <- function(ftp, latitude, longitude, grid,
                              proj = c("longlat", "utm"), percents) {
  # Probably going to deprecate this - much easier to work w/ rasters then
  # convert to other spatial features if needed

  ## Create lat/long coordinates for each footprint grid box
  if (fp_is_empty(ftp)) {return(NA)}
  fp_p <- percent_footprint(ftp)
  footprint_c <- fp_p %>%
    dplyr::mutate(
      d = sqrt(x ^ 2 + y ^ 2),
      # calculate theta, correct for quadrants
      theta = ifelse(
        x > 0 & y > 0, atan(y / x) * 180 / pi,
        ifelse(
          x > 0 & y < 0,
          atan(y / x) * 180 / pi + 360, atan(y / x) * 180 / pi + 180
        )
      ),
      # convert theta to bearing using modulus
      b = (450 - theta) %% 360,
      # compute lat/long coordinates
      lat = geosphere::destPoint(p = c(longitude, latitude), b = b, d = d)[, 2],
      lon = geosphere::destPoint(p = c(longitude, latitude), b = b, d = d)[, 1]
    ) %>%
    dplyr::select(x = lon, y = lat, z = z)

  if (proj == "utm") {
    x <- footprint_c$x
    y <- footprint_c$y
    z <- footprint_c$z
    xy <- data.frame(ID = 1:length(x), X = x, Y = y, Z = z)
    sp::coordinates(xy) <- c("X", "Y")
    sp::proj4string(xy) <- sp::CRS("+proj=longlat +datum=WGS84")
    xy_proj <- sp::spTransform(xy, sp::CRS("+proj=utm +zone=18 ellps=WGS84"))
    footprint_r <- raster::rasterize(
      xy_proj,
      raster::raster(
        raster::extent(xy_proj), ncols = grid, nrows = grid,
        sp::CRS("+proj=utm +zone=18 ellps=WGS84")
      ),
      "Z", fun = mean
    )
  }

  # Convert footprint to projected raster
  if (proj == "longlat") {
    footprint_r <- raster::rasterize(
      footprint_c[, 1:2],
      raster::raster(
        raster::extent(footprint_c[, 1:2]), ncol = grid, nrow = grid,
        sp::CRS("+proj=longlat +datum=NAD83")
      ),
      footprint_c[,3], fun = mean
    )
  }
  # Convert raster to contour lines, return SpatialLinesDataFrame
  raster::rasterToContour(footprint_r, levels = percents)
}

percent_footprint <- function(x, percent = 100, weights = FALSE) {
  # Might deprecate this - way too clunky, and fp_percent should do the job
  rotate <- function(x) t(apply(x, 2, rev))
  fpx <- x$x
  fpy <- x$y
  fpz <- x$z
  if (weights == TRUE) fpw <- x$fpw else fpw <- matrix(NA, nrow(fpx), ncol(fpx))
  df <- data.frame(
    x = numeric(length(fpx)),
    y = numeric(length(fpy)),
    z = numeric(length(fpz)),
    zp = numeric(length(cumsum(fpz))),
    w = numeric(length(fpw)),
    wp = numeric(length(cumsum(fpw)))
  )
  df$z <- as.vector(fpz)
  df$zp <- as.vector(df$zp)
  df$x <- as.vector(fpx)
  df$y <- as.vector(fpy)
  df$w <- as.vector(df$w)
  df$wp <- as.vector(df$wp)
  t <- df[order(-df$z), ]
  t$cumz <- cumsum(t$z) * 100
  m <- xtabs(cumz ~ x + y, data = t, addNA = TRUE)
  msub <- ifelse(m <= percent, m, NA)
  f <- x
  f$z <- rotate(rotate(rotate(unclass(msub))))
  f <- as.data.frame(as.table(f$z))
  f$y <- as.numeric(as.character(f$y))
  f$x <- as.numeric(as.character(f$x))
  names(f)[3] <- "z"
  f
}

calculate_cover <- function(x, image, type) {
  if (fp_is_empty(x)) return(NA)
  fp_sub <- image %>%
    raster::extract(x, df = TRUE, factors = TRUE)
  tot <- nrow(fp_sub)
  patch <- fp_sub %>%
    tidyr::gather(key, value, -ID) %>%
    dplyr::filter(value == rlang::as_name(type)) %>%
    # Count number of cells in each patch type
    dplyr::tally() %>%
    dplyr::pull(n)
  # Return proportion
  patch / tot
}

calculate_overlap <- function(area, ftp) {
  # 0-buffers to prevent topology issues
  area <- rgeos::gBuffer(area, width = 0)
  ftp <- rgeos::gBuffer(ftp, width = 0)
  # calculate footprint area
  ftp_area <- rgeos::gArea(ftp)
  # create polygon of overlap
  overlap <- raster::intersect(ftp, area)
  # calculate overlap area
  overlap_area <- rgeos::gArea(overlap)
  # calculate proportion of ftp within site area
  p_overlap <- overlap_area / ftp_area
  return(p_overlap)
}

ndwi <- function(x) {
  # Takes RGB+NI image and returns normalized difference water index
  num <- raster::subset(x, 2) - raster::subset(x, 4) # g - ni
  den <- raster::subset(x, 2) + raster::subset(x, 4) # g + ni
  num / den
}

ndvi <- function(x) {
  # Takes RGB+NI image and returns normalized difference vegetation index
  num <- raster::subset(x, 4) - raster::subset(x, 1) # ni - r
  den <- raster::subset(x, 4) + raster::subset(x, 1) # ni + r
  num / den
}