Skip to content
Snippets Groups Projects
Commit 04b68d44 authored by Bock, Andy's avatar Bock, Andy
Browse files

Stylistic Changes

parent b311c2eb
No related branches found
No related tags found
1 merge request!36Style changes/ non-contributing area work
Expand all Hide whitespace changes
Inline Side-by-side
workspace/CoastalSegs.Rmd deleted 100644 → 0
+ 0
− 133
View file @ b311c2eb
---
title: "Coastal_Segs10"
output: html_document
editor_options:
chunk_output_type: console
---
title: CoastalSegs10.Rmd
Project: GFv2.0
Date: 7/2020
Author: [Marilyn Santiago](msant@usgs.gov)
Script purpose:: To extract nhd v2 flowlines with terminal at Ftype coastal
and which terminal DA is <=10,
will check if COMIDs already in the agg_fline,
Will add huc12 and check if GagesIII at those lines are not in the POIs dataset
Generates 0 order catchments
Need to tackle large DA ( > 10 sqkm terminal coastal outlets not on a POI)
```{r}
knitr::opts_chunk$set(error = TRUE)
```
```{r setup}
library(nhdplusTools)
library(dplyr)
library(sf)
library(rgdal)
library(tidyverse)
# Load user-defined data paths
data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
# Load NHDPlus Data if precprocessed to RDS format
nhdplus_path(data_paths$nhdplus_gdb)
staged_nhd <- stage_national_data(nhdplus_data = data_paths$nhdplus_gdb,
output_path = data_paths$nhdplus_dir)
# Load custom functions
source("R/utils.R")
source("R/NHD_navigate.R")
# NHDPlus V2 Region
RPU <-"01a"
VPU <-substr(RPU, 1, 2)
# Source Data
huc12_xWalk <- file.path(data_paths$nhdplus_dir, "CrosswalkTable_NHDplus_HU12.csv")
# Inputs
nav_gpkg <- file.path("cache", paste0("GF_", VPU, ".gpkg"))
nhdflowline <- "NHDFlowline"
coastFlowlines <- "coast_NCA"
pois_all <- paste0("POIs_", hydReg)
nseg_raw <- paste0("nsegment_raw_", VPU)
poi_res_coast <- paste0("gage_rescoast_", hydReg)
# Output Geopackage where we are writing coast output layers
out_gpkg <- file.path("cache", paste0(RPU,"_coast.gpkg"))
# Hy_Refactor output gpkg
href_gpkg <- file.path("cache", paste0(RPU, ".gpkg"))
# Output layers
segments_coast_lt10 <- "coast_lt10" #flowlines that terminate at the coast w DA<=10
segments_coast_all <- "coast_all"
dup_line <-"duplicate_fline"
nhdCoast_net <- paste0("nhdCoast_network_", hydReg)
nhdCoast_ND <- paste0("nhdCoast_ND_", hydReg)
cat0 <- "catchment_coast"
```
```{r Coastal flowlines}
if(needs_layer(out_gpkg, nhdCoast_net)) {
# nhd features by the RPU
nhd <- read_sf(href_gpkg, "nhd_flowline")
# NHD flowlines which are coastlines
nhd_Coastline <- read_rds(staged_nhd$flowline) %>% filter(FTYPE == "Coastline")
# Coastal < 10 sqkm network generated in NHD_Navigate
coastNHD <- read_sf(nav_gpkg, coastFlowlines)
#########*******move to NHD_NAV
# Non-dendritic NHD flowlines
non_dendritic <- nhd %>% filter(StreamOrde != StreamCalc)
# Retrieve terminal paths
terminal_fl <- nhd %>% filter(TerminalFl == 1) %>% inner_join(st_drop_geometry(nhd_Coastline), by = c("DnHydroseq" = "Hydroseq"))
nhdCoast_COMs <- unique(unlist(lapply(unique(terminal_fl$COMID.x), NetworkNav, st_drop_geometry(nhd))))
coastNHD <- nhd %>% bind_rows(nhd %>% filter(COMID %in% nhdCoast_COMs))
#********************
# crosswalk to HU12, filter b VPu
nhd_to_HU12 <- read.csv(huc12_xWalk, colClasses = c("character", "integer", "character")) %>% filter(grepl(paste0("^", VPU), .data$HUC_12))
nhd_to_HU12_seg <- nhd_to_HU12 %>% filter(FEATUREID %in% coastNHD$COMID)
# Get the coastal flowlines
nhdCoast_network <- nhd %>% filter(COMID %in% coastNHD$COMID) %>% left_join(nhd_to_HU12, by = c("COMID" = "FEATUREID"))
# Write out geopackage layer representing coastal segments
# Note this is NOT correct and includes terminal non-coastal segments
write_sf(nhdCoast_network, out_gpkg, nhdCoast_net)
} else {
# Load HUC12 POIs as the tmpPOIs if they already exist
nhdCoast_network <- read_sf(out_gpkg, nhdCoast_net)
nhd_Coastline <- read_rds(staged_nhd$flowline) %>% filter(FTYPE == "Coastline")
nhd_to_HU12 <- read.csv(huc12_xWalk) %>% filter(grepl(paste0("^", VPU), .data$HUC_12))
}
```
```{r Zero-order Coastal Catchments}
if(needs_layer(out_gpkg, cat0)) {
# catchment features
cat <- read_rds(staged_nhd$catchment)
# set same as the agg_fline crs here
cat0 <- cat %>% dplyr::filter(FEATUREID %in% nhd_Coastline$COMID) %>% st_transform(5070)
# get RPUID from nhd attributes
nhdAtt <- readRDS(staged_nhd$attributes) %>% select(COMID, RPUID)
# crosswalk to HU12
nhd_to_HU12_segCats <- nhd_to_HU12 %>% dplyr::select(FEATUREID, HUC_12) %>% dplyr::filter(FEATUREID %in% cat0$FEATUREID)
#nhd_to_HU12_Cats <- nhd_to_HU12 %>% filter(FEATUREID %in% nhd_to_HU12_segCats$FEATUREID) #HUC_12 %in% nhd_to_HU12_segCats$HUC_12 &
cat0 <- cat0 %>% inner_join(nhd_to_HU12, by = c("FEATUREID")) %>% inner_join(nhdAtt, by = c("FEATUREID" = "COMID"))
# Do something here like Group/aggregate to HUC12
write_sf(cat0, out_gpkg, cat0_lyr)
} else {
print("Output gpkg already exists, need to delete to reprocess:")
print(paste0(out_gpkg))
}
```
workspace/NHD_navigate.Rmd
+ 205
− 130
View file @ 04b68d44
This diff is collapsed.
workspace/Non_contributing.Rmd deleted 100644 → 0
+ 0
− 203
View file @ b311c2eb
title: non_contributing_cat_v2.Rmd
Project: GFv2.0
Date: 7/2020
Author: [Marilyn Santiago](msant@usgs.gov)
Script purpose: To extract Sink Catchments from V2, assign a HUC12
by the crosswalk table (cwHUC_12), then dissolve by the cw_HUC12 then assign HUCs by location (HUC_12) to the dissolved sinks, then classify. classify dissolved sinks based on area ratio of the dissolved sinkarea/cw_HUC12 area
small (0,0.95) full [0.95,1.3] large > 1.3
Process all RPUs
Required:
catchments and the crosswalk (table) to HU12
RPU data from
NHDPlus attributes
```{r}
knitr::opts_chunk$set(error = TRUE)
```
```{r setup}
library(nhdplusTools)
library(dplyr)
library(sf)
library(rgdal)
library(tidyverse)
# Load user-defined data paths
data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
# Load NHDPlus Data if precprocessed to RDS format
nhdplus_path(data_paths$nhdplus_gdb)
staged_nhd <- stage_national_data(nhdplus_data = data_paths$nhdplus_gdb,
output_path = data_paths$nhdplus_dir)
hydReg <- "13"
# Load custom functions
source("R/utils.R")
source("R/NHD_navigate.R")
# Input layers
themPOIs <- st_read("cache/CONUS.gpkg", "POIs")
segs <- st_read("cache/CONUS.gpkg", "nsegment_raw")
HUC12_out <- read_sf(data_paths$hu12_points_path, "hu_outlets")
nhd <- readRDS(staged_nhd$flowline)
POIs_att <- nhd %>% filter(COMID %in% themPOIs$COMID)
nhd_to_HU12 <- read.csv(file.path(data_paths$nhdplus_dir, "CrosswalkTable_NHDplus_HU12.csv"),
colClasses = c("character", "integer", "character"), header = T)
HUC12_lyr <- read_rds(file.path(data_paths$nhdplus_dir, "HUC12.rds")) %>%
select(HUC_10, HUC_12, AreaHUC12, HU_12_TYPE, HU_10_TYPE, HU_12_DS) %>% st_transform(5070) %>%
mutate(HUC_10 = as.character(HUC_10), HUC_12 = as.character(HUC_12),
HU_10_DS = substr(HU_12_DS, 1, 10))
sink_lyr <- readOGR(dsn = data_paths$nhdplus_gdb, layer = "Sink") %>% st_as_sf() %>% select(SINKID, FEATUREID, SOURCEFC, InRPU) %>% rename(RPUID = InRPU, SOURCEFC_Sink = SOURCEFC) %>%
mutate(RPUID = as.character(RPUID)) %>% mutate(FEATUREID = ifelse(FEATUREID == -9999, SINKID, FEATUREID))
# Bring in catchments and identify sinks, set crs same as the agg_fline, set precision
nhd_cats <- read_rds(staged_nhd$catchment) %>% st_transform(5070) %>% st_set_precision(10) %>% st_make_valid()
# Output Geopackage where we are writing coast output layers
out_gpkg <- file.path("cache", paste0("NC_R13.gpkg"))
# Output layers
sink_Net <- "Terminal_Segs"
sink_cat <- "Terminal_Cats"
sinkFinal <- "sink_Final"
```
```{r filter by region for now}
# subset layers initially for hydroregion
nhd <- nhd %>% filter(VPUID == hydReg)
HUC12_out <- HUC12_out %>% filter(grepl(paste0("^", hydReg,".*"), .data$HUC12)) %>% st_drop_geometry()
POIs_att <- nhd %>% filter(COMID %in% themPOIs$COMID)
nhd_to_HU12 <- nhd_to_HU12 %>% filter(grepl(paste0("^", hydReg,".*"),.data$HUC_12))
HUC12_lyr <- HUC12_lyr %>% filter(grepl(paste0("^", hydReg,".*"), .data$HUC_12)) %>%
st_cast("MULTIPOLYGON") %>% st_cast("POLYGON") %>% st_make_valid()
sink_lyr <- sink_lyr %>% filter(grepl(paste0("^", hydReg,".*"), .data$RPUID))
reg_cats <- nhd_cats %>% inner_join(nhd_to_HU12 %>% select(FEATUREID, HUC_12), by = "FEATUREID")
st_write(reg_cats, out_gpkg, "reg_cats", delete_layer = TRUE)
```
```{r Initial NCA screening, concentrate on closed HUC12s not in thematic POIs, and 1:1 correspondence with NHD catchments}
# Closed catchments from HUC12 dataset
HUC10_closed <- HUC12_lyr %>% filter(HU_10_TYPE == "C")
HUC12_closed <- HUC12_lyr %>%
filter(HU_12_TYPE == "C", !HUC_12 %in% HUC10_closed$HUC_12) %>% bind_rows(HUC10_closed)
# ALL HUC12s that are closed and upstream contributions
while(nrow(st_drop_geometry(HUC12_lyr) %>%
filter(!HUC_12 %in% HUC12_closed$HUC_12, HU_12_DS %in% HUC12_closed$HUC_12)) > 0){
HUC12_sub <- HUC12_lyr %>% filter(HU_12_DS %in% HUC12_closed$HUC_12) %>% st_drop_geometry()
HUC12_closed <- HUC12_closed %>%
bind_rows(HUC12_lyr %>% filter(!HUC_12 %in% HUC12_closed$HUC_12, HU_12_DS %in% HUC12_closed$HUC_12))
print(nrow(st_drop_geometry(HUC12_lyr) %>%
filter(!HUC_12 %in% HUC12_closed$HUC_12, HU_12_DS %in% HUC12_closed$HUC_12)))
}
# HUC12s not represented in thematic POI dataset (may remove later)
NCA_WBD <- HUC12_closed %>% filter(!HUC_12 %in% themPOIs$Type_HUC12)
st_write(NCA_WBD, out_gpkg, "NCA_WBD", delete_layer = T)
# NHDPlus catchments within these HUC12s as identified by the crosswalk
st_write(reg_cats %>% filter(HUC_12 %in% NCA_WBD$HUC_12), out_gpkg, "nhd_HUC12", delete_layer = T)
# Get closed basins from WBD and retrieve associated NHD+ catches
ClosedBasins <- reg_cats %>% filter(HUC_12 %in% HUC12_closed$HUC_12) %>%
st_as_sf() %>% group_by(HUC_12) %>% summarise() %>%
st_cast("MULTIPOLYGON") %>% st_cast("POLYGON") %>% st_make_valid() %>%
mutate(nhdCatch_sqKM = as.numeric(st_area(.)) * 0.000001)
# Intersect the two without constraining for HUC12
sink_byHUC_sum_int <- sf::st_intersection(select(ClosedBasins, HUC_12),
select(NCA_WBD, HUC_12)) %>%
st_cast("MULTIPOLYGON") %>% st_cast("POLYGON") %>% st_make_valid() %>%
dplyr::mutate(intArea = as.numeric(st_area(.)) * 0.000001)
# Calculate the Coefficent of areal correspondence (amount of overlap)
Final_NC_CAC <- sink_byHUC_sum_int %>% filter(HUC_12 == HUC_12.1) %>%
inner_join(ClosedBasins %>% st_drop_geometry(), by = "HUC_12") %>%
inner_join(HUC12_lyr %>% st_drop_geometry(), by = "HUC_12") %>%
mutate(intArea = if_else(is.na(intArea), 0, intArea)) %>%
mutate(CAC = (intArea / ((nhdCatch_sqKM - intArea) + (AreaHUC12 - intArea) + intArea)))
# Not final spatial layers, just means the mapping is correct
Final_CAC_HUC12 <- Final_NC_CAC %>% group_by(HUC_12) %>% summarize(sumCAC = sum(CAC), n = n()) %>%
st_cast("MULTIPOLYGON") %>% st_cast("POLYGON") %>% st_make_valid() %>% filter(sumCAC > 0.90)
# Write out HUC12s for which meet the criteria
st_write(sink_byHUC_sum_int, out_gpkg, "HUC12_int", delete_layer = T)
st_write(Final_CAC_HUC12, out_gpkg, "CAC_HUC12_Fit", delete_layer = T)
##### End of step 1
```
```{r Find spatial correpsondence at scales larger than HUC12}
# Terminal segments, ensure not coastal terminal
sink_seg <- nhd %>% filter(TerminalFl == 1)
# cAtchments that are considered sinks
sink_cats <- nhd_cats %>% filter(SOURCEFC == "Sink" | FEATUREID %in% sink_lyr$FEATUREID) %>%
left_join(nhd_to_HU12 %>% select(FEATUREID, HUC_12), by = "FEATUREID") %>%
bind_rows(reg_cats %>% filter(FEATUREID %in% sink_seg$COMID))
# Dendritic flowlines that terminate at catchment sinks
if(needs_layer(out_gpkg, sink_Net)){
sink_onNet <- nhd %>%
filter(COMID %in% c(sink_cats$FEATUREID, sink_lyr$FEATUREID), !LevelPathI %in% POIs_att$LevelPathI)
upNet <- unique(unlist(lapply(unique(sink_onNet$COMID), NetworkNav, st_drop_geometry(nhd), withTrib = "TRUE")))
# Write out all
sinkNet <- nhd %>% filter(COMID %in% c(upNet, sink_onNet$COMID))
st_write(sinkNet, out_gpkg, sink_Net, delete_layer = TRUE)
} else {
sinkNet <- st_read(out_gpkg, sink_Net)
}
sink_cats <- sink_cats %>% bind_rows(reg_cats %>% filter(FEATUREID %in% sinkNet$COMID)) %>%
distinct(., .keep_all = T)
st_write(sink_cats, out_gpkg, sink_cat)
# Intersect the two without constraining for HUC12 area
sink_byNHDCatch <- sf::st_intersection(select(sink_cats, FEATUREID),
select(HUC12_lyr, HUC_12)) %>%
st_cast("MULTIPOLYGON") %>% st_cast("POLYGON") %>% st_make_valid() %>%
group_by(FEATUREID, HUC_12) %>% summarize() %>% st_as_sf() %>% ungroup() %>%
mutate(intArea = as.numeric(st_area(.)) * .000001) %>%
inner_join(st_drop_geometry(sink_cats) %>% select(FEATUREID, AreaSqKM), by = "FEATUREID") %>%
inner_join(st_drop_geometry(HUC12_lyr) %>% select(HUC_12, AreaHUC12), by = "HUC_12") %>%
st_cast("MULTIPOLYGON") %>% st_cast("POLYGON") %>% st_make_valid()
#st_write(sink_byNHDCatch, out_gpkg, "nhd_wbd_int", delete_layer = TRUE)
# Find NHD catch that are 100% within existing HUC12
HUC_12_addendum <- st_drop_geometry(sink_byNHDCatch) %>%
#group_by(FEATUREID) %>% mutate(propAreaNHD = intArea / nhdArea) %>%
#group_by(HUC_12) %>% mutate(propAreaWBD = intArea / wbdArea) %>%
group_by(FEATUREID, HUC_12) %>% summarize(sum = sum(intArea),
nhd = sum(AreaSqKM), wbd = mean(AreaHUC12), count=n()) %>%
mutate(CAC = (sum / ((nhd - sum) + (wbd - sum) + sum)), PercInNHD = sum/nhd, PercInWBD = sum/wbd)
# most of the NHD catchment within intersected area, can be assigned
HUC12_Fit <- HUC_12_addendum %>% filter(PercInNHD > 0.95)
wbd_in_nhd <- sink_cats %>% filter(FEATUREID %in% HUC12_Fit$FEATUREID) %>%
st_cast("MULTIPOLYGON") %>% st_cast("POLYGON") %>% st_make_valid()
st_write(wbd_in_nhd, out_gpkg, "wbd_in_nhd", delete_layer = T)
NHD_Fit <- HUC_12_addendum %>% filter(PercInWBD > 0.95) %>% group_by(FEATUREID) %>%
summarize(intArea = sum(sum), wbdArea = (sum(wbd)))
nhd_in_wbd <- sink_cats %>% filter(FEATUREID %in% NHD_Fit$FEATUREID) %>%
st_cast("MULTIPOLYGON") %>% st_cast("POLYGON") %>% st_make_valid()
st_write(nhd_in_wbd, out_gpkg, "nhd_in_wbd", delete_layer = T)
} else {
# load sink_cats_HUC12, nhd_to_HU12, sink_lyr, sink_Final
print("Output gpkg already exists, need to delete to reprocess:")
}
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment