Basic Use
The goal of hydrofabric3D is to generate DEM-based cross sections for hydrographic networks.
Installation
You can install the development version of hydrofabric3D from GitHub with:
Example
This is a basic example which shows you how to cut cross sections for a network.
Define Network
library(hydrofabric3D)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
(net = linestring %>%
mutate(bf_width = exp(0.700 + 0.365* log(totdasqkm))))
#> Simple feature collection with 325 features and 5 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 77487.09 ymin: 890726.5 xmax: 130307.4 ymax: 939129.8
#> Projected CRS: NAD83 / Conus Albers
#> # A tibble: 325 × 6
#> nhdplus_comid geometry comid totdasqkm dist_m bf_width
#> * <chr> <LINESTRING [m]> <dbl> <dbl> <dbl> <dbl>
#> 1 101 (128525.6 892408.3, 128565.7 … 1.01e2 7254. 3.25e3 51.7
#> 2 24599575 (128084.7 892952.4, 128525.6 … 2.46e7 7249. 7.00e2 51.6
#> 3 1078635 (127687.6 893270.4, 127799.7 … 1.08e6 7248. 5.22e2 51.6
#> 4 1078637 (124942.8 893959.6, 124948.2 … 1.08e6 68.2 4.17e3 9.41
#> 5 1078639 (125523.1 892528, 125657.3 89… 1.08e6 7180. 2.76e3 51.5
#> 6 1078577 (123219.9 902292.8, 123233.5 … 1.08e6 19.8 9.91e3 5.99
#> 7 1078575 (121975.5 909050.8, 122028.9 … 1.08e6 41.3 1.87e4 7.83
#> 8 1078657 (124263.8 892410.4, 124420.6 … 1.08e6 7179. 1.66e3 51.5
#> 9 1078663 (125628.9 892216, 125555.7 89… 1.08e6 0.099 7.54e2 0.866
#> 10 1078643 (124248.1 892440.7, 124263.8 … 1.08e6 7178. 3.41e1 51.5
#> # ℹ 315 more rows
plot(net$geometry)
Cut cross sections
(transects = cut_cross_sections(net = net,
crosswalk_id = "comid",
cs_widths = pmax(50, net$bf_width * 7),
num = 10) )
#> Densifying
#> Smoothing
#> Cutting transects
#> Formatting transects
#> Simple feature collection with 2372 features and 6 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 77486.95 ymin: 890563.3 xmax: 130368.6 ymax: 939130
#> Projected CRS: NAD83 / Conus Albers
#> # A tibble: 2,372 × 7
#> comid cs_id cs_lengthm cs_measure ds_distance sinuosity
#> <dbl> <int> <dbl> <dbl> <dbl> <dbl>
#> 1 101 1 362. 0.393 12.8 1.05
#> 2 101 2 362. 11.3 369. 1.00
#> 3 101 3 362. 25.2 822. 1.01
#> 4 101 4 362. 38.1 1242. 1.05
#> 5 101 5 362. 50.1 1632. 1.08
#> 6 101 6 362. 60.2 1961. 1.19
#> 7 101 7 362. 68.7 2238. 1.34
#> 8 101 8 362. 76.1 2480. 1.34
#> 9 101 9 362. 86.1 2805. 1.15
#> 10 101 10 362. 99.7 3250. 1.00
#> # ℹ 2,362 more rows
#> # ℹ 1 more variable: geometry <LINESTRING [m]>
plot(transects$geometry)
Define Cross section points
(pts = cross_section_pts(transects,
crosswalk_id = "comid",
dem = "/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/USGS_Seamless_DEM_1.vrt"))
#> Simple feature collection with 24386 features and 7 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 77487.45 ymin: 890577.2 xmax: 130362.9 ymax: 939128.2
#> Projected CRS: NAD83 / Conus Albers
#> # A tibble: 24,386 × 8
#> comid cs_id pt_id cs_lengthm relative_distance Z points_per_cs
#> <dbl> <int> <int> <dbl> <dbl> <dbl> <dbl>
#> 1 101 1 1 362. 0 41.9 12
#> 2 101 1 2 362. 32.9 41.1 12
#> 3 101 1 3 362. 65.7 42.1 12
#> 4 101 1 4 362. 98.6 40.6 12
#> 5 101 1 5 362. 131. 38.8 12
#> 6 101 1 6 362. 164. 36.2 12
#> 7 101 1 7 362. 197. 37.9 12
#> 8 101 1 8 362. 230. 41.6 12
#> 9 101 1 9 362. 263. 44.2 12
#> 10 101 1 10 362. 296. 44.3 12
#> # ℹ 24,376 more rows
#> # ℹ 1 more variable: geometry <POINT [m]>Classify Cross section points
(classified_pts = classify_points(pts, crosswalk_id = "comid"))
#> Simple feature collection with 24386 features and 14 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 77487.45 ymin: 890577.2 xmax: 130362.9 ymax: 939128.2
#> Projected CRS: NAD83 / Conus Albers
#> # A tibble: 24,386 × 15
#> comid cs_id pt_id cs_lengthm relative_distance Z points_per_cs class
#> <dbl> <int> <int> <dbl> <dbl> <dbl> <int> <chr>
#> 1 101 1 1 362. 0 41.9 12 left_bank
#> 2 101 1 2 362. 32.9 41.7 12 channel
#> 3 101 1 3 362. 65.7 41.2 12 channel
#> 4 101 1 4 362. 98.6 40.5 12 channel
#> 5 101 1 5 362. 131. 38.5 12 channel
#> 6 101 1 6 362. 164. 36.2 12 bottom
#> 7 101 1 7 362. 197. 38.6 12 channel
#> 8 101 1 8 362. 230. 41.2 12 channel
#> 9 101 1 9 362. 263. 43.4 12 channel
#> 10 101 1 10 362. 296. 44.1 12 right_bank
#> # ℹ 24,376 more rows
#> # ℹ 7 more variables: point_type <chr>, bottom <dbl>, left_bank <dbl>,
#> # right_bank <dbl>, valid_banks <lgl>, has_relief <lgl>, geometry <POINT [m]>Explore!
library(ggplot2)
ggplot(data = filter(classified_pts, comid == 101) ) +
geom_line(aes(x = relative_distance, y = Z)) +
geom_point(aes(x = relative_distance, y = Z, color = class)) +
facet_wrap(~cs_id, scales = "free") +
theme_minimal() +
theme(legend.position = "bottom")
Time to get 2372 transects and 24386 classified points …
system.time({
cs = net %>%
cut_cross_sections(crosswalk_id = "comid",
cs_widths = pmax(50, net$bf_width * 7),
num = 10) %>%
cross_section_pts(
crosswalk_id = "comid",
dem = '/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/USGS_Seamless_DEM_1.vrt') %>%
classify_points(crosswalk_id = "comid")
})
#> Densifying
#> Smoothing
#> Cutting transects
#> Formatting transects
#> user system elapsed
#> 18.391 0.076 20.199