Querying data from the USGS National Water Information System#
The USGS National Water Information System (NWIS) provides access to surface-water data at approximately 8,500 sites. The data are served online—most in near realtime—to meet many diverse needs. For additional information regarding NWIS, visit the USGS website.
from HydroGenerate.hydropower_potential import calculate_hp_potential
import pandas as pd
import matplotlib.pyplot as plt
Use streamflow data from a USGS stream gauge.#
The USGS has a Python package (dataretrieval) designed to simplify the process of loading hydrologic data into the Python environment. dataretrieval can be pip installed into a Python environment. The example code shows how to donwload data using dataretrieval and format it for hydropower potential evaluation in HydroGenerate. The example used, and additional information, is avlaiable as a HydroShare Resource. Information about dataretrieval is avaliable directly on the USGS GitHUb Repository.
# Install dataretrieval, if it doesn't exist in Python environment, using:
# !pip install dataretrieval
# import data retrieval
from dataretrieval import nwis
# Set the parameters needed to retrieve data
siteNumber = "10109000" # LOGAN RIVER ABOVE STATE DAM, NEAR LOGAN, UT
parameterCode = "00060" # Discharge
startDate = "2010-01-01"
endDate = "2020-12-31"
# Retrieve the data
dailyStreamflow = nwis.get_dv(sites=siteNumber, parameterCd=parameterCode, start=startDate, end=endDate)
# convert to pandas dataframe
flow = pd.DataFrame(dailyStreamflow[0])
flow.head()
| site_no | 00060_Mean | 00060_Mean_cd | |
|---|---|---|---|
| datetime | |||
| 2010-01-01 00:00:00+00:00 | 10109000 | 108.0 | A |
| 2010-01-02 00:00:00+00:00 | 10109000 | 107.0 | A |
| 2010-01-03 00:00:00+00:00 | 10109000 | 106.0 | A |
| 2010-01-04 00:00:00+00:00 | 10109000 | 102.0 | A |
| 2010-01-05 00:00:00+00:00 | 10109000 | 105.0 | A |
# Use the streamflow data from a USGS stream gauge - Cont.
# The flow data for this example was created in the previous step.
head = 100 # ft
power = None
penstock_length = 120 # ft
hp_type = 'Diversion'
flow_column = '00060_Mean' # name of the column containing the flow data
hp = calculate_hp_potential(flow= flow, rated_power= power, head= head,
penstock_headloss_calculation= True,
hydropower_type= hp_type, penstock_length= penstock_length,
flow_column= flow_column, annual_caclulation= True)
# Explore output
print('Design flow (cfs):', hp.design_flow)
print('Head_loss at design flow (ft):', round(hp.penstock_design_headloss, 2))
print('Turbine type:', hp.turbine_type)
print('Rated Power (Kw):', round(hp.rated_power, 2))
print('Net head (ft):', round(hp.net_head, 2))
print('Generator Efficiency:',hp.generator_efficiency)
print('Head Loss method:',hp.penstock_headloss_method)
print('Penstock length (ft):', hp.penstock_length)
print('Penstock diameter (ft):', round(hp.penstock_diameter,2))
print('\nPandas dataframe output: \n', hp.dataframe_output)
print('Annual output: \n', hp.annual_dataframe_output)
Design flow (cfs): 200.0
Head_loss at design flow (ft): 8.94
Turbine type: Crossflow
Rated Power (Kw): 1193.76
Net head (ft): 91.06
Generator Efficiency: 0.98
Head Loss method: Darcy-Weisbach
Penstock length (ft): 120.0
Penstock diameter (ft): 3.31
Pandas dataframe output:
site_no 00060_Mean 00060_Mean_cd power_kW \
datetime
2010-01-01 00:00:00+00:00 10109000 108.0 A 502.584462
2010-01-02 00:00:00+00:00 10109000 107.0 A 456.688144
2010-01-03 00:00:00+00:00 10109000 106.0 A 398.442172
2010-01-04 00:00:00+00:00 10109000 102.0 A 0.000000
2010-01-05 00:00:00+00:00 10109000 105.0 A 324.023491
... ... ... ... ...
2020-12-27 00:00:00+00:00 10109000 91.0 A 0.000000
2020-12-28 00:00:00+00:00 10109000 85.9 A 0.000000
2020-12-29 00:00:00+00:00 10109000 83.0 A 0.000000
2020-12-30 00:00:00+00:00 10109000 82.3 A 0.000000
2020-12-31 00:00:00+00:00 10109000 90.1 A 0.000000
turbine_flow_cfs efficiency energy_kWh
datetime
2010-01-01 00:00:00+00:00 108.0 0.575850 NaN
2010-01-02 00:00:00+00:00 107.0 0.527893 10960.515465
2010-01-03 00:00:00+00:00 106.0 0.464683 9562.612123
2010-01-04 00:00:00+00:00 102.0 0.000000 0.000000
2010-01-05 00:00:00+00:00 105.0 0.381307 7776.563783
... ... ... ...
2020-12-27 00:00:00+00:00 91.0 0.000000 0.000000
2020-12-28 00:00:00+00:00 85.9 0.000000 0.000000
2020-12-29 00:00:00+00:00 83.0 0.000000 0.000000
2020-12-30 00:00:00+00:00 82.3 0.000000 0.000000
2020-12-31 00:00:00+00:00 90.1 0.000000 0.000000
[4018 rows x 7 columns]
Annual output:
annual_turbinedvolume_ft3 mean_annual_effienciency \
datetime
2010 1319.865870 0.398728
2011 1835.676204 0.687777
2012 1442.086842 0.559378
2013 1087.183892 0.213385
2014 1290.223844 0.401570
2015 1208.405282 0.285440
2016 1374.823115 0.478596
2017 1917.789260 0.772958
2018 1416.658356 0.546933
2019 1458.516249 0.568326
2020 1333.259716 0.402846
total_annual_energy_KWh revenue_M$ capacity_factor
datetime
2010 4.266764e+06 0.248326 0.408017
2011 8.650045e+06 0.503433 0.827177
2012 6.069435e+06 0.353241 0.580401
2013 2.392200e+06 0.139226 0.228759
2014 4.319024e+06 0.251367 0.413015
2015 3.219217e+06 0.187358 0.307844
2016 5.142542e+06 0.299296 0.491765
2017 9.640457e+06 0.561075 0.921887
2018 5.757931e+06 0.335112 0.550613
2019 6.103963e+06 0.355251 0.583703
2020 4.347384e+06 0.253018 0.415727
# Plot results
# Columns: discharge_cfs site_id power_kW efficiency energy_kWh
plt.rcParams['figure.figsize'] = [14, 7]
df = hp.dataframe_output.copy()
fig, ax1 = plt.subplots()
color_plot = 'tab:red'
ax1.set_xlabel('Days')
ax1.set_ylabel('Flow rate (cfs)', color=color_plot)
ax1.plot(df['turbine_flow_cfs'], label="Flow rate", color=color_plot)
ax1.tick_params(axis='y', labelcolor=color_plot)
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
color_plot2 = 'tab:blue'
ax2.set_ylabel('Power (kW)', color=color_plot2) # we already handled the x-label with ax1
ax2.plot(df['power_kW'], label="Power", color=color_plot2)
ax2.tick_params(axis='y', labelcolor=color_plot2)
ax1.grid(True, axis='both', color='k',linestyle='--',alpha=0.4)
plt.title("Yearly flow data from USGS and potential power")
fig.tight_layout() # otherwise the right y-label is slightly clipped
#plt.savefig(os.path.join('..','fig','usgs_twin_falls_flow_power.jpg'))
plt.show()
