summarize-layer-by-huc12
Source:vignettes/summarize-layer-by-huc12.Rmd
summarize-layer-by-huc12.Rmd
library(middlesnake)
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
library(sf)
#> Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUEGet watersheds by Conservation District
automatically clips watersheds to district boundary
columbia_huc12 <- get_huc_by_cd("Columbia")
#> ℹ Downloading huc12 boundaries from Watershed Boundary Dataset...
#> ℹ Assigning HUCs to Columbia...
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> ✔ Downloaded and processed 50 HUC polygons.Show geoserver layers
show_geoserver_layers()
#> Querying GeoServer WMS capabilities...
#> ✔ Retrieved 18 WMS layers.
#> [1] "all_columbia_huc_12s"
#> [2] "columbia-critical-aquifer"
#> [3] "columbia-critical-conservation-areas"
#> [4] "columbia-frequently-flooded"
#> [5] "columbia-geologic-hazard"
#> [6] "columbia-wetlands"
#> [7] "columbia_change_detect"
#> [8] "columbia_county_boundary"
#> [9] "columbia_fema"
#> [10] "dnr_forest_priority"
#> [11] "full_lidar_hillshade"
#> [12] "salmon-recovery-portal"
#> [13] "salmon-recovery-projects"
#> [14] "soil_sample_locations_orig"
#> [15] "usda_2011_colored2"
#> [16] "usda_2024_colored2"
#> [17] "washington-wria"
#> [18] "wildfire_perimiters_2023"Get the layer that need summarized
Here we use frequently flooded areas
flooded <- get_geoserver_layer("columbia-frequently-flooded")
#> ℹ Attempting to download vector layer "columbia-frequently-flooded" via WFS
#> ✔ Layer downloaded: 131 features.Summarize
How many frequently flooded acres are inside each huc 12 watershed in Columbia County?
columbia_flooded_huc12 <- summarize_acres_by_huc12(huc12 = columbia_huc12,
overlay = flooded)
#> ℹ Validating geometries...
#> ℹ Performing spatial intersection...
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> ℹ Summarizing point or line data by feature count...
#> ✔ Summary complete.returns huc12 (numeric code) and acres. If
you want the watershed name add it
columbia_flooded_huc12 %>%
dplyr::left_join(columbia_huc12 %>%
select(name,huc12)
)
#> Joining with `by = join_by(huc12)`
#> # A tibble: 35 × 4
#> huc12 count name geometry
#> <chr> <int> <chr> <POLYGON [°]>
#> 1 170601070704 14 Tucannon River ((-118.0833 46.53686, -11…
#> 2 170701020307 14 Patit Creek ((-117.8352 46.35517, -11…
#> 3 170701020401 14 Payne Hollow-Touchet River ((-117.9852 46.31925, -11…
#> 4 170701020308 11 Lower North Fork Touchet River ((-117.8228 46.2928, -117…
#> 5 170601070703 10 Town of Starbuck-Tucannon River ((-117.9791 46.55071, -11…
#> 6 170701020304 8 South Fork Touchet River ((-117.9051 46.27083, -11…
#> 7 170701020402 8 Whiskey Creek ((-117.9633 46.21827, -11…
#> 8 170601070701 7 Smith Hollow-Tucannon River ((-117.9958 46.50787, -11…
#> 9 170701020501 7 Upper Whetstone Hollow ((-117.7978 46.37677, -11…
#> 10 170601070608 6 Town of Marengo-Tucannon River ((-117.7862 46.46699, -11…
#> # ℹ 25 more rowsLet’s try wetlands
wetlands <- get_geoserver_layer("columbia-wetlands")
#> ℹ Attempting to download vector layer "columbia-wetlands" via WFS
#> ✔ Layer downloaded: 618 features.This time we want to stratify by wetland type within each huc 12
We’ll set the group_vars argument
columbia_wetlands_huc12 <- summarize_acres_by_huc12(huc12 = columbia_huc12,overlay = wetlands,group_vars = "WETLAND_TY")
#> ℹ Validating geometries...
#> ℹ Performing spatial intersection...
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> ℹ Summarizing point or line data by feature count...
#> ✔ Summary complete.
columbia_wetlands_huc12
#> # A tibble: 90 × 3
#> huc12 WETLAND_TY count
#> <chr> <chr> <int>
#> 1 170601070608 Freshwater Forested/Shrub Wetland 44
#> 2 170701020304 Riverine 40
#> 3 170601070608 Freshwater Emergent Wetland 34
#> 4 170701020401 Riverine 24
#> 5 170701020303 Freshwater Pond 23
#> 6 170701020303 Riverine 19
#> 7 170701020304 Freshwater Pond 19
#> 8 170701020302 Freshwater Pond 18
#> 9 170701020401 Freshwater Forested/Shrub Wetland 18
#> 10 170601070704 Freshwater Forested/Shrub Wetland 17
#> # ℹ 80 more rowsUse a layer not on Geoserver
We already have Columbia watersheds to use as aoi (i.e. Area of Interest) to automatically subset the layer we are grabbing. In this example is using WA DNR Priority areas.
dnr_forest_priority <- get_gis_from_url(url = "https://services.arcgis.com/4x406oNViizbGo13/arcgis/rest/services/Eastern_WA_Forest_Health_Priority_HUC6_Watersheds/FeatureServer/0/query?where=1%3D1&outFields=*&outSR=4326&f=geojson",aoi = columbia_huc12)
#> ℹ Requesting GIS data from <https://services.arcgis.com/4x406oNViizbGo13/arcgis/rest/services/Eastern_WA_Forest_Health_Priority_HUC6_Watersheds/FeatureServer/0/query?where=1%3D1&outFields=*&outSR=4326&f=geojson&geometry=-118.2422963,45.9990865,-117.6032963,46.6246534&geometryType=esriGeometryEnvelope&inSR=4326&spatialRel=esriSpatialRelIntersects>
#> ✔ Parsed 38 features from GeoJSON.
columbia_forest_priority_huc12 <- summarize_acres_by_huc12(huc12 = columbia_huc12,
overlay = dnr_forest_priority)
#> ℹ Validating geometries...
#> ℹ Performing spatial intersection...
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> ℹ Reprojecting to EPSG:5070 for accurate area calculation...
#> ℹ Calculating area in acres...
#> ✔ Summary complete.