R/aggregate_OCN.R
In lucarraro/OCNet: Optimal Channel Networks
Defines functions
aggregate_OCN
Documented in
aggregate_OCN
aggregate_OCN <- function(OCN,
thrA=0.002*OCN$FD$nNodes*OCN$cellsize^2,
streamOrderType="Strahler",
maxReachLength=Inf,
equalizeLengths=FALSE,
breakpoints=NULL,
displayUpdates=FALSE){
if (!("slope" %in% names(OCN$FD))){
stop('Missing fields in OCN. You should run landscape_OCN prior to aggregate_OCN.')
}
if (maxReachLength < OCN$cellsize*sqrt(2)){
stop("maxReachLength cannot be smaller than OCN$cellsize*sqrt(2).")
}
#t1 <- Sys.time()
if (thrA==0) maxReachLength <- OCN$cellsize*sqrt(2)
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# BUILD NETWORK AT RN LEVEL ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
if(displayUpdates){message("Calculating network at RN level... \r", appendLF = FALSE)}
#print('Crop data at FD level to RN level...',quote=FALSE);
RN_mask <- as.vector(OCN$FD$A >= thrA)# RN_mask allows to sample RN-level values from matrices/vectors at FD level
RN_to_FD <- which(RN_mask) # RN_to_FD[i] is the pixel ID at the FD level of the pixel whose ID at the RN level is i
FD_to_RN <- RN_mask*cumsum(as.numeric(RN_mask)) # FD_to_RN[i] is the pixel ID at the RN level of the pixel whose ID at the FD level is i
# if pixel i at FD level doesn't belong to RN, then FD_to_RN[i]=0
Nnodes_RN <- length(RN_to_FD)
W_RN <- OCN$FD$W[RN_mask,,drop=FALSE]
W_RN <- W_RN[,RN_mask,drop=FALSE]
Outlet_RN <- FD_to_RN[OCN$FD$outlet]
Outlet_RN <- Outlet_RN[Outlet_RN!=0] # remove outlets if the corresponding catchment size is lower than threshold
DownNode_RN <- numeric(Nnodes_RN)
# for (i in 1:Nnodes_RN){
# if (!(i %in% Outlet_RN)){
# DownNode_RN[i] <- which(W_RN[i,]==1)
# }}
tmp <- W_RN@rowpointers
NotOutlet <- which((tmp[-1] - tmp[-length(tmp)])==1)
DownNode_RN[NotOutlet] <- W_RN@colindices
# reverse downNode_RN
DownNode_RN_rev <- vector("list",Nnodes_RN)
for (i in 1:Nnodes_RN){
d <- DownNode_RN[i]
if (d!=0){DownNode_RN_rev[[d]] <- c(DownNode_RN_rev[[d]],i) }}
A_RN <- OCN$FD$A[RN_mask]
X_RN <- OCN$FD$X[RN_mask]
Y_RN <- OCN$FD$Y[RN_mask]
Z_RN <- OCN$FD$Z[RN_mask]
Length_RN <- OCN$FD$leng[RN_mask]
# Drainage density
DrainageDensity_RN <- sum(Length_RN)/(OCN$dimX*OCN$dimY*OCN$cellsize^2)
# Connectivity indices at pixel level
DegreeIn <- colSums(W_RN)
DegreeOut <- rowSums(W_RN)
Confluence <- DegreeIn>1
Source <- DegreeIn==0
SourceOrConfluence <- Source|Confluence
ConfluenceNotOutlet <- Confluence&(DownNode_RN!=0)
ChannelHeads <- SourceOrConfluence #Source|ConfluenceNotOutlet
OutletNotChannelHead <- (DownNode_RN==0)&(!ChannelHeads)
IsNodeAG <- SourceOrConfluence|OutletNotChannelHead
IsNodeAG[breakpoints] <- TRUE
whichNodeAG <- which(IsNodeAG)
# Calculate slope for each pixel of the river network
Slope_RN <- OCN$FD$slope[RN_mask]
# sort nodes in downstream direction
ind_sort <- sort(A_RN, index.return=TRUE)
ind_sort <- ind_sort$ix
Upstream_RN <- vector("list",Nnodes_RN)
Nupstream_RN <- numeric(Nnodes_RN)
for (i in 1:Nnodes_RN){
ups <- as.numeric(DownNode_RN_rev[[ind_sort[i]]])
nodes <- numeric(0)
for (u in ups){ nodes <- c(nodes, Upstream_RN[[u]])}
Upstream_RN[[ind_sort[i]]] <- c(nodes, ind_sort[i])
Nupstream_RN[ind_sort[i]] <- length(Upstream_RN[[ind_sort[i]]])
}
# RN_to_CM[i] indicates outlet to which reach i drains
RN_to_CM <- numeric(Nnodes_RN)
for (i in 1:OCN$nOutlet){
RN_to_CM[Upstream_RN[[Outlet_RN[i]]]] <- i
}
if (displayUpdates){message("Calculating network at RN level... 100.0%\n", appendLF = FALSE)}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# BUILD NETWORK AT AG LEVEL ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
if(displayUpdates){message("Calculating network at AG level... \r", appendLF = FALSE)}
# Vector that attributes reach ID to all river network pixels
#print('Define nodes of aggregated network...',quote=FALSE);
Nnodes_AG <- sum(IsNodeAG)
Length_AG <- numeric(Nnodes_AG)
RN_to_AG <- numeric(Nnodes_RN)
AG_to_RN <- vector("list", Nnodes_AG)
reachID <- 1
X_AG <- NaN*numeric(Nnodes_AG)
Y_AG <- NaN*numeric(Nnodes_AG)
Z_AG <- NaN*numeric(Nnodes_AG)
A_AG <- NaN*numeric(Nnodes_AG)
# new version
new_maxLength <- maxReachLength
while (length(whichNodeAG) != 0){ # explore all AG Nodes
i <- whichNodeAG[1] # select the first
RN_to_AG[i] <- reachID
AG_to_RN[[reachID]] <- i
X_AG[reachID] <- X_RN[i]
Y_AG[reachID] <- Y_RN[i]
Z_AG[reachID] <- Z_RN[i]
A_AG[reachID] <- A_RN[i]
if (equalizeLengths){
if (IsNodeAG[i]){ # if the node is pre-established AG node, re-estimate new_maxReachLength; otherwise, use previously calculated value
j <- DownNode_RN[i]
tmp_length <- Length_RN[i]
while (!IsNodeAG[j] && j!=0) {
tmp_length <- tmp_length + Length_RN[j]
j_old <- j
j <- DownNode_RN[j]}
n_splits <- max(1, ceiling((tmp_length+1.5*OCN$cellsize)/maxReachLength)) # it works even if tmp_length < maxReachLength. min is 1, to prevent issues when i=outlet
# +1.5 cellsize to avoid unsplittable reaches
new_maxLength <- min(tmp_length/n_splits + 1.5*OCN$cellsize, maxReachLength) # sum +1.5 cellsize to avoid creating additional reaches due to rounding
new_maxLength <- max(new_maxLength, 1.5*OCN$cellsize)
}
} else {new_maxLength <- maxReachLength}
lengReach <- Length_RN[i]
j <- DownNode_RN[i]
j0 <- j
tmp <- NULL
while (!IsNodeAG[j] && j!=0 && lengReach <= new_maxLength) {
tmp <- c(tmp, j)
lengReach <- lengReach + Length_RN[j]
j_old <- j
j <- DownNode_RN[j]}
if (lengReach > new_maxLength){
j <- j_old
whichNodeAG <- c(whichNodeAG[1], j, whichNodeAG[-1]) # new channel head is added at the front, so that same new_maxLength is applied
# second position (first is to be deleted)
ChannelHeads[j] <- 1
lengReach <- lengReach - Length_RN[j]
tmp <- tmp[-length(tmp)]
}
Length_AG[reachID] <- lengReach
RN_to_AG[tmp] <- reachID
AG_to_RN[[reachID]] <- c(AG_to_RN[[reachID]], tmp)
reachID <- reachID + 1
whichNodeAG <- whichNodeAG[-1]
}
# # old version
# if (equalizeLengths){
# while (length(whichNodeAG) != 0){ # explore all AG Nodes
# i <- whichNodeAG[1] # select the first
# RN_to_AG[i] <- reachID
# AG_to_RN[[reachID]] <- i
# j <- DownNode_RN[i]
# X_AG[reachID] <- X_RN[i]
# Y_AG[reachID] <- Y_RN[i]
# Z_AG[reachID] <- Z_RN[i]
# A_AG[reachID] <- A_RN[i]
# Length_AG[reachID] <- Length_RN[i]
# tmp_length <- Length_RN[i]
# tmp <- NULL
# j0 <- j
# while (!IsNodeAG[j] && j!=0) {
# tmp_length <- tmp_length + Length_RN[j]
# j_old <- j
# j <- DownNode_RN[j]}
#
# if (tmp_length > maxReachLength){
# n_splits <- ceiling(tmp_length/maxReachLength)
# new_maxLength <- tmp_length/n_splits
# new_maxLength <- max(new_maxLength, 1.5*OCN$cellsize)
# j <- j0
# while (!IsNodeAG[j] && j!=0 && Length_AG[reachID] <= new_maxLength) {
# tmp <- c(tmp, j)
# RN_to_AG[j] <- reachID
# AG_to_RN[[reachID]] <- c(AG_to_RN[[reachID]], tmp)
# Length_AG[reachID] <- Length_AG[reachID] + Length_RN[j]
# j_old <- j
# j <- DownNode_RN[j]}
# if (Length_AG[reachID] > new_maxLength){
# j <- j_old
# Length_AG[reachID] <- Length_AG[reachID] - Length_RN[j]
# ChannelHeads[j] <- 1
# whichNodeAG <- c(whichNodeAG,j)}
#
# } else {
# RN_to_AG[tmp] <- reachID
# AG_to_RN[[reachID]] <- c(AG_to_RN[[reachID]], tmp)
# Length_AG[reachID] <- tmp_length
# }
#
# reachID <- reachID + 1
# whichNodeAG <- whichNodeAG[-1]
# }
# } else {
# while (length(whichNodeAG) != 0){ # explore all AG Nodes
# i <- whichNodeAG[1] # select the first
# RN_to_AG[i] <- reachID
# AG_to_RN[[reachID]] <- i
# j <- DownNode_RN[i]
# X_AG[reachID] <- X_RN[i]
# Y_AG[reachID] <- Y_RN[i]
# Z_AG[reachID] <- Z_RN[i]
# A_AG[reachID] <- A_RN[i]
# #Length_AG[reachID] <- Length_RN[i]
# tmp_length <- Length_RN[i]
# tmp <- NULL
# j0 <- j # useless?
# while (!IsNodeAG[j] && j!=0 && tmp_length <= maxReachLength) {
# tmp <- c(tmp, j)
# tmp_length <- tmp_length + Length_RN[j]
# j_old <- j
# j <- DownNode_RN[j]}
# if (tmp_length > maxReachLength){
# j <- j_old
# whichNodeAG <- c(whichNodeAG, j)
# ChannelHeads[j] <- 1
# tmp_length <- tmp_length - Length_RN[j]
# tmp <- tmp[-length(tmp)]
# }
# Length_AG[reachID] <- tmp_length
# RN_to_AG[tmp] <- reachID
# AG_to_RN[[reachID]] <- c(AG_to_RN[[reachID]], tmp)
# reachID <- reachID + 1
# whichNodeAG <- whichNodeAG[-1]
# }
# }
Nnodes_AG <- length(X_AG)
# FD_to_SC: vector of length OCN$FD$nNodes containing subcatchmentID for every pixel of the catchment
# AG_to_FD: list containing FD indices of pixels belonging to a given reach
# SC_to_FD: list containing FD indices of pixels belonging to a given subcatchment
FD_to_SC <- numeric(OCN$FD$nNodes)
# initialize FD_to_SC by attributing SC values to pixels belonging to AG level
FD_to_SC[RN_mask] <- RN_to_AG
# attribute new SC values to pixels corresponding to outlets of catchments without reaches (because the drained area of the catchment is < thrA)
Nnodes_SC <- Nnodes_AG + sum(OCN$FD$A[OCN$FD$outlet]<thrA)
FD_to_SC[OCN$FD$outlet[OCN$FD$A[OCN$FD$outlet] < thrA]] <- (Nnodes_AG+1):Nnodes_SC
IndexHeadpixel <- which(OCN$FD$A==OCN$cellsize^2) # find FD pixels corresponding to headwaters
AG_to_FD <- vector("list", Nnodes_AG)
for(i in 1:Nnodes_AG) { # attribute river network pixels to fields of the AG_to_FD list
AG_to_FD[[i]] <- RN_to_FD[AG_to_RN[[i]]]
}
SC_to_FD <- AG_to_FD[1:Nnodes_AG] # initialize SC_to_FD by attributing the pixels that belong to reaches
# add pixels corresponding to outlets of catchments without reaches
if (Nnodes_SC > Nnodes_AG){
for (i in (Nnodes_AG+1):Nnodes_SC){
SC_to_FD[[i]] <- OCN$FD$outlet[OCN$FD$A[OCN$FD$outlet]<thrA][i-Nnodes_AG]
}}
# for (i in 1:length(IndexHeadpixel)){ # i: index that spans all headwater pixels
# p <- IndexHeadpixel[i] # p: ID of headwater pixel
# pNew <- p; # pNew: pixel downstream of p
# k <- 0; # k: SC value of pixel pNew
# sub_p <- integer(0) # sub_p is the subset of pixels downstream of pixel p
# while (k==0){ # continue downstream movement until a pixel to which the SC has already been attributed is found
# k <- FD_to_SC[pNew]
# if (k==0){
# sub_p <- c(sub_p,pNew)
# pNew <- OCN$FD$downNode[pNew]
# }}
# FD_to_SC[sub_p] <- k
# SC_to_FD[[k]] <- c(SC_to_FD[[k]],sub_p)
# }
ll <- continue_FD_SC(IndexHeadpixel, FD_to_SC, SC_to_FD, OCN$FD$downNode)
FD_to_SC <- ll$FD_to_SC
SC_to_FD <- ll$SC_to_FD
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# CALCULATE PROPERTIES AT AG LEVEL ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
#print('W matrix at AG level...',quote=FALSE);
# Adjacency matrix at reach level
DownNode_AG <- numeric(Nnodes_AG)
W_AG <- spam(0,Nnodes_AG,Nnodes_AG)
ind <- matrix(0,Nnodes_AG,2)
reachID <- sum(ChannelHeads) + 1
for (i in 1:Nnodes_RN){
if (DownNode_RN[i] != 0 && RN_to_AG[DownNode_RN[i]] != RN_to_AG[i]) {
DownNode_AG[RN_to_AG[i]] <- RN_to_AG[DownNode_RN[i]]
ind[RN_to_AG[i],] <- c(RN_to_AG[i],DownNode_AG[RN_to_AG[i]])
}
}
ind <- ind[-which(ind[,1]==0),]
W_AG[ind] <- 1
Outlet_AG <- RN_to_AG[Outlet_RN]
# reverse downNode_AG
DownNode_AG_rev <- vector("list",Nnodes_AG)
for (i in 1:Nnodes_AG){
d <- DownNode_AG[i]
if (d!=0){DownNode_AG_rev[[d]] <- c(DownNode_AG_rev[[d]],i) }}
# Upstream_AG : list containing IDs of all reaches upstream of each reach (plus reach itself)
# sort nodes in downstream direction
ind_sort <- sort(A_AG, index.return=TRUE)
ind_sort <- ind_sort$ix
Upstream_AG <- vector("list",Nnodes_AG)
Nupstream_AG <- numeric(Nnodes_AG)
for (i in 1:Nnodes_AG){
ups <- as.numeric(DownNode_AG_rev[[ind_sort[i]]])
nodes <- numeric(0)
for (u in ups){ nodes <- c(nodes, Upstream_AG[[u]])}
Upstream_AG[[ind_sort[i]]] <- c(nodes, ind_sort[i])
Nupstream_AG[ind_sort[i]] <- length(Upstream_AG[[ind_sort[i]]])
}
# AG_to_CM[i] indicates outlet to which reach i drains
AG_to_CM <- numeric(Nnodes_AG)
for (i in 1:OCN$nOutlet){
AG_to_CM[Upstream_AG[[Outlet_AG[i]]]] <- i
}
ind_sort <- sort(A_AG, index.return=T)
ind_sort <- ind_sort$ix
if (streamOrderType=="Strahler"){
# calculate Strahler stream order
StreamOrder_AG <- numeric(Nnodes_AG)
for (j in ind_sort){
tmp <- DownNode_AG_rev[[j]] # set of reaches draining into j
if (length(tmp)>0){
IncreaseOrder <- sum(StreamOrder_AG[tmp]==max(StreamOrder_AG[tmp])) # check whether tmp reaches have the same stream order
if (IncreaseOrder > 1) {
StreamOrder_AG[j] <- 1 + max(StreamOrder_AG[tmp]) # if so, increase stream order
} else {StreamOrder_AG[j] <- max(StreamOrder_AG[tmp])} # otherwise, keep previous stream order
} else {StreamOrder_AG[j] <- 1} # if j is an headwater, impose StreamOrder = 1
}
} else if (streamOrderType=="Shreve"){
# calculate Shreve stream order
StreamOrder_AG <- numeric(Nnodes_AG)
for (j in ind_sort){
tmp <- DownNode_AG_rev[[j]] # set of reaches draining into j
if (length(tmp)>0){
StreamOrder_AG[j] <- sum(StreamOrder_AG[tmp])
} else {StreamOrder_AG[j] <- 1} # if j is an headwater, impose StreamOrder = 1
}
}
# Calculate slopes of reaches
Slope_AG <- numeric(Nnodes_AG)
for (i in 1:Nnodes_AG){
if (!(i %in% Outlet_AG))
Slope_AG[i] <- (Z_AG[i] - Z_AG[DownNode_AG[i]])/Length_AG[i]
}
if(displayUpdates){message("Calculating network at AG level... 100.0%\n", appendLF = FALSE)}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# CALCULATE PROPERTIES AT SC LEVEL ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
if(displayUpdates){message("Calculating network at SC level... \r", appendLF = FALSE)}
#print(sprintf('Elapsed time %.2f s',difftime(Sys.time(),t1,units='secs')),quote=FALSE)
#t1 <- Sys.time()
#print('Subcatchment properties...',quote=FALSE)
# calculate subcatchment properties: Local Elevation, Local Drained Area, Upstream Area
Z_SC <- numeric(Nnodes_SC)
Alocal_SC <- numeric(Nnodes_SC)
for (i in 1:Nnodes_SC) {
Z_SC[i] <- mean(OCN$FD$Z[SC_to_FD[[i]]])
Alocal_SC[i] <- length(SC_to_FD[[i]])*OCN$cellsize^2
}
# drained area at AG level: note that the first Nnodes_AG elements of Alocal_SC correspond to subcatchments with reaches
# Areach_AG: includes the areas drained by the reaches
Areach_AG <- numeric(Nnodes_AG)
for (i in 1:Nnodes_AG) {
Areach_AG[i] <- sum(Alocal_SC[Upstream_AG[[i]]])
}
# coordinates of AG nodes considered at the downstream end of the respective edge
XReach <- numeric(Nnodes_AG)
YReach <- numeric(Nnodes_AG)
ZReach <- numeric(Nnodes_AG)
for (i in 1:Nnodes_AG){
tmp <- AG_to_RN[[i]]
ind <- which(A_RN[tmp]==max(A_RN[tmp]))
node <- tmp[ind]
XReach[i] <- X_RN[node]
YReach[i] <- Y_RN[node]
ZReach[i] <- Z_RN[node]
}
XReach[Outlet_AG] <- NaN
YReach[Outlet_AG] <- NaN
ZReach[Outlet_AG] <- NaN
# build neighbouring nodes at FD level
# find list of possible neighbouring pixels
movement <- matrix(c(0,-1,-1,-1,0,1,1,1,1,1,0,-1,-1,-1,0,1),nrow=2,byrow=TRUE)
if (length(OCN$typeInitialState)!=0 | OCN$FD$nNodes==OCN$dimX*OCN$dimY){ # all OCNs
NeighbouringNodes <- NN_OCN(OCN$dimX, OCN$dimY, OCN$periodicBoundaries, movement)
# NeighbouringNodes <- vector("list", OCN$dimX*OCN$dimY)
# cont_node <- 0
# for (cc in 1:OCN$dimX) {
# for (rr in 1:OCN$dimY) {
# cont_node <- cont_node + 1
# neigh_r <- rep(rr,8)+movement[1,]
# neigh_c <- rep(cc,8)+movement[2,]
# if (OCN$periodicBoundaries == TRUE){
# neigh_r[neigh_r==0] <- OCN$dimY
# neigh_c[neigh_c==0] <- OCN$dimX
# neigh_r[neigh_r>OCN$dimY] <- 1
# neigh_c[neigh_c>OCN$dimX] <- 1
# }
# NotAboundary <- neigh_r>0 & neigh_r<=OCN$dimY & neigh_c>0 & neigh_c<=OCN$dimX # only effective when periodicBoundaries=FALSE
# NeighbouringNodes[[cont_node]] <- neigh_r[NotAboundary] + (neigh_c[NotAboundary]-1)*OCN$dimY
# }}
} else {
NeighbouringNodes <- NN_river(OCN$dimX, OCN$dimY, OCN$periodicBoundaries, movement, OCN$FD$toDEM, OCN$FD$nNodes)
# NeighbouringNodes <- vector("list", OCN$dimX*OCN$dimY)
# for (i in 1:OCN$FD$nNodes){
# nodeDEM <- OCN$FD$toDEM[i]
# cc <- (nodeDEM %% OCN$dimX); if (cc==0) cc <- OCN$dimX
# rr <- (nodeDEM - cc)/OCN$dimX + 1
# neigh_r <- rep(rr,8)+movement[1,]
# neigh_c <- rep(cc,8)+movement[2,]
# if (OCN$periodicBoundaries == TRUE){
# neigh_r[neigh_r==0] <- OCN$dimY
# neigh_c[neigh_c==0] <- OCN$dimX
# neigh_r[neigh_r>OCN$dimY] <- 1
# neigh_c[neigh_c>OCN$dimX] <- 1
# }
# NotAboundary <- neigh_r>0 & neigh_r<=OCN$dimY & neigh_c>0 & neigh_c<=OCN$dimX # only effective when periodicBoundaries=FALSE
# NeighbouringNodes[[nodeDEM]] <- (neigh_r[NotAboundary]-1)*OCN$dimX + neigh_c[NotAboundary]
# }
}
if (OCN$FD$nNodes < OCN$dimX*OCN$dimY){ # general contour OCNs and real rivers
NeighbouringNodes <- NN_FD(OCN$FD$nNodes, OCN$dimX, OCN$dimY, NeighbouringNodes, OCN$FD$toDEM)
# NeighbouringNodes_FD <- vector("list", OCN$FD$nNodes)
# DEM_to_FD <- numeric(OCN$dimX*OCN$dimY)
# DEM_to_FD[OCN$FD$toDEM] <- 1:OCN$FD$nNodes
# for (i in 1:OCN$FD$nNodes){
# indDEM <- OCN$FD$toDEM[i]
# tmp <- DEM_to_FD[NeighbouringNodes[[indDEM]]]
# NeighbouringNodes_FD[[i]] <- tmp[tmp != 0]
# }
# NeighbouringNodes <- NeighbouringNodes_FD
}
# Subcatchment adjacency matrix: find which subcatchments have borders in common
# W_SC <- spam(0,Nnodes_SC,Nnodes_SC)
# indices <- matrix(0,Nnodes_SC*20,2)
# k <- 1
# for (i in 1:Nnodes_SC){
# set <- SC_to_FD[[i]]
# nodes <- numeric(0)
# for (s in set){ nodes <- union(nodes, FD_to_SC[NeighbouringNodes[[s]]])}
# NeighSubcatch <- setdiff(nodes, i)
# indices[k:(k+length(NeighSubcatch)-1),1] <- i
# indices[k:(k+length(NeighSubcatch)-1),2] <- NeighSubcatch
# k <- k + length(NeighSubcatch)
# if (displayUpdates){
# if ((i %% max(1,round(Nnodes_SC*0.01)))==0){
# message(sprintf("Calculating network at SC level... %.1f%%\r",i/Nnodes_AG*100), appendLF = FALSE)}}
# }
# indices <- indices[1:(k-1),]
# W_SC[indices] <- 1
ll <- WSC(Nnodes_SC,SC_to_FD,FD_to_SC,NeighbouringNodes)
W_SC <- spam(0,Nnodes_SC,Nnodes_SC)
W_SC[cbind(ll[[1]],ll[[2]])] <- 1
# X,Y of subcatchment centroids
X_SC <- numeric(Nnodes_SC)
Y_SC <- numeric(Nnodes_SC)
for (i in 1:Nnodes_SC){
X_SC[i] <- mean(OCN$FD$X[SC_to_FD[[i]]])
Y_SC[i] <- mean(OCN$FD$Y[SC_to_FD[[i]]])
}
if(displayUpdates){message("Calculating network at SC level... 100.0%\n", appendLF = FALSE)}
#%%%%%%%%%%%%%%%%%%%%%#
# EXPORT VARIABLES ####
#%%%%%%%%%%%%%%%%%%%%%#
#FD level
OCN$FD[["toRN"]] <- FD_to_RN
OCN$FD[["toSC"]] <- FD_to_SC
# RN level
OCN$RN[["A"]] <- A_RN
OCN$RN[["W"]] <- W_RN
OCN$RN[["downNode"]] <- DownNode_RN
OCN$RN[["drainageDensity"]] <- DrainageDensity_RN
OCN$RN[["leng"]] <- Length_RN
OCN$RN[["nNodes"]] <- Nnodes_RN
OCN$RN[["nUpstream"]] <- Nupstream_RN
OCN$RN[["outlet"]] <- Outlet_RN
OCN$RN[["slope"]] <- Slope_RN
OCN$RN[["toFD"]] <- RN_to_FD
OCN$RN[["toAGReach"]] <- RN_to_AG
OCN$RN[["toCM"]] <- RN_to_CM
OCN$RN[["upstream"]] <- Upstream_RN
OCN$RN[["X"]] <- X_RN
OCN$RN[["Y"]] <- Y_RN
OCN$RN[["Z"]] <- Z_RN
# AG level
OCN$AG[["A"]] <- A_AG
OCN$AG[["AReach"]] <- Areach_AG
OCN$AG[["W"]] <- W_AG
OCN$AG[["downNode"]] <- DownNode_AG
OCN$AG[["leng"]] <- Length_AG
OCN$AG[["nNodes"]] <- Nnodes_AG
OCN$AG[["nUpstream"]] <- Nupstream_AG
OCN$AG[["outlet"]] <- Outlet_AG
OCN$AG[["slope"]] <- Slope_AG
OCN$AG[["streamOrder"]] <- StreamOrder_AG
OCN$AG[["ReachToFD"]] <- AG_to_FD
OCN$AG[["ReachToRN"]] <- AG_to_RN
OCN$AG[["toCM"]] <- AG_to_CM
OCN$AG[["upstream"]] <- Upstream_AG
OCN$AG[["X"]] <- X_AG
OCN$AG[["XReach"]] <- XReach
OCN$AG[["Y"]] <- Y_AG
OCN$AG[["YReach"]] <- YReach
OCN$AG[["Z"]] <- Z_AG
OCN$AG[["ZReach"]] <- ZReach
# SC level
OCN$SC[["ALocal"]] <- Alocal_SC
OCN$SC[["W"]] <- W_SC
OCN$SC[["nNodes"]] <- Nnodes_SC
OCN$SC[["toFD"]] <- SC_to_FD
OCN$SC[["X"]] <- X_SC
OCN$SC[["Y"]] <- Y_SC
OCN$SC[["Z"]] <- Z_SC
# other
OCN$thrA <- thrA
OCN$streamOrderType <- streamOrderType
OCN$maxReachLength <- maxReachLength
# re-define AG_to_RN, AG_to_FD, RN_to_AG considering AG nodes as pixels and not reaches
AG_to_FDnode <- numeric(Nnodes_AG)
AG_to_RNnode <- numeric(Nnodes_AG)
for (i in 1:Nnodes_AG){
tmpFD <- AG_to_FD[[i]]
AG_to_FDnode[i] <- tmpFD[OCN$FD$A[tmpFD]==min(OCN$FD$A[tmpFD])]
tmpRN <- AG_to_RN[[i]]
AG_to_RNnode[i] <- tmpRN[OCN$RN$A[tmpRN]==min(OCN$RN$A[tmpRN])]
}
RN_to_AGnode <- numeric(Nnodes_RN)
for (i in 1:Nnodes_AG){
RN_to_AGnode[AG_to_RNnode[i]] <- i
}
OCN$RN[["toAG"]] <- RN_to_AGnode
OCN$AG[["toFD"]] <- AG_to_FDnode
OCN$AG[["toRN"]] <- AG_to_RNnode
invisible(OCN)
}
lucarraro/OCNet documentation built on May 1, 2024, 5:24 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions aggregate_OCN
Documented in aggregate_OCN
aggregate_OCN <- function(OCN,
thrA=0.002*OCN$FD$nNodes*OCN$cellsize^2,
streamOrderType="Strahler",
maxReachLength=Inf,
equalizeLengths=FALSE,
breakpoints=NULL,
displayUpdates=FALSE){
if (!("slope" %in% names(OCN$FD))){
stop('Missing fields in OCN. You should run landscape_OCN prior to aggregate_OCN.')
}
if (maxReachLength < OCN$cellsize*sqrt(2)){
stop("maxReachLength cannot be smaller than OCN$cellsize*sqrt(2).")
}
#t1 <- Sys.time()
if (thrA==0) maxReachLength <- OCN$cellsize*sqrt(2)
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# BUILD NETWORK AT RN LEVEL ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
if(displayUpdates){message("Calculating network at RN level... \r", appendLF = FALSE)}
#print('Crop data at FD level to RN level...',quote=FALSE);
RN_mask <- as.vector(OCN$FD$A >= thrA)# RN_mask allows to sample RN-level values from matrices/vectors at FD level
RN_to_FD <- which(RN_mask) # RN_to_FD[i] is the pixel ID at the FD level of the pixel whose ID at the RN level is i
FD_to_RN <- RN_mask*cumsum(as.numeric(RN_mask)) # FD_to_RN[i] is the pixel ID at the RN level of the pixel whose ID at the FD level is i
# if pixel i at FD level doesn't belong to RN, then FD_to_RN[i]=0
Nnodes_RN <- length(RN_to_FD)
W_RN <- OCN$FD$W[RN_mask,,drop=FALSE]
W_RN <- W_RN[,RN_mask,drop=FALSE]
Outlet_RN <- FD_to_RN[OCN$FD$outlet]
Outlet_RN <- Outlet_RN[Outlet_RN!=0] # remove outlets if the corresponding catchment size is lower than threshold
DownNode_RN <- numeric(Nnodes_RN)
# for (i in 1:Nnodes_RN){
# if (!(i %in% Outlet_RN)){
# DownNode_RN[i] <- which(W_RN[i,]==1)
# }}
tmp <- W_RN@rowpointers
NotOutlet <- which((tmp[-1] - tmp[-length(tmp)])==1)
DownNode_RN[NotOutlet] <- W_RN@colindices
# reverse downNode_RN
DownNode_RN_rev <- vector("list",Nnodes_RN)
for (i in 1:Nnodes_RN){
d <- DownNode_RN[i]
if (d!=0){DownNode_RN_rev[[d]] <- c(DownNode_RN_rev[[d]],i) }}
A_RN <- OCN$FD$A[RN_mask]
X_RN <- OCN$FD$X[RN_mask]
Y_RN <- OCN$FD$Y[RN_mask]
Z_RN <- OCN$FD$Z[RN_mask]
Length_RN <- OCN$FD$leng[RN_mask]
# Drainage density
DrainageDensity_RN <- sum(Length_RN)/(OCN$dimX*OCN$dimY*OCN$cellsize^2)
# Connectivity indices at pixel level
DegreeIn <- colSums(W_RN)
DegreeOut <- rowSums(W_RN)
Confluence <- DegreeIn>1
Source <- DegreeIn==0
SourceOrConfluence <- Source|Confluence
ConfluenceNotOutlet <- Confluence&(DownNode_RN!=0)
ChannelHeads <- SourceOrConfluence #Source|ConfluenceNotOutlet
OutletNotChannelHead <- (DownNode_RN==0)&(!ChannelHeads)
IsNodeAG <- SourceOrConfluence|OutletNotChannelHead
IsNodeAG[breakpoints] <- TRUE
whichNodeAG <- which(IsNodeAG)
# Calculate slope for each pixel of the river network
Slope_RN <- OCN$FD$slope[RN_mask]
# sort nodes in downstream direction
ind_sort <- sort(A_RN, index.return=TRUE)
ind_sort <- ind_sort$ix
Upstream_RN <- vector("list",Nnodes_RN)
Nupstream_RN <- numeric(Nnodes_RN)
for (i in 1:Nnodes_RN){
ups <- as.numeric(DownNode_RN_rev[[ind_sort[i]]])
nodes <- numeric(0)
for (u in ups){ nodes <- c(nodes, Upstream_RN[[u]])}
Upstream_RN[[ind_sort[i]]] <- c(nodes, ind_sort[i])
Nupstream_RN[ind_sort[i]] <- length(Upstream_RN[[ind_sort[i]]])
}
# RN_to_CM[i] indicates outlet to which reach i drains
RN_to_CM <- numeric(Nnodes_RN)
for (i in 1:OCN$nOutlet){
RN_to_CM[Upstream_RN[[Outlet_RN[i]]]] <- i
}
if (displayUpdates){message("Calculating network at RN level... 100.0%\n", appendLF = FALSE)}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# BUILD NETWORK AT AG LEVEL ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
if(displayUpdates){message("Calculating network at AG level... \r", appendLF = FALSE)}
# Vector that attributes reach ID to all river network pixels
#print('Define nodes of aggregated network...',quote=FALSE);
Nnodes_AG <- sum(IsNodeAG)
Length_AG <- numeric(Nnodes_AG)
RN_to_AG <- numeric(Nnodes_RN)
AG_to_RN <- vector("list", Nnodes_AG)
reachID <- 1
X_AG <- NaN*numeric(Nnodes_AG)
Y_AG <- NaN*numeric(Nnodes_AG)
Z_AG <- NaN*numeric(Nnodes_AG)
A_AG <- NaN*numeric(Nnodes_AG)
# new version
new_maxLength <- maxReachLength
while (length(whichNodeAG) != 0){ # explore all AG Nodes
i <- whichNodeAG[1] # select the first
RN_to_AG[i] <- reachID
AG_to_RN[[reachID]] <- i
X_AG[reachID] <- X_RN[i]
Y_AG[reachID] <- Y_RN[i]
Z_AG[reachID] <- Z_RN[i]
A_AG[reachID] <- A_RN[i]
if (equalizeLengths){
if (IsNodeAG[i]){ # if the node is pre-established AG node, re-estimate new_maxReachLength; otherwise, use previously calculated value
j <- DownNode_RN[i]
tmp_length <- Length_RN[i]
while (!IsNodeAG[j] && j!=0) {
tmp_length <- tmp_length + Length_RN[j]
j_old <- j
j <- DownNode_RN[j]}
n_splits <- max(1, ceiling((tmp_length+1.5*OCN$cellsize)/maxReachLength)) # it works even if tmp_length < maxReachLength. min is 1, to prevent issues when i=outlet
# +1.5 cellsize to avoid unsplittable reaches
new_maxLength <- min(tmp_length/n_splits + 1.5*OCN$cellsize, maxReachLength) # sum +1.5 cellsize to avoid creating additional reaches due to rounding
new_maxLength <- max(new_maxLength, 1.5*OCN$cellsize)
}
} else {new_maxLength <- maxReachLength}
lengReach <- Length_RN[i]
j <- DownNode_RN[i]
j0 <- j
tmp <- NULL
while (!IsNodeAG[j] && j!=0 && lengReach <= new_maxLength) {
tmp <- c(tmp, j)
lengReach <- lengReach + Length_RN[j]
j_old <- j
j <- DownNode_RN[j]}
if (lengReach > new_maxLength){
j <- j_old
whichNodeAG <- c(whichNodeAG[1], j, whichNodeAG[-1]) # new channel head is added at the front, so that same new_maxLength is applied
# second position (first is to be deleted)
ChannelHeads[j] <- 1
lengReach <- lengReach - Length_RN[j]
tmp <- tmp[-length(tmp)]
}
Length_AG[reachID] <- lengReach
RN_to_AG[tmp] <- reachID
AG_to_RN[[reachID]] <- c(AG_to_RN[[reachID]], tmp)
reachID <- reachID + 1
whichNodeAG <- whichNodeAG[-1]
}
# # old version
# if (equalizeLengths){
# while (length(whichNodeAG) != 0){ # explore all AG Nodes
# i <- whichNodeAG[1] # select the first
# RN_to_AG[i] <- reachID
# AG_to_RN[[reachID]] <- i
# j <- DownNode_RN[i]
# X_AG[reachID] <- X_RN[i]
# Y_AG[reachID] <- Y_RN[i]
# Z_AG[reachID] <- Z_RN[i]
# A_AG[reachID] <- A_RN[i]
# Length_AG[reachID] <- Length_RN[i]
# tmp_length <- Length_RN[i]
# tmp <- NULL
# j0 <- j
# while (!IsNodeAG[j] && j!=0) {
# tmp_length <- tmp_length + Length_RN[j]
# j_old <- j
# j <- DownNode_RN[j]}
#
# if (tmp_length > maxReachLength){
# n_splits <- ceiling(tmp_length/maxReachLength)
# new_maxLength <- tmp_length/n_splits
# new_maxLength <- max(new_maxLength, 1.5*OCN$cellsize)
# j <- j0
# while (!IsNodeAG[j] && j!=0 && Length_AG[reachID] <= new_maxLength) {
# tmp <- c(tmp, j)
# RN_to_AG[j] <- reachID
# AG_to_RN[[reachID]] <- c(AG_to_RN[[reachID]], tmp)
# Length_AG[reachID] <- Length_AG[reachID] + Length_RN[j]
# j_old <- j
# j <- DownNode_RN[j]}
# if (Length_AG[reachID] > new_maxLength){
# j <- j_old
# Length_AG[reachID] <- Length_AG[reachID] - Length_RN[j]
# ChannelHeads[j] <- 1
# whichNodeAG <- c(whichNodeAG,j)}
#
# } else {
# RN_to_AG[tmp] <- reachID
# AG_to_RN[[reachID]] <- c(AG_to_RN[[reachID]], tmp)
# Length_AG[reachID] <- tmp_length
# }
#
# reachID <- reachID + 1
# whichNodeAG <- whichNodeAG[-1]
# }
# } else {
# while (length(whichNodeAG) != 0){ # explore all AG Nodes
# i <- whichNodeAG[1] # select the first
# RN_to_AG[i] <- reachID
# AG_to_RN[[reachID]] <- i
# j <- DownNode_RN[i]
# X_AG[reachID] <- X_RN[i]
# Y_AG[reachID] <- Y_RN[i]
# Z_AG[reachID] <- Z_RN[i]
# A_AG[reachID] <- A_RN[i]
# #Length_AG[reachID] <- Length_RN[i]
# tmp_length <- Length_RN[i]
# tmp <- NULL
# j0 <- j # useless?
# while (!IsNodeAG[j] && j!=0 && tmp_length <= maxReachLength) {
# tmp <- c(tmp, j)
# tmp_length <- tmp_length + Length_RN[j]
# j_old <- j
# j <- DownNode_RN[j]}
# if (tmp_length > maxReachLength){
# j <- j_old
# whichNodeAG <- c(whichNodeAG, j)
# ChannelHeads[j] <- 1
# tmp_length <- tmp_length - Length_RN[j]
# tmp <- tmp[-length(tmp)]
# }
# Length_AG[reachID] <- tmp_length
# RN_to_AG[tmp] <- reachID
# AG_to_RN[[reachID]] <- c(AG_to_RN[[reachID]], tmp)
# reachID <- reachID + 1
# whichNodeAG <- whichNodeAG[-1]
# }
# }
Nnodes_AG <- length(X_AG)
# FD_to_SC: vector of length OCN$FD$nNodes containing subcatchmentID for every pixel of the catchment
# AG_to_FD: list containing FD indices of pixels belonging to a given reach
# SC_to_FD: list containing FD indices of pixels belonging to a given subcatchment
FD_to_SC <- numeric(OCN$FD$nNodes)
# initialize FD_to_SC by attributing SC values to pixels belonging to AG level
FD_to_SC[RN_mask] <- RN_to_AG
# attribute new SC values to pixels corresponding to outlets of catchments without reaches (because the drained area of the catchment is < thrA)
Nnodes_SC <- Nnodes_AG + sum(OCN$FD$A[OCN$FD$outlet]<thrA)
FD_to_SC[OCN$FD$outlet[OCN$FD$A[OCN$FD$outlet] < thrA]] <- (Nnodes_AG+1):Nnodes_SC
IndexHeadpixel <- which(OCN$FD$A==OCN$cellsize^2) # find FD pixels corresponding to headwaters
AG_to_FD <- vector("list", Nnodes_AG)
for(i in 1:Nnodes_AG) { # attribute river network pixels to fields of the AG_to_FD list
AG_to_FD[[i]] <- RN_to_FD[AG_to_RN[[i]]]
}
SC_to_FD <- AG_to_FD[1:Nnodes_AG] # initialize SC_to_FD by attributing the pixels that belong to reaches
# add pixels corresponding to outlets of catchments without reaches
if (Nnodes_SC > Nnodes_AG){
for (i in (Nnodes_AG+1):Nnodes_SC){
SC_to_FD[[i]] <- OCN$FD$outlet[OCN$FD$A[OCN$FD$outlet]<thrA][i-Nnodes_AG]
}}
# for (i in 1:length(IndexHeadpixel)){ # i: index that spans all headwater pixels
# p <- IndexHeadpixel[i] # p: ID of headwater pixel
# pNew <- p; # pNew: pixel downstream of p
# k <- 0; # k: SC value of pixel pNew
# sub_p <- integer(0) # sub_p is the subset of pixels downstream of pixel p
# while (k==0){ # continue downstream movement until a pixel to which the SC has already been attributed is found
# k <- FD_to_SC[pNew]
# if (k==0){
# sub_p <- c(sub_p,pNew)
# pNew <- OCN$FD$downNode[pNew]
# }}
# FD_to_SC[sub_p] <- k
# SC_to_FD[[k]] <- c(SC_to_FD[[k]],sub_p)
# }
ll <- continue_FD_SC(IndexHeadpixel, FD_to_SC, SC_to_FD, OCN$FD$downNode)
FD_to_SC <- ll$FD_to_SC
SC_to_FD <- ll$SC_to_FD
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# CALCULATE PROPERTIES AT AG LEVEL ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
#print('W matrix at AG level...',quote=FALSE);
# Adjacency matrix at reach level
DownNode_AG <- numeric(Nnodes_AG)
W_AG <- spam(0,Nnodes_AG,Nnodes_AG)
ind <- matrix(0,Nnodes_AG,2)
reachID <- sum(ChannelHeads) + 1
for (i in 1:Nnodes_RN){
if (DownNode_RN[i] != 0 && RN_to_AG[DownNode_RN[i]] != RN_to_AG[i]) {
DownNode_AG[RN_to_AG[i]] <- RN_to_AG[DownNode_RN[i]]
ind[RN_to_AG[i],] <- c(RN_to_AG[i],DownNode_AG[RN_to_AG[i]])
}
}
ind <- ind[-which(ind[,1]==0),]
W_AG[ind] <- 1
Outlet_AG <- RN_to_AG[Outlet_RN]
# reverse downNode_AG
DownNode_AG_rev <- vector("list",Nnodes_AG)
for (i in 1:Nnodes_AG){
d <- DownNode_AG[i]
if (d!=0){DownNode_AG_rev[[d]] <- c(DownNode_AG_rev[[d]],i) }}
# Upstream_AG : list containing IDs of all reaches upstream of each reach (plus reach itself)
# sort nodes in downstream direction
ind_sort <- sort(A_AG, index.return=TRUE)
ind_sort <- ind_sort$ix
Upstream_AG <- vector("list",Nnodes_AG)
Nupstream_AG <- numeric(Nnodes_AG)
for (i in 1:Nnodes_AG){
ups <- as.numeric(DownNode_AG_rev[[ind_sort[i]]])
nodes <- numeric(0)
for (u in ups){ nodes <- c(nodes, Upstream_AG[[u]])}
Upstream_AG[[ind_sort[i]]] <- c(nodes, ind_sort[i])
Nupstream_AG[ind_sort[i]] <- length(Upstream_AG[[ind_sort[i]]])
}
# AG_to_CM[i] indicates outlet to which reach i drains
AG_to_CM <- numeric(Nnodes_AG)
for (i in 1:OCN$nOutlet){
AG_to_CM[Upstream_AG[[Outlet_AG[i]]]] <- i
}
ind_sort <- sort(A_AG, index.return=T)
ind_sort <- ind_sort$ix
if (streamOrderType=="Strahler"){
# calculate Strahler stream order
StreamOrder_AG <- numeric(Nnodes_AG)
for (j in ind_sort){
tmp <- DownNode_AG_rev[[j]] # set of reaches draining into j
if (length(tmp)>0){
IncreaseOrder <- sum(StreamOrder_AG[tmp]==max(StreamOrder_AG[tmp])) # check whether tmp reaches have the same stream order
if (IncreaseOrder > 1) {
StreamOrder_AG[j] <- 1 + max(StreamOrder_AG[tmp]) # if so, increase stream order
} else {StreamOrder_AG[j] <- max(StreamOrder_AG[tmp])} # otherwise, keep previous stream order
} else {StreamOrder_AG[j] <- 1} # if j is an headwater, impose StreamOrder = 1
}
} else if (streamOrderType=="Shreve"){
# calculate Shreve stream order
StreamOrder_AG <- numeric(Nnodes_AG)
for (j in ind_sort){
tmp <- DownNode_AG_rev[[j]] # set of reaches draining into j
if (length(tmp)>0){
StreamOrder_AG[j] <- sum(StreamOrder_AG[tmp])
} else {StreamOrder_AG[j] <- 1} # if j is an headwater, impose StreamOrder = 1
}
}
# Calculate slopes of reaches
Slope_AG <- numeric(Nnodes_AG)
for (i in 1:Nnodes_AG){
if (!(i %in% Outlet_AG))
Slope_AG[i] <- (Z_AG[i] - Z_AG[DownNode_AG[i]])/Length_AG[i]
}
if(displayUpdates){message("Calculating network at AG level... 100.0%\n", appendLF = FALSE)}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# CALCULATE PROPERTIES AT SC LEVEL ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
if(displayUpdates){message("Calculating network at SC level... \r", appendLF = FALSE)}
#print(sprintf('Elapsed time %.2f s',difftime(Sys.time(),t1,units='secs')),quote=FALSE)
#t1 <- Sys.time()
#print('Subcatchment properties...',quote=FALSE)
# calculate subcatchment properties: Local Elevation, Local Drained Area, Upstream Area
Z_SC <- numeric(Nnodes_SC)
Alocal_SC <- numeric(Nnodes_SC)
for (i in 1:Nnodes_SC) {
Z_SC[i] <- mean(OCN$FD$Z[SC_to_FD[[i]]])
Alocal_SC[i] <- length(SC_to_FD[[i]])*OCN$cellsize^2
}
# drained area at AG level: note that the first Nnodes_AG elements of Alocal_SC correspond to subcatchments with reaches
# Areach_AG: includes the areas drained by the reaches
Areach_AG <- numeric(Nnodes_AG)
for (i in 1:Nnodes_AG) {
Areach_AG[i] <- sum(Alocal_SC[Upstream_AG[[i]]])
}
# coordinates of AG nodes considered at the downstream end of the respective edge
XReach <- numeric(Nnodes_AG)
YReach <- numeric(Nnodes_AG)
ZReach <- numeric(Nnodes_AG)
for (i in 1:Nnodes_AG){
tmp <- AG_to_RN[[i]]
ind <- which(A_RN[tmp]==max(A_RN[tmp]))
node <- tmp[ind]
XReach[i] <- X_RN[node]
YReach[i] <- Y_RN[node]
ZReach[i] <- Z_RN[node]
}
XReach[Outlet_AG] <- NaN
YReach[Outlet_AG] <- NaN
ZReach[Outlet_AG] <- NaN
# build neighbouring nodes at FD level
# find list of possible neighbouring pixels
movement <- matrix(c(0,-1,-1,-1,0,1,1,1,1,1,0,-1,-1,-1,0,1),nrow=2,byrow=TRUE)
if (length(OCN$typeInitialState)!=0 | OCN$FD$nNodes==OCN$dimX*OCN$dimY){ # all OCNs
NeighbouringNodes <- NN_OCN(OCN$dimX, OCN$dimY, OCN$periodicBoundaries, movement)
# NeighbouringNodes <- vector("list", OCN$dimX*OCN$dimY)
# cont_node <- 0
# for (cc in 1:OCN$dimX) {
# for (rr in 1:OCN$dimY) {
# cont_node <- cont_node + 1
# neigh_r <- rep(rr,8)+movement[1,]
# neigh_c <- rep(cc,8)+movement[2,]
# if (OCN$periodicBoundaries == TRUE){
# neigh_r[neigh_r==0] <- OCN$dimY
# neigh_c[neigh_c==0] <- OCN$dimX
# neigh_r[neigh_r>OCN$dimY] <- 1
# neigh_c[neigh_c>OCN$dimX] <- 1
# }
# NotAboundary <- neigh_r>0 & neigh_r<=OCN$dimY & neigh_c>0 & neigh_c<=OCN$dimX # only effective when periodicBoundaries=FALSE
# NeighbouringNodes[[cont_node]] <- neigh_r[NotAboundary] + (neigh_c[NotAboundary]-1)*OCN$dimY
# }}
} else {
NeighbouringNodes <- NN_river(OCN$dimX, OCN$dimY, OCN$periodicBoundaries, movement, OCN$FD$toDEM, OCN$FD$nNodes)
# NeighbouringNodes <- vector("list", OCN$dimX*OCN$dimY)
# for (i in 1:OCN$FD$nNodes){
# nodeDEM <- OCN$FD$toDEM[i]
# cc <- (nodeDEM %% OCN$dimX); if (cc==0) cc <- OCN$dimX
# rr <- (nodeDEM - cc)/OCN$dimX + 1
# neigh_r <- rep(rr,8)+movement[1,]
# neigh_c <- rep(cc,8)+movement[2,]
# if (OCN$periodicBoundaries == TRUE){
# neigh_r[neigh_r==0] <- OCN$dimY
# neigh_c[neigh_c==0] <- OCN$dimX
# neigh_r[neigh_r>OCN$dimY] <- 1
# neigh_c[neigh_c>OCN$dimX] <- 1
# }
# NotAboundary <- neigh_r>0 & neigh_r<=OCN$dimY & neigh_c>0 & neigh_c<=OCN$dimX # only effective when periodicBoundaries=FALSE
# NeighbouringNodes[[nodeDEM]] <- (neigh_r[NotAboundary]-1)*OCN$dimX + neigh_c[NotAboundary]
# }
}
if (OCN$FD$nNodes < OCN$dimX*OCN$dimY){ # general contour OCNs and real rivers
NeighbouringNodes <- NN_FD(OCN$FD$nNodes, OCN$dimX, OCN$dimY, NeighbouringNodes, OCN$FD$toDEM)
# NeighbouringNodes_FD <- vector("list", OCN$FD$nNodes)
# DEM_to_FD <- numeric(OCN$dimX*OCN$dimY)
# DEM_to_FD[OCN$FD$toDEM] <- 1:OCN$FD$nNodes
# for (i in 1:OCN$FD$nNodes){
# indDEM <- OCN$FD$toDEM[i]
# tmp <- DEM_to_FD[NeighbouringNodes[[indDEM]]]
# NeighbouringNodes_FD[[i]] <- tmp[tmp != 0]
# }
# NeighbouringNodes <- NeighbouringNodes_FD
}
# Subcatchment adjacency matrix: find which subcatchments have borders in common
# W_SC <- spam(0,Nnodes_SC,Nnodes_SC)
# indices <- matrix(0,Nnodes_SC*20,2)
# k <- 1
# for (i in 1:Nnodes_SC){
# set <- SC_to_FD[[i]]
# nodes <- numeric(0)
# for (s in set){ nodes <- union(nodes, FD_to_SC[NeighbouringNodes[[s]]])}
# NeighSubcatch <- setdiff(nodes, i)
# indices[k:(k+length(NeighSubcatch)-1),1] <- i
# indices[k:(k+length(NeighSubcatch)-1),2] <- NeighSubcatch
# k <- k + length(NeighSubcatch)
# if (displayUpdates){
# if ((i %% max(1,round(Nnodes_SC*0.01)))==0){
# message(sprintf("Calculating network at SC level... %.1f%%\r",i/Nnodes_AG*100), appendLF = FALSE)}}
# }
# indices <- indices[1:(k-1),]
# W_SC[indices] <- 1
ll <- WSC(Nnodes_SC,SC_to_FD,FD_to_SC,NeighbouringNodes)
W_SC <- spam(0,Nnodes_SC,Nnodes_SC)
W_SC[cbind(ll[[1]],ll[[2]])] <- 1
# X,Y of subcatchment centroids
X_SC <- numeric(Nnodes_SC)
Y_SC <- numeric(Nnodes_SC)
for (i in 1:Nnodes_SC){
X_SC[i] <- mean(OCN$FD$X[SC_to_FD[[i]]])
Y_SC[i] <- mean(OCN$FD$Y[SC_to_FD[[i]]])
}
if(displayUpdates){message("Calculating network at SC level... 100.0%\n", appendLF = FALSE)}
#%%%%%%%%%%%%%%%%%%%%%#
# EXPORT VARIABLES ####
#%%%%%%%%%%%%%%%%%%%%%#
#FD level
OCN$FD[["toRN"]] <- FD_to_RN
OCN$FD[["toSC"]] <- FD_to_SC
# RN level
OCN$RN[["A"]] <- A_RN
OCN$RN[["W"]] <- W_RN
OCN$RN[["downNode"]] <- DownNode_RN
OCN$RN[["drainageDensity"]] <- DrainageDensity_RN
OCN$RN[["leng"]] <- Length_RN
OCN$RN[["nNodes"]] <- Nnodes_RN
OCN$RN[["nUpstream"]] <- Nupstream_RN
OCN$RN[["outlet"]] <- Outlet_RN
OCN$RN[["slope"]] <- Slope_RN
OCN$RN[["toFD"]] <- RN_to_FD
OCN$RN[["toAGReach"]] <- RN_to_AG
OCN$RN[["toCM"]] <- RN_to_CM
OCN$RN[["upstream"]] <- Upstream_RN
OCN$RN[["X"]] <- X_RN
OCN$RN[["Y"]] <- Y_RN
OCN$RN[["Z"]] <- Z_RN
# AG level
OCN$AG[["A"]] <- A_AG
OCN$AG[["AReach"]] <- Areach_AG
OCN$AG[["W"]] <- W_AG
OCN$AG[["downNode"]] <- DownNode_AG
OCN$AG[["leng"]] <- Length_AG
OCN$AG[["nNodes"]] <- Nnodes_AG
OCN$AG[["nUpstream"]] <- Nupstream_AG
OCN$AG[["outlet"]] <- Outlet_AG
OCN$AG[["slope"]] <- Slope_AG
OCN$AG[["streamOrder"]] <- StreamOrder_AG
OCN$AG[["ReachToFD"]] <- AG_to_FD
OCN$AG[["ReachToRN"]] <- AG_to_RN
OCN$AG[["toCM"]] <- AG_to_CM
OCN$AG[["upstream"]] <- Upstream_AG
OCN$AG[["X"]] <- X_AG
OCN$AG[["XReach"]] <- XReach
OCN$AG[["Y"]] <- Y_AG
OCN$AG[["YReach"]] <- YReach
OCN$AG[["Z"]] <- Z_AG
OCN$AG[["ZReach"]] <- ZReach
# SC level
OCN$SC[["ALocal"]] <- Alocal_SC
OCN$SC[["W"]] <- W_SC
OCN$SC[["nNodes"]] <- Nnodes_SC
OCN$SC[["toFD"]] <- SC_to_FD
OCN$SC[["X"]] <- X_SC
OCN$SC[["Y"]] <- Y_SC
OCN$SC[["Z"]] <- Z_SC
# other
OCN$thrA <- thrA
OCN$streamOrderType <- streamOrderType
OCN$maxReachLength <- maxReachLength
# re-define AG_to_RN, AG_to_FD, RN_to_AG considering AG nodes as pixels and not reaches
AG_to_FDnode <- numeric(Nnodes_AG)
AG_to_RNnode <- numeric(Nnodes_AG)
for (i in 1:Nnodes_AG){
tmpFD <- AG_to_FD[[i]]
AG_to_FDnode[i] <- tmpFD[OCN$FD$A[tmpFD]==min(OCN$FD$A[tmpFD])]
tmpRN <- AG_to_RN[[i]]
AG_to_RNnode[i] <- tmpRN[OCN$RN$A[tmpRN]==min(OCN$RN$A[tmpRN])]
}
RN_to_AGnode <- numeric(Nnodes_RN)
for (i in 1:Nnodes_AG){
RN_to_AGnode[AG_to_RNnode[i]] <- i
}
OCN$RN[["toAG"]] <- RN_to_AGnode
OCN$AG[["toFD"]] <- AG_to_FDnode
OCN$AG[["toRN"]] <- AG_to_RNnode
invisible(OCN)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.