In [34]:
#This notebook can post-process a QUALTX output file
#and plot water quality profiles along all the modeled reaches
#ET 20160514
#Enter the location of the QUALTX output file
#This program can parse both urls and local paths
Model = 'Doe Branch'
#Model = 'Rhodairs Gully'
if Model == 'Doe Branch':
#Doe Branch
QUALTXoutfile = r'C:\Users\Ernest\Ernest_Sandbox\TurboQUALTX\TurboQUALTX\QUALTX_outfiles\DOEOCT4.OUT'
permit_nums = ['10698-003'] #permit number
DOstd = 5 #mg/L #DO standard
length_km = 7.5 #distance downstream from permit discharger location (for mapping)
if Model == 'Rhodairs Gully':
#QUALTXoutfile = r'https://raw.githubusercontent.com/Harefoot/TurboQUALTX/master/RHOTRUNC.OUT'
QUALTXoutfile = r'C:\Users\Ernest\Ernest_Sandbox\TurboQUALTX\TurboQUALTX\QUALTX_outfiles\RHOTRUNC.OUT'
permit_nums = ['10838-003']
DOstd = 4 #mg/L #DO standard
length_km = 10 #distance downstream from permit discharger location (for mapping)
In [35]:
#Import python packages
import pandas as pd
pd.options.mode.chained_assignment = None # default='warn'
import os
import sys
#sys.path.insert(0,r'M:\Library\Python\Packages')
import ET_Utils.QUALTX_Utils
import numpy as np
import bokeh
from bokeh.plotting import figure, show
from bokeh.models import HoverTool, BoxSelectTool
from bokeh.io import output_notebook
output_notebook()
In [36]:
#Parse the QUALTX output file
Scenario = QUALTXoutfile
#WQC_of_interest = ['DO_MG/L','BOD_MG/L','NH3_MG/L','NO3+2_MG/L','PHOS_MG/L']
WQC_of_interest = ['DO_MG/L','BOD_MG/L','NH3_MG/L']#,'NO3+2_MG/L','PHOS_MG/L']
StreamNames, Rdf, Hdf, AllWQdf = ET_Utils.QUALTX_Utils.Process_QUALTX(QUALTXoutfile,
Scenario,DOstd = DOstd,
WQC_of_interest = WQC_of_interest,
loc = 1,plot_pdf = 0)
In [37]:
#Plot a map
import folium
import numpy as np
from folium.element import IFrame
from ET_Utils.WATERS_Utils import point_indexing_client
from ET_Utils.WATERS_Utils import event_indexing_client
#bbox = [-87.40, 24.25, -74.70, 36.70]
tiles = ('http://services.arcgisonline.com/arcgis/rest/'
'services/Ocean/World_Ocean_Base/MapServer/tile/{z}/{y}/{x}')
attr = 'ESRI'
tiles_esriRelief = ('http://server.arcgisonline.com/ArcGIS/rest/'
'services/World_Shaded_Relief/MapServer/tile/{z}/{y}/{x}')
attr_esriRelief ='ESRI shaded relief'
tiles_stamen='Stamen Terrain'
attr_stamen = 'Stamen Terrain'
outfall_csv = r'C:\Users\Ernest\Ernest_Sandbox\TurboQUALTX\TurboQUALTX\GIS\tceq_wastewater_outfalls.csv'
odf = pd.read_csv(outfall_csv)
odf = odf[odf['PERMIT_NUM'].isin(permit_nums)].reset_index()
#df.isin([1, 3, 12, 'a']
#bbox = [-106.6460, 25.8371, -93.5083, 36.5007]
bbox = [min(odf['X']),min(odf['Y']),max(odf['X']),max(odf['Y'])]
location = np.array(bbox).reshape(2, 2).mean(axis=0).tolist()[::-1]
mapa = folium.Map(location=location, zoom_start=13,max_zoom=13)
#mapa = folium.Map(location=location, zoom_start=13,
# tiles=tiles_esriRelief, attr=attr_esriRelief,max_zoom = 13)
for i in odf.index:#[0:10]:
popup ='GIS number: '+odf['GIS_NUMBER'][i]+'\n'+' Permittee: '+odf['PERMITTEE'][i]
folium.Marker([odf['Y'][i], odf['X'][i]], popup=popup).add_to(mapa)
#print location
#mock_events = event_indexing_client.trace_downstream(odf['Y'][i], odf['X'][i], length_km)
lon, lat = (odf['X'][i], odf['Y'][i])
#lon, lat = (-94.024207547, 29.99177026900003)
#lon, lat = (-96.7970489739,32.77673479287502)
print lon,lat
try:
#mock_events = event_indexing_client.trace_downstream(lat, lon, length_km)
#folium.GeoJson(mock_events).add_to(mapa)
path = point_indexing_client.trace_to_stream(lat, lon)
folium.GeoJson(path).add_to(mapa)
end_point = path['features'][1]
lon, lat = end_point['geometry']['coordinates']
length_km = length_km - path['features'][0]['properties']['length_km']
print('Length (km) of trace along NHD streams: {}'.format(length_km))
events = event_indexing_client.trace_downstream(lat, lon, length_km)
folium.GeoJson(events).add_to(mapa)
except TypeError:
print "No NHDPlus flowlines to trace. Note to self: plot circle with radius = length km centered around discharge point"
#folium.Marker([45.3288, -121.6625], popup='Mt. Hood Meadows').add_to(mapa)
#folium.Marker([45.3311, -121.7113], popup='Timberline Lodge').add_to(mapa)
mapa
Out[37]:
In [38]:
#Set plotting parameters
plot_width=900
plot_height=600
y_label = 'Concentration (mg/L)'
x_label = 'Distance (km)'
#Set line width
lws = np.zeros(len(WQC_of_interest))+2
#Set line colors
plot_colors = ['blue','black','green','cyan','orange']
#Set line symbols
plot_symbols = ['-','-','-','-','-']
#Set marker sizes
markersizes = np.zeros(len(WQC_of_interest))+1
#Set label sizes
x_labelsize = 15
y_labelsize = 15
In [39]:
# Loop through each stream and plot the WQ profile plot
TOOLS="hover,box_zoom,wheel_zoom,pan,reset"
for StreamName in StreamNames:
WQdf = AllWQdf[AllWQdf['Stream']==StreamName].copy()
SRdf = Rdf[Rdf['StreamName']==StreamName].reset_index()
x_range = [max(WQdf['ENDING_DIST'])*1.025,min(WQdf['ENDING_DIST'])]
#Determine range of y-axis
maxy = []
for tempWQC in WQC_of_interest:
maxy.append(max(WQdf[tempWQC]))
y_range = [0,round(max(maxy))]#*1.5)]
p = figure(title=StreamName, tools=TOOLS,x_range = x_range,y_range=y_range,plot_width=plot_width, plot_height=plot_height)
if StreamName == StreamNames[0]:
hover = p.select(dict(type=HoverTool))
hover.tooltips = [
# format text in tooltip
("Conc ", "$y mg/L"),
]
p.xaxis.axis_label = x_label
p.yaxis.axis_label = y_label
for k in range(0,len(WQC_of_interest)):
# add a line renderer
p.line(np.asarray(WQdf['ENDING_DIST']),np.asarray(WQdf[WQC_of_interest[k]]),
legend=WQC_of_interest[k], line_width= lws[k],
color = plot_colors[k],name = WQC_of_interest[k])
#Overplot reach start and end points of each reach
if len(SRdf) > 0:
Points1 = list(SRdf['BEGIN NAME'])
Rkms1 = list(SRdf['BEGIN REACH KM'])
if len(SRdf) > 1:
Points2 = [SRdf['END NAME'][SRdf.index[-1]]]
Rkms2 = [SRdf['END REACH KM'][SRdf.index[-1]]]
else:
Points2 = list(SRdf['END NAME'])
Rkms2 = list(SRdf['END REACH KM'])
Points = Points1 + Points2
#Rkms = Rkms1.append(Rkms2)
Rkms = Rkms1+Rkms2
#p1.text([65,65,65],[65,65,65], text=[ str(i) for i in x], alpha=0.5, text_font_size="5pt", text_baseline="middle",
#text_align="center")
for i in range(0,len(Points)):
#print Points[i]
p.line([Rkms[i],Rkms[i]],[-999,999],'..',line_dash = 'dotted',color = 'grey')
p.text([Rkms[i]],[np.mean(y_range)],text = [Points[i]],text_align = 'center',
angle = 3.14159265/2,text_font_size = '9pt',text_color = 'grey')
if DOstd != np.nan:
DOstds = [DOstd,DOstd-0.2]
vas = ['bottom','top']
for i in range(len(DOstds)):
p.line([-999,999],[DOstds[i],DOstds[i]],'..',line_dash = 'dotted',color = 'red')
p.text([np.mean(x_range)],[DOstds[i]],["DO std = "+"{:4.1f}".format(DOstds[i])+" mg/L"],
text_baseline = vas[i],text_font_size = '9pt',text_color = 'red')
pass
p.xgrid.grid_line_color = None
p.ygrid.grid_line_color = None
show(p)