R/splitreach.rivernet.r

In rivernet: Read, Analyze and Plot River Networks

splitreach <- function(net,reachind,fract, ...) UseMethod("splitreach")


splitreach.rivernet <- function(net,reachind,fract,...)
{
  # check input:
  
  if ( fract <= 0 | fract >= 1 )
  {
    print("*** splitreach.rivernet: argument fract is not between zero and one (exclusive zero and one)")
    return(net)
  }
  if ( reachind < 1 | reachind > length(net$reaches) )
  {
    print("*** splitreach.rivernet: argument reachind must be between 1 and the number of reaches")
  }
  if ( is.factor(net$attrib.reach$Reach_ID) ) net$attrib.reach$Reach_ID <- as.character(net$attrib.reach$Reach_ID)
  
  # calculate split:
  
  x <- net$reaches[[reachind]]$x
  y <- net$reaches[[reachind]]$y
  z <- net$reaches[[reachind]]$z
  length.inc <- rep(NA,length(x))
  length.inc[1] <- 0
  for ( i in 2:length(length.inc) ) length.inc[i] <- length.inc[i-1] + sqrt((x[i]-x[i-1])^2+(y[i]-y[i-1])^2) 
  length1 <- fract*length.inc[length(length.inc)]
  length2 <- (1-fract)*length.inc[length(length.inc)]
  ind <- max(which(length.inc<length1))
  x.split <- ((length.inc[ind+1]-length1)*x[ind]+(length1-length.inc[ind])*x[ind+1])/(length.inc[ind+1]-length.inc[ind])
  y.split <- ((length.inc[ind+1]-length1)*y[ind]+(length1-length.inc[ind])*y[ind+1])/(length.inc[ind+1]-length.inc[ind])
  z.split <- ((length.inc[ind+1]-length1)*z[ind]+(length1-length.inc[ind])*z[ind+1])/(length.inc[ind+1]-length.inc[ind])
  x1 <- x[1:(ind+1)]; x1[ind+1] <- x.split
  y1 <- y[1:(ind+1)]; y1[ind+1] <- y.split
  z1 <- z[1:(ind+1)]; z1[ind+1] <- z.split
  x2 <- x[ind:length(x)]; x2[1] <- x.split
  y2 <- y[ind:length(x)]; y2[1] <- y.split
  z2 <- z[ind:length(x)]; z2[1] <- z.split
  
  # create new node:
  
  n.nodes <- length(net$nodes)
  net$nodes[[n.nodes+1]] <- list()
  net$nodes[[n.nodes+1]]$x <- x.split
  net$nodes[[n.nodes+1]]$y <- y.split
  net$attrib.node <- rbind(net$attrib.node,NA)
  net$attrib.node[n.nodes+1,"Node"] <- max(net$attrib.node$Node,na.rm=TRUE) + 1
  net$attrib.node[n.nodes+1,"Node_ID"] <- paste("Split",
                                                net$attrib.node$Node_ID[net$attrib.reach$node_start[reachind]],
                                                net$attrib.node$Node_ID[net$attrib.reach$node_end[reachind]],
                                                sep="_")
  net$attrib.node[n.nodes+1,"x"] <- x.split
  net$attrib.node[n.nodes+1,"y"] <- y.split
  
  # modify reach and create new reach:
  
  n.reaches <- length(net$reaches)
  net$attrib.reach <- rbind(net$attrib.reach,net$attrib.reach[reachind,])
  ID <- net$attrib.reach$Reach_ID[reachind]
  
  net$reaches[[reachind]]$n <- length(x1)
  net$reaches[[reachind]]$x <- x1
  net$reaches[[reachind]]$y <- y1
  net$reaches[[reachind]]$z <- z1
  net$reaches[[reachind]]$length <- length1
  net$attrib.reach$x_end[reachind] <- x.split
  net$attrib.reach$y_end[reachind] <- y.split
  net$attrib.reach$z_end[reachind] <- z.split
  net$attrib.reach$node_end[reachind] <- n.nodes+1
  net$attrib.reach$length[reachind] <- length1
  net$attrib.reach$Reach_ID[reachind] <- paste(ID,"1",sep="_")
  
  net$reaches[[n.reaches+1]] <- list()
  net$reaches[[n.reaches+1]]$n <- length(x2)
  net$reaches[[n.reaches+1]]$x <- x2
  net$reaches[[n.reaches+1]]$y <- y2
  net$reaches[[n.reaches+1]]$z <- z2
  net$reaches[[n.reaches+1]]$length <- length2
  net$attrib.reach$x_start[n.reaches+1] <- x.split
  net$attrib.reach$y_start[n.reaches+1] <- y.split
  net$attrib.reach$z_start[n.reaches+1] <- z.split
  net$attrib.reach$node_start[n.reaches+1] <- n.nodes+1
  net$attrib.reach$length[n.reaches+1] <- length2
  net$attrib.reach$Reach[n.reaches+1] <- max(net$attrib.reach$Reach) + 1
  net$attrib.reach$Reach_ID[n.reaches+1] <- paste(ID,"2",sep="_")
  
  # update variables generated by analyze.rivernet if this function was run already:
  
  if ( !is.null(net$paths) )
  {
    net$attrib.reach$n_end[reachind] <- 1
    net$attrib.reach$n_start[n.reaches+1] <- 1
    net$attrib.reach$endreach[reachind] <- net$attrib.reach$n_start[reachind]==0 | net$attrib.reach$n_end[reachind]==0
    net$attrib.reach$endreach[n.reaches+1] <- net$attrib.reach$n_start[n.reaches+1]==0 | net$attrib.reach$n_end[n.reaches+1]==0
    if ( net$attrib.reach$downstream[reachind] )
    {
      net$nodes[[n.nodes+1]]$from_reach <- reachind
      net$nodes[[n.nodes+1]]$to_reach   <- n.reaches+1
      net$nodes[[net$attrib.reach$node_end[n.reaches+1]]]$from_reach <- n.reaches+1
      net$reaches[[n.reaches+1]]$from_node <- n.nodes+1
      net$reaches[[n.reaches+1]]$to_node <- net$reaches[[reachind]]$to_node
      net$reaches[[reachind]]$to_node <- n.nodes+1
      net$attrib.reach$outlet[reachind] <- FALSE
      net$attrib.reach$headwater[n.reaches+1] <- FALSE
      net$attrib.reach$reach_down[reachind] <- n.reaches+1
      for ( i in 1:length(net$paths) )
      {
        ind <- match(reachind,net$paths[[i]])
        {
          if ( !is.na(ind) ) 
          {
            if ( ind == length(net$paths[[i]]) )
            {
              net$paths[[i]] <- c(net$paths[[i]],n.reaches+1)
            }
            else
            {
              net$paths[[i]] <- c(net$paths[[i]][1:ind],n.reaches+1,net$paths[[i]][(ind+1):length(net$paths[[i]])])
            }
          }
        }
      }
    }
    else
    {
      net$nodes[[n.nodes+1]]$from_reach <- n.reaches+1
      net$nodes[[n.nodes+1]]$to_reach   <- reachind
      net$nodes[[net$attrib.reach$node_end[n.reaches+1]]]$to_reach <- n.reaches+1
      net$reaches[[n.reaches+1]]$to_node <- n.nodes+1
      net$reaches[[n.reaches+1]]$from_node <- net$reaches[[reachind]]$from_node
      net$reaches[[reachind]]$from_node <- n.nodes+1
      net$attrib.reach$outlet[n.reaches+1] <- FALSE
      net$attrib.reach$headwater[reachind] <- FALSE
      net$attrib.reach$reach_down[n.reaches+1] <- reachind
      for ( i in 1:length(net$paths) )
      {
        ind <- match(reachind,net$paths[[i]])
        {
          if ( !is.na(ind) ) 
          {
            if ( ind == 1 )
            {
              net$paths[[i]] <- c(n.reaches+1,net$paths[[i]])
            }
            else
            {
              net$paths[[i]] <- c(net$paths[[i]][1:(ind-1)],n.reaches+1,net$paths[[i]][ind:length(net$paths[[i]])])
            }
          }
        }
      }
    }
  }
  
  return(net)
}