#' determine index of downstream
#' \code{sp.RiverOrder}
#' @param sp SpatialLines*
#' @param coord coordinates of \code{sp}.
#' @return Index of downstream for each segements
#' @export
sp.RiverDown <- function(sp, coord = extractCoords(sp)){
ft = FromToNode(sp, coord = coord)
nsp = length(sp)
idown = rep(-3, nsp)
for(i in 1:nsp){
pto = ft[i,2]
id = which(ft[,1] == pto)
if(length(id) == 1){
idown[i] = id
}
}
idown
}
#' determine River path, disolve the river network based on downstream-relationship.
#' \code{sp.RiverPath}
#' @param sp SpatialLines*
#' @param idown downstream index of \code{sp}.
#' @param coord coordinates of \code{sp}.
#' @return List of River Path(SegIDs, PointIDs, sp)
#' @export
#' @examples
# data(sh)
# riv=sh$riv
# x=sp.CutSptialLines(riv, tol=200)
# px(x)
# xx=sp.RiverPath(x)
# px(xx$sp)
sp.RiverPath <- function(sp,
coord = extractCoords(sp, unique = TRUE),
idown = sp.RiverDown(sp, coord = coord)){
pt.list <- unlist(coordinates(sp), recursive = FALSE)
id.list<-lapply(pt.list, function(x) { xy2ID(xy = x, coord = coord)})
ft = FromToNode(sp, coord = coord)
nsp = length(sp)
p.frq = data.frame(table(unlist(id.list)) )
p0 = ft[!(ft[, 1] %in% ft[,2]), 1]
tmp = p.frq[p.frq[,2]>2, 1] # more than twice.
p1 = as.numeric(levels(tmp))[tmp]
p2 = ft[!(ft[, 2] %in% ft[,1]), 2] # TO that not in FROM
# jid = which(table(as.matrix(ft)) > 2) #id of joint points
riv.keyid = c(which(idown<0), which(ft[,2] %in% p1) )
updown = cbind(1:nsp, idown)
goUp <- function(updown, id0){
uid = which(idown == id0) # id of my upstream
if(length(uid) == 1){ #up stream exist
ret = c(goUp(updown, uid), id0)
}else{
ret = id0
}
ret
}
StreamPath = lapply(riv.keyid, function(x) goUp(updown, x))
#========Point ID==========
nstr = length(StreamPath)
pl = list()
for(i in 1:nstr){
sid = StreamPath[[i]]
pid = unique(unlist(lapply(1:length(sid), function(x) id.list[[sid[x]]])))
# print(sid)
# print(pid)
pl[[i]] = pid
}
#=========Spatial Lines===============
ll = list()
for(i in 1:nstr){
sid = StreamPath[[i]]
pid = unique(unlist(lapply(1:length(sid), function(x) id.list[[sid[x]]])))
ll[[i]] = sp::Lines( sp::Line(coord[pid, ] ), ID = i)
}
spx = sp::SpatialLines(ll, proj4string = raster::crs(sp))
ret <- list(SegIDs = StreamPath,
PointIDs= pl,
sp = spx)
}
#' calculate river order
#' \code{sp.RiverOrder}
#' @param sp SpatialLines*
#' @param coord coordinates of \code{sp}.
#' @return Stream Order of SpatialLines*
#' @export
sp.RiverOrder <- function(sp, coord = extractCoords(sp)){
msg='sp.RiverOrder::'
get1st <- function(x){
fr = x[,2]
to = x[,3]
tb= table(x[,2:3])
pid = as.numeric(names(tb))
ncount = as.numeric(tb)
p.out = to[which(! to %in% fr)]
p.jnt = pid[which(ncount> 2)]
p.key = c(p.out, p.jnt)
# print(p.key)
tid=fr[ !fr %in% to] #target from-id
# plot(riv.simp)
# points(coord[ ])
# points(coord[p.key,], col=2)
nd = length(tid)
ret = NULL
# message(msg, 'Number of Streams = ', nd)
for(i in 1:nd){
cr = tid[i]
for(j in 1:100000){
# message(msg, 'Stream ',i,'\tSeg ', j, '\tFrom Node ', cr)
rid = which(fr %in% cr)
ret=c(ret, rid)
sid= to[rid] #ID of to-point;
# if(length(sid) < 1 | length(sid)>1){
# print(sid)
# }
if( any(sid %in% p.key) ){ #if sid IS IN key-points(outlets or joint point).
break;
}else{
cr = sid;
}
}
}
ret;
}
# ext=raster::extent(sp)
# sp.tmp = rgeos::gSimplify(sp, tol=max(diff(ext[1:2] ), diff(ext[3:4])) )
ft = unique(FromToNode(sp, coord))
x = cbind(1:length(sp), ft)
y =x
x.ord =x[,1]*0
for(i in 1:10000){
ids = y[,1]
message(msg, 'Order = ', i)
id = get1st(y)
x.ord[ ids[id] ] = i
y=rbind(y[-1 * id ,])
ny = length(y)
# message(msg, 'Left Segments =', ny)
if(ny<=0){ break }
}
x.ord
}
#' return the Index of nodes in SpatialData
#' \code{NodeIDList}
#' @param sp SpatialLines*
#' @param coord Coordinate of vertex in \code{sp}
#' @return list of Index of each SpatialData, line or polygon
#' @export
NodeIDList <- function(sp, coord = extractCoords(sp, unique = TRUE) ){
pt.list <- unlist(coordinates(sp), recursive = FALSE)
id.list<-lapply(pt.list, function(x) { xy2ID(x, coord = coord)})
}
#' return the FROM and TO nodes index of the SpatialLines
#' \code{FromToNode}
#' @param sp SpatialLines*
#' @param coord Coordinate of vertex in \code{sp}
#' @return FROM and TO nodes index of the SpatialLines
#' @export
FromToNode <- function(sp, coord = extractCoords(sp, unique = TRUE) ){
id.list = NodeIDList(sp, coord=coord)
frto = cbind(
unlist(lapply(id.list, function(x) x[1])),
unlist(lapply(id.list, function(x) x[length(x)] ) )
)
colnames(frto)=c('FrNode', 'ToNode')
frto
}
#' parameters for river types
#' \code{RiverType}
#' @param n number of types
#' @return data.frame of parameters
#' @export
RiverType <- function(n){
cn = c('Index', 'Depth', 'BankSlope',
'Width', 'Sinuosity', 'Manning',
'Cwr', 'KsatH', 'BedThick')
nc = length(cn)
rtype = cbind(1:n,
5.0 + 0.5 * 1:n, #Depth
0, #bankslope
2*1:n, #Width
1.0, #Sinuosity
4.63e-07, #manning's n, day/m^(1/3)
0.6, #CWR
0.1, #KsatH
.1 # Bed Sediment Thickness.
)
colnames(rtype) = cn
rtype
}
#' Cut the river by PIHM.mesh
#' \code{sp.RiverSeg}
#' @param sp.mesh PIHM mesh
#' @param sp.riv river
#' @return SpatialLinesDataFrame, the \code{PIHM.river} was cut by \code{PIHM.mesh}
#' @export
sp.RiverSeg <- function(sp.mesh, sp.riv){
sp = sp::SpatialPolygonsDataFrame(sp.mesh,
data = data.frame('iEle'=1:length(sp.mesh) ),
match.ID = FALSE)
sr =sp::SpatialLinesDataFrame(sp.riv,
data = data.frame('iRiv' = 1:length(sp.riv)),
match.ID = FALSE )
seg = raster::intersect(sr, sp)
seg
}
#'
#' #' Dissolve River segments
#' #' \code{sp.RiverDissolve}
#' #' @param sp SpatialLines
#' #' @return SpatialLinesDataFrame
#' #' @export
#' sp.DissolveLines <- function(sp){
#' msg='sp.DissolveLines::'
#' id2Line <- function(id, key){
#' lp=key
#' ct = count(as.numeric(id))
#' cn = as.numeric( names(ct))
#' # ctf = count(id[,1]) #count of frNode
#' # ctt = count(id[,2])
#' lfr = which(id[,1] == key)
#' pto = which( cn == id[lfr,2] )
#' if(length(pto)>0){
#' if( ct[pto] == 2 ) {
#' # 1 - header/outlet, >2 - joint point
#' # message(msg, key, '\t', pto)
#' ret = id2Line(id, key = pto)
#' return (c(lp, ret) )
#' }else{
#' return(c(key, cn[pto]) )
#' }
#' }else{
#' return(c(key, cn[pto]) )
#' }
#' }
#' cdu = extractCoords(sp, unique = T)
#' ft=FromToNode(sp, coord = cdu)
#' ct= count(ft)
#' pnode = as.numeric(names(ct))[which(ct != 2)]
#'
#' # plot(sp)
#' sl=list()
#' k=0
#' nnode = length(pnode)
#' for(i in 1:nnode ){
#' key = pnode[i]
#' # points(cdu[key, 1], cdu[key, 2], pch=20)
#' if( length(which(ft[,1] == key)) > 0){
#' k = k +1
#' # message(msg, k, '\t', i, '/', nnode )
#' ids = id2Line(ft, key)
#' # cdf[ids] = k
#' # points(cdu[ids, ], col=k, pch=k)
#' # text(cdu[ids[1], 1], cdu[ids[1], 2], paste(k), col=2 )
#' ln = sp::Line( cdu[ids,] )
#' sl[[k]] = sp::Lines(ln, ID=paste(k) )
#' }else{
#' # print(key)
#' }
#' }
#' sll = sp::SpatialLines(sl)
#' # ord= sp.RiverOrder(sll)
#' df = data.frame('INDEX'=1:length(sl))
#' sld = sp::SpatialLinesDataFrame(sll, data = df)
#' sld
#' }
#' Find Shared points in SpatialLine*
#' \code{SharedPoints}
#' @param sp SpatialLines
#' @return Topologic relationship between components of SpatialLines.
#' @export
SharedPoints <- function(sp){
msg='SharedPoints::'
# sp=riv
pl = extractCoords(sp, aslist = TRUE)
nsp = length(sp)
ret = matrix(0, nrow = nsp, ncol=nsp)
for(i in 1:(nsp-1) ){
for(j in (i+1):nsp ){
plot(rbind(pl[[i]], pl[[j]]), asp=1);
points(pl[[i]], col=2)
points(pl[[j]], col=3)
cp = CommonPoints(pl[[i]], pl[[j]])
if( is.null(cp) ){
}else{
message(msg, i, '/', j)
# print(ret)
# print(cp)
# readline("go")
ret[i,j] = cp[,1]
ret[j,i] = cp[,2]
}
}
}
ret
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.