HARPgroup / om

Object-oriented Meta-model
0 stars 1 forks source link

elfgen for State Plan & VWP #216

Open jdkleiner opened 4 years ago

jdkleiner commented 4 years ago

Step 0.0 - Test recent elfgen package updates

library(elfgen)

# Example input parameters for retrieving data from VAHydro
watershed.code <- 'nhd_huc8_02070005'
watershed.bundle <- 'watershed'
watershed.ftype <- 'nhd_huc8'
x.metric <- 'erom_q0001e_mean'
y.metric <- 'aqbio_nt_total'
y.sampres <- 'species'
datasite <- 'RESTRICTED - CONTACT elfgen MAINTAINER'

# elfdata_vahydro() function for retrieving data from VAHydro
watershed.df <- elfdata_vahydro(watershed.code,watershed.bundle,watershed.ftype,x.metric,y.metric,y.sampres,datasite)
# clean_vahydro() function for cleaning data by removing any stations where the ratio of DA:Q is greater than 1000, also aggregates to the maximum richness value at each flow value
watershed.df <- clean_vahydro(watershed.df)

elf_quantile <- 0.80
breakpt <- bkpt_pwit("watershed.df" = watershed.df, "quantile" = elf_quantile, "blo" = 100, "bhi" = 1000)  

elf <- elfgen("watershed.df" = watershed.df,
              "quantile" = elf_quantile,
              "breakpt" = breakpt,
              "yaxis_thresh" = 53, 
              "xlabel" = "Mean Annual Flow (ft3/s)",
              "ylabel" = "Fish Species Richness")

elf$plot
elf$stats

Step 1.0 - Build automated code for generating ELF analysis for VWP scenarios

site <- "http://deq2.bse.vt.edu/d.dh"

inputs <- list( varkey = "erom_q0001e_mean", featureid = 320509, entity_type = "dh_feature" )

property dataframe returned

dataframe <- getProperty(inputs, site)


- [x] Generate elfgen plot & overlay analysis point on ggplot object
      - Example below of what I'm thinking. I want to be able to see where on the ELF curve an intake is located (using `[MEAN FLOW @ INTAKE]` on the x-axis), and the species richness associated with that mean flow
![elfgen_vwp_1](https://user-images.githubusercontent.com/29379385/88217929-4a474200-cc2d-11ea-85a3-55e88ff47b46.png)
- [x] Calculate percent richness change at the analysis point 
`richness_change(elf$stats, "pctchg" = [FLOW REDUCTION @ INTAKE], "xval" = [MEAN FLOW @ INTAKE])`
   - Calculate equations from ELF confidence intervals and evaluate for richness change +/- X
- [x] HARP latest code: https://github.com/HARPgroup/om/blob/master/R/summarize/waterSupply_elfgen.R
- [x] Additional HARP work: https://github.com/HARPgroup/hydro-tools/tree/master/HARP-2020/Elfgen
![fig elfgen nhd_huc8_02070005 erom_q0001e_mean aqbio_nt_total](https://user-images.githubusercontent.com/29379385/94742498-c3f15f80-0343-11eb-9164-817284a0d42d.png)
- `flow_reduction_pct <- 10`
  - **"The absolute richness change is 0.23 +/-0.12"** 
  - **"The percent richness change is 1.33 +0.81/-0.71"**

### Step 2.0 - Calculate Percent Richness Change by River Segment & runid (in R)
- [x] Iterate through each VAHydro river segment:

DOWNLOAD RSEG LAYER DIRECT FROM VAHYDRO

localpath <- tempdir() filename <- paste("vahydro_riversegs_export.csv",sep="") destfile <- paste(localpath,filename,sep="\") download.file(paste(site,"vahydro_riversegs_export",sep=""), destfile = destfile, method = "libcurl") RSeg.csv <- read.csv(file=paste(localpath , filename,sep="\"), header=TRUE, sep=",")


- [x] Use a single rseg for testing purposes
  - [x] Find the river segment outlet nhdplus segment
    - [x] Load nhdplus segments contained inside the VAHydro rseg
    - [x] Retrieve drainage area for nhdplus segments
      - [x] `nhdp_drainage_sqmi` property
      - [x] The following view should work (replacing the hydroid for the particuar rseg of interest)
`http://deq2.bse.vt.edu/d.dh/dh-feature-containing-export/68370/watershed/nhdplus/nhdp_drainage_sqmi`
    - [x] Select nhdplus segment with largest drainage area, this can be considered  the *rseg outlet*
  - [x] Retrieve `consumptive_use_frac` for rseg for each runid of interest
  - [x] Perform the following using Fixed Method `breakpt <- 530 cfs`
    - [x] Calculate percent richness change +/- X for this `consumptive_use_frac` at the *rseg outlet* for containing **HUC8**
    - [x] Calculate percent richness change +/- X for this `consumptive_use_frac` at the *rseg outlet*  for containing **HUC10**
- [x] POST the richness change via REST to VAHydro rseg runid scenario - use propname "elfgen_richness_change_[HUC Level]"
  - [x] Data model _(may get revised as we go)_:
![elfgen Application_3 2](https://user-images.githubusercontent.com/29379385/98324273-d6e9f600-1fb9-11eb-8f2d-3a296a0758fe.png)
  - [x] Each rseg with have unique values for species loss +/- X for each of the model scenarios
  - [x] May also want to POST the `elf$stats` for the elf used for calculating the richness change (would allow us to filter on n or rsquared if we choose to exclude those with very few points or poor fitting elfs etc.)
- Notes:
  - [x] Scotty: is the `consumptive_use_frac` metric cumulative or only representing withdrawals within a particular rseg? - Ans: Yes, it is cumulative
  - [x] Only calculate species loss for those rsegs with mean annual flow below the breakpoint? - Ans: Generate all, but only use those below the breakpoint for analysis 
  - [x] Are we interested in mean annual flow only, or a particular month, or all months? (is `consumptive_use_frac` strictly an annual metric?) - Ans: Focus only on MAF for now

### Step 2.1 - Running Analysis Script From DEQ server
#### Step 2.1.0 - Queries & Commands
see:https://github.com/HARPgroup/vahydro/wiki/2020-Model-Runs-&-Analysis-Scripts

Input arguments: run_rseg_elfgen.R [rseg model pid] [run id] [huc level] [datasource]

cd /var/www/R Rscript /opt/model/om/R/summarize/run_rseg_elfgen.R 4711908 11 huc8 VAHydro-EDAS

Create tmp file of rseg model pids

deq2.bse.vt.edu ssh dbase2 psql drupal.dh03

copy ( SELECT model.pid AS model_pid FROM dh_feature AS rseg INNER JOIN dh_properties as model ON ( model.featureid = rseg.hydroid
AND model.varid IN (SELECT hydroid FROM dh_variabledefinition WHERE varkey = 'om_water_model_node' AND propcode = 'vahydro-1.0') )
WHERE rseg.bundle = 'watershed' AND rseg.ftype = 'vahydro' GROUP BY model.pid ORDER BY model.pid ASC ) to '/tmp/rseg_models.txt' ; COPY 686 Ctrl + D to exit psql console Ctrl + D to close dbase2 connection

Copy file to vahydro directory

cd /var/www/html/files/vahydro sftp dbase2 get /tmp/rseg_models.txt bye

##########

TEST RUNS

##########

Use the "head" command to take the first 5 records for a test run

head -n 5 rseg_models.txt > rseg_models_test5.txt

Bash for loop

cd /var/www/R

!/bin/sh

for r in 11 13; do echo $r for i in cat /var/www/html/files/vahydro/rseg_models_test5.txt; do echo $i Rscript /opt/model/om/R/summarize/run_rseg_elfgen.R $i $r huc8 VAHydro-EDAS
done done

##########

FULL RUN

##########

Bash for loop

cd /var/www/R

Execute using 'nohup' to hide output

!/bin/sh

for r in 13; do echo $r for i in cat /var/www/html/files/vahydro/rseg_models.txt; do echo $i nohup Rscript /opt/model/om/R/summarize/run_rseg_elfgen.R $i $r huc8 VAHydro-EDAS
done done

tail -f /home/jkleiner/nohup.out


#### Step 2.1.1 - Batches run
- 686 Riverseg Models
  - huc 8 
  - VAHydro-EDAS
  - runid
    - [x] 11 (2020) - _Kicked off 2:00pm 11/30/20_
    - [ ] 12 (2030)
    - [x] 13 (2040) - _Kicked off 11:00am 12/1/20 -  Finished 5:00pm_
    - [x] 18 (Exempt Users) - _Kicked off 5:00pm 12/1/20_
    - [x] ~~17 (Dry CC)~~ - _CC uses that 10yr sim period, which introduces odd consumptive_use_frac results_
    - [x] ~~19 (Med CC)~~ - _CC uses that 10yr sim period, which introduces odd consumptive_use_frac results_
    - [x] ~~20 (Wet CC)~~ - _CC uses that 10yr sim period, which introduces odd consumptive_use_frac results_

#### Step 2.1.2 - Rseg re-runs

!/bin/sh

for r in 11; do echo $r for i in cat /var/www/html/files/vahydro/rseg_models_reruns_list_12.4.20.txt; do echo $i nohup Rscript /opt/model/om/R/summarize/run_rseg_elfgen.R $i $r huc8 VAHydro-EDAS
done done


- 12 Riverseg Models
  - huc 8 
  - VAHydro-EDAS
  - runid
    - [x] 11 (2020) - _Kicked off 9:00am 12/4/20_
    - [x] 13 (2040)
    - [x] 18 (Exempt Users)

### Step 3.0 - Options for presenting this information 
#### Step 3.1 - Tabular:
- [ ] Simplest option, present data in tables with rseg `consumptive_use_frac` and estimates for percent richness change
#### Step 3.2 - Maps: 
- [ ] One option: 
![VA_runid_13_consumptive_use_frac_map](https://user-images.githubusercontent.com/29379385/95114451-fb7c5500-0711-11eb-9738-b649076641ec.png)

#### Step 3.3 - Calculate an "eco deficit" 
- [ ] Calculate an "eco deficit" as the difference in area under the black line (natural ELF) and the red line (post-withdrawal)
  - This is a "next step" item, not as critical right now as getting the data stored in VAHydro
  - Rob's rugged example image below: 
![image](https://user-images.githubusercontent.com/4571170/94030993-69c52d00-fd8c-11ea-90e5-bc1ba3ff1418.png) 
- [ ] Try generating a similar image by retrieving the rseg species loss +/- X from VAHydro 
  - Use a single minor basin for testing, PS?

### Step 4.0 - _Future Item_ - Push data to VAHydro (Specific to VWP modeling scenarios)
-  Thinking this script should run automatically for VWP model scenarios, and at a minimum should produce the items below
   1. Plot image(s) - post images to the run scenario as dh_image_file(s), similar to waterSupplyElement.R
   2. Post properties 
      - post a property for the potential taxa loss determined using the richness_change() function above. (This is the primary piece of info we're interested in, what is the potential loss of richness at this location for the withdrawal they are proposing?)
         - Use confidence interval from the elf plot to determine the +/- for the richness loss value?
           - elf_dnt_05, elf_dnt_95, elf_dnt_50; elf dnt = change in nt
      - additional properties we may want: 
         - post a property for the maximum richness value at this location as seen in the example plot above (where the dotted line intersects the y-axis)
         - post a property for the `[MEAN FLOW @ INTAKE]`
         - whats the maximum observed richness value in this watershed?
jdkleiner commented 3 years ago

Re-run elfgen analysis for rsegs that have been rerun since 11/30/20

# Create tmp file of rseg model pids
copy (
SELECT model.pid AS model_pid   
  FROM dh_feature AS rseg
  INNER JOIN dh_properties as model
    ON (
      model.featureid = rseg.hydroid  
      AND model.varid IN (SELECT hydroid 
                          FROM dh_variabledefinition 
                          WHERE varkey = 'om_water_model_node' AND
                                propcode = 'vahydro-1.0')
    )    

  INNER JOIN dh_properties as run_object
    ON (
      run_object.featureid = model.pid  
      AND run_object.varid IN (SELECT hydroid 
                          FROM dh_variabledefinition 
                          WHERE varkey = 'om_scenario' AND
                                run_object.propname = 'runid_11')
    )      

  WHERE rseg.bundle = 'watershed' AND
        rseg.ftype = 'vahydro' AND
        TO_CHAR(TO_TIMESTAMP(run_object.modified), 'MM/DD/YYYY HH24:MI:SS') > '11/30/2020'
  GROUP BY model.pid, run_object.pid
  ORDER BY TO_CHAR(TO_TIMESTAMP(run_object.modified), 'MM/DD/YYYY HH24:MI:SS') DESC  
) to '/tmp/rseg_models_runid_11_12.16.20.txt'
;
COPY 41
Ctrl + D to exit psql console
Ctrl + D to close dbase2 connection

# Copy file to vahydro directory
cd /var/www/html/files/vahydro
sftp dbase2
get /tmp/rseg_models_runid_11_12.16.20.txt
bye

###################
# runid_11 FULL RUN
###################
#Bash for loop
cd /var/www/R  

#Execute using 'nohup' to hide output

#!/bin/sh
for r in 11; do
  echo $r
  for i in `cat /var/www/html/files/vahydro/rseg_models_runid_11_12.16.20.txt`; do 
    echo $i
    nohup Rscript /opt/model/om/R/summarize/run_rseg_elfgen.R $i $r huc8 VAHydro-EDAS    
  done
done

tail -f /home/jkleiner/nohup.out