-
Notifications
You must be signed in to change notification settings - Fork 26
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
[1pt] Tool to compare synthetic rating curve with benchmark rating cu…
…rve (sierra test) (#332) * `usgs_gage_crosswalk.py`: generates `usgs_elev_table.csv` in run_by_unit.py with elevation at USGS gages * `rating_curve_comparison.py`: post-processing script to plot and calculate metrics between FIM and USGS rating curve data * updates `aggregate_fim_outputs.py` call argument in `fim_run.sh` from 4 jobs to 6 jobs (optimizing API performance) * reroutes median elevation data from `add_crosswalk.py` and `rem.py` to new file (depreciating `hand_ref_elev_table.csv`) *- adds new files to `viz_whitelist` in `output_cleanup.py`
- Loading branch information
1 parent
4df6d1d
commit 44b6333
Showing
8 changed files
with
590 additions
and
31 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,124 @@ | ||
| #!/usr/bin/env python3 | ||
|
|
||
| import os | ||
| import geopandas as gpd | ||
| import pandas as pd | ||
| import numpy as np | ||
| import rasterio | ||
| import argparse | ||
| import pygeos | ||
| from shapely.wkb import dumps, loads | ||
| import warnings | ||
| warnings.simplefilter(action='ignore', category=FutureWarning) | ||
|
|
||
| ''' Get elevation at adjusted USGS gages locations | ||
| Parameters | ||
| ---------- | ||
| usgs_gages_filename : str | ||
| File name of USGS stations layer. | ||
| dem_filename : str | ||
| File name of original DEM. | ||
| input_flows_filename : str | ||
| File name of FIM streams layer. | ||
| input_catchment_filename : str | ||
| File name of FIM catchment layer. | ||
| wbd_buffer_filename : str | ||
| File name of buffered wbd. | ||
| dem_adj_filename : str | ||
| File name of thalweg adjusted DEM. | ||
| output_table_filename : str | ||
| File name of output table. | ||
| ''' | ||
|
|
||
|
|
||
| def crosswalk_usgs_gage(usgs_gages_filename,dem_filename,input_flows_filename,input_catchment_filename,wbd_buffer_filename,dem_adj_filename,output_table_filename): | ||
|
|
||
| wbd_buffer = gpd.read_file(wbd_buffer_filename) | ||
| usgs_gages = gpd.read_file(usgs_gages_filename, mask=wbd_buffer) | ||
| dem_m = rasterio.open(dem_filename,'r') | ||
| input_flows = gpd.read_file(input_flows_filename) | ||
| input_catchment = gpd.read_file(input_catchment_filename) | ||
| dem_adj = rasterio.open(dem_adj_filename,'r') | ||
|
|
||
| if input_flows.HydroID.dtype != 'int': input_flows.HydroID = input_flows.HydroID.astype(int) | ||
|
|
||
| # Identify closest HydroID | ||
| closest_catchment = gpd.sjoin(usgs_gages, input_catchment, how='left', op='within').reset_index(drop=True) | ||
| closest_hydro_id = closest_catchment.filter(items=['site_no','HydroID','min_thal_elev','med_thal_elev','max_thal_elev', 'order_']) | ||
| closest_hydro_id = closest_hydro_id.dropna() | ||
|
|
||
| # Get USGS gages that are within catchment boundaries | ||
| usgs_gages = usgs_gages.loc[usgs_gages.site_no.isin(list(closest_hydro_id.site_no))] | ||
|
|
||
| columns = ['location_id','HydroID','dem_elevation','dem_adj_elevation','min_thal_elev', 'med_thal_elev','max_thal_elev','str_order'] | ||
| gage_data = [] | ||
|
|
||
| # Move USGS gage to stream | ||
| for index, gage in usgs_gages.iterrows(): | ||
|
|
||
| # Get stream attributes | ||
| hydro_id = closest_hydro_id.loc[closest_hydro_id.site_no==gage.site_no].HydroID.item() | ||
| str_order = str(int(closest_hydro_id.loc[closest_hydro_id.site_no==gage.site_no].order_.item())) | ||
| min_thal_elev = round(closest_hydro_id.loc[closest_hydro_id.site_no==gage.site_no].min_thal_elev.item(),2) | ||
| med_thal_elev = round(closest_hydro_id.loc[closest_hydro_id.site_no==gage.site_no].med_thal_elev.item(),2) | ||
| max_thal_elev = round(closest_hydro_id.loc[closest_hydro_id.site_no==gage.site_no].max_thal_elev.item(),2) | ||
|
|
||
| # Convert headwater point geometries to WKB representation | ||
| wkb_gages = dumps(gage.geometry) | ||
|
|
||
| # Create pygeos headwater point geometries from WKB representation | ||
| gage_bin_geom = pygeos.io.from_wkb(wkb_gages) | ||
|
|
||
| # Closest segment to headwater | ||
| closest_stream = input_flows.loc[input_flows.HydroID==hydro_id] | ||
| wkb_closest_stream = dumps(closest_stream.geometry.item()) | ||
| stream_bin_geom = pygeos.io.from_wkb(wkb_closest_stream) | ||
|
|
||
| # Linear reference headwater to closest stream segment | ||
| gage_distance_to_line = pygeos.linear.line_locate_point(stream_bin_geom, gage_bin_geom) | ||
| referenced_gage = pygeos.linear.line_interpolate_point(stream_bin_geom, gage_distance_to_line) | ||
|
|
||
| # Convert geometries to wkb representation | ||
| bin_referenced_gage = pygeos.io.to_wkb(referenced_gage) | ||
|
|
||
| # Convert to shapely geometries | ||
| shply_referenced_gage = loads(bin_referenced_gage) | ||
|
|
||
| # Sample rasters at adjusted gage | ||
| dem_m_elev = round(list(rasterio.sample.sample_gen(dem_m,shply_referenced_gage.coords))[0].item(),2) | ||
| dem_adj_elev = round(list(rasterio.sample.sample_gen(dem_adj,shply_referenced_gage.coords))[0].item(),2) | ||
|
|
||
| # Append dem_m_elev, dem_adj_elev, hydro_id, and gage number to table | ||
| site_elevations = [str(gage.site_no), str(hydro_id), dem_m_elev, dem_adj_elev, min_thal_elev, med_thal_elev, max_thal_elev,str(str_order)] | ||
| gage_data.append(site_elevations) | ||
|
|
||
|
|
||
| elev_table = pd.DataFrame(gage_data, columns=columns) | ||
|
|
||
| if not elev_table.empty: | ||
| elev_table.to_csv(output_table_filename,index=False) | ||
|
|
||
|
|
||
| if __name__ == '__main__': | ||
|
|
||
| parser = argparse.ArgumentParser(description='Crosswalk USGS sites to HydroID and get elevations') | ||
| parser.add_argument('-gages','--usgs-gages-filename', help='USGS gages', required=True) | ||
| parser.add_argument('-dem','--dem-filename',help='DEM',required=True) | ||
| parser.add_argument('-flows','--input-flows-filename', help='DEM derived streams', required=True) | ||
| parser.add_argument('-cat','--input-catchment-filename', help='DEM derived catchments', required=True) | ||
| 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) | ||
|
|
||
| args = vars(parser.parse_args()) | ||
|
|
||
| usgs_gages_filename = args['usgs_gages_filename'] | ||
| dem_filename = args['dem_filename'] | ||
| input_flows_filename = args['input_flows_filename'] | ||
| input_catchment_filename = args['input_catchment_filename'] | ||
| wbd_buffer_filename = args['wbd_buffer_filename'] | ||
| dem_adj_filename = args['dem_adj_filename'] | ||
| output_table_filename = args['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) |
Oops, something went wrong.