search

bcgov / watershedBC

Summarize watershed attributes in BC
https://bcgov-env.shinyapps.io/watershedBC/
Apache License 2.0
6 stars 0 forks source link

Roll up FWA watersheds by stream order #18

Closed bevingtona closed 8 months ago

bevingtona commented 8 months ago

Roll up FWA watersheds by stream order as an option

image

library(stringr)
library(dplyr)
library(bcdata)
library(sf)

# OPTION 1: GET WATERSHEDS FROM BC DATA ####

  ws <- bcdc_query_geodata("freshwater-atlas-watersheds") %>% 
    filter(WATERSHED_GROUP_CODE == "BOWR") %>% 
    collect()

  ws2 <- ws %>% select(WATERSHED_ORDER,FWA_WATERSHED_CODE) %>% 
    mutate(iFWA_1 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,1],
           iFWA_2 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,2],
           iFWA_3 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,3],
           iFWA_4 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,4],
           iFWA_5 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,5],
           iFWA_6 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,6],
           iFWA_7 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,7],
           iFWA_8 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,8],
           iFWA_9 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,9],
           iFWA_10 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,10],
           iFWA_11 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,11],
           iFWA_12 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,12],
           iFWA_13 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,13],
           iFWA_14 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,14],
           iFWA_15 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,15],
           iFWA_16 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,16],
           iFWA_17 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,17],
           iFWA_18 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,18],
           iFWA_19 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,19],
           iFWA_20 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,20],
           iFWA_21 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,21]) 

  ws_all <- bind_rows(lapply(paste0("iFWA_",1:21), function(i=paste0("iFWA_",1:21)[1]){
    print(i)
    ws3 <- ws2 %>% rename(iFWA = i)
    ws3 %>% pull(iFWA) %>% unique() %>% length()
    ws3 %>% filter(as.numeric(iFWA) > 0) %>% group_by(iFWA) %>% summarize(WATERSHED_ORDER = max(WATERSHED_ORDER)) 
    }))

  st_write(ws_all, "C:/Users/bevin/Desktop/BOWR.gpkg")
# OPTION 2: GET WATERSHEDS FROM ####

  l <- st_layers("F:/FWA_WATERSHEDS_POLY/FWA_WATERSHEDS_POLY.gdb")[1] %>% as.data.frame() %>% 
                filter(!name %in% c("_max_watershed_feature_id",
                                    "_fwa_watersheds_poly_updates_BACKUP",
                                    "_fwa_watersheds_poly_updates"))
  lapply(l[,1], function(id = "PARS"){

    if(!file.exists(paste0("C:/Users/bevin/Desktop/",id,".gpkg"))){
    ws <- st_read("F:/FWA_WATERSHEDS_POLY/FWA_WATERSHEDS_POLY.gdb", layer = id)

    ws2 <- ws %>% select(WATERSHED_ORDER,FWA_WATERSHED_CODE) %>% 
      mutate(iFWA_1 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,1],
             iFWA_2 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,2],
             iFWA_3 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,3],
             iFWA_4 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,4],
             iFWA_5 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,5],
             iFWA_6 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,6],
             iFWA_7 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,7],
             iFWA_8 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,8],
             iFWA_9 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,9],
             iFWA_10 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,10],
             iFWA_11 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,11],
             iFWA_12 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,12],
             iFWA_13 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,13],
             iFWA_14 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,14],
             iFWA_15 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,15],
             iFWA_16 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,16],
             iFWA_17 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,17],
             iFWA_18 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,18],
             iFWA_19 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,19],
             iFWA_20 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,20],
             iFWA_21 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,21]) 

    ws_all <- bind_rows(lapply(paste0("iFWA_",1:21), function(i=paste0("iFWA_",1:21)[1]){
      print(i)
      ws2 %>% dplyr::rename(iFWA = i )%>% filter(as.numeric(iFWA) > 0) %>% group_by(iFWA) %>% summarize(WATERSHED_ORDER = max(WATERSHED_ORDER, na.rm = T))
      }))

    ws_all <- ws_all %>% 
      st_make_valid() %>% 
      mutate(id = "PARS", area_m = round(as.numeric(st_area(.)), 0)) %>% 
      arrange(-WATERSHED_ORDER) 

    st_write(ws_all, paste0("C:/Users/bevin/Desktop/",id,".gpkg"))}
    })
bevingtona commented 8 months ago

Related to #2

Jhydromet commented 8 months ago

I didn't look super closely at this code but I've noticed rolling up FWA_WATERSHED_CODE upstream runs into problems in some areas where two line rivers occur.

Stuart lake has same FWA_WATERSHED_CODE as Stuart river, so if you attempt to roll up codes from the outlet of the lake, it will include the entire Stuart river, and therefore all tribes to the river below the lake.

On Sat, Jan 13, 2024, 11:19 a.m. Alexandre Bevington < @.***> wrote:

Related to #2 https://github.com/bcgov/watershedBC/issues/2

— Reply to this email directly, view it on GitHub https://github.com/bcgov/watershedBC/issues/18#issuecomment-1890744782, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMVJ7EZAU3UVVUDLGPZFRPTYOLM55AVCNFSM6AAAAABBZQSRR6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQOJQG42DINZYGI . You are receiving this because you are subscribed to this thread.Message ID: @.***>

bevingtona commented 8 months ago

Good point!

I also noticed the way I showed above only computes by WATERSHED_GROUP_CODE... so the Fraser looks like this:

image

I'm testing a way to combine all of the WATERSHED_GROUP_CODE outputs such that I can then dissolve the entire Fraser into a single polygon, and so forth. I'm not certain this will fix the Stuart Lake issue, but I think it will. Seeing as the Stuart Lake watersheds will have the initial Fraser code of 100, then the Nechako, then Stuart River, etc. Standby!

bevingtona commented 8 months ago

This works OK, it builds nested watersheds for each stream segment by stream order for the freshwater atlas. Good for nested basins, bad for point of interest basins.

image

l <- st_layers("FWA_WATERSHEDS_POLY.gdb")[1] %>% as.data.frame() %>% 
  filter(!name %in% c("_max_watershed_feature_id",
                      "_fwa_watersheds_poly_updates_BACKUP",
                      "_fwa_watersheds_poly_updates"))

ws_d <- lapply(l[,1], function(id = "PARS"){

  if(!file.exists(paste0(id,".gpkg"))){

    ws <- st_read("FWA_WATERSHEDS_POLY.gdb", layer = id)

    ws2 <- ws %>% select(WATERSHED_ORDER,FWA_WATERSHED_CODE) %>% 
      mutate(iFWA_1 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,1],
             iFWA_2 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,2],
             iFWA_3 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,3],
             iFWA_4 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,4],
             iFWA_5 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,5],
             iFWA_6 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,6],
             iFWA_7 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,7],
             iFWA_8 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,8],
             iFWA_9 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,9],
             iFWA_10 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,10],
             iFWA_11 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,11],
             iFWA_12 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,12],
             iFWA_13 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,13],
             iFWA_14 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,14],
             iFWA_15 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,15],
             iFWA_16 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,16],
             iFWA_17 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,17],
             iFWA_18 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,18],
             iFWA_19 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,19],
             iFWA_20 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,20],
             iFWA_21 = str_split(FWA_WATERSHED_CODE, "-", simplify = T)[,21]) 

    ws_all <- bind_rows(lapply(paste0("iFWA_",1:21), function(i=paste0("iFWA_",1:21)[1]){
      print(i)
      ws2 %>% dplyr::rename(iFWA = i )%>% filter(as.numeric(iFWA) > 0) %>% group_by(iFWA) %>% summarize(WATERSHED_ORDER = max(WATERSHED_ORDER, na.rm = T)) %>% mutate(class = i)
    }))

    ws_all <- ws_all %>% 
      st_make_valid() %>% 
      mutate(id = id, area_m = round(as.numeric(st_area(.)), 0)) %>% 
      arrange(-WATERSHED_ORDER) 

    st_write(ws_all, paste0(id,".gpkg"))
  }
  ws_all
  })
l <- list.files(full.names = T, pattern = "*.gpkg$")

fwa_list <- paste0("iFWA_",1:21)

t <- lapply(fwa_list, function(fwa="iFWA_1"){

  print(fwa)

  all <- bind_rows(lapply(l, function(i=l[1]){

    print(paste(fwa, i))

    name <- st_layers(i)[[1]]

    st_read(i, query = paste0("SELECT * FROM ",name," WHERE class = '",fwa,"'"), quiet = T)

  }))

  all <- all %>% group_by(iFWA) %>% summarize(WATERSHED_ORDER= as.character(max(WATERSHED_ORDER)))
  st_write(all, paste0(fwa,".gpkg"))
  all
  })

tt <- bind_rows(t)
st_write(tt, "iFWA.gpkg")
Jhydromet commented 8 months ago

I found an error in the script I'd emailed, adding the corrected here for your reference. Also realized we were talking about two different methods! This code fails at lake outlets with two line rivers because the lake will have the same FWA_WATERSHED_CODE as its draining river.. It also fails if you want to map to a specific point along a two line river.

I think this works for small (non-two line river) watersheds.. but looking at LOCAL_WATERSHED_CODE doesn't present an obvious solution to these issues either.

All this is to say - Using the FWA Watersheds dataset alone will NOT allow you to reliably map watersheds to a point, but maybe a logic statement could determine whether the point is along a simple stream segment - then apply this method.

Be forewarned the query is slow for large watersheds.

image

`library(tidyverse) library(bcmaps) library(bcdata)

pp <- sf::st_sf(st_sfc(st_point(c(-124.270991, 54.417258))), crs =4326) %>% st_transform(src = 4326, crs = 3005)

fwa_watershed_code <- bcdata::bcdc_query_geodata("https://catalogue.data.gov.bc.ca/dataset/3ee497c4-57d7-47f8-b030-2e0c03f8462a") %>% filter(INTERSECTS(pp)) %>% dplyr::select(FWA_WATERSHED_CODE) %>% collect()

code <- str_remove_all(fwa_watershed_code$FWA_WATERSHED_CODE, "-000000")

cql <- paste0("FWA_WATERSHED_CODE like '%", code, "%'")

fwa_watershed_queried <- bcdata::bcdc_query_geodata("https://catalogue.data.gov.bc.ca/dataset/3ee497c4-57d7-47f8-b030-2e0c03f8462a") %>% filter(CQL(cql)) %>% collect()

fwa_watershed.m <- st_transform(fwa_watershed_queried, src = 3005, crs = 4326) %>% as_Spatial()

pp.m <- st_transform(pp, src = 3005, crs = 4326) %>% as_Spatial()`

Jhydromet commented 8 months ago

Just a follow up - Looks like from the FWA manual that its worth exploring the "DOWNSTREAM_ROUTE_MEASURE" to solve this double line river issue I'm having

bevingtona commented 8 months ago

Nice, thanks for the code! Agreed. The FWA watersheds aren't great for point of interest watershed mapping because of the inability to delineate watersheds smaller than the polygons provided in the FWA watersheds layer.

But I think we are trying to do different things. I'm not trying to do POI watersheds. I'm trying to roll up the FWA by watershed order such that I can pre-compute whole watersheds by stream order using the FWA. So if you click on a 3rd order stream that has no name, it will show you the roll up of that 3rd order stream with the upstream polygons.

So for your Stuart River example below Stuart Lake... if you select polygons that intersect Stuart River, then you'd get Stuart River, Nechako River and Fraser River:

image

I think I'm ready to add the "roll up by FWA watershed order" dataset to the app and I'll close the issue once the app is updated.

For your question of using the FWA for POI watersheds... maybe make a new issue? I think the only way is to do a hybrid raster-vector watershed delineation. So use your method to map the watershed to the closest watershed polygon, then use a raster based watershed delineation script in only the lowest watershed polygon to resolve the sub-basin watershed.... not sure that is well explained. In any case, I think it is a separate question.

Jhydromet commented 8 months ago

This makes sense - feel free to close the issue! I'm nervous that the raster method will lead to errors related to snapping points.. but surely your suggestion on querying closest FWA watershed would minimize that risk.. I'm giving one final crack at doing this solely with FWA queries, otherwise I"ll open an issue on your tool!