library(sf)

## Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3

library(igraph)

## 
## Attaching package: 'igraph'

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

# Read
dat = st_read("water system.shp", stringsAsFactors = FALSE)

## Reading layer `water system' from data source `/home/michael/Dropbox/BGU/1_Other/Oren_Yosef/p_01_dissolve_and_network_cluster/water system.shp' using driver `ESRI Shapefile'
## Simple feature collection with 419 features and 3 fields
## geometry type:  MULTILINESTRING
## dimension:      XY
## bbox:           xmin: 179180.7 ymin: 645988.5 xmax: 182479.4 ymax: 649840.6
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=31.73439361111111 +lon_0=35.20451694444445 +k=1.0000067 +x_0=219529.584 +y_0=626907.39 +ellps=GRS80 +units=m +no_defs

plot(dat)

# To MULTILINESTRING to LINESTRING
dat = st_cast(dat, "LINESTRING")

## Warning in st_cast.sf(dat, "LINESTRING"): repeating attributes for all sub-
## geometries for which they may not be constant

for(i in unique(dat$KAV_DIAM)) {

  s = dat$KAV_DIAM == i
  
  # Create graph based on intersection
  n = st_intersects(dat[s, ], dat[s, ], sparse = FALSE)
  g = graph_from_adjacency_matrix(n, mode = "undirected")
  
  # Find network components and set IDs
  cluster = clusters(g)
  dat$id[s] = paste(cluster$membership, dat$KAV_DIAM[s], sep = "+")

}

plot(dat[, "id"])

# "Dissolve"
dat = aggregate(
  dat[, "KAV_DIAM"],
  list(id = dat$id),
  FUN = mean
)

# Write
st_write(dat, "result.shp", delete_dsn = TRUE)

## Deleting source `/home/michael/Dropbox/BGU/1_Other/Oren_Yosef/p_01_dissolve_and_network_cluster/result.shp' using driver `ESRI Shapefile'
## Writing layer `result' to data source `/home/michael/Dropbox/BGU/1_Other/Oren_Yosef/p_01_dissolve_and_network_cluster/result.shp' using driver `ESRI Shapefile'
## features:       983
## fields:         2
## geometry type:  Unknown (any)