search

HARPgroup / cbp_wsm

1 stars 0 forks source link

Subwatershed Creation in cbp/hspf #85

Open rburghol opened 1 year ago

rburghol commented 1 year ago

Overview

This will be used to create new sub-watersheds by sub-dividing old watersheds in the base CBP model domain.
This process is completed, with the final operational script detailed here: https://github.com/HARPgroup/cbp_wsm/issues/98

Tasks

Routing

Hydro parameters:

FTABLES

Land Use

WDM files

Testing

rburghol commented 1 year ago

Making river scenario vahydro_2022

mkdir input/param/river/vahydro_2022
mkdir input/param/river/vahydro_2022/ftables
mkdir input/param/river/vahydro_2022/variable_ftables
# copy all segments for the James River
cbp copy_river_parameters.csh P620171001WQg vahydro_2022 J

This results in:

ls input/param/river/vahydro_2022/
NUTRX.csv  OXRX.csv  PLANK.csv  RQUAL.csv  SEDTRN.csv  ftables  variable_ftables

And FTABLES:

ls input/param/river/vahydro_2022/ftables/
JA0_7291_7290.ftable  JL1_7080_7190.ftable  JL6_6990_6960.ftable  JU2_7140_7330.ftable
JA1_7600_7570.ftable  JL1_7170_6800.ftable  JL6_7150_6890.ftable  JU2_7180_7380.ftable
JA1_7640_7280.ftable  JL1_7190_7250.ftable  JL6_7160_7440.ftable  JU2_7360_7000.ftable
JA2_7290_0001.ftable  JL1_7200_7250.ftable  JL6_7320_7150.ftable  JU2_7450_7360.ftable
JA2_7410_7470.ftable  JL1_7530_7430.ftable  JL6_7430_7320.ftable  JU3_6380_6900.ftable
JA2_7550_7280.ftable  JL2_6240_6520.ftable  JL6_7440_7430.ftable  JU3_6640_6790.ftable
JA2_7570_7480.ftable  JL2_6440_6441.ftable  JL7_6800_7070.ftable  JU3_6650_7300.ftable
JA4_7280_7340.ftable  JL2_6441_6520.ftable  JL7_7030_6800.ftable  JU3_6790_7260.ftable
JA4_7340_7470.ftable  JL2_6850_6890.ftable  JL7_7070_0001.ftable  JU3_6950_7330.ftable
JA4_7470_7480.ftable  JL2_7110_7120.ftable  JL7_7100_7030.ftable  JU3_7400_7510.ftable
JA5_7480_0001.ftable  JL2_7120_6970.ftable  JU1_6290_6590.ftable  JU3_7490_7400.ftable
JB0_7051_0001.ftable  JL2_7240_7350.ftable  JU1_6300_6650.ftable  JU4_7000_7300.ftable
JB0_7052_0001.ftable  JL2_7250_7090.ftable  JU1_6340_6650.ftable  JU4_7260_7380.ftable
JB1_8090_0001.ftable  JL2_7350_7090.ftable  JU1_6590_6600.ftable  JU4_7330_7000.ftable
JB2_7800_0001.ftable  JL3_7020_7100.ftable  JU1_6880_7260.ftable  JU4_7380_7160.ftable
JB3_6820_7053.ftable  JL3_7090_7150.ftable  JU1_7560_7500.ftable  JU5_7300_7510.ftable
JB3_7053_0001.ftable  JL4_6520_6710.ftable  JU1_7630_7490.ftable  JU5_7420_7160.ftable
JL1_6560_6440.ftable  JL4_6710_6740.ftable  JU1_7690_7490.ftable  JU5_7500_7420.ftable
JL1_6760_6910.ftable  JL6_6740_7100.ftable  JU1_7750_7560.ftable  JU5_7510_7500.ftable
JL1_6770_6850.ftable  JL6_6890_6990.ftable  JU2_6410_6640.ftable
JL1_6910_6960.ftable  JL6_6960_6970.ftable  JU2_6600_6810.ftable
JL1_6940_7200.ftable  JL6_6970_6740.ftable  JU2_6810_6900.ftable

Add additional files left out by copy_river_parameters.csh:

rburghol commented 1 year ago

Verify/Test FTABLE location.

cbp copy_river_parameters_all.csh P620171001WQg vahydro_2022
# edit config/control/river/vadeq_2021.con
# and set PARAMETERS = vahydro_2022
# now try to run_rug
cbp run_rug.csh vadeq_2021 PS2_5550_5560
# runs great! See result in http://deq1.bse.vt.edu:81/p6/vadeq/tmp/uci/river/vadeq_2021/PS2_5550_5560.uci
rburghol commented 1 year ago

Create a Test Sub-watershed

rburghol commented 1 year ago

Code snippet to create a model/load an existing model:

riverseg<- RomFeature$new(
  ds,
  list(
    hydrocode=paste('vahydrosw_wshed_',rseg_name, sep = ''),
    ftype=rseg_ftype,
    bundle='watershed'
  ),
  TRUE
)

model <- RomProperty$new(
  ds,
  list(
    varkey="om_model_element", 
    propname=riverseg$name,
    featureid=riverseg$hydroid, 
    entity_type="dh_feature", 
    propcode=model_version
  ), 
  TRUE
)
model$save(TRUE)
megpritch commented 1 year ago

PS2_5560_5100_linville_creek

Downstream: 5560 New ID: ? (1) choose random # , not pre-existing (2) +1 to downstream: 5561

Generate new ID = standalone script

juliabruneau commented 1 year ago

subsheds_naming.R

Link to script: https://github.com/HARPgroup/HARParchive/blob/master/HARP-2022-Summer/AutomatedScripts/subsheds_naming.R

Use: .../subsheds_naming.R [subshed name] [path to save new name] Ex: Rscript ~/HARParchive/HARP-2022-Summer/AutomatedScripts/subsheds_naming.R PS2_5560_5100_linville_creek ~/HARParchive/HARP-2022-Summer/AutomatedScripts/Subsheds

Result: PS2_5561_5560.txt with the content:

PS2_5561_5560

EDITS:

megpritch commented 1 year ago

Revisions to subsheds_naming.R

The repeat loop mentioned before now looks like this:

repeat {
  if(length(as.numeric(unique_ids)) >= 9998) {
    print('No unique names available.')
    break #script stops here and will not write a new name
  }
  if(sub_id %in% unique_ids) {
    sub_id <- sub_id + 1
  } 
  if(sub_id > 9999) { 
    sub_id <- 0002
  }
  else {
    sub_id <- formatC(sub_id, width = 4, format = "d", flag = "0") #make sure new id stays 4 digits
    new_name <- paste(splits[[1]][[1]],sub_id, downstr, sep = "_") 
    row <- data.frame(new_name, wordname)
    colnames(row) <- colnames(full_list)
    new_list <- rbind(head(full_list,-1), row)
    new_list <- new_list[order(new_list$river), ]
    end <- data.frame(tail(full_list,1))
    new_list <- rbind(new_list, end)

    write.table(new_list,
                file=list,
                append = TRUE,
                quote = FALSE,
                sep = ",",
                row.names = FALSE,
                col.names = TRUE)

    break
  }
}
rburghol commented 1 year ago

Note @megpritch @juliabruneau : line 1 (the header) must end with "***" not "...", as it causes the routine to bomb :)

juliabruneau commented 1 year ago

Note @megpritch @juliabruneau : line 1 (the header) must end with "***" not "...", as it causes the routine to bomb :)

@rburghol there was difficulties in matching what was interpreted as in R. It reads in R as ... . It seems that the header remained as after this morning's testing, so can we assume that ... is still saved as ?

rburghol commented 1 year ago

@megpritch @juliabruneau -- sorry, I changec that manually so that I could test some things. Gotta figure out a way to make that work :). I suggest setting the value of that first line, 2nd column to "name***", right before saving and see what happens.

juliabruneau commented 1 year ago

Got it - we will include that!

juliabruneau commented 1 year ago

/opt/model/meta_model/run_model hsp2_cbp6 subsheds PS2_5568_5560 auto river 04_ps_sep_div

juliabruneau commented 1 year ago

subsheds_naming.bat

Use: bash /subsheds_naming.bat [path to subshed list] [master csv] [model version] Ex: bash HARParchive/HARP-2022-Summer/AutomatedScripts/SubshedsCreation/subsheds_naming.bat HARParchive/HARP-2022-Summer/AutomatedScripts/SubshedsCreation/subshed_riversegs.txt rivernames.csv cbp-6.0

Ex: OD3_8720_8900_ut_leatherwood from text file Result: OD3_8723_8720

In VAHydro: image

  • [ ] set path name of routing file for R routine

I'm not exactly sure what this means/if it is included in the batch script?

megpritch commented 1 year ago

Re-run of get_subshed_riversegs.R

When trying to run the script in terminal to update our list of subshed names, I get Error: C stack usage 9428040 is too close to the limit. Whenever this script is run in R, Rstudio freezes and the REST retrieval always takes 5-10 minutes because it is such a large dataset before we subset it into just subwatersheds. I think this may be related to the terminal error, but I'm not sure what I can change to fix it. For now, I am generating the subshed list in Rstudio and pushing to Github.

image

For reference, this is the whole script minus the setup:

rsegs <- RomFeature$new(ds,list(ftype='vahydro',bundle='watershed'),TRUE)
subshed_hydrocodes <- subset(rsegs$hydrocode, nchar(rsegs$hydrocode)>29)
subshed_hydrocodes <- data.frame(subshed_hydrocodes) #collect hydrocodes into data frame

#make hydrocodes into riverseg's only (for ftable_creation purposes)
subshed_riversegs <- gsub("vahydrosw_wshed_","",as.character(subshed_hydrocodes$subshed_hydrocodes))
subshed_riversegs <- data.frame(subshed_riversegs)

#write to csv
file <- paste(path, 'subshed_riversegs.txt', sep='')
write.table(subshed_riversegs, file = file, quote = FALSE, row.names = FALSE, col.names = FALSE)
megpritch commented 1 year ago

Diagram & Comments for subsheds_naming.R

@durelles Clearer comments have been added to the naming script and here is the jamboard we made while organizing the script:

image
juliabruneau commented 1 year ago

proportion_landuse.R

The main steps in this Rscript:

  1. Pull the subshed area value from VAHydro
  2. Sum together the land segments of the main watershed
  3. Calculate the proportion: subshed area/watershed area
  4. Checking if proportioning was performed already in the list: 4.1 Look for existing subsheds; sum together; calculate the proportion 4.2 if the proportion from VAHydro equals the one calculated from the list == end the script
  5. If no subsheds exists; subtract the proportion from main watershed & add that proportion to the new subshed (making sure this is done to ALL land segments within a river segment)
  6. Re-write the table with the new proportioned land use values

I had a question about using the sqldf; I'm receiving an error when trying to use our argument for the river segment, instead of 'PS2_5560_5100' typed out. Can you notice an error in this:

receiving_landuses <- sqldf("select * from landuse_full where riverseg = 'main_seg'")

(main_seg = PS2_5560_5100)

rburghol commented 1 year ago

Find watershed areas for all models of a given type: 


library("sqldf")
library("stringr") #for str_remove()
library("hydrotools")
library("openmi.om")

# Load Libraries
basepath='/var/www/R';
site <- "http://deq1.bse.vt.edu/d.dh"    #Specify the site of interest, either d.bet OR d.dh
source("/var/www/R/config.local.private");
source(paste(basepath,'config.R',sep='/'))

# get the DA, need to grab a model output first in order to ensure segments with a channel subcomp
# are included
#

# GET DA
df_area <- data.frame(
  'model_version' = c('vahydro-1.0',  'vahydro-1.0',  'vahydro-1.0'),
  'runid' = c('runid_11', '0.%20River%20Channel', 'local_channel'),
  'runlabel' = c('QBaseline_2020', 'comp_da', 'subcomp_da'),
  'metric' = c('Qbaseline', 'drainage_area', 'drainage_area')
)
da_data <- om_vahydro_metric_grid(
  metric = metric, runids = df_area,
  base_url = paste(site,'entity-model-prop-level-export',sep="/")
)
da_data <- sqldf(
  "select pid, comp_da, subcomp_da, riverseg,
   CASE
    WHEN comp_da is null then subcomp_da
    ELSE comp_da
    END as da
   from da_data
  ")
megpritch commented 1 year ago

@rburghol When trying to run the code you gave us we receive the error Error in httr::add_headers(HTTP_X_CSRF_TOKEN = token) : object 'token' not found :

image
rburghol commented 1 year ago

FYI @juliabruneau - 2 follow-ups:

rburghol commented 1 year ago

@megpritch -- Question are you running this from the deq2 machine or your local machine?

Maybe it is as simple as using a slightly newer syntax to load the datasource stuff. Could you try the following and let me know?

library("sqldf")
library("stringr") #for str_remove()
site <- "http://deq1.bse.vt.edu/d.dh"    #Specify the site of interest, either d.bet OR d.dh
#----------------------------------------------
# Load Libraries
basepath='/var/www/R';
source(paste(basepath,'config.R',sep='/'))

instead of:

library("sqldf")
library("stringr") #for str_remove()
library("hydrotools")
library("openmi.om")

# Load Libraries
basepath='/var/www/R';
site <- "http://deq1.bse.vt.edu/d.dh"    #Specify the site of interest, either d.bet OR d.dh
source("/var/www/R/config.local.private");
source(paste(basepath,'config.R',sep='/'))
megpritch commented 1 year ago

@rburghol - I am running from my local machine. The chunk you just sent did not work either. Neither did our usual method: 

# Setup
suppressPackageStartupMessages(library("hydrotools")) #needed to pull values from VAHydro 
# Link data source (ds) to VAHydro
basepath='/var/www/R';
source("/var/www/R/config.R") #will need file in same folder/directory
ds <- RomDataSource$new(site, rest_uname = rest_uname)
ds$get_token(rest_pw)

However, every time we run we still get this: 

reading from http://deq1.bse.vt.edu:81/d.dh
[1] "REST AUTH INFO HAS BEEN SUPPLIED"
[1] "RETRIEVING REST TOKEN"
[1] "Login attempt successful"
[1] "token = _3vBjqU1cGUXZSQWjl1RKwShbxbNZcPhqktSlrImLMQ"
>

So I am not sure why it cannot find the token.

Julia ran what you sent in the terminal and it works, but uses a different token : token = 'wrRjwZK6jtr975TbkEmwxBtEbW3ff_ICJ5PzMEFsrzM' If we manually put this token in R then it works on the local machine as well.

rburghol commented 1 year ago

Ahh @megpritch good debugging -- the function om_vahydro_metric_grid() assumes that token is a global variable, that is, token is NOT passed in as an argument. That is poor coding practice, and honestly, I can't believe I haven't realized this by now.
The message that you showed:

[1] "Login attempt successful"
[1] "token = _3vBjqU1cGUXZSQWjl1RKwShbxbNZcPhqktSlrImLMQ"

Comes from the config.R script, and sets the token variable to allow secure conversations with the REST service. So it should be set. And also, whatever is the value of that token should be the value that would work for you when you type your token = 'wrRj...' statement.

@jdkleiner -- I am drafting a new version of this function, instead of modifying the one that we have, We should send the ds in as a param, and use the auth retrieval metric there. If the ds is not sent in, it uses the old method (or tries) and tells the user that it is deprecated. It will require reinstalling hydrotools however... which isn't too painful. Steps:

Script 1: Update Hydrotools. 

library("devtools")
install_github("HARPgroup/hydro-tools")

Then add Script 2 into your code at the beginning.

Script 2: Load the datasource into variable ds at beginning of scripts. Finally, try this code

# add this after source("/var/www/R/config.R")
ds = RomDataSource$new(site, rest_uname)
ds$get_token(rest_pw)

Script 3: Add ds = ds to the argument that you pass into om_vahydro_metric_grid().


da_data <- om_vahydro_metric_grid(
  metric = metric, runids = df_area,
  base_url = paste(site,'entity-model-prop-level-export',sep="/"),
  ds = ds
)

Are you all able to move forward with the data, and you said it does work on deq2, so, we can at least tackle the next step?

juliabruneau commented 1 year ago

@rburghol I did all the steps above, and received this error:

[1] "Found info"
  model_version    runid       runlabel    metric
1   vahydro-1.0 runid_11 QBaseline_2020 Qbaseline
[1] "retrieving  http://deq1.bse.vt.edu:81/d.dh/entity-model-prop-level-export/all/dh_feature/watershed/vahydro/vahydro-1.0/runid_11/Qbaseline"
 Error in vahydro_auth_read(url, token, ctype = "text/csv", delim = ",") : 
object 'token' not found
rburghol commented 1 year ago

Doing some testing here. First time thru works like a champ. Second time through has some trouble. One issue is that the old subwatershed does not get deleted, and the area calcs get squirrely. I am investigating.

juliabruneau commented 1 year ago

success with proportion_landuse.R (now area_propor.R)

Doing some testing here. First time thru works like a champ. Second time through has some trouble. One issue is that the old subwatershed does not get deleted, and the area calcs get squirrely. I am investigating.

Next:

megpritch commented 1 year ago

Two New Generic-Use Scripts and their Usage

copy_parent.R

Takes the main riverseg values and duplicates them for the subshed. Use: Rscript /copy_parent.R [file to edit] [subshed] [main_seg] Example: Rscript HARP-2022-Summer/AutomatedScripts/SubshedsCreation/copy_parent.R 'HARP-2022-Summer/AutomatedScripts/SubshedsCreation/wF180615RXAPXXXW_l2w.csv' 'PS2_5568_5560' 'PS2_5560_5100'

Files this works with:

area_propor.R

Takes the main riverseg values and proportions them for the subshed. Use: Rscript /area_propor.R [file to edit] [subshed] [main_seg] [subshed_da] [columns_to_proportion] Example: Rscript HARP-2022-Summer/AutomatedScripts/SubshedsCreation/area_propor.R 'HARP-2022-Summer/AutomatedScripts/SubshedsCreation/land_water_area.csv' 'PS2_5568_5560' 'PS2_5560_5100' '46.04' '-1:-2'

Files this works with:

rburghol commented 1 year ago

@megpritch @juliabruneau -- great work pushing this forward. Thanks also for putting up some example use/testing scripts. I will put mine below, as well. Let's get together on the batch script at your earliest convenience. To your questions: I think that SCRORG should also be done with the area weighting technique -- as noted, we don't actually use that (it is for organic matter scouring), it is just necessary to use some of the CBP routines that require it when analyzing hydrology. Also, septic should be handled specially wit ha copy statement. Both are described in the list at the top of the issue -- they are not perfect solutions, but they will work given that septic and sediment don't influence our hydrology, and these will be overwritten when the CBP simulates this in the future, so they are just a place holder to allow us to run hydrology.

My tests:

# copy the land_water_area file to the local directory for testing purposes
cp config/catalog/geo/vahydro/land_water_area.csv lwa.csv

# run the script for area 46.04
Rscript /opt/model/HARParchive/HARP-2022-Summer/AutomatedScripts/SubshedsCreation/proportion_landuse.R PS2_5560_5100 PS2_5568_5560 46.04 lwa.csv
# look for the receiving segment areas, to make sure no duplicates: 
fgrep 5560_ lwa.csv
PS2_5560_5100,N51165,103656.746042785
PS2_5560_5100,N51171,57427.0449502788
PS2_5560_5100,N51660,497.759821436728

# look for the subwatershed segment areas, to make sure no duplicates: 
fgrep 5568_ lwa.csv
PS2_5568_5560,N51165,18902.5801572155
PS2_5568_5560,N51171,10472.2496297213
PS2_5568_5560,N51660,90.7702130632719

# run the script for a new area:
Rscript /opt/model/HARParchive/HARP-2022-Summer/AutomatedScripts/SubshedsCreation/proportion_landuse.R PS2_5560_5100 PS2_5568_5560 25.1 lwa.csv

# look for the receiving segment areas, to make sure no duplicates: 
fgrep 5560_ lwa.csv
PS2_5560_5100,N51165,112254.053351475
PS2_5560_5100,N51171,62190.0533613641
PS2_5560_5100,N51660,539.044101661422

# look for the subwatershed segment areas, to make sure no duplicates: 
fgrep 5568_ lwa.csv
PS2_5568_5560,N51165,10305.2728485254
PS2_5568_5560,N51171,5709.24121863603
PS2_5568_5560,N51660,49.4859328385777

# go back to the original subwatershed area to insure that it is reversible:
Rscript /opt/model/HARParchive/HARP-2022-Summer/AutomatedScripts/SubshedsCreation/proportion_landuse.R PS2_5560_5100 PS2_5568_5560 46.04 lwa.csv

fgrep 5560_ lwa.csv
PS2_5560_5100,N51165,103656.746042785
PS2_5560_5100,N51171,57427.0449502789
PS2_5560_5100,N51660,497.759821436728

fgrep 5568_ lwa.csv
PS2_5568_5560,N51165,18902.5801572155
PS2_5568_5560,N51171,10472.2496297212
PS2_5568_5560,N51660,90.7702130632716

# now try with land use areas
subshed="PS2_5568_5560"
downstream="PS2_5560_5100"
darea=46.04
fgrep 5560_ $lu_file
PS2_5560_5100,N51165,1756.573852539,1543.499145508,3520.169189453,960.945739746,3001.164550781,2712.096435547,478.605255127,312.409790039,3272.076416016,93.118621826,232.714904785,953.464050293,13167.125976563,31907.837890625,103.151504517,212.606628418,1928.662066817,3122.345136404,170.526861668,960.980539322,8342.058493316,112.452880859,38.694198608,0.944061875,8.76594162,46.049324036,53.806312561,0,0,0,0,0,0,0,0,254.137008667,39939.290039063,11.329862595,43.119621277,2424.030373573,874.579755783,0,0,0
PS2_5560_5100,N51171,649.578063965,1515.272338867,2004.469360352,941.467163086,626.188842773,1383.432739258,244.13520813,574.1328125,335.051696777,73.389328003,169.168136597,973.683044434,6561.417480469,13924.020507813,19.79283905,46.106086731,1093.889055014,1431.119960785,146.281178117,505.34324193,2904.128371701,0,0,0,0,0,2.084458113,0,0,0,0,0,0,0,0,107.269371033,30199.81640625,17.853826523,19.131643295,864.975316048,566.09667635,0,0,0
PS2_5560_5100,N51660,1.07519424,3.665212154,1.136198282,5.175999641,1.657707453,6.216020584,1.096945405,1.806896925,0,0.131079838,0.327385128,1.336273193,8.666325569,16.483577728,0.084208697,0.174154788,29.797811508,3.61521244,0,0.829137743,18.563562393,41.381149292,92.106185913,0.349092931,12.327147484,176.085693359,4.279939651,0,0,0,0,0,0,0,0,0,129.676681519,0,0,21.864974976,8.620264053,0,0,0

Rscript /opt/model/p6/vadeq/run/resegment/area_propor.R $lu_file $subshed $downstream $darea '-1:-2'
# breaks.
rburghol commented 1 year ago

@juliabruneau @megpritch Here is what I'm thinking for the batch script:

Error in Proportion

subshed="PS2_5568_5560"
downstream="PS2_5560_5100"
darea=46.04

# using the scripts I copied
Rscript /opt/model/p6/vadeq/run/resegment/area_propor.R $lu_file $subshed $downstream $darea '-1:-2'
Error in `[.data.frame`(subsheds, cols) : undefined columns selected
Calls: area_propor -> [ -> [.data.frame
Execution halted

# using the original syntax 
Rscript /opt/model/HARParchive/HARP-2022-Summer/AutomatedScripts/SubshedsCreation/area_propor.R '/opt/model/HARParchive/HARP-2022-Summer/AutomatedScripts/SubshedsCreation/land_water_area.csv' 'PS2_5568_5560' 'PS2_5560_5100' '46.04' ',-1:-2'
Error in `[.data.frame`(subsheds, cols) : undefined columns selected
Calls: area_propor -> [ -> [.data.frame
Execution halted

Proposed code for add_sub_watershed.sh

# we give our script:
# - hydrocode
# - receiving stream segment
# - drainage area of subshed
# - model version (cbp-6.0)
# - cbp scenario (subsheds)
# load name, rseg and area info
hydrocode=$1
downstream=$2
darea=$3
model_version=$4
scenario=$5

# get info
GEO=`cbp get_config subsheds river GEO`

# name it or retrieve the name if it already exists
read -r subshed downstream <<< "$(Rscript /opt/model/p6/vadeq/run/resegment/subsheds_naming.R PS2_5560_5100_linville_creek /opt/model/p6/vadeq/config/catalog/${GEO}/vahydro/rivernames.csv cbp-6.0)"

# set watershed area
Rscript /opt/model/p6/vadeq/run/resegment/area_propor.R /opt/model/p6/vadeq/config/catalog/${GEO}/vahydro/land_water_area.csv $subshed $downstream $darea

# set land use area
# iterate thru multiple lines if they exist (subsheds has one but others have more)
# see next comment for my bash work to iterate through multiple potential land use data sets

# copy transport stuff
Rscript /opt/model/p6/vadeq/input/param/transport/wF180615RXAPXXXW_l2w.csv $subshed $downstream
# ... repeat for other transport things
rburghol commented 1 year ago

Iterate through multiple land use datasets

LANDUSE=`cbp get_config subsheds river 'LAND USE'`
# echo $LANDUSE
# 1984 07 01 2013VAHYDRO2018615 2055 07 01 2013VAHYDRO2018615


- To iterate through, we know that each line has 4 values, year, month, day, and dataset name
- so we iterate through and know that every time we reach the 4th, 8th, etc element we have our dataset

cnt=0 for i in $LANDUSE; do ((cnt++)) if [[ $cnt -eq 1 ]]; then yr=$i fi if [[ $cnt -eq 4 ]]; then lu_file="input/scenario/river/land_use/landuse${i}.csv" echo "LU file for $yr = " $lu_file Rscript /opt/model/p6/vadeq/run/resegment/area_propor.R $lu_file $subshed $downstream $darea

# reset our counter
cnt=0

fi done

juliabruneau commented 1 year ago

Thank you @rburghol for the work on the batch. Right off the bat, I noticed that we input drainage area manually as an argument; isn't there a way that we could extract it from get_subsheds_riverseg , since we modified it to create the table of subsheds with da?

As below:

read -r subshed downstream <<< "$(Rscript ...

I'm now using the same strategy to extract drainage area from get_subsheds_riverseg.R

read -r darea <<< "$(Rscript /opt/model/p6/vadeq/run/resegment/get_subsheds_riverseg.R /opt/model/p6/vadeq/config/catalog/${GEO}/vahydro/)"

Please correct me if I understood that code wrong. I added cat(paste(darea)) into the Rscript in order to extract the correct value.

juliabruneau commented 1 year ago

debugging area_propor.R

(@rburghol ) I created a copy of this file (area_propor_copy.R) in order to isolate where we are receiving the error, and based on my messages:

:~/HARParchive$ Rscript HARP-2022-Summer/AutomatedScripts/SubshedsCreation/area_propor_copy.R 'HARP-2022-Summer/AutomatedScripts/SubshedsCreation/land_water_area.csv' 'PS2_5568_5560' 'PS2_5560_5100' '46.04' '-1:-2'
Warning message:
no DISPLAY variable so Tk is not available
main segs and subsheds subset
row of zeros created
Error in `[.data.frame`(subsheds, , cols) : undefined columns selected
Calls: area_propor -> [ -> [.data.frame
Execution halted

We are getting an error only when run in the terminal, when calculating the sum of areas: sub_area <- sum(subsheds[cols])/640 . I have tested on my local R and in R on the terminal, and errors do not occur at any point.

rburghol commented 1 year ago

@juliabruneau I am taking a look at get_subsheds_riverseg.R now, first time I have cracked it open and I will respond shortly!

rburghol commented 1 year ago

@juliabruneau OK - we only want to run the command to get_subshed_riversegs.R get the list of VAHydro-subshed models once saved as subshed_list.csv, to get our list of sub-watersheds from VAHydro, but we will not be using that later on. So, that is why we want to input DA to our script, as we won't want to make that REST call every time we run the script, and also ,we want to have the ability to change drainage area on these models manually, as we move forward there won't be a vahydro-1.0 model for new segments that we create.

juliabruneau commented 1 year ago

Success with area_propor.R

(@rburghol) After testing, the Rscript is working for both land_water_area.R and the land use files.

Ex:

juliasb@deq2:~/HARParchive$ Rscript HARP-2022-Summer/AutomatedScripts/SubshedsCreation/area_propor.R 'HARP-2022-Summer/AutomatedScripts/SubshedsCreation/land_water_area.csv' 'PS2_5568_5560' 'PS2_5560_5100' '46.04'
Warning message:
no DISPLAY variable so Tk is not available
proportioned the river segment and subshed

I will implement the changes in the batch now.

rburghol commented 1 year ago

@juliabruneau good stuff. In order to verify that the thing does work, rather than just no errors, I would like to use the fgrep statement from above in my tests with a clean file which can ensure no duplicates. I also want to replicate the tests where I proportion with one drainage area, then change that drainage area, then changed it back to ensure that it's reversible. They are back in this comment: https://github.com/HARPgroup/cbp_wsm/issues/85#issuecomment-1385672385

Per your question, there are at least two transport files that need to be dealt with, check the entries above. Every file that I listed above will ultimately have to be handled with one of our function routines. That means the transport files, the "hydro parameters" files, "land use" etc. Some of those WDM files will need to just be a copy from parent… I will try to put a checklist together when I have a minute, but they should be all up there in the bulleted list.

rburghol commented 1 year ago

Testing the component scripts on all the files before adding to the bat:

darea=9.52
subshed="JL1_6562_6560"
downstream="JL1_6560_6440"

Rscript $CBP_ROOT/run/resegment/copy_parent.R $CBP_ROOT/input/param/river/${PARAMS}/HYDR.csv $subshed $downstream
created subshed containing parent values
fgrep JL1_6562_6560 input/param/river/${PARAMS}/HYDR.csv
# JL1_6562_6560,1,17.01,249.28,NA,0.5,NA,293.16,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA
juliabruneau commented 1 year ago

@rburghol I'm trying to perform testing on the batch script - script-by-script, and I ran into this hydrotools error:

Error: package or namespace load failed for ‘hydrotools’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]):
 namespace ‘rlang’ 1.0.4 is already loaded, but >= 1.0.6 is required

How can I update hydrotools in my terminal? I tried in R in my terminal like we did before but it's not working.

juliabruneau commented 1 year ago

TESTING with different drainage areas

Goal with fgrep:

  1. New subshed created in rivernames.csv
  2. area_propor.R - Proportioning done correctly (no duplicates), ex csv: land_water_area.csv
  3. copy_parent.R - Copying the parent done correctly (no duplicates), ex csv: gen_info_rseg.csv

Does proportioning work correctly if da changes?

Scenario 1. YM2_6120_6430_ni_reservoir with da = 25.268455795755

/opt/model/p6/vadeq$ ./run/resegment/add_sub_watershed.bat YM2_6120_6430_ni_reservoir YM2_6120_6430 cbp-6.0 subsheds 25.268455795755

/opt/model/p6/vadeq$ fgrep YM2_6122_6120 config/catalog/geo/vahydro/rivernames.csv
YM2_6122_6120,Ni River Reservoir

/opt/model/p6/vadeq$ fgrep YM2_6122_6120 config/catalog/geo/vahydro/land_water_area.csv
YM2_6122_6120,N51033,3341.43704423947
YM2_6122_6120,N51137,318.983672820896
YM2_6122_6120,N51177,12511.3909922228
opt/model/p6/vadeq$ fgrep YM2_6120_6430 config/catalog/geo/vahydro/land_water_area.csv
YM2_6120_6430,N51033,30538.2575057605
YM2_6120_6430,N51137,2915.27430017911
YM2_6120_6430,N51177,114344.838707777

/opt/model/p6/vadeq$ fgrep YM2_6122_6120 input/param/river/vahydro_2022/gen_info_rseg.csv
YM2_6122_6120,3,0,S,60
/opt/model/p6/vadeq$ fgrep YM2_6120_6430 input/param/river/vahydro_2022/gen_info_rseg.csv
YM2_6120_6430,3,0,S,60

Scenario 2. YM2_6120_6430_ni_reservoir with da = 20

/opt/model/p6/vadeq$ fgrep YM2_6122_6120 config/catalog/geo/vahydro/land_water_area.csv
YM2_6122_6120,N51033,2644.74969998034
YM2_6122_6120,N51137,252.475794642333
YM2_6122_6120,N51177,9902.77450537733
/opt/model/p6/vadeq$ fgrep YM2_6120_6430 config/catalog/geo/vahydro/land_water_area.csv
YM2_6120_6430,N51033,31234.9448500196
YM2_6120_6430,N51137,2981.78217835767
YM2_6120_6430,N51177,116953.455194622

/opt/model/p6/vadeq$ fgrep YM2_6122_6120 input/param/river/vahydro_2022/gen_info_rseg.csv
YM2_6122_6120,3,0,S,60
/opt/model/p6/vadeq$ fgrep YM2_6120_6430 input/param/river/vahydro_2022/gen_info_rseg.csv
YM2_6120_6430,3,0,S,60

Scenario 3. YM2_6120_6430_ni_reservoir with da = 25.268455795755 same as in scenario 1.

/opt/model/p6/vadeq$ fgrep YM2_6122_6120 config/catalog/geo/vahydro/land_water_area.csv
YM2_6122_6120,N51033,3341.43704423948
YM2_6122_6120,N51137,318.983672820897
YM2_6122_6120,N51177,12511.3909922228
/opt/model/p6/vadeq$ fgrep YM2_6120_6430 config/catalog/geo/vahydro/land_water_area.csv
YM2_6120_6430,N51033,30538.2575057605
YM2_6120_6430,N51137,2915.27430017911
YM2_6120_6430,N51177,114344.838707777

/opt/model/p6/vadeq$ fgrep YM2_6122_6120 input/param/river/vahydro_2022/gen_info_rseg.csv
YM2_6122_6120,3,0,S,60
/opt/model/p6/vadeq$ fgrep YM2_6120_6430 input/param/river/vahydro_2022/gen_info_rseg.csv
YM2_6120_6430,3,0,S,60

Conclusion:

juliabruneau commented 1 year ago

TESTING - Northern River

Goal with fgrep:

  1. New subshed created in rivernames.csv
  2. area_propor.R - Proportioning done correctly (no duplicates), ex csv: land_water_area.csv
  3. copy_parent.R - Copying the parent done correctly (no duplicates), ex csv: gen_info_rseg.csv

Norther River Segment: RU4_5640_6030_lake_pelham with da = 26.15

/opt/model/p6/vadeq$ ./run/resegment/add_sub_watershed.bat RU4_5640_6030_lake_pelham RU4_5640_6030 cbp-6.0 subsheds 26.15

/opt/model/p6/vadeq$ fgrep RU4_5642_5640 config/catalog/geo/vahydro/rivernames.csv
RU4_5642_5640,Lake Pelham: Mountain Run

/opt/model/p6/vadeq$ fgrep RU4_5642_5640 config/catalog/geo/vahydro/land_water_area.csv
RU4_5642_5640,N51047,8209.21335628797
RU4_5642_5640,N51061,7183.42999536222
RU4_5642_5640,N51179,1343.35664834981
/opt/model/p6/vadeq$ fgrep RU4_5640_6030 config/catalog/geo/vahydro/land_water_area.csv
RU4_5640_6030,N51047,72415.502403712
RU4_5640_6030,N51061,63366.8135446378
RU4_5640_6030,N51179,11850.0814116502

/opt/model/p6/vadeq$ fgrep RU4_5642_5640 input/param/river/vahydro_2022/gen_info_rseg.csv
RU4_5642_5640,3,0,S,60
/opt/model/p6/vadeq$ fgrep RU4_5640_6030 input/param/river/vahydro_2022/gen_info_rseg.csv
RU4_5640_6030,3,0,S,60
megpritch commented 1 year ago

Runs That Result in Errors

./run/resegment/add_sub_watershed.bat [subshed hydrocode] [downstream rseg] cbp-6.0 subsheds [da]

Debugging: fgrep for the downstream river segment in the first file after generating the new subshed name: (Since that's where we wee the first error.)

/opt/model/p6/vadeq$ fgrep OR3_7740_8271 config/catalog/geo/vahydro/land_water_area.csv
megpritch commented 1 year ago

Runs That Have Been Successful (with fgrep too)

fgrep Routine: (re-stated just so it's easy to copy)

  1. Run .bat in /opt/model/p6/vadeq$: 
    ./run/resegment/add_sub_watershed.bat [subshed hydrocode] [downstream rseg] cbp-6.0 subsheds [da]

  2. Check for proportions/duplicates:

    fgrep [new subshed rseg] config/catalog/geo/vahydro/rivernames.csv

fgrep [new subshed rseg] config/catalog/geo/vahydro/land_water_area.csv
fgrep [downstr rseg] config/catalog/geo/vahydro/land_water_area.csv

fgrep [new subshed rseg] input/param/river/vahydro_2022/gen_info_rseg.csv
fgrep [downstr rseg] input/param/river/vahydro_2022/gen_info_rseg.csv

  3. Repeat: re-run .bat and check for duplicates