R/riverdist_1.R

#' Pythagorean Distance
#' @description Pythagorean distance between two points.  Called internally.
#' @param p1 X-Y coordinates of point 1
#' @param p2 X-Y coordinates of point 2
#' @return Distance (numeric)
#' @author Matt Tyers
#' @examples
#' point1 <- c(1,3)
#' point2 <- c(4,7)
#' 
#' pdist(point1,point2)
#' @export
pdist <- function(p1,p2) {
  dist <- sqrt((p1[1]-p2[1])^2 + (p1[2]-p2[2])^2)
  return(dist)
}

#' Total Pythagorean Distance
#' @description Total Pythagorean distance of a sequence of points.  Called internally.
#' @param xy A matrix of X-Y coordinates of the sequence of points.
#' @return Distance (numeric)
#' @author Matt Tyers
#' @examples
#' points <- matrix(c(1:10), nrow=5, ncol=2, byrow=FALSE)
#' 
#' pdisttot(xy=points)
#' @export
pdisttot <- function(xy) {
  n <- dim(xy)[1]
  if(n==1) dist <- 0
  if(n==2) dist <- pdist(xy[1,],xy[2,])
  if(n>2) dist <- sqrt(((xy[1:(n-1),1] - xy[2:n,1])^2) + ((xy[1:(n-1),2] - xy[2:n,2])^2))
  return(sum(dist))
}


#' Add Cumulative Distance to a River Network
#' @description Adds a vector of cumulative distances to a river network.  Called internally.
#' @param rivers The river network object to use.
#' @return Returns an object of class \code{"rivernetwork"} containing all
#'   spatial and topological information.  See \link{rivernetwork-class}.
#' @author Matt Tyers
#' @examples
#' Gulk1 <- addcumuldist(rivers=Gulk)
#' @export
addcumuldist <- function(rivers) {
  cumuldist <- list()
  for(i in 1:length(rivers$lines)) {
    xy <- rivers$lines[[i]]
    n <- dim(xy)[1]
    cumuldist[[i]] <- c(0,cumsum(sqrt(((xy[1:(n-1),1] - xy[2:n,1])^2) + ((xy[1:(n-1),2] - xy[2:n,2])^2))))
  }
  rivers$cumuldist <- cumuldist
  return(rivers)
}


#' Calculate the Connectivity Matrix for a River Network
#' @description Calculates the connectivity matrix for a river network, during import and editing.  Called internally.
#' @param lines A list of coordinate matrices, each corresponding to a line segment.
#' @param tolerance The spatial tolerance for establishing connectivity.
#' @return A matrix with topological information.  See the \code{$connections} element of the \link{rivernetwork-class}.
#' @author Matt Tyers
#' @examples
#' Gulk_connections <- calculateconnections(lines=Gulk$lines, tolerance=Gulk$tolerance)
#' @export
calculateconnections <- function(lines,tolerance) {
  length <- length(lines)
  # defining a connectivity matrix...
  # connection type 1: beginning - beginning
  # connection type 2: beginning - end
  # connection type 3: end - beginning
  # connection type 4: end - end
  # connection type 5: beginning - beginning and end - end
  # connection type 6: beginning - end and end - beginning
  connections <- matrix(NA,nrow=length,ncol=length)
  begmat <- matrix(unlist(sapply(lines,function(xy) xy[1,],simplify=F)),ncol=2,byrow=T)
  endmat <- matrix(unlist(sapply(lines,function(xy) xy[nrow(xy),],simplify=F)),ncol=2,byrow=T)
  pdist2 <- function(p1,p2mat) {
    dist <- sqrt((p1[1]-p2mat[,1])^2 + (p1[2]-p2mat[,2])^2)
    return(dist)
  }  
  
  for(i in 1:length) {
    mat1 <- pdist2(begmat[i,],begmat)
    mat2 <- pdist2(begmat[i,],endmat)
    mat3 <- pdist2(endmat[i,],begmat)
    mat4 <- pdist2(endmat[i,],endmat)
    connections[i,mat1<tolerance] <- 1
    connections[i,mat2<tolerance] <- 2
    connections[i,mat3<tolerance] <- 3
    connections[i,mat4<tolerance] <- 4
    connections[i,(mat1<tolerance & mat4<tolerance)] <- 5
    connections[i,(mat2<tolerance & mat3<tolerance)] <- 6
  }
  diag(connections) <- NA
  return(connections)
}


#' Create a River Network Object from a Shapefile
#' @description Uses \link[rgdal]{readOGR} in package 'rgdal' to read a river 
#'   shapefile, and establishes connectivity of segment endpoints based on 
#'   spatial proximity.
#' @param sp SpatialLinesDataFrame object. optional.
#' @param path File path, default is the current working directory.
#' @param layer Name of the shapefile, without the .shp extension.
#' @param tolerance Snapping tolerance of segment endpoints to determine 
#'   connectivity.  Default is 100, therefore care should be exercised when 
#'   working with larger units of distance, such as km.
#' @param reproject A valid Proj.4 projection string, if the shapefile is to be 
#'   re-projected.  Re-projection is done using \link[sp]{spTransform} in 
#'   package 'sp'.
#' @param supplyprojection A valid Proj.4 projection string, if the input shapefile does not have the projection information attached.
#' @return Returns an object of class \code{"rivernetwork"} containing all
#'   spatial and topological information.  See \link{rivernetwork-class}.
#' @note Since distance can only be calculated using projected coordinates, 
#'   \code{line2network()} will generate an error if a non-projected input 
#'   shapefile is detected.  To resolve this, the shapefile can be re-projected 
#'   in a GIS environment, or using \code{reproject=}, shown in the second 
#'   example below.
#' @author Matt Tyers, Joseph Stachelek
#' @importFrom rgdal readOGR
#' @importFrom sp is.projected
#' @importFrom sp CRS
#' @importFrom sp spTransform
#' @examples 
#' filepath <- system.file("extdata", package="riverdist")
#' 
#' Gulk_UTM5 <- line2network(path=filepath, layer="Gulk_UTM5")
#' plot(Gulk_UTM5)
#' 
#' ## Reading directly from a SpatialLinesDataFrame object
#' 
#' sp <- rgdal::readOGR(dsn = filepath, layer = "Gulk_UTM5", verbose = FALSE)
#' Gulk_UTM5 <- line2network(sp)
#' plot(Gulk_UTM5)
#' 
#' ## Re-projecting in Alaska Albers Equal Area projection:
#' 
#' AKalbers <- "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 
#'     +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
#'     
#' Gulk_AKalbers <- line2network(path=filepath, layer="Gulk_UTM5", reproject=AKalbers)
#' plot(Gulk_AKalbers)
#' 
#' @export
line2network <- function(sp = NA, path=".", layer = NA, tolerance=100, 
                         reproject=NULL, supplyprojection=NULL) {
  
  if(suppressWarnings(is.na(sp))){
    sp <- suppressWarnings(rgdal::readOGR(dsn = path, layer = layer, verbose = F))   }
  
  if(class(sp)!="SpatialLinesDataFrame"){ 
    stop("Specified shapefile is not a linear feature.")
  }
  
  if(is.na(sp@proj4string@projargs) & !is.null(supplyprojection)){ 
    sp@proj4string@projargs <- supplyprojection
  }
    
  if(is.na(sp@proj4string@projargs)){ 
    stop("Shapefile projection information is missing.  Use supplyprojection= to specify a Proj.4 projection to use.  If the input shapefile is in WGS84 geographic (long-lat) coordinates, this will be +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 (in double-quotes).  If so, it must also be reprojected using reproject=.")
  }
    
  proj4 <- strsplit(sp@proj4string@projargs,split=" ")
  projected <- sp::is.projected(sp)
  if(is.null(reproject) & !projected) stop("Distances can only be computed from a projected coordinate system.  Use reproject= to specify a Proj.4 projection to use.")
  
  if(!is.null(reproject)) {
    sp <- sp::spTransform(sp,sp::CRS(reproject))
    proj4 <- strsplit(sp@proj4string@projargs,split=" ")
  }
  
  units <- "unknown"
  for(i in 1:length(proj4[[1]])) {
    if(proj4[[1]][i]!="") {
      proj4arg <- strsplit(proj4[[1]][i],split="=")
      if(proj4arg[[1]][1]=="+units") {
        units <- proj4arg[[1]][2]
        cat('\n',"Units:",proj4arg[[1]][2],'\n')
      }
    }
  }
  
  if(length(sp@lines) > 1) {
    sp_line <- NA
    sp_seg <- NA
    lines <- list()
    j<-1
    for(i in 1:length(sp@lines)) {
      for(k in 1:length(sp@lines[i][[1]]@Lines)) {
        lines[[j]] <- sp@lines[i][[1]]@Lines[[k]]@coords
        sp_line[j] <- i
        sp_seg[j] <- k
        j<-j+1
      }
    }
  }
  if(length(sp@lines) == 1) {
    lines <- sp@lines[1][[1]]@Lines    # extracting just a list of lines and coordinates
    
    length <- length(lines) # number of line segments
    
    lines.new <- list()
    for(i in 1:length) {
      lines.new[[i]] <- lines[[i]]@coords
    }
    lines <- lines.new 
    sp_line <- rep(1,length)
    sp_seg <- 1:length
  }
  length <- length(lines)
  
  rivID <- 1:length
  lineID <- data.frame(rivID,sp_line,sp_seg)
  
  connections <- calculateconnections(lines=lines, tolerance=tolerance)
  
  if(any(connections %in% 5:6)) braided <- TRUE
  
  # making a vector of total segment lengths
  lengths <- rep(NA,length)
  for(i in 1:length) {
    lengths[i] <- pdisttot(lines[[i]])
  }
  
  names <- rep(NA,length)
  mouth.seg <- NA
  mouth.vert <- NA
  mouth <- list(mouth.seg,mouth.vert)
  names(mouth) <- c("mouth.seg","mouth.vert")
  sequenced <- FALSE
  braided <- NA
  
  cumuldist <- list()
  for(i in 1:length) {
    xy <- lines[[i]]
    n <- dim(xy)[1]
    cumuldist[[i]] <- c(0,cumsum(sqrt(((xy[1:(n-1),1] - xy[2:n,1])^2) + ((xy[1:(n-1),2] - xy[2:n,2])^2))))
  }
  
  out.names <- c("sp","lineID","lines","connections","lengths","names","mouth","sequenced","tolerance","units","braided","cumuldist")
  out <- list(sp,lineID,lines,connections,lengths,names,mouth,sequenced,tolerance,units,braided,cumuldist)
  names(out) <- out.names
  class(out) <- "rivernetwork"
  
  length1 <- length(out$lengths)
  suppressMessages(out <- removeduplicates(out))
  length2 <- length(out$lengths)
  if(length2<length1) cat('\n',"Removed",length1-length2,"duplicate segments.",'\n')
  suppressMessages(out <- removemicrosegs(out))
  length3 <- length(out$lengths)
  if(length3<length2) cat('\n',"Removed",length2-length3,"segments with lengths shorter than the connectivity tolerance.",'\n')
  
  return(out)
}


#' Convert a Point Shapefile to River Locations
#' @description This function reads a point shapefile and determines the closest
#'   vertex in the river network to each point of XY data, returning a data
#'   frame with river locations, defined as segment numbers and vertex
#'   numbers, along with the data table read from the input shapefile.
#' @param path File path, default is the current working directory.
#' @param layer Name of the shapefile, without the .shp extension.
#' @param rivers The river network object to use.
#' @return A data frame of river locations, with segment numbers in
#'   \code{$seg}, vertex numbers in \code{$vert}, snapping distances in \code{$snapdist}, and the remaining columns
#'   corresponding to the data table in the input point shapefile.
#' @author Matt Tyers
#' @note If the input shapefile is detected to be in a different projection than
#'   the river network, the input shapefile will be re-projected before
#'   conversion to river locations.
#' @importFrom rgdal readOGR
#' @importFrom sp proj4string
#' @importFrom sp CRS
#' @importFrom sp spTransform
#' @examples 
#' filepath <- system.file("extdata", package="riverdist")
#' 
#' fakefish_UTM5 <- pointshp2segvert(path=filepath, layer="fakefish_UTM5", rivers=Gulk)
#' head(fakefish_UTM5)
#' 
#' plot(x=Gulk)
#' points(fakefish_UTM5$x, fakefish_UTM5$y)
#' riverpoints(seg=fakefish_UTM5$seg, vert=fakefish_UTM5$vert, rivers=Gulk, pch=16, col=2)
#' 
#' @export
pointshp2segvert <- function(path=".",layer,rivers) {
  shp <- rgdal::readOGR(dsn=path,layer=layer,pointDropZ=TRUE)
  if(sp::proj4string(shp) != sp::proj4string(rivers$sp)) {
    cat('\n',"Point projection detected as different from river network.  Re-projecting points before snapping to river network...")
    projection <- sp::CRS(sp::proj4string(rivers$sp))
    shp <- sp::spTransform(shp,projection) ##
  }
  segvert <- xy2segvert(x=shp@coords[,1],y=shp@coords[,2],rivers=rivers)
  outdf <- cbind(segvert,shp@data)
  return(outdf)
}


#' Plotting a River Network
#' @description S3 plotting method for the \link{rivernetwork-class}.  Produces a map of all river segments of a river network object.
#' @aliases mapriver
#' @param x The river network object to plot
#' @param segmentnum Whether or not to plot segment numbers (defaults to TRUE)
#' @param offset Whether to offset segment numbers from lines (defaults to TRUE)
#' @param lwd Line width
#' @param cex Global character expansion factor for plotting
#' @param scale Whether or not to give x- and y-axes the same scale
#' @param color How to differentiate segments.  If \code{color==TRUE} (default),
#'   segments will be drawn in solid lines with differing colors.  If
#'   \code{color==FALSE}, segments will be drawn in the same color with differing line
#'   types.
#' @param empty Creates an empty plot if set to \code{TRUE}.  Suppresses differentiation by line type if \code{color==FALSE}, and suppresses segment number labels.  Defaults to \code{FALSE}.
#' @param linecol Line color to use if \code{empty} is \code{TRUE} or \code{color} is \code{FALSE}.  Defaults to black.
#' @param xlab Label for X-axis (defaults to "")
#' @param ylab Label for Y-axis (defaults to "")
#' @param ... Additional plotting arguments (see \link[graphics]{par})
#' @author Matt Tyers
#' @note This function is intended to provide basic visual checks for the user,
#'   not for any real mapping.
#' @examples 
#' data(Gulk)
#' plot(x=Gulk)
#' @method plot rivernetwork
#' @importFrom grDevices rgb
#' @importFrom graphics plot
#' @importFrom graphics par
#' @importFrom graphics text
#' @importFrom stats cor
#' @importFrom graphics axTicks
#' @export
plot.rivernetwork <- function(x,segmentnum=TRUE,offset=TRUE,lwd=1,cex=.6,scale=TRUE,color=TRUE,empty=FALSE,linecol=1,xlab="",ylab="",...) {
  if(class(x)!="rivernetwork") stop("Argument 'x' must be of class 'rivernetwork'.  See help(line2network) for more information.")
  lines <- x$lines
  length <- length(lines)
  allx <- unlist(lapply(lines,FUN='[',TRUE,1))
  ally <- unlist(lapply(lines,FUN='[',TRUE,2))
  xmax <- max(allx, na.rm=T)
  xmin <- min(allx, na.rm=T)
  ymax <- max(ally, na.rm=T)
  ymin <- min(ally, na.rm=T)
  plot(c(xmin,xmax),c(ymin,ymax),col="white",cex.axis=.6,asp=1,xlab=xlab,ylab=ylab,...=...)
  parusr <- par("usr")
  midptx <- midpty <- NA
  if(color&!empty) {
    linecolors <- rgb((sin(1:length)+1)/2.3,(cos(7*1:length)+1)/2.3,(sin(3*(length-(1:length)))+1)/2.3)
    ltys <- rep(1,length)
  }
  if(!color&!empty) {
    linecolors <- rep(linecol,length)
    ltys <- 1:length
  }
  if(empty) {
    linecolors <- rep(linecol,length)
    ltys <- rep(1,length)
    segmentnum <- F
  }
  
  xtext <- ytext <- direction <- rep(NA,length)
  for(j in 1:length) {
    lines(lines[[j]],col=linecolors[j],lty=ltys[j],lwd=lwd)
    
    xplot <- lines[[j]][,1][lines[[j]][,1]>parusr[1] & lines[[j]][,1]<parusr[2] & lines[[j]][,2]>parusr[3] & lines[[j]][,2]<parusr[4]]
    yplot <- lines[[j]][,2][lines[[j]][,1]>parusr[1] & lines[[j]][,1]<parusr[2] & lines[[j]][,2]>parusr[3] & lines[[j]][,2]<parusr[4]]
    if(length(xplot)>0) midptx[j] <- mean(xplot)
    if(length(yplot)>0) midpty[j] <- mean(yplot)
    if(length(xplot)==0 | length(yplot)==0) midptx[j] <- midpty[j] <-NA
    middle <- floor(length(xplot)/2)
    directionj <- 3+(abs(xplot[floor(length(xplot)*.25)]-xplot[floor(length(xplot)*.75)]) < abs(yplot[floor(length(yplot)*.25)]-yplot[floor(length(yplot)*.75)]))
    direction[j] <- ifelse(length(directionj)==0,1,directionj)
    xtext[j] <- ifelse(length(xplot)>0,xplot[middle],NA)
    ytext[j] <- ifelse(length(yplot)>0,yplot[middle],NA)
  }
  if(segmentnum) {
    if(offset) text(xtext,ytext,labels=1:length,pos=direction,cex=cex,col=linecolors)
    else text(xtext,ytext,labels=1:length,cex=cex,col=linecolors)
  }
  if(scale) {
    if(length>1) corthing <- cor(midptx,midpty,use="complete.obs")
    if(length==1) corthing <- (lines[[1]][1,1]-lines[[1]][dim(lines[[1]])[1],1])*(lines[[1]][1,2]-lines[[1]][dim(lines[[1]])[1],2])
    axticks1 <- axTicks(1)
    if(is.na(corthing)) legendleft <- F   ######
    else legendleft <- corthing<=0   ######
    if(legendleft) {   ######
      if(length(axticks1)>5) {
        scalex <- axticks1[c(1,3)]
        labx <- axticks1[2]
      }
      if(length(axticks1)<=5) {
        scalex <- axticks1[c(1,2)]
        labx <- mean(axticks1[c(1,2)])
      }
      scaley <- axTicks(2)[1]
    }
    if(!legendleft) {   ######
      if(length(axticks1)>5) {
        scalex <- axticks1[c(length(axticks1)-2,length(axticks1))]
        labx <- axticks1[length(axticks1)-1]
      }
      if(length(axticks1)<=5) {
        scalex <- axticks1[c(length(axticks1)-1,length(axticks1))]
        labx <- mean(axticks1[c(length(axticks1)-1,length(axticks1))])
      }
      scaley <- axTicks(2)[1]
    }
    lines(scalex,rep(scaley,2))
    if(x$units=="m") text(labx,scaley,labels=paste((scalex[2]-scalex[1])/1000,"km"),pos=3,cex=cex)
    if(x$units!="m") text(labx,scaley,labels=paste((scalex[2]-scalex[1]),x$units),pos=3,cex=cex)
  }
}


scalebar <- function(rivers,cex=.6) {
  lines <- rivers$lines
  length <- length(lines)
  
  midptx <- midpty <- rep(NA,length)
  parusr <- par("usr")
  for(j in 1:length) {
    xplot <- lines[[j]][,1][lines[[j]][,1]>parusr[1] & lines[[j]][,1]<parusr[2] & lines[[j]][,2]>parusr[3] & lines[[j]][,2]<parusr[4]]
    yplot <- lines[[j]][,2][lines[[j]][,1]>parusr[1] & lines[[j]][,1]<parusr[2] & lines[[j]][,2]>parusr[3] & lines[[j]][,2]<parusr[4]]
    if(length(xplot)>0) midptx[j] <- mean(xplot)
    if(length(yplot)>0) midpty[j] <- mean(yplot)
  }
  
  if(length>1) corthing <- cor(midptx,midpty,use="complete.obs")
  if(length==1) corthing <- (lines[[1]][1,1]-lines[[1]][dim(lines[[1]])[1],1])*(lines[[1]][1,2]-lines[[1]][dim(lines[[1]])[1],2])
  axticks1 <- axTicks(1)
  if(corthing<=0) {
    if(length(axticks1)>5) {
      scalex <- axticks1[c(1,3)]
      labx <- axticks1[2]
    }
    if(length(axticks1)<=5) {
      scalex <- axticks1[c(1,2)]
      labx <- mean(axticks1[c(1,2)])
    }
    scaley <- axTicks(2)[1]
  }
  if(corthing>0) {
    if(length(axticks1)>5) {
      scalex <- axticks1[c(length(axticks1)-2,length(axticks1))]
      labx <- axticks1[length(axticks1)-1]
    }
    if(length(axticks1)<=5) {
      scalex <- axticks1[c(length(axticks1)-1,length(axticks1))]
      labx <- mean(axticks1[c(length(axticks1)-1,length(axticks1))])
    }
    scaley <- axTicks(2)[1]
  }
  lines(scalex,rep(scaley,2))
  if(rivers$units=="m") text(labx,scaley,labels=paste((scalex[2]-scalex[1])/1000,"km"),pos=3,cex=cex)
  if(rivers$units!="m") text(labx,scaley,labels=paste((scalex[2]-scalex[1]),rivers$units),pos=3,cex=cex)
}


# plotrivernetwork_OLD <- function(x,segmentnum=TRUE,offset=TRUE,lwd=1,cex=.6,scale=TRUE,color=TRUE,empty=FALSE,xlab="",ylab="",...) {
#   rivers <- x
#   if(class(rivers)!="rivernetwork") stop("Argument 'rivers' must be of class 'rivernetwork'.  See help(line2network) for more information.")
#   lines <- rivers$lines
#   length <- length(lines)
#   xmin <- min(lines[[1]][,1])
#   xmax <- max(lines[[1]][,1])
#   ymin <- min(lines[[1]][,2])
#   ymax <- max(lines[[1]][,2])
#   if(length>1) {
#     for(j in 2:length) {
#       if(min(lines[[j]][,1])<xmin) xmin <- min(lines[[j]][,1])
#       if(max(lines[[j]][,1])>xmax) xmax <- max(lines[[j]][,1])
#       if(min(lines[[j]][,2])<ymin) ymin <- min(lines[[j]][,2])
#       if(max(lines[[j]][,2])>ymax) ymax <- max(lines[[j]][,2])
#       # print(c(j,xmin,xmax,ymin,ymax))
#     }
#   }
#   plot(c(xmin,xmax),c(ymin,ymax),col="white",cex.axis=.6,asp=1,xlab=xlab,ylab=ylab,...=...)
#   midptx <- midpty <- NA
#   for(j in 1:length) {
#     if(!color&!empty) lines(lines[[j]],col=1,lty=j,lwd=lwd)
#     linecolor <- rgb((sin(j)+1)/2.3,(cos(7*j)+1)/2.3,(sin(3*(length-j))+1)/2.3)
#     if(color&!empty) lines(lines[[j]],col=linecolor,lty=1,lwd=lwd)
#     if(color&empty) lines(lines[[j]],col=1,lty=1,lwd=lwd)
#     linelength <- dim(lines[[j]])[1]
#     
#     xplot <- lines[[j]][,1][lines[[j]][,1]>par("usr")[1] & lines[[j]][,1]<par("usr")[2] & lines[[j]][,2]>par("usr")[3] & lines[[j]][,2]<par("usr")[4]]
#     yplot <- lines[[j]][,2][lines[[j]][,1]>par("usr")[1] & lines[[j]][,1]<par("usr")[2] & lines[[j]][,2]>par("usr")[3] & lines[[j]][,2]<par("usr")[4]]
#     if(length(xplot)>0) midptx[j] <- mean(xplot)
#     if(length(yplot)>0) midpty[j] <- mean(yplot)
#     if(length(xplot)==0 | length(yplot)==0) midptx[j] <- midpty[j] <-NA
#     middle <- floor(length(xplot)/2)
#     direction <- 3+(abs(xplot[floor(length(xplot)*.25)]-xplot[floor(length(xplot)*.75)]) < abs(yplot[floor(length(yplot)*.25)]-yplot[floor(length(yplot)*.75)]))
#     if(!offset) direction <- NULL
#     xtext <- ifelse(length(xplot)>0,xplot[middle],NA)
#     ytext <- ifelse(length(yplot)>0,yplot[middle],NA)
#     if(segmentnum&!empty) text(x=xtext,y=ytext,labels=j,pos=direction,cex=cex,col=ifelse(color,linecolor,1))
#   }
#   if(scale) {
#     if(length>1) corthing <- cor(midptx,midpty,use="complete.obs")
#     if(length==1) corthing <- (lines[[1]][1,1]-lines[[1]][dim(lines[[1]])[1],1])*(lines[[1]][1,2]-lines[[1]][dim(lines[[1]])[1],2])
#     if(corthing<=0) {
#       if(length(axTicks(1))>5) {
#         scalex <- axTicks(1)[c(1,3)]
#         labx <- axTicks(1)[2]
#       }
#       if(length(axTicks(1))<=5) {
#         scalex <- axTicks(1)[c(1,2)]
#         labx <- mean(axTicks(1)[c(1,2)])
#       }
#       scaley <- axTicks(2)[1]
#     }
#     if(corthing>0) {
#       if(length(axTicks(1))>5) {
#         scalex <- axTicks(1)[c(length(axTicks(1))-2,length(axTicks(1)))]
#         labx <- axTicks(1)[length(axTicks(1))-1]
#       }
#       if(length(axTicks(1))<=5) {
#         scalex <- axTicks(1)[c(length(axTicks(1))-1,length(axTicks(1)))]
#         labx <- mean(axTicks(1)[c(length(axTicks(1))-1,length(axTicks(1)))])
#       }
#       scaley <- axTicks(2)[1]
#     }
#     lines(scalex,rep(scaley,2))
#     if(rivers$units=="m") text(labx,scaley,labels=paste((scalex[2]-scalex[1])/1000,"km"),pos=3,cex=cex)
#     if(rivers$units!="m") text(labx,scaley,labels=paste((scalex[2]-scalex[1]),rivers$units),pos=3,cex=cex)
#   }
# }


plotrivernetwork2 <- function(x,segmentnum=TRUE,offset=TRUE,lwd=1,cex=.6,scale=TRUE,color=TRUE,empty=FALSE,xlab="",ylab="",...) {
  if(class(x)!="rivernetwork") stop("Argument 'x' must be of class 'rivernetwork'.  See help(line2network) for more information.")
  lines <- x$lines
  length <- length(lines)
  allx <- unlist(lapply(lines,FUN='[',TRUE,1))
  ally <- unlist(lapply(lines,FUN='[',TRUE,2))
  xmax <- max(allx, na.rm=T)
  xmin <- min(allx, na.rm=T)
  ymax <- max(ally, na.rm=T)
  ymin <- min(ally, na.rm=T)
  plot(c(xmin,xmax),c(ymin,ymax),col="white",cex.axis=.6,asp=1,xlab=xlab,ylab=ylab,...=...)
  parusr <- par("usr")
  midptx <- midpty <- NA
  if(color&!empty) {
    linecolors <- rgb((sin(1:length)+1)/2.3,(cos(7*1:length)+1)/2.3,(sin(3*(length-(1:length)))+1)/2.3)
    ltys <- rep(1,length)
  }
  if(!color&!empty) {
    linecolors <- rep(1,length)
    ltys <- 1:length
  }
  if(empty) {
    ltys <- linecolors <- rep(1,length)
  }
  
  xtext <- ytext <- direction <- rep(NA,length)
  for(j in 1:length) {
    lines(lines[[j]],col=linecolors[j],lty=ltys[j],lwd=lwd)
    
    xplot <- lines[[j]][,1][lines[[j]][,1]>parusr[1] & lines[[j]][,1]<parusr[2] & lines[[j]][,2]>parusr[3] & lines[[j]][,2]<parusr[4]]
    yplot <- lines[[j]][,2][lines[[j]][,1]>parusr[1] & lines[[j]][,1]<parusr[2] & lines[[j]][,2]>parusr[3] & lines[[j]][,2]<parusr[4]]
    if(length(xplot)>0) midptx[j] <- mean(xplot)
    if(length(yplot)>0) midpty[j] <- mean(yplot)
    if(length(xplot)==0 | length(yplot)==0) midptx[j] <- midpty[j] <-NA
    middle <- floor(length(xplot)/2)
    direction[j] <- 3+(abs(xplot[floor(length(xplot)*.25)]-xplot[floor(length(xplot)*.75)]) < abs(yplot[floor(length(yplot)*.25)]-yplot[floor(length(yplot)*.75)]))
    xtext[j] <- ifelse(length(xplot)>0,xplot[middle],NA)
    ytext[j] <- ifelse(length(yplot)>0,yplot[middle],NA)
  }
  if(offset) text(xtext,ytext,labels=1:length,pos=direction,cex=cex,col=linecolors)
  else text(xtext,ytext,labels=1:length,cex=cex,col=linecolors)
  if(scale) {
    if(length>1) corthing <- cor(midptx,midpty,use="complete.obs")
    if(length==1) corthing <- (lines[[1]][1,1]-lines[[1]][dim(lines[[1]])[1],1])*(lines[[1]][1,2]-lines[[1]][dim(lines[[1]])[1],2])
    axticks1 <- axTicks(1)
    if(corthing<=0) {
      if(length(axticks1)>5) {
        scalex <- axticks1[c(1,3)]
        labx <- axticks1[2]
      }
      if(length(axticks1)<=5) {
        scalex <- axticks1[c(1,2)]
        labx <- mean(axticks1[c(1,2)])
      }
      scaley <- axTicks(2)[1]
    }
    if(corthing>0) {
      if(length(axticks1)>5) {
        scalex <- axticks1[c(length(axticks1)-2,length(axticks1))]
        labx <- axticks1[length(axticks1)-1]
      }
      if(length(axticks1)<=5) {
        scalex <- axticks1[c(length(axticks1)-1,length(axticks1))]
        labx <- mean(axticks1[c(length(axticks1)-1,length(axticks1))])
      }
      scaley <- axTicks(2)[1]
    }
    lines(scalex,rep(scaley,2))
    if(x$units=="m") text(labx,scaley,labels=paste((scalex[2]-scalex[1])/1000,"km"),pos=3,cex=cex)
    if(x$units!="m") text(labx,scaley,labels=paste((scalex[2]-scalex[1]),x$units),pos=3,cex=cex)
  }
}

# plotrivernetwork2(abstreams)
# 
# par(mfrow=c(5,5))
# system.time(plot(abstreams))
# system.time(plotrivernetwork2(abstreams))
# microbenchmark(plot(abstreams))
# microbenchmark(plotrivernetwork2(abstreams))


#' Highlight Segments
#' @description Plots a river network object and displays specified segments in bold, for easy identification.
#' @param seg A vector of segments to highlight
#' @param rivers The river network object to use
#' @param cex The character expansion factor to use for segment labels
#' @param lwd The line width to use for highlighted segments
#' @param add Whether to add the highlighted segments to an existing plot (\code{TRUE}) or call a new plot (\code{FALSE}).  Defaults to \code{FALSE}.
#' @param color Whether to display segment labels as the same color as the segments.  Defaults to \code{FALSE}.
#' @param ... Additional plotting arguments (see \link[graphics]{par})
#' @author Matt Tyers
#' @examples
#' data(Kenai3)
#' plot(Kenai3)
#' highlightseg(seg=c(10,30,68),rivers=Kenai3)
#' @importFrom grDevices rgb
#' @importFrom graphics plot
#' @importFrom graphics par
#' @importFrom graphics text
#' @export
highlightseg <- function(seg,rivers,cex=0.8,lwd=3,add=FALSE,color=FALSE,...) {
  length<-length(rivers$lines)
  if(!add) plot(rivers,color=FALSE,segmentnum=FALSE,...=...)
  lines<-rivers$lines
  midptx <- midpty <- NA
  for(j in seg) {
    lines(lines[[j]],col=rgb((sin(j)+1)/2.3,(cos(7*j)+1)/2.3,(sin(3*(length-j))+1)/2.3),lty=1,lwd=lwd)
    linelength <- dim(lines[[j]])[1]
    xplot <- lines[[j]][,1][lines[[j]][,1]>par("usr")[1] & lines[[j]][,1]<par("usr")[2] & lines[[j]][,2]>par("usr")[3] & lines[[j]][,2]<par("usr")[4]]
    yplot <- lines[[j]][,2][lines[[j]][,1]>par("usr")[1] & lines[[j]][,1]<par("usr")[2] & lines[[j]][,2]>par("usr")[3] & lines[[j]][,2]<par("usr")[4]]
    if(length(xplot)>0) midptx[j] <- mean(xplot)
    if(length(yplot)>0) midpty[j] <- mean(yplot)
    if(length(xplot)==0 | length(yplot)==0) midptx[j] <- midpty[j] <-NA
    text(midptx[j],midpty[j],j,cex=cex,col=ifelse(color,rgb((sin(j)+1)/2.3,(cos(7*j)+1)/2.3,(sin(3*(length-j))+1)/2.3),1))
  }
}


#' Check Connectivity of a River Network Object
#' @description Produces a graphical check of the connectivity of a river 
#'   network object.  It produces a \link{plot} of the river network object,
#'   and overlays red dots at non-connected endpoints and green dots at 
#'   connected endpoints.
#' @param rivers The river network object to check
#' @param add Whether call a new plot (\code{FALSE}) or add dots to an existing
#'   plot (\code{TRUE}).  Defaults to \code{FALSE}.
#' @param ... Additional plotting arguments (see \link[graphics]{par})
#' @author Matt Tyers
#' @seealso \link{line2network}
#' @examples 
#' data(Gulk)
#' topologydots(rivers=Gulk)
#' @importFrom graphics plot
#' @importFrom graphics points
#' @export
topologydots <- function(rivers,add=FALSE,...) {
  # if(class(rivers)!="rivernetwork") stop("Argument 'rivers' must be of class 'rivernetwork'.  See help(line2network) for more information.")
  # if(!add) plot(rivers,color=F,...=...)
  # connections <- rivers$connections
  # lines <- rivers$lines
  # connections[is.na(connections)==T]<-0
  # for(i in 1:dim(connections)[2]) {
  #   low <- 2
  #   high <- 2
  #   if(length(connections[i,][connections[i,]==1 | connections[i,]==2 | connections[i,]==5 | connections[i,]==6])>0) low <- 3
  #   if(length(connections[i,][connections[i,]==3 | connections[i,]==4 | connections[i,]==5 | connections[i,]==6])>0) high <- 3
  #   points(lines[[i]][1,1],lines[[i]][1,2],pch=16,col=low)
  #   points(lines[[i]][dim(lines[[i]])[1],1],
  #          lines[[i]][dim(lines[[i]])[1],2],pch=16,col=high)
  # }
  if(class(rivers)!="rivernetwork") stop("Argument 'rivers' must be of class 'rivernetwork'.  See help(line2network) for more information.")
  if(!add) plot(rivers,color=F,...=...)
  connections <- rivers$connections
  lines <- rivers$lines
  connections[is.na(connections)==T]<-0
  low <- high <- rep(2,dim(connections)[2])
  coordslow <- coordshigh <- matrix(NA,nrow=length(low),ncol=2)
  con1 <- rowSums(connections==1)
  con2 <- rowSums(connections==2)
  con3 <- rowSums(connections==3)
  con4 <- rowSums(connections==4)
  con5 <- rowSums(connections==5)
  con6 <- rowSums(connections==6)
  low[con1+con2+con5+con6>0] <- 3
  high[con3+con4+con5+con6>0] <- 3
  for(i in 1:dim(connections)[2]) {
    coordslow[i,] <- lines[[i]][1,]
    coordshigh[i,] <- lines[[i]][dim(lines[[i]])[1],]
  }
  points(coordslow, pch=16, col=low)
  points(coordshigh, pch=16, col=high)
}


#' Check Which Segments are Connected to a Given Segment.
#' @description Returns which segments are connected to a specified segment within a river network.  It may be useful for error checking.
#' @param seg The segment to check
#' @param rivers The river network object it belongs to
#' @return A vector of segment numbers
#' @author Matt Tyers
#' @examples 
#' data(Gulk)
#' plot(Gulk)
#' whoconnected(seg=4, rivers=Gulk)
#' @export
whoconnected <- function(seg,rivers) {
  connections1 <- !is.na(rivers$connections)
  connected <- which(connections1[seg,])
  return(connected)
}


#' Convert XY Coordinates to River Locations
#' @description This function determines the closest vertex in the river network
#'   to each point of XY data and returns a list of river locations, defined 
#'   as segment numbers and vertex numbers.
#' @param x A vector of x-coordinates to transform
#' @param y A vector of y-coordinates to transform
#' @param rivers The river network object to use
#' @return A data frame of river locations, with segment numbers in \code{$seg}, vertex numbers in \code{$vert}, and the snapping distance for each point in \code{$snapdist}.
#' @author Matt Tyers
#' @note Conversion to river locations is only valid if the input XY 
#'   coordinates and river network are in the same projected coordinate system. 
#'   Point data in geographic coordinates can be projected using 
#'   \link[rgdal]{project} in package 'rgdal', and an example is shown below.
#' @examples 
#' data(Gulk,fakefish)
#' head(fakefish)
#' 
#' fakefish.riv <- xy2segvert(x=fakefish$x, y=fakefish$y, rivers=Gulk)
#' head(fakefish.riv)
#' 
#' plot(x=Gulk, xlim=c(862000,882000), ylim=c(6978000,6993000))
#' points(fakefish$x, fakefish$y, pch=16, col=2)
#' riverpoints(seg=fakefish.riv$seg, vert=fakefish.riv$vert, rivers=Gulk, pch=15, col=4)
#' 
#' 
#' ## converting a matrix of points stored in long-lat to Alaska Albers Equal Area:
#' data(line98, Kenai1)
#' head(line98)  # note that coordinates are stored in long-lat, NOT lat-long
#' 
#' library(rgdal)
#' line98albers <- project(line98,proj="+proj=aea +lat_1=55 +lat_2=65 
#'     +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 
#'     +ellps=GRS80 +towgs84=0,0,0")
#' head(line98albers)
#' 
#' zoomtoseg(seg=c(162,19), rivers=Kenai1)
#' points(line98albers)
#' @export
xy2segvert <- function(x,y,rivers) {
  if(class(rivers)!="rivernetwork") stop("Argument 'rivers' must be of class 'rivernetwork'.  See help(line2network) for more information.")
  if(any(is.na(x))|any(is.na(y))|!is.numeric(x)|!is.numeric(y)) stop("Missing or non-numeric coordinates.")
  
  lengthlength <- length(unlist(lines))/2
  
  whichseg <- whichvert <- allx <- ally <- rep(NA,lengthlength)
  istart <- 1
  for(i in 1:length(rivers$lines)) {
    linelength <- dim(rivers$lines[[i]])[1]
    whichseg[istart:(istart+linelength-1)] <- i
    whichvert[istart:(istart+linelength-1)] <- 1:linelength
    allx[istart:(istart+linelength-1)] <- rivers$lines[[i]][,1]
    ally[istart:(istart+linelength-1)] <- rivers$lines[[i]][,2]
    istart <- istart+linelength
  }
  
  pdist2 <- function(p1,p2x,p2y) {
    dist <- sqrt((p1[1]-p2x)^2 + (p1[2]-p2y)^2)
    return(dist)
  }
  
  min.i <- mapply(function(x,y) which.min(pdist2(c(x,y),allx,ally))[1],x,y)  # if there's a tie, accept the first vertex
  seg <- whichseg[min.i]
  vert <- whichvert[min.i]
  snapdist <- sqrt((x-allx[min.i])^2 + (y-ally[min.i])^2)
  out <- data.frame(seg,vert,snapdist)
  return(out)
}

# 
# xy2segvert_OLD <- function(x,y,rivers) {
#   if(class(rivers)!="rivernetwork") stop("Argument 'rivers' must be of class 'rivernetwork'.  See help(line2network) for more information.")
#   if(any(is.na(x))|any(is.na(y))|!is.numeric(x)|!is.numeric(y)) stop("Missing or non-numeric coordinates.")
#   seg <- rep(NA,length(x))
#   vert <- rep(NA,length(x))
#   snapdist <- rep(NA,length(x))
#   lines <- rivers$lines
#   
#   pdist2 <- function(p1,p2mat) {
#     dist <- sqrt((p1[1]-p2mat[,1])^2 + (p1[2]-p2mat[,2])^2)
#     return(dist)
#   }
#   
#   for(i in 1:length(x)) {
#     min.dist1 <- 1000000
#     for(river.i in 1:length(lines)) {
#       dists <- pdist2(c(x[i],y[i]),lines[[river.i]])
#       min.i <- min(dists)
#       if(min.i < min.dist1) {
#         min.dist1 <- min.i
#         seg[i] <- river.i
#         vert[i] <- which(dists==min.i)[1]  # if there's a tie, accept the first vertex
#         snapdist[i] <- pdist(c(x[i],y[i]),lines[[river.i]][vert[i],])
#       }
#     }
#   }
#   out <- data.frame(seg,vert,snapdist)
#   return(out)
# }
# 
# rivtest <- Kenai3
# abpts <- xysim(5000,rivtest)
# asdf <- xy2segvert2(x=abpts$x,y=abpts$y,rivers=rivtest)
# 
# system.time(asdf1 <- xy2segvert(x=abpts$x,y=abpts$y,rivers=rivtest))
# system.time(asdf2 <- xy2segvert2(x=abpts$x,y=abpts$y,rivers=rivtest))
# 
# microbenchmark(asdf1 <- xy2segvert(x=abpts$x,y=abpts$y,rivers=rivtest))
# microbenchmark(asdf2 <- xy2segvert2(x=abpts$x,y=abpts$y,rivers=rivtest))
# all.equal(asdf1,asdf2)
# 



#' Draw Points from River Locations
#' @description Adds points to an active plot.  Works like \link[graphics]{points} but with river locations (segments and vertices) rather than xy coordinates.
#' @param seg A vector of segments
#' @param vert A vector of vertices
#' @param rivers The river network object to use
#' @param pch Point character, as a vector or single value
#' @param col Point color, as a vector or single value
#' @param jitter Maximum amount of random noise to add to "jitter" points if desired, so points do not overlap one another
#' @param ... Additional arguments for \link{points}
#' @author Matt Tyers
#' @examples
#' data(fakefish,Gulk)
#' 
#' plot(x=Gulk, xlim=c(862000,882000), ylim=c(6978000,6993000))
#' points(x=fakefish$x, y=fakefish$y, pch=16, col=2)
#' riverpoints(seg=fakefish$seg, vert=fakefish$vert, rivers=Gulk, pch=15, col=4)
#' 
#' plot(x=Gulk, empty=TRUE)
#' with(fakefish, riverpoints(seg=seg, vert=vert, rivers=Gulk, 
#'        pch=16, col=flight, jitter=1000))
#' @importFrom graphics points
#' @importFrom stats runif
#' @export
riverpoints <- function(seg,vert,rivers,pch=1,col=1,jitter=0,...) {
  if(class(rivers)!="rivernetwork") stop("Argument 'rivers' must be of class 'rivernetwork'.  See help(line2network) for more information.")
  lines <- rivers$lines
  if(max(seg,na.rm=T)>length(lines) | min(seg,na.rm=T)<1) stop("Invalid segment numbers specified.")
  
  x <- y <- rep(NA,length(seg))
  for(i in 1:length(seg)) {
    x[i] <- lines[[seg[i]]][vert[i],1]
    y[i] <- lines[[seg[i]]][vert[i],2]
  }
  if(jitter!=0) {
    x <- x+runif(length(x),min=-jitter,max=jitter)
    y <- y+runif(length(x),min=-jitter,max=jitter)
  }
  points(x,y,pch=pch,col=col,...=...)
}

# riverpoints_OLD <- function(seg,vert,rivers,pch=1,col=1,jitter=0,...) {
#   if(class(rivers)!="rivernetwork") stop("Argument 'rivers' must be of class 'rivernetwork'.  See help(line2network) for more information.")
#   if(length(pch)==1) pch <- rep(pch,length(seg))
#   if(length(col)==1) col <- rep(col,length(seg))
#   lines <- rivers$lines
#   if(max(seg,na.rm=T)>length(lines) | min(seg,na.rm=T)<1) stop("Invalid segment numbers specified.")
#   
#   for(i in 1:length(seg)) {
#     if(jitter!=0) {
#       jitterx <- runif(length(lines[[seg[i]]][vert[i],1]),min=-jitter,max=jitter)
#       jittery <- runif(length(lines[[seg[i]]][vert[i],2]),min=-jitter,max=jitter)
#     }
#     else {
#       jitterx <- 0
#       jittery <- 0
#     }
#     if(vert[i]>dim(lines[[seg[i]]])[1] | vert[i]<1) stop("Invalid vertex numbers specified.")
#     points(lines[[seg[i]]][vert[i],1] + jitterx,
#            lines[[seg[i]]][vert[i],2] + jittery, pch=pch[i],col=col[i],...=...)
#   }
# }


# 
# plot(rivtest)
# system.time(riverpoints(asdf2$seg,asdf2$vert,rivtest))
# system.time(riverpoints2(asdf2$seg,asdf2$vert,rivtest,pch=16,col=4))
# 
# microbenchmark(riverpoints(asdf2$seg,asdf2$vert,rivtest))
# microbenchmark(riverpoints2(asdf2$seg,asdf2$vert,rivtest))



#' Trim a River Network Object to Specified Segments
#' @description Removes line segments from a river network object.  User can specify which segments to remove (\code{trim}) or which segments to keep (\code{trimto}).
#' @param trim Vector of line segments to remove
#' @param trimto Vector of line segments to keep
#' @param rivers The river network object
#' @return A new river network object
#' @seealso \link{line2network}
#' @note Specifying segments in both trim and trimto arguments will result in an error.
#' @author Matt Tyers
#' @examples
#' data(Kenai1)
#' plot(x=Kenai1)
#' 
#' Kenai1.trim <- trimriver(trim=c(46,32,115,174,169,114,124,142,80), rivers=Kenai1)
#' plot(x=Kenai1.trim)
#'  
#' Kenai1.trim.2 <- trimriver(trimto=c(20,57,118,183,45,162,39,98,19), rivers=Kenai1)
#' plot(x=Kenai1.trim.2)
#' @export
trimriver <- function(trim=NULL,trimto=NULL,rivers) {
  if(class(rivers)!="rivernetwork") stop("Argument 'rivers' must be of class 'rivernetwork'.  See help(line2network) for more information.")
  segs <- 1:length(rivers$lines)
  if(!is.null(trim) & !is.null(trimto)) {
    stop("Error - cannot use both trim and trimto arguments")
  }
  if(!is.null(trim)) {
    if(max(trim,na.rm=T)>length(rivers$lines) | min(trim,na.rm=T)<1) stop("Invalid segment numbers specified.")
    segs <- segs[-trim]
  }
  if(!is.null(trimto)) {
    if(max(trimto,na.rm=T)>length(rivers$lines) | min(trimto,na.rm=T)<1) stop("Invalid segment numbers specified.")
    segs <- segs[trimto]
  }
  trimmed.rivers <- rivers
  trimmed.rivers$lines <- trimmed.rivers$lines[segs]
  trimmed.rivers$connections <- as.matrix(trimmed.rivers$connections[segs,segs])
  if(!is.na(trimmed.rivers$braided)) {
    if(trimmed.rivers$braided & !any(trimmed.rivers$connections %in% 5:6)) trimmed.rivers$braided <- NA
  }
  trimmed.rivers$names <- rivers$names[segs]
  trimmed.rivers$lengths <- rivers$lengths[segs]
  if(length(segs)==0) stop("Error - resulting river network has no remaining line segments")
  if(!is.na(trimmed.rivers$mouth$mouth.seg) & !is.na(trimmed.rivers$mouth$mouth.vert)) {
    if(any(segs==rivers$mouth$mouth.seg)) {
      trimmed.rivers$mouth$mouth.seg <- which(segs==trimmed.rivers$mouth$mouth.seg)
    }
    if(!any(segs==rivers$mouth$mouth.seg)) {
      trimmed.rivers$mouth$mouth.seg <- NA
      trimmed.rivers$mouth$mouth.vert <- NA
      message("River mouth must be redefined - see help(setmouth)")
    }
  }
  
  if(!is.null(rivers$segroutes) | !is.null(rivers$distlookup)) {
    trimmed.rivers$segroutes <- NULL
    trimmed.rivers$distlookup <- NULL
    warning("Segment routes and/or distance lookup must be rebuilt - see help(buildsegroutes).")
  }
  trimmed.rivers <- addcumuldist(trimmed.rivers)
  
  # updating sp object
  id <- rivers$lineID
  sp_lines1 <- rivers$sp@lines[unique(id[segs,2])]   
  j<-1
  for(i in unique(id[segs,2])) {
    sp_lines1[[j]]@Lines <- sp_lines1[[j]]@Lines[id[segs,3][id[segs,2]==i]]
    j<-j+1
  }
  rivID <- NA
  sp_line <- NA
  sp_seg <- NA
  k<-1
  for(i in 1:length(sp_lines1)) {
    for(j in 1:length(sp_lines1[[i]]@Lines)) {
      sp_line[k] <- i
      sp_seg[k] <- j
      k<-k+1
    }
  }
  rivID <- 1:(k-1)
  lineID <- data.frame(rivID,sp_line,sp_seg)
  trimmed.rivers$lineID <- lineID
  trimmed.rivers$sp@lines <- sp_lines1
  if(dim(rivers$sp@data)[1]==max(rivers$lineID[,2]) & dim(rivers$sp@data)[1]>1) {
    trimmed.rivers$sp@data <- rivers$sp@data[unique(rivers$lineID[segs,2]),]
  }
  message("Note: any point data already using the input river network must be re-transformed to river coordinates using xy2segvert() or ptshp2segvert().")
  return(trimmed.rivers)
}


#' Trim a River Network to a Set of X-Y Coordinates
#' @description Removes line segments from a river network object that are not
#'   adjacent to a set of point data, given in X-Y coordinates.
#' @param x Vector of x-coordinates of point data to buffer around
#' @param y Vector of y-coordinates of point data to buffer around
#' @param rivers The river network object to use
#' @param method Three methods are available.  If \code{"snap"} is specified
#'   (the default), only the closest segment to each point is retained.   If
#'   \code{"snaproute"} is specified, segments are also retained that will
#'   maintain total connectivity in the resulting river network.  If
#'   \code{"buffer"} is specified, all segments with endpoints or midpoints
#'   within \code{dist} units of the input locations are retained.
#' @param dist Distance to use for buffering, if \code{method="buffer"}.  If
#'   this is not specified, the maximum spread in the x- and y- direction will
#'   be used.
#' @return A new river network object (see \link{rivernetwork})
#' @seealso \link{line2network}
#' @note If \code{method=="buffer"}, only distances to segment endpoints 
#'   and midpoints are checked, and still only whole segments are removed.
#' @author Matt Tyers
#' @examples
#' data(Koyukuk2)
#' x <- c(139241.0, 139416.1, 124600.1, 122226.8)
#' y <- c(1917577, 1913864, 1898723, 1898792)
#' 
#' plot(x=Koyukuk2)
#' points(x, y, pch=15, col=4)
#' legend(par("usr")[1], par("usr")[4], legend="points to buffer around", pch=15, col=4, cex=.6)
#' 
#' Koyukuk2.buf1 <- trimtopoints(x, y, rivers=Koyukuk2, method="snap")
#' plot(x=Koyukuk2.buf1)
#' points(x, y, pch=15, col=4)
#' 
#' Koyukuk2.buf2 <- trimtopoints(x, y, rivers=Koyukuk2, method="snaproute")
#' plot(x=Koyukuk2.buf2)
#' points(x, y, pch=15, col=4)
#' 
#' Koyukuk2.buf3 <- trimtopoints(x, y, rivers=Koyukuk2, method="buffer", dist=1000)
#' plot(x=Koyukuk2.buf3)
#' points(x, y, pch=15, col=4)
#' @export
trimtopoints <- function(x,y,rivers,method="snap",dist=NULL) {
  if(class(rivers)!="rivernetwork") stop("Argument 'rivers' must be of class 'rivernetwork'.  See help(line2network) for more information.")
  if(is.null(dist)) dist <- max(c(max(x)-min(x)),c(max(y)-min(y)))
  
  if(method=="buffer") {
    keep <- 0
    k <- 1
    for(i in 1:length(x)) {
      for(j in 1:length(rivers$lines)) {
        line.len <- dim(rivers$lines[[j]])[1]
        if((pdist(c(x[i],y[i]),rivers$lines[[j]][1,]) <= dist) | 
           (pdist(c(x[i],y[i]),rivers$lines[[j]][line.len,]) <= dist) |
           (pdist(c(x[i],y[i]),rivers$lines[[j]][floor(line.len/2),]) <= dist)) {
          if(length(keep[keep==j])==0) {
            keep[k] <- j
            k <- k+1
          }
        }
      }
    }
  }
  
  if(method=="bufferroute") {
    keep <- 0
    k <- 1
    for(i in 1:length(x)) {
      for(j in 1:length(rivers$lines)) {
        line.len <- dim(rivers$lines[[j]])[1]
        if((pdist(c(x[i],y[i]),rivers$lines[[j]][1,]) <= dist) | 
           (pdist(c(x[i],y[i]),rivers$lines[[j]][line.len,]) <= dist) |
           (pdist(c(x[i],y[i]),rivers$lines[[j]][floor(line.len/2),]) <= dist)) {
          if(length(keep[keep==j])==0) {
            keep[k] <- j
            k <- k+1
          }
        }
      }
    }
    
    keep.bool <- rep(FALSE,length(rivers$lines))
    keep.bool[keep] <- TRUE
    
    for(i in keep) {
      for(j in keep) {
        if(i!=j) {
          route1 <- detectroute(start=i,end=j,rivers=rivers)
          keep.bool[route1] <- TRUE
        }
      }
    }
    keep <- (1:length(rivers$lines))[keep.bool]
  }
  
  if(method=="snap") {
    segvert <- xy2segvert(x,y,rivers=rivers)
    keep.bool <- rep(FALSE,length(rivers$lines))
    keep.bool[unique(segvert$seg)] <- TRUE
    keep <- (1:length(rivers$lines))[keep.bool]
  }
  
  if(method=="snaproute") {
    segvert <- xy2segvert(x=x,y=y,rivers)
    keep.bool <- rep(FALSE,length(rivers$lines))
    keep.bool[unique(segvert$seg)] <- TRUE
    
    for(i in segvert$seg) {
      for(j in segvert$seg) {
        if(i!=j) {
          route1 <- detectroute(start=i,end=j,rivers=rivers)
          keep.bool[route1] <- TRUE
        }
      }
    }
    keep <- (1:length(rivers$lines))[keep.bool]
  }
  
  rivers1 <- rivers
  rivers1$lines <- rivers1$lines[keep]
  rivers1$connections <- as.matrix(rivers1$connections[keep,keep])
  if(!is.na(rivers1$braided)) {
    if(rivers1$braided & !any(rivers1$connections %in% 5:6)) rivers1$braided <- NA
  }
  rivers1$names <- rivers$names[keep]
  rivers1$lengths <- rivers$lengths[keep]
  if(keep[1]==0) stop("Error - resulting river network has no remaining line segments")
  
  if(!is.na(rivers1$mouth$mouth.seg) & !is.na(rivers1$mouth$mouth.vert)) {
    if(any(keep==rivers$mouth$mouth.seg)) {
      rivers1$mouth$mouth.seg <- which(keep==rivers1$mouth$mouth.seg)
    }
    if(!any(keep==rivers$mouth$mouth.seg)) {
      rivers1$mouth$mouth.seg <- NA
      rivers1$mouth$mouth.vert <- NA
    }
  }
  
  if(!is.null(rivers1$segroutes) | !is.null(rivers1$distlookup)) {
    rivers1$segroutes <- NULL
    rivers1$distlookup <- NULL
    warning("Segment routes and/or distance lookup must be rebuilt - see help(buildsegroutes).")
  }
  rivers1 <- addcumuldist(rivers1)
  
  # updating sp object
  id <- rivers$lineID
  sp_lines1 <- rivers$sp@lines[unique(id[keep,2])]  
  j<-1
  for(i in unique(id[keep,2])) {
    sp_lines1[[j]]@Lines <- sp_lines1[[j]]@Lines[id[keep,3][id[keep,2]==i]]
    j<-j+1
  }
  rivID <- NA
  sp_line <- NA
  sp_seg <- NA
  k<-1
  for(i in 1:length(sp_lines1)) {
    for(j in 1:length(sp_lines1[[i]]@Lines)) {
      sp_line[k] <- i
      sp_seg[k] <- j
      k<-k+1
    }
  }
  rivID <- 1:(k-1)
  lineID <- data.frame(rivID,sp_line,sp_seg)
  rivers1$lineID <- lineID
  rivers1$sp@lines <- sp_lines1
  if(dim(rivers$sp@data)[1]==max(rivers$lineID[,2]) & dim(rivers$sp@data)[1]>1) {
    rivers1$sp@data <- rivers$sp@data[unique(rivers$lineID[keep,2]),]
  }
  message("Note: any point data already using the input river network must be re-transformed to river coordinates using xy2segvert() or ptshp2segvert().")
  return(rivers1)
}


#' Remove Unconnected Segments
#' @description Detects and removes segments that are not connected to the river mouth.
#' @param rivers The river network object to use.
#' @author Matt Tyers
#' @note This function is called within \link{cleanup}, which is recommended in most cases.
#' @examples 
#' data(Koyukuk2)
#' Koy_subset <- trimriver(trimto=c(30,28,29,3,19,27,4),rivers=Koyukuk2)
#' Koy_subset <- setmouth(seg=1,vert=427,rivers=Koy_subset)
#' plot(Koy_subset)
#' 
#' Koy_subset_trim <- removeunconnected(Koy_subset)
#' plot(Koy_subset_trim)
#' @export
removeunconnected <- function(rivers) {
  if(is.na(rivers$mouth$mouth.seg)) stop("River mouth must be specified.")
  origin <- rivers$mouth$mouth.seg
  
  takeout <- NULL
  k <- 1
  
  dists <- rep(NA,length(rivers$lines))
  for(i in 1:length(rivers$lines)) {
    dists[i] <- pdist(rivers$lines[[rivers$mouth$mouth.seg]][rivers$mouth$mouth.vert,],
                      c(mean(rivers$lines[[i]][,1]),mean(rivers$lines[[i]][,2])))
  }
  order <- order(dists)
  rivers1 <- rivers
  rivers1$connections <- rivers$connections[order,order]
  rivers1$lines <- rivers$lines[order]
  origin1 <-which(order==origin)
  
  
  for(i in 1:length(rivers1$lines)) {
    if(is.na(detectroute(end=origin1,start=i,rivers=rivers1,stopiferror=FALSE,algorithm="sequential")[1])) {
      takeout[k] <- i
      k <- k+1
    }
  }
  
  if(!is.null(takeout)) takeout <- order[takeout]
  
  suppressMessages(rivers2 <- trimriver(trim=takeout,rivers=rivers))
  message("Note: any point data already using the input river network must be re-transformed to river coordinates using xy2segvert() or ptshp2segvert().")
  return(rivers2)
}
mbtyers/riverdist documentation built on May 22, 2019, 12:58 p.m.