hydrogeo.R

library(sf)

library(deldir)
require(igraph)
# require(gstat)
require(smoothr)
require(tripack)
require(spatstat)
library(fields)
library(gstat)
require(lwgeom)
require(terra)




#' Densify the polyline knots
#'
#' Produce almost equispaced knots on a polyline.
#' @param x [\code{matrix(m, 2)}|\code{sf}] Either a mx2 matrix whose rows 
#'          contain the knot coordinates of the polyline or a \code{sf} object.
#' @param dx [\code{dx(1)}] Distance between the knots.
#' @param method [\code{character(1)}] Interpolation method used to densify 
#'        the polyline.
#' @return [\code{matrix(n, 2)}|\code{sf}] Either a matrix of coordinates or a 
#'          \code{sf} object of the n equidistant knots.
#' @export
densify <- function(x, dx, method = c("linear", "nearest", 
                                      "pchip", "cubic", "spline")){
  if(inherits(x, "sf")){
    x <- .densifySf(x = x, dx = dx, method = method)
  }else if(is.matrix(x) || is.data.frame(x)){
    x <- .densifyMatrix(x = x, dx = dx, method = method)
  }
  return(x)
}

.densifyMatrix <- function(x, dx, method = c("linear", "nearest", 
                                             "pchip", "cubic", "spline")){
  # tst <- duplicated(x)
  dxy <- posLine(x)
  tst <- diff(dxy) == 0
  x <- x[!tst, ]
  dxy <- dxy[!tst]
  dxyi <- seq(min(dxy), max(dxy), by = dx)
  xint <- matrix(nrow = length(dxyi), ncol = ncol(x))
  # iterate over the dimensions of x (x, y, z, ...)
  for(i in 1:ncol(x)){
    xint[, i] <- signal::interp1(dxy, x[, i], dxyi, method = method )
  }
  return(xint)    
}

# posLine <- function(x){
#   cumsum(c(0, sqrt(apply(diff(x)^2, 1, sum))))
# }

.densifySf <- function(xsf, dx, method = c("linear", "nearest", 
                                           "pchip", "cubic", "spline")){
  x <- st_coordinates(xsf)
  xdens <- .densifyMatrix(x, method = method, dx = dx)
  xdens_pts <- st_as_sf(as.data.frame(xdens[, 1:2]), crs = st_crs(xsf), coords = 1:2)
  xdens_pts <- xdens_pts[!duplicated(xdens_pts),]
  # isNa <- apply(xdens_pts, 1, anyNA)
  xy <- st_cast(st_combine(xdens_pts), 
                as.character(st_geometry_type(xsf)))
  return(xy)
}

# return coordinates of mid-of-arcs of delaunay triangle arcs
# x = output from tripack::tri.mesh()
.getMidArcs <- function(tti, d){
  xyi <- cbind(d$x[tti], d$y[tti])
  (xyi + xyi[c(2,3,1),])/2
}

getMidArcs <- function(x){
  tt <- tripack::triangles(x)
  tst <- t(apply(tt[,1:3], 1, .getMidArcs, x))
  md <- rbind(tst[, c(1, 4)], tst[, c(2, 5)], tst[, c(3, 6)])
  md[ !duplicated(md), ]
}


#' Centerline of polygon
#' 
#' Return the centerline of a polygon. The algorithm is based on this post:
#' https://stackoverflow.com/questions/9595117/identify-a-linear-feature-on-a-raster-map-and-return-a-linear-shape-object-using/9643004#9643004
#' The polygon can be optionnally densified (introduce equidistant knots).
#' @param x [\code{sf}] Polygon
#' @param method [\code{character(1)}] Interpolation method used to densify 
#'        the polyline. Only used if \code{dx} is not \code{NULL}.
#' @param dx [\code{dx(1)}|\code{NULL}] Distance between the knots. If 
#'           \code{dx = NULL}, the polygon \code{x} will not be densified.
#' @return [\code{sf}] A polyline (class \code{sf})
# if dx != NULL, density the sf object x
# compute delaunay triangulation
# compute mid points of arcs (and remove duplicates)
# compute a graph tree to sort the points
# source: https://stackoverflow.com/questions/9595117/identify-a-linear-feature-on-a-raster-map-and-return-a-linear-shape-object-using/9643004#9643004
centerline <- function(x, method = c("linear", "nearest", 
                                     "pchip", "cubic", "spline"), dx = NULL){
  if(!is.null(dx)){
    x <- .densifySf(x, method = method, dx = dx)
  }
  xy <- st_coordinates(x)
  # xy <- xy[- which(duplicated(xy)), ]
  
  ### compute triangulation
  d <- tripack::tri.mesh(xy[,1],xy[,2], duplicate = "remove")
  #  not working -> always throwing an error d=deldir(xy[,1],xy[,2], eps = 1e-5)
  
  ### find midpoints of triangle sides
  mids <- getMidArcs(d)
  
  ### get points that are inside the polygon 
  mids_sf <- st_as_sf(data.frame(mids), coords = 1:2, crs = st_crs(x))
  # small buffer to remove the points that lay on the polyogn
  x_buf <- st_buffer(x, -1)
  sel <- st_intersection(mids_sf, x_buf)
  
  ### now build a minimum spanning tree weighted on the distance
  G <- igraph::graph.adjacency(as.matrix(dist(st_coordinates(sel))), weighted = TRUE, mode = "upper")
  minSpaTree <- igraph::minimum.spanning.tree(G, weighted = TRUE)
  
  ### get a diameter
  path <- igraph::get.diameter(minSpaTree)
  
  if(length(path)!= igraph::vcount(minSpaTree)){
    stop("Path not linear - try increasing dens parameter")
  }
  
  ### path should be the sequence of points in order
  # AA <- list(sel=sel[path+1,], tree = minSpaTree)
  # 
  # centerline <- st_coordinates(AA$sel)
  # centerline <- centerline[!apply(centerline, 1, anyNA), ]
  
  # selp <- sel[path + 1,]
  selp <- sel[path ,]
  tst <- apply((st_coordinates(selp)), 1, anyNA)
  selpath <- st_cast(st_combine( selp[!tst, ]), "LINESTRING")
  return(selpath)
}

#' Relative position of a knots along a polyline
#' 
#' Relative position of a knots along a polyline
#' @param [\code{matrix(m,n)}] A matrix whose rows store the coordinates. If
#'        \code{n = 2}, 2D distances are compute. If \code{n = 3}, 3D distances
#'        are computed. Etc.
#' @return [\numeric(m)] The relative positions of the m knots along the 
#'          polyline.
posLine <- function(x){
  if(inherits(x, "sfc")) x <- st_coordinates(x)[, 1:2]
  x <- as.matrix(x)
  c(0, cumsum(sqrt(rowSums(diff(x)^2))))
}


#' Points perpendicular to a polyline
#'
#' Compute oriented transects parallel to a polyline (one transect per points), 
#' the knots on the transects lie on the perpendicular to the polyline knots.
#' @param xy [\code{matrix(n, 2)}] Coordinates matrix with rows corresponding to the
#'           2D point coordinates (x and y positions).
#' @param d [\code{numeric(m)}] \code{m}-length vector defining the distance 
#'          between the transects (e.g., the perpendicular points) and the 
#'          polyline. If m > 1 several transects are returned.
#' @return [\code{list}] A list with elements x and y of dimension (n, m).
#' @export
# TODO: use projective geometry
perpPoints <- function(xy, d){
  xy2 <- xy[c(1,1:nrow(xy), nrow(xy)),]
  xlat <- matrix(nrow = nrow(xy), ncol = length(d))
  ylat <- matrix(nrow = nrow(xy), ncol = length(d))
  for(i in 2:(nrow(xy2) - 1)){
    if( xy2[i - 1, 2] == xy2[i + 1, 2]){
      xlat[i-1, ] <- xy2[i, 1]
      ylat[i-1, ] <- xy2[i, 2] + d
    }else if(xy2[i - 1, 1] == xy2[i + 1, 1] ){
      xlat[i-1, ] = xy2[i, 1] + d
      ylat[i-1, ] = xy2[i, 2]
    }else{
      #get the slope of the line
      m <- ((xy2[i - 1, 2] - xy2[i + 1, 2])/(xy2[i - 1, 1] - xy2[i + 1, 1]))
      #get the negative reciprocal, 
      n_m <- -1/m
      sng <- sign(xy2[i + 1, 1] - xy2[i - 1, 1])
      DD <- d / sqrt( n_m^2 + 1)
      if( m < 0){
        DD <- -DD
      }
      if(sng > 0){
        xlat[i-1, ] = xy2[i, 1] +  DD
        ylat[i-1, ] = xy2[i, 2] + n_m * DD
      }else{
        xlat[i-1, ] = xy2[i, 1] - DD
        ylat[i-1, ] = xy2[i, 2] - n_m * DD
      }
    }  
  }
  return(list(x = xlat, y = ylat))
}


smoothThisChannel <- function(u, d, zmin){
  umin <- u[1]
  if(!is.null(zmin)) u[ u < zmin] <- NA
  for(i in seq_along(u)){
    if(!is.na(u[i])){
      if(u[i] > umin){
        u[i] <- NA
      }else{
        umin <- u[i]
      }
    }
  }
  una <- is.na(u)
  u[una] <- signal::interp1(x  = d[!una],
                            y  = u[!una],
                            xi = d[una],
                            method = "linear",
                            extrap = TRUE)
  return(u)
}



rasterizePoints <- function(x, res, att, method = c("tps", "idw"), 
                            xmin, xmax, ymin, ymax, clipch = FALSE){
  if(missing(att)) stop("You must specify 'att'!!!")
  method <- match.arg(method, c("tps", "idw"))
  if(!inherits(x, "sf")){
    if(is.matrix(x)) x <- data.frame(x)
    x <- st_as_sf(x, coords = 1:2)
  }
  
  if(missing(xmin) && missing(xmax) && missing(ymin) && missing(ymax)){
    xmin = st_bbox(x)[1]
    ymin = st_bbox(x)[2]
    xmax = st_bbox(x)[3]
    ymax = st_bbox(x)[4]
  }else{
    x <- st_crop(x, xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax)
  }
  if(missing(res)){
    res <- max(c(xmax - xmin, ymax - ymin))/100
  }
  r <- raster::raster(xmn = xmin,
                      xmx = xmax,
                      ymn = ymin,
                      ymx = ymax, 
                      resolution = res)
  
  x_df <- data.frame(st_coordinates(x)[, 1:2], v = x[[att]])
  names(x_df)[1] <- "x"
  names(x_df)[2] <- "y"
  i <- !is.na(x_df$v)
  x_df <- x_df[i, ]
  
  if(method == "idw"){
    m <- gstat(formula=v~1, locations= ~x+y, data=x_df, set = list(idp = 1), nmin = 5, nmax = 10)
    
  }else if(method == "tps"){
    m <- Tps(x_df[, c("x", "y")], x_df$v)
  }
  r <- interpolate(r, m)
  if(isTRUE(clipch)){
    ch <- st_convex_hull(x)
    r <- crop(r, as(ch, "Spatial"))
  }
  return(r)
  
}