R/riverdist_1.R

Defines functions removeunconnected trimtopoints trimriver riverpoints xy2segvert whoconnected topologydots highlightseg plotrivernetwork2 scalebar plot.rivernetwork pointshp2segvert line2network calculateconnections addcumuldist pdisttot pdist

Documented in addcumuldist calculateconnections highlightseg line2network pdist pdisttot plot.rivernetwork pointshp2segvert removeunconnected riverpoints topologydots trimriver trimtopoints whoconnected xy2segvert

#' 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"
#'     
#' 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 (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")
#' 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)
}

Try the riverdist package in your browser

Any scripts or data that you put into this service are public.

riverdist documentation built on Sept. 10, 2021, 9:06 a.m.