Add mainstems attribute for sierra test

Sierra test considered all USGS gage locations to be mainstems even though many actually occurred with tributaries. This resulted in unrealistic comparisons as incorrect gages were assigned to mainstems segments. This feature branch identifies gages that are on mainstems via attribute field.

- Modifies usgs_gage_crosswalk.py to filter out gages from the usgs_gages.gpkg layer such that for a "MS" run, only consider gages that contain rating curve information (via curve attribute) and are also mainstems gages (via mainstems attribute).
- Modifies usgs_gage_crosswalk.py to filter out gages from the usgs_gages.gpkg layer such that for a "FR" run, only consider gages that contain rating curve information (via curve attribute) and are not mainstems gages (via mainstems attribute).
- Modifies how mainstems segments are determined by using the nwm_flows_ms.gpkg as a lookup to determine if the NWM segment specified by WRDS for a gage site is a mainstems gage.
- Adds a mainstem attribute field to usgs_gages.gpkg that indicates whether a gage is located on a mainstems river.
- Adds NWM_FLOWS_MS variable to the .env and .env.template files.
- Adds the extent argument specified by user when running fim_run.sh to usgs_gage_crosswalk.py.

CHANGELOG.md

@@ -1,5 +1,23 @@
 All notable changes to this project will be documented in this file.
 We follow the [Semantic Versioning 2.0.0](http://semver.org/) format.
+
+## v3.0.15.1 - 2021-04-13 - [PR #355](https://github.com/NOAA-OWP/cahaba/pull/355)
+
+Sierra test considered all USGS gage locations to be mainstems even though many actually occurred with tributaries. This resulted in unrealistic comparisons as incorrect gages were assigned to mainstems segments. This feature branch identifies gages that are on mainstems via attribute field. 
+## Changes
+
+- Modifies `usgs_gage_crosswalk.py` to filter out gages from the `usgs_gages.gpkg` layer such that for a "MS" run, only consider gages that contain rating curve information (via `curve` attribute) and are also mainstems gages (via `mainstems` attribute). 
+- Modifies `usgs_gage_crosswalk.py` to filter out gages from the `usgs_gages.gpkg` layer such that for a "FR" run, only consider gages that contain rating curve information (via `curve` attribute) and are not mainstems gages (via `mainstems` attribute). 
+- Modifies how mainstems segments are determined by using the `nwm_flows_ms.gpkg` as a lookup to determine if the NWM segment specified by WRDS for a gage site is a mainstems gage. 
+
+## Additions
+
+- Adds a `mainstem` attribute field to `usgs_gages.gpkg` that indicates whether a gage is located on a mainstems river. 
+- Adds `NWM_FLOWS_MS` variable to the `.env` and `.env.template` files.
+- Adds the `extent` argument specified by user when running `fim_run.sh` to `usgs_gage_crosswalk.py`.
+
+<br/><br/>
+
 ## v3.0.15.0 - 2021-04-08 - [PR #340](https://github.com/NOAA-OWP/cahaba/pull/340)
 
 Implementing a prototype technique to estimate the missing bathymetric component in the HAND-derived synthetic rating curves. The new Bathymetric Adjusted Rating Curve (BARC) function is built within the `fim_run.sh` workflow and will ingest bankfull geometry estimates provided by the user to modify the cross section area used in the synthetic rating curve generation.

src/run_by_unit.sh

@@ -436,7 +436,7 @@ Tcount
 echo -e $startDiv"USGS Crosswalk $hucNumber"$stopDiv
 date -u
 Tstart
-$srcDir/usgs_gage_crosswalk.py -gages $inputDataDir/usgs_gages/usgs_gages.gpkg -dem $outputHucDataDir/dem_meters.tif -flows $outputHucDataDir/demDerived_reaches_split_filtered_addedAttributes_crosswalked.gpkg -cat $outputHucDataDir/gw_catchments_reaches_filtered_addedAttributes_crosswalked.gpkg -wbd $outputHucDataDir/wbd_buffered.gpkg -dem_adj $dem_thalwegCond -outtable $outputHucDataDir/usgs_elev_table.csv
+$srcDir/usgs_gage_crosswalk.py -gages $inputDataDir/usgs_gages/usgs_gages.gpkg -dem $outputHucDataDir/dem_meters.tif -flows $outputHucDataDir/demDerived_reaches_split_filtered_addedAttributes_crosswalked.gpkg -cat $outputHucDataDir/gw_catchments_reaches_filtered_addedAttributes_crosswalked.gpkg -wbd $outputHucDataDir/wbd_buffered.gpkg -dem_adj $dem_thalwegCond -outtable $outputHucDataDir/usgs_elev_table.csv -e $extent
 Tcount
 
 ## CLEANUP OUTPUTS ##

src/usgs_gage_crosswalk.py

@@ -32,7 +32,7 @@
 '''
 
 
-def crosswalk_usgs_gage(usgs_gages_filename,dem_filename,input_flows_filename,input_catchment_filename,wbd_buffer_filename,dem_adj_filename,output_table_filename):
+def crosswalk_usgs_gage(usgs_gages_filename,dem_filename,input_flows_filename,input_catchment_filename,wbd_buffer_filename,dem_adj_filename,output_table_filename,extent):
 
     wbd_buffer = gpd.read_file(wbd_buffer_filename)
     usgs_gages = gpd.read_file(usgs_gages_filename, mask=wbd_buffer)
@@ -41,8 +41,12 @@ def crosswalk_usgs_gage(usgs_gages_filename,dem_filename,input_flows_filename,in
     input_catchment = gpd.read_file(input_catchment_filename)
     dem_adj = rasterio.open(dem_adj_filename,'r')
 
-    #Query out usgs_gages that don't have rating curve data
-    usgs_gages = usgs_gages.query('curve == "yes"')
+    #MS extent use gages that are mainstem
+    if extent == "MS":
+        usgs_gages = usgs_gages.query('curve == "yes" & mainstem == "yes"')
+    #FR extent use gages that are not mainstem
+    if extent == "FR":
+        usgs_gages = usgs_gages.query('curve == "yes" & mainstem == "no"')
 
     if input_flows.HydroID.dtype != 'int': input_flows.HydroID = input_flows.HydroID.astype(int)
 
@@ -113,7 +117,7 @@ def crosswalk_usgs_gage(usgs_gages_filename,dem_filename,input_flows_filename,in
     parser.add_argument('-wbd','--wbd-buffer-filename', help='WBD buffer', required=True)
     parser.add_argument('-dem_adj','--dem-adj-filename', help='Thalweg adjusted DEM', required=True)
     parser.add_argument('-outtable','--output-table-filename', help='Table to append data', required=True)
-
+    parser.add_argument('-e', '--extent', help="extent configuration entered by user when running fim_run.sh", required = True)
     args = vars(parser.parse_args())
 
     usgs_gages_filename = args['usgs_gages_filename']
@@ -123,5 +127,6 @@ def crosswalk_usgs_gage(usgs_gages_filename,dem_filename,input_flows_filename,in
     wbd_buffer_filename = args['wbd_buffer_filename']
     dem_adj_filename = args['dem_adj_filename']
     output_table_filename = args['output_table_filename']
+    extent = args['extent']
 
-    crosswalk_usgs_gage(usgs_gages_filename,dem_filename,input_flows_filename,input_catchment_filename,wbd_buffer_filename, dem_adj_filename,output_table_filename)
+    crosswalk_usgs_gage(usgs_gages_filename,dem_filename,input_flows_filename,input_catchment_filename,wbd_buffer_filename, dem_adj_filename,output_table_filename, extent)

tools/.env.template

@@ -1,3 +1,4 @@
 API_BASE_URL=
 EVALUATED_SITES_CSV=
 WBD_LAYER=
+NWM_FLOWS_MS=

tools/rating_curve_get_usgs_curves.py

@@ -23,6 +23,8 @@
 load_dotenv()
 API_BASE_URL = os.getenv("API_BASE_URL")
 WBD_LAYER = os.getenv("WBD_LAYER")
+EVALUATED_SITES_CSV = os.getenv("EVALUATED_SITES_CSV")
+NWM_FLOWS_MS = os.getenv("NWM_FLOWS_MS")
 
 def get_all_active_usgs_sites():
     '''
@@ -245,7 +247,16 @@ def usgs_rating_to_elev(list_of_gage_sites, workspace=False, sleep_time = 1.0):
     sites_with_data = pd.DataFrame({'location_id':all_rating_curves['location_id'].unique(),'curve':'yes'})
     acceptable_sites_gdf = acceptable_sites_gdf.merge(sites_with_data, on = 'location_id', how = 'left')
     acceptable_sites_gdf.fillna({'curve':'no'},inplace = True)    
-
+    #Add mainstems attribute to acceptable sites
+    print('Attributing mainstems sites')
+    #Import mainstems segments used in run_by_unit.sh
+    ms_df = gpd.read_file(NWM_FLOWS_MS)
+    ms_segs = ms_df.ID.astype(str).to_list()
+    #Populate mainstems attribute field
+    acceptable_sites_gdf['mainstem'] = 'no'
+    acceptable_sites_gdf.loc[acceptable_sites_gdf.eval('feature_id in @ms_segs'),'mainstem'] = 'yes' 
+
+
     #If workspace is specified, write data to file.
     if workspace:
         #Write rating curve dataframe to file