#Stream Network Analysis by Narayan Adhikari
#Importing stream network in r
Stream_Network <- st_read("D:/Final_Seminar/Stream_Network.shp")#read the data## Reading layer `Stream_Network' from data source `D:\Final_Seminar\Stream_Network.shp' using driver `ESRI Shapefile'
## Simple feature collection with 687 features and 4 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: -85.00014 ymin: 34.99986 xmax: -83.99986 ymax: 36.00014
## geographic CRS: WGS 84class(Stream_Network)#cheking class## [1] "sf" "data.frame"str(Stream_Network)#structure of Stream_Network## Classes 'sf' and 'data.frame': 687 obs. of 5 variables:
## $ arcid : num 1 2 3 4 5 6 7 8 9 10 ...
## $ grid_code: num 1 3 1 1 3 1 1 1 1 2 ...
## $ from_node: num 8 9 10 12 11 13 19 20 22 23 ...
## $ to_node : num 3 1 9 11 6 7 2 21 23 4 ...
## $ geometry :sfc_LINESTRING of length 687; first list element: 'XY' num [1:4, 1:2] -84.8 -84.8 -84.8 -84.8 36 ...
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA
## ..- attr(*, "names")= chr [1:4] "arcid" "grid_code" "from_node" "to_node"st_crs(Stream_Network)#Coordinate Reference System## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]nrow(Stream_Network)#no. of rows## [1] 687ncol(Stream_Network)#no. of column## [1] 5colnames(Stream_Network)## [1] "arcid" "grid_code" "from_node" "to_node" "geometry"plot(Stream_Network)
summary(Stream_Network)#getting summary ## arcid grid_code from_node to_node
## Min. : 1.0 Min. :1.00 Min. : 8.0 Min. : 1.0
## 1st Qu.:172.5 1st Qu.:1.00 1st Qu.:182.5 1st Qu.:180.5
## Median :344.0 Median :1.00 Median :356.0 Median :361.0
## Mean :344.0 Mean :1.83 Mean :357.5 Mean :356.6
## 3rd Qu.:515.5 3rd Qu.:2.00 3rd Qu.:533.5 3rd Qu.:535.0
## Max. :687.0 Max. :5.00 Max. :705.0 Max. :715.0
## geometry
## LINESTRING :687
## epsg:4326 : 0
## +proj=long...: 0
##
##
## NACharacters <- Stream_Network [is.na(Stream_Network), ]#checking if there is na character in the data set
NACharacters## Simple feature collection with 0 features and 4 fields
## bbox: xmin: NA ymin: NA xmax: NA ymax: NA
## geographic CRS: WGS 84
## [1] arcid grid_code from_node to_node geometry
## <0 rows> (or 0-length row.names)#Adding the watershed polygon of interest
watershed <- st_read("D:/Final_Seminar/watershed.shp")#read the data## Reading layer `watershed' from data source `D:\Final_Seminar\watershed.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -85.00014 ymin: 35.40605 xmax: -84.19249 ymax: 36.00014
## geographic CRS: WGS 84class(watershed)## [1] "sf" "data.frame"st_crs(watershed)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]#Chainging CRS of the shapefiles
watershed_prj <- st_transform(watershed, crs=32136)#Chaiging into NAD 83/ TN meter usinf epsg code
watershed_prj## Simple feature collection with 1 feature and 2 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 690444.6 ymin: 119543.9 xmax: 763300.1 ymax: 186294.2
## projected CRS: NAD83 / Tennessee
## Id gridcode geometry
## 1 1315 202 POLYGON ((693552.5 121042.2...strm_network_prj <- st_transform(Stream_Network, st_crs(watershed_prj))
strm_network_prj## Simple feature collection with 687 features and 4 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 690688.8 ymin: 74477.87 xmax: 781811.4 ymax: 186634.7
## projected CRS: NAD83 / Tennessee
## First 10 features:
## arcid grid_code from_node to_node geometry
## 1 1 1 8 3 LINESTRING (709919.4 185368...
## 2 2 3 9 1 LINESTRING (706238.2 185308...
## 3 3 1 10 9 LINESTRING (706394.6 184801...
## 4 4 1 12 11 LINESTRING (762256.7 185513...
## 5 5 3 11 6 LINESTRING (763318.1 185687...
## 6 6 1 13 7 LINESTRING (774933.2 185493...
## 7 7 1 19 2 LINESTRING (708972 182983.8...
## 8 8 1 20 21 LINESTRING (779924.9 183960...
## 9 9 1 22 23 LINESTRING (719535.4 182750...
## 10 10 2 23 4 LINESTRING (720177.8 182497...plot(watershed_prj)
#Clip geometry using st_intersection
clip_basin <- st_intersection(strm_network_prj, watershed_prj)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesclip_basin## Simple feature collection with 221 features and 6 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 690720.6 ymin: 122535.5 xmax: 760842.6 ymax: 184474.6
## projected CRS: NAD83 / Tennessee
## First 10 features:
## arcid grid_code from_node to_node Id gridcode
## 17 17 1 26 32 1315 202
## 18 18 1 14 32 1315 202
## 21 21 1 15 35 1315 202
## 23 23 1 29 37 1315 202
## 24 24 2 32 38 1315 202
## 26 26 1 39 35 1315 202
## 29 29 2 35 46 1315 202
## 31 31 1 45 48 1315 202
## 32 32 2 48 38 1315 202
## 36 36 3 38 37 1315 202
## geometry
## 17 LINESTRING (728177.8 182008...
## 18 LINESTRING (730386 184337.6...
## 21 LINESTRING (751502.8 184474...
## 23 LINESTRING (734692 181242.9...
## 24 LINESTRING (728520.3 180858...
## 26 LINESTRING (750952.5 180026...
## 29 LINESTRING (746770.3 180556...
## 31 LINESTRING (722348.2 178659...
## 32 LINESTRING (723145.1 178130...
## 36 LINESTRING (728535.1 179841...plot(clip_basin)
#Plotting the stream orders within the watershed boundary
clip_basin <- clip_basin %>% mutate(Stream_order=clip_basin$grid_code)
colnames(clip_basin)## [1] "arcid" "grid_code" "from_node" "to_node" "Id"
## [6] "gridcode" "geometry" "Stream_order"StreamColors <- c("purple","green","grey","red", "blue")[unique((clip_basin$Stream_order))]
StreamColors## [1] "purple" "green" "grey" "red" "blue"plot(st_geometry(watershed_prj), main="Stream Order map of the Cumberland Plateau")
plot(clip_basin["Stream_order"], add=TRUE,pal=StreamColors)
legend("bottomright", legend=levels(as.factor(clip_basin$Stream_order)), fill=(StreamColors))
#Stream_order distribution in Cumberland
stream_order <- factor((clip_basin$Stream_order))#converted stream order into factor
par(mar=c(5, 4, 4, 2))
plot(stream_order, col="blue", xlab="Stream order", ylab="No. of streams")
#Calculating the drainage density
Area_basin <- watershed_prj %>%
mutate(basin_area=st_area(geometry))
Total_area <- sum(Area_basin$basin_area)
Total_area## 2918942453 [m^2]class(Total_area)## [1] "units"Total_area <- as.numeric(Total_area/1000000)
Total_area## [1] 2918.942Tot_len_stream <- clip_basin %>%
mutate(Streams_length=st_length(geometry))
class(Tot_len_stream$Streams_length)## [1] "units"Total_stream <- as.numeric(Tot_len_stream$Streams_length)
Total_stream_length <- (sum(Total_stream)/1000)
Total_stream_length## [1] 860.367class(Total_stream_length)## [1] "numeric"Drainage_density <- Total_stream_length/Total_area
Drainage_density## [1] 0.294753#Sampling area demarcation
Fifth_order <- clip_basin %>% filter(clip_basin$Stream_order==5)
Fifth_order## Simple feature collection with 16 features and 7 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 690720.6 ymin: 123174.7 xmax: 716286.2 ymax: 154489.5
## projected CRS: NAD83 / Tennessee
## First 10 features:
## arcid grid_code from_node to_node Id gridcode
## 1 207 5 214 224 1315 202
## 2 217 5 224 234 1315 202
## 3 263 5 234 279 1315 202
## 4 270 5 279 287 1315 202
## 5 275 5 287 292 1315 202
## 6 286 5 292 302 1315 202
## 7 303 5 302 323 1315 202
## 8 310 5 321 330 1315 202
## 9 312 5 323 321 1315 202
## 10 328 5 330 342 1315 202
## geometry Stream_order
## 1 LINESTRING (716286.2 154489... 5
## 2 LINESTRING (715808.5 152572... 5
## 3 LINESTRING (713851.3 150327... 5
## 4 LINESTRING (710157.3 144055... 5
## 5 LINESTRING (710447.5 142980... 5
## 6 LINESTRING (710669.5 141318... 5
## 7 LINESTRING (708448 139780.9... 5
## 8 LINESTRING (702601.2 135951... 5
## 9 LINESTRING (709174.7 135968... 5
## 10 LINESTRING (702104.6 135329... 5plot(Fifth_order["Stream_order"])
Sampling_site <- st_buffer(Fifth_order, 2000)
Sampling_site## Simple feature collection with 16 features and 7 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 688720.7 ymin: 121174.8 xmax: 718286.1 ymax: 156489.4
## projected CRS: NAD83 / Tennessee
## First 10 features:
## arcid grid_code from_node to_node Id gridcode
## 1 207 5 214 224 1315 202
## 2 217 5 224 234 1315 202
## 3 263 5 234 279 1315 202
## 4 270 5 279 287 1315 202
## 5 275 5 287 292 1315 202
## 6 286 5 292 302 1315 202
## 7 303 5 302 323 1315 202
## 8 310 5 321 330 1315 202
## 9 312 5 323 321 1315 202
## 10 328 5 330 342 1315 202
## geometry Stream_order
## 1 POLYGON ((713949.6 153308.3... 5
## 2 POLYGON ((712843 152054.7, ... 5
## 3 POLYGON ((708967.3 147011.7... 5
## 4 POLYGON ((708446.8 143019.1... 5
## 5 POLYGON ((708691.5 142023.7... 5
## 6 POLYGON ((708438.4 142090.4... 5
## 7 POLYGON ((705887.9 139195, ... 5
## 8 POLYGON ((703667.7 134081.4... 5
## 9 POLYGON ((705555.2 133360, ... 5
## 10 POLYGON ((699715.2 134467.1... 5plot(st_geometry(Sampling_site), main="Proposed rock sampling Sites in the drainage basin")
plot(st_geometry(watershed_prj), add=TRUE)
#reading dem in r
dem <- raster("D:/Final_Seminar/Cumberland.tif")
crs(dem)## CRS arguments: +proj=longlat +datum=WGS84 +no_defspar(mar=c(5, 4, 4, 2))
plot(dem)
#changing crs of raster
dem_prj <- projectRaster(dem, crs = "+init=epsg:32136")
crs(dem_prj)## CRS arguments:
## +proj=lcc +lat_0=34.3333333333333 +lon_0=-86 +lat_1=36.4166666666667
## +lat_2=35.25 +x_0=600000 +y_0=0 +datum=NAD83 +units=m +no_defs#Minimum and Maximum elevation of the area
cellStats(dem_prj, min)## [1] 175.5502cellStats(dem_prj, max)## [1] 1661.282#masking and cropping the dem
dem_mask <- mask(dem_prj, watershed_prj)
dem_mask## class : RasterLayer
## dimensions : 3658, 3680, 13461440 (nrow, ncol, ncell)
## resolution : 25.2, 30.8 (x, y)
## extent : 689994, 782730, 74261.11, 186927.5 (xmin, xmax, ymin, ymax)
## crs : +proj=lcc +lat_0=34.3333333333333 +lon_0=-86 +lat_1=36.4166666666667 +lat_2=35.25 +x_0=600000 +y_0=0 +datum=NAD83 +units=m +no_defs
## source : memory
## names : Cumberland
## values : 175.5502, 921.7294 (min, max)plot(dem_mask)
dem_cropped <- crop(dem_mask, watershed_prj)
par(mar=c(5, 4, 4, 2))
plot(dem_cropped)
#Plotting the digital elevation model
brk1 <- c(200, 300, 400, 500, 600, 700, 800, 900)
par(mar=c(5, 4, 4, 2))
plot(dem_cropped, main="Elevation model of the watershed", breaks=brk1, col=brewer.pal(n = 8, name = "RdBu"))
plot(st_geometry(watershed_prj), add=TRUE)
plot(clip_basin, add=TRUE)## Warning in plot.sf(clip_basin, add = TRUE): ignoring all but the first attribute
par(mar=c(5, 4, 4, 2))
hist(dem_cropped, main="Hisotgram of the elevation of the watershed", col="blue", xlab="Elevation (meters)", ylab="Frequency(pixels)")
#Plotting the slope model of the watershed
slope <- terrain(dem_cropped, opt="slope", unit="degrees", neighbors=8, filename="slope", overwrite=TRUE)
slope## class : RasterLayer
## dimensions : 2167, 2891, 6264797 (nrow, ncol, ncell)
## resolution : 25.2, 30.8 (x, y)
## extent : 690447.6, 763300.8, 119537.1, 186280.7 (xmin, xmax, ymin, ymax)
## crs : +proj=lcc +lat_0=34.3333333333333 +lon_0=-86 +lat_1=36.4166666666667 +lat_2=35.25 +x_0=600000 +y_0=0 +datum=NAD83 +units=m +no_defs
## source : D:/Final_Seminar/slope.grd
## names : slope
## values : 0, 56.72872 (min, max)brk2 <- c(10, 20, 30, 40, 50, 60)
par(mar=c(5, 4, 4, 2))
plot(slope, main="Slope map of the watershed", breaks=brk2, col=brewer.pal(n = 6, name = "Set1"))
plot(st_geometry(watershed_prj), add=TRUE)
plot(clip_basin, add=TRUE)## Warning in plot.sf(clip_basin, add = TRUE): ignoring all but the first attribute
par(mar=c(5, 4, 4, 2))
hist(slope, main="Histogram of the slope of the watershed", col="blue", xlab="Slope (degrees)", ylab="Frequency(pixels)")## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2%
## of the raster cells were used. 100000 values used.
#Plotting aspect map of the area
brk3 <- c(45, 90, 135, 180, 225, 270, 315, 360)
aspect <- terrain(dem_cropped, opt="aspect", unit="degrees", neighbors=8, filename="aspect", overwrite=TRUE)
par(mar=c(5, 4, 4, 2))
plot(aspect,main="Aspect map of the watershed", breaks=brk3, col=brewer.pal(n = 8, name = "Set1"))
plot(st_geometry(watershed_prj), add=TRUE)
plot(clip_basin, add=TRUE)## Warning in plot.sf(clip_basin, add = TRUE): ignoring all but the first attribute
par(mar=c(5, 4, 4, 2))
hist(aspect, main= "Histogram of the aspect of the watershed", col= "blue", xlab="Aspect (degrees)", ylab="Frequency(pixels)")## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2%
## of the raster cells were used. 100000 values used.