R/actualSARN.R

In craigbrinkerhoff/SARNr: SARN algorithm in R for any resolution imagery

##############################
## FUNCTIONS FOR SARNr
## Winter 2022
## Craig Brinkerhoff
##############################


#' SARN network joining
#' 
#' Joins DEM and RS river networks
#' 
#' @param  data: sarndata object generated by sarn_data()
#' @param bufferSize: Lateral buffer distance for joining river networks (meters)
#' 
#' @import dplyr
#' @import terra
#' @import sf
#' 
#' @return sarnnetwork object with two vector river networks: DEM/RS 'coincident' rivers and DEM/RS 'not coincident' rivers
#' 
#' @export sarn_joinNetworks
sarn_joinNetworks <- function(data, bufferSize) {
  print("Joining DEM and RS networks...")
  
  ## PREPROCESSING
  rivnet_dem_buf <- st_buffer(data$dem_network, dist=bufferSize, endCapStyle='ROUND')
  #rivnet_rs <- st_collection_extract(st_intersection(data$rs_network, rivnet_dem_buf), type='LINESTRING')
  
  #SNAP DEM TO RIVERS (used in orig. SARN algorithm but not used here)
  # This is necessary to that DEM river pixels at now adjacent to RS rive r pixels when we do raster algebra
#  rivnet_dem_snap <- st_snap(data$dem_network, data$rs_network, bufferSize)
#  rivnet_dem_snap_buf <- st_buffer(rivnet_dem_snap, dist=bufferSize, endCapStyle = 'ROUND')
  
  v <- vect(rivnet_dem_buf) #necessary and dumb conversion from sf vector to terra vector...
  rivnet_dem_buf_rast <- rast(extent = ext(data$riverMask), crs = crs(v), resolution=res(data$riverMask)) #use extent from RS network
  rivnet_dem_buf_rast <- rasterize(x=v, y=rivnet_dem_buf_rast)
  rivnet_combo_rast <- (rivnet_dem_buf_rast+1) + data$riverMask #3 equals dem+rs synonymous rivers
  
  #only keep coincident pixels
  rivnet_combo_rast[rivnet_combo_rast<3] <- NA

  rivnet_combo <- as.polygons(rivnet_combo_rast)
  rivnet_combo <- st_as_sf(rivnet_combo) #back to SF...
  rivnet_combo <- st_make_valid(rivnet_combo) #repair geometries
  
  #GET 'COINCIDENT NETWORK': DEM WITH A CORRESPONDING RS RIVER
  coincident_network <- st_intersection(rivnet_combo, data$dem_network)
  coincident_network <- st_collection_extract(coincident_network, "LINESTRING") #can sometimes return points if line only intersects polygon at a point, so we just remove those
  coincident_network <- dplyr::select(coincident_network, c('cat', 'geometry'))
  colnames(coincident_network) <- c('reachID', 'geometry')
  coincident_network <- do.call(rbind,lapply(1:nrow(coincident_network),function(i){st_cast(coincident_network[i,],"LINESTRING")}))
  
  #GET 'NON-COINCIDENT NETWORK': DEM WITH NO CORRESPONDING RS RIVER
  non_coincident_network <- st_erase(data$dem_network, rivnet_combo)
  non_coincident_network <- st_collection_extract(non_coincident_network, "LINESTRING") #can sometimes return points if line only intersects polygon at a point, so we just remove those
  non_coincident_network <- dplyr::select(non_coincident_network, c('cat', 'geometry'))
  colnames(non_coincident_network) <- c('reachID', 'geometry')
  non_coincident_network <- do.call(rbind,lapply(1:nrow(non_coincident_network),function(i){st_cast(non_coincident_network[i,],"LINESTRING")}))
  
  datalist <- list(coincidentNetwork_init = coincident_network,
              not_coincidentNetwork_init = non_coincident_network,
              dem = data$dem,
              trimmedNetwork_fin = NA,
              trimmedFlag = 'No')
  
  out <- structure(c(datalist),
                   class = c("sarnnetwork"))
  
  return(out)
}


#' SARN network clipping
#' 
#' Iteratively trims river network to the actively-flowing network
#' 
#' @param networks: sarnnetwork object produced by sarn_joinNetworks()
#' @param crs_code: epsg code for the projection you are using (for the trimmed nodes)
#' @param printOutput: Do you want the results of the iterative trimming to print to console? Default= Yes
#' 
#' @import sf
#' @import dplyr
#' 
#' @return sarnnetwork object with trimmed river network (output$trimmedNetwork)
#'
#' @export sarn_trimNetworks
sarn_trimNetworks <- function(networks, crs_code, printOutput="Yes") {
  print('Giving the river network a much needed hair cut!!')
  
  printOutput <- ifelse(printOutput == 'yes' || printOutput == 'Yes', 1, 0)
  
  #use temporary networks to do the iterative trimming
  coincidentNetwork <- networks$coincidentNetwork_init
  not_coincidentNetwork <- networks$not_coincidentNetwork_init

  #TRIM NETWORK USING NODE TOPOLOGY
  #Identify initial RS end points
  RS_points1 <- st_line_sample(coincidentNetwork, sample = 0)
  RS_points2 <- st_line_sample(coincidentNetwork, sample = 1)
  
  RS_points1 <- st_as_sfc(do.call(rbind,lapply(1:length(RS_points1),function(i){st_cast(RS_points1[i],"POINT")})), crs=crs_code)
  RS_points2 <- st_as_sfc(do.call(rbind,lapply(1:length(RS_points2),function(i){st_cast(RS_points2[i],"POINT")})), crs=crs_code)
  RS_points <- st_as_sfc(rbind(RS_points1, RS_points2), crs=crs_code) %>% 
    st_sf %>% 
    st_cast
  
  #iteratively trim the network back until the RS end points are 1st order
  flag <- 1
  while (flag != 0) {
    #get DEM reach nodes (terminal start nodes, terminal end nodes, and intermediate nodes which are those that are simultaneously start/end nodes)
    DEM_points_start <- st_line_sample(not_coincidentNetwork, sample=0)
    DEM_points_end <- st_line_sample(not_coincidentNetwork, sample=1)
    DEM_points_intermediate <- st_intersection(DEM_points_start, DEM_points_end)
    
    DEM_points_intermediate <- st_as_sfc(do.call(rbind,lapply(1:length(DEM_points_intermediate),function(i){st_cast(DEM_points_intermediate[i],"POINT")})), crs=crs_code) %>% 
      st_sf %>% 
      st_cast
    DEM_points_start <- st_as_sfc(do.call(rbind,lapply(1:length(DEM_points_start),function(i){st_cast(DEM_points_start[i],"POINT")})), crs=crs_code) %>% 
      st_sf %>% 
      st_cast
    DEM_points_end <- st_as_sfc(do.call(rbind,lapply(1:length(DEM_points_end),function(i){st_cast(DEM_points_end[i],"POINT")})), crs=crs_code) %>% 
      st_sf %>%
      st_cast
    DEM_points_bothEnds <- rbind(DEM_points_start, DEM_points_end)
    
    #identify nodes to trim to. These correspond to terminal start nodes that are neither 1) RS-visible nor 2) intermediate DEM nodes
    #(remember that the terminal/intermediate node distinction is relative to the current algorithm iteration)
    trimmingNodes <- st_erase(DEM_points_bothEnds, DEM_points_intermediate)
    trimmingNodes <- st_erase(trimmingNodes, RS_points)
    
    #Check if there are still nodes left to remove (i.e. does this file exist). If not, break the algorithm
    flag <- nrow(trimmingNodes)
    
    #trim not-coincident network by removing reaches associated with the terminal nodes just identified
    trimmedReaches <- (st_intersects(not_coincidentNetwork, trimmingNodes, sparse=F) *1)
    trimmedReaches <- rowSums(trimmedReaches, na.rm=T) #using some boolean math to invert this since i can't 'erase' lines using nodes apparently
    not_coincidentNetwork <- not_coincidentNetwork[which(trimmedReaches == 0),]
    
    if(printOutput == 1){
      print(paste0('trimmed ', flag, ' reaches'))
    }
  }
  
  coincidentNetwork$type <- 1 #RS
  not_coincidentNetwork$type <- 0 #DEM
  
  #MERGE COINCIDENT AND NON-COINCIDENT NETWORKS INTO SINGLE HYDROGRAPHY PRODUCT
  rivnet_fin <- rbind(coincidentNetwork, not_coincidentNetwork)
  
  #FINALIZE
  networks$trimmedNetwork_fin <- rivnet_fin
  networks$trimmedFlag <- 'Yes'
  
  return(networks)
}