This notebook creates a figure for the SciPy 2014 John Hunter Excellence in Plotting contest (https://conference.scipy.org/scipy2014/plotting_contest/). A version of this figure appears in a 2013 Environmental Chemistry article (E. Couzo, H.E. Jeffries, W. Vizuete. Houston's rapid ozone increases: preconditions and geographic origins. Environ. Chem. 2013, 10, 260-268. doi:10.1071/EN13040)
Houston, you have an ozone problem.¶
Ozone (O$_3$) is a harmful air pollutant that forms when nitrogen oxides and hydrocarbons react in the presence of sunlight. Houston, Texas, has some of the worst O$_3$ pollution in the U.S. because (1) Houston often suffers from stagnant meteorological conditions that lead to O$_3$ accumulation, (2) nitrogen oxide emissions are large due to the vast number of automobiles, and (3) the petrochemical industry emits highly reactive compounds, which accelerate the rate of O$_3$ production. The last reason - industrial petrochemical emissions - has been the focus of intense study and regulatory action.
Many of Houston's highest O$_3$ concentrations are characterized by rapid increases in concentrations of at least 40 ppb in one hour, or 60 ppb in two hours. These rapid increases are called non-typical O$_3$ changes (NTOCs) and are unique to Houston. It is believed that many NTOCs are caused by sudden releases of reactive hydrocarbons from petrochemical facilities. Using publicly available air quality monitoring data (http://www.tceq.state.tx.us/) from 2000-2011, this analysis maps the geographic origins of NTOCs at six air quality monitors across Houston.
This figure shows that NTOCs are measured at air quality monitors when they are downwind of Houston's petrochemical industry, the heart of which is marked by the black star in the figure. Each of the six radar plots is centered on the location of an air quality monitoring state. The theta-axis shows wind direction, and the r-axis shows wind speed. Each red marker shows the wind speed and direction when a NTOC was measured, i.e., when the air quality monitor detected an extremely rapid O$_3$ increase. Although there is some scatter in each plot, there is strong preference for a narrow range of wind directions. The red markers are clustered in the direction of Houston's petrochemical industry (black star), which means NTOCs tend to occur at monitors when they are downwind of the pterochemical industry. The blue markers in each radar plot show the wind speed and direction that were measured when daily peak O$_3$ was measured. These data do not cluster in a particualr direction meaning that peak O$_3$ comes from a diversity of directions.
This figure links geographically the NTOC phenomenon to Houston's petrochemical industry. It shows that, while peak O$_3$ has no preference for wind direction, NTOCs tend to happen only when downwind of the petrochemical industry. The figure also shows that NTOCs are measured when wind speeds are relatively low (usually under 10 km/hr, always under 20 km/hr). If wind speeds are too high, industrial petrochemical emissions are transported out of the region before they have a chance to form O$_3$.
%pylab inline
import numpy as np
import csv
import pylab
from mpl_toolkits.basemap import Basemap
#get a list of all NTOCs by monitor
NTOC_filepath = 'observed_NTOC_any_crit_2000_2011.txt'
NTOC_file = open(NTOC_filepath, 'rU')
NTOC_reader = csv.reader(NTOC_file)
NTOC_data = [line for line in NTOC_reader][2:-1]
def TCEQwd_to_POLARwd(tceq_wd):
# converts WD from TCEQ coordinates (O degree North, clockwise rotation) to
# more standard coordinates (0 degree East, counterclockwise rotation)
polar_wd = (360 - tceq_wd + 90) % 360
return polar_wd
def get_peak_O3_wswd(mon='mon'):
O3_filepath = '{0}-O3.txt'.format(mon)
WS_filepath = '{0}-ws.txt'.format(mon)
WD_filepath = '{0}-wd.txt'.format(mon)
O3_file = open(O3_filepath, 'rU')
WS_file = open(WS_filepath, 'rU')
WD_file = open(WD_filepath, 'rU')
O3_reader = csv.reader(O3_file)
WS_reader = csv.reader(WS_file)
WD_reader = csv.reader(WD_file)
O3_datalines = [line for line in O3_reader][6:]
WS_datalines = [line for line in WS_reader][5:]
WD_datalines = [line for line in WD_reader][5:]
peak_O3_wswd = [[], []]
for day in range(len(WS_datalines)):
O3_line = O3_datalines[day]
O3 = O3_line[1:]
O3_array = np.array(O3, dtype='float')
hr = O3_array.argmax() #time of max 1-hr O3
WS_line = WS_datalines[day]
WS = WS_line[1:]
WS_array = np.array(WS, dtype='float') * 1.60934 # convert miles to kilometers
WD_line = WD_datalines[day]
WD = WD_line[1:]
WD_array = np.array(WD, dtype='float')
WD_array_polar = TCEQwd_to_POLARwd(WD_array)
peak_O3_wswd[0].append(WS_array[hr])
peak_O3_wswd[1].append(WD_array_polar[hr])
return np.array(peak_O3_wswd)
def get_NTOC_wswd(mon='mon'):
mon_NTOC_list = [ntoc for ntoc in NTOC_data if ' '+mon in ntoc] # ' '+mon required because of leading space in front of mon name in NTOC_file
WS_filepath = '{0}-ws.txt'.format(mon)
WD_filepath = '{0}-wd.txt'.format(mon)
WS_file = open(WS_filepath, 'rU')
WD_file = open(WD_filepath, 'rU')
WS_reader = csv.reader(WS_file)
WD_reader = csv.reader(WD_file)
WS_datalines = [line for line in WS_reader][5:]
WD_datalines = [line for line in WD_reader][5:]
NTOC_wswd = [[], []]
for day in mon_NTOC_list:
ntocdatestr = day[0]
ntochr = int(day[1])
ntocWS_lines = [wsline for wsline in WS_datalines if ntocdatestr in wsline][0]
ntocWS = ntocWS_lines[1:]
ntocWS_array = np.array(ntocWS, dtype='float') * 1.60934 # convert miles to kilometers
ntocWD_lines = [wdline for wdline in WD_datalines if ntocdatestr in wdline][0]
ntocWD = ntocWD_lines[1:]
ntocWD_array = np.array(ntocWD, dtype='float')
ntocWD_array_polar = TCEQwd_to_POLARwd(ntocWD_array)
NTOC_wswd[0].append(ntocWS_array[ntochr])
NTOC_wswd[1].append(ntocWD_array_polar[ntochr])
return np.array(NTOC_wswd)
WALV_peakO3_wswd = get_peak_O3_wswd('WALV')
HALC_peakO3_wswd = get_peak_O3_wswd('HALC')
CLIN_peakO3_wswd = get_peak_O3_wswd('CLIN')
HCQA_peakO3_wswd = get_peak_O3_wswd('HCQA')
MACP_peakO3_wswd = get_peak_O3_wswd('MACP')
TXCT_peakO3_wswd = get_peak_O3_wswd('TXCT')
WALV_NTOC_wswd = get_NTOC_wswd('WALV')
HALC_NTOC_wswd = get_NTOC_wswd('HALC')
CLIN_NTOC_wswd = get_NTOC_wswd('CLIN')
HCQA_NTOC_wswd = get_NTOC_wswd('HCQA')
MACP_NTOC_wswd = get_NTOC_wswd('MACP')
TXCT_NTOC_wswd = get_NTOC_wswd('TXCT')
map_dict = dict(projection = 'lcc', lat_1 = 45., lat_2 = 33., lat_0 = 40., lon_0 = -97.,
llcrnrlon=-95.6, urcrnrlon=-94.85, llcrnrlat=29.35, urcrnrlat=30.05)
fig = pylab.figure(figsize=(12,8))
fig.suptitle(' Wind speed/direction at time of peak O$_3$ and NTOC', fontsize=20)
map = Basemap(**map_dict)
c = map.drawcounties(zorder=1)
c.set_facecolor('0.9')
map.drawmapboundary(fill_color='aqua')
map.readshapefile('Highways', 'hwy', color='black', linewidth=0.75)
map.drawparallels(np.arange(0,360,0.5), labels=[1,0,0,1], dashes=[1,5])
map.drawmeridians(np.arange(0,360,0.5), labels=[1,0,0,1], dashes=[1,5])
#get lat/lon of Sam Houston Tollway bridge (proxy for center of Houston's petrochemical industry) from Google maps
bridge_lat = 29.736236
bridge_lon = -95.145945
s = map.scatter(bridge_lon, bridge_lat, latlon=True, marker='*', color='black', edgecolor='black', s=1000, zorder=2)
WALVfig = fig.add_axes([.59,.58,0.15,0.15], polar=True, frameon=False)
WALVfig.scatter(WALV_peakO3_wswd[1] * np.pi/180, WALV_peakO3_wswd[0], color='blue', marker='x', s=75, lw=1)
WALVfig.scatter(WALV_NTOC_wswd[1] * np.pi/180, WALV_NTOC_wswd[0], color='red', marker='x', s=75, lw=1)
HALCfig = fig.add_axes([.39,.66,0.15,0.15], polar=True, frameon=False)
HALCfig.scatter(HALC_peakO3_wswd[1] * np.pi/180, HALC_peakO3_wswd[0], color='blue', marker='x', s=75, lw=1)
HALCfig.scatter(HALC_NTOC_wswd[1] * np.pi/180, HALC_NTOC_wswd[0], color='red', marker='x', s=75, lw=1)
CLINfig = fig.add_axes([.42,.47,0.15,0.15], polar=True, frameon=False)
CLINfig.scatter(CLIN_peakO3_wswd[1] * np.pi/180, CLIN_peakO3_wswd[0], color='blue', marker='x', s=75, lw=1)
CLINfig.scatter(CLIN_NTOC_wswd[1] * np.pi/180, CLIN_NTOC_wswd[0], color='red', marker='x', s=75, lw=1)
HCQAfig = fig.add_axes([.28,.35,0.15,0.15], polar=True, frameon=False)
HCQAfig.scatter(HCQA_peakO3_wswd[1] * np.pi/180, HCQA_peakO3_wswd[0], color='blue', marker='x', s=75, lw=1)
HCQAfig.scatter(HCQA_NTOC_wswd[1] * np.pi/180, HCQA_NTOC_wswd[0], color='red', marker='x', s=75, lw=1)
MACPfig = fig.add_axes([.33,.23,0.15,0.15], polar=True, frameon=False)
MACPfig.scatter(MACP_peakO3_wswd[1] * np.pi/180, MACP_peakO3_wswd[0], color='blue', marker='x', s=75, lw=1)
MACPfig.scatter(MACP_NTOC_wswd[1] * np.pi/180, MACP_NTOC_wswd[0], color='red', marker='x', s=75, lw=1)
TXCTfig = fig.add_axes([.61,.13,0.15,0.15], polar=True, frameon=False)
TXCTfig.scatter(TXCT_peakO3_wswd[1] * np.pi/180, TXCT_peakO3_wswd[0], color='blue', marker='x', s=75, lw=1)
TXCTfig.scatter(TXCT_NTOC_wswd[1] * np.pi/180, TXCT_NTOC_wswd[0], color='red', marker='x', s=75, lw=1)
for radar in [WALVfig, HALCfig, CLINfig, HCQAfig, MACPfig, TXCTfig]:
theta_axis = radar.get_xaxis()
r_axis = radar.get_yaxis()
theta_axis.set_ticks(np.arange(0,360,90) *np.pi/180 +2*np.pi/180)
theta_axis.set_ticklabels([])
r_axis.set_ticks([0,10,20,30])
r_axis.set_ticklabels([''])
radar.set_ylim(0,35)
theta_axis.grid(b=True, which='both', color='black', lw=1.5, ls='-')
r_axis.grid(b=True, which='both', color='black', lw=.75, ls='-')
key = fig.add_axes([.75,0.4,.25,.25], polar=True, frameon=False)
key.set_title('Key\n\n', fontsize=16)
ktheta_axis = key.get_xaxis()
kr_axis = key.get_yaxis()
ktheta_axis.set_ticks(np.arange(0,360,90) *np.pi/180 +2*np.pi/180)
ktheta_axis.set_ticklabels(['E','N','W','S'])
kr_axis.set_ticks([0,10,20,30])
kr_axis.set_ticklabels([''])
key.set_ylim(0,35)
ktheta_axis.grid(b=True, which='both', color='black', lw=1.5, ls='-')
kr_axis.grid(b=True, which='both', color='black', lw=.75, ls='-')
for r in ['0','10','20','30']:
key.text(60 *np.pi/180, int(r), r+' km/hr', size=8, ha='center', va='center', bbox=dict(facecolor='white'))
key.scatter([], [], color='blue', marker='x', s=75, lw=1, label='Peak O$_3$')
key.scatter([], [], color='red', marker='x', s=75, lw=1, label='NTOC')
key.legend(scatterpoints=1, loc=(.18,-.5))
fig.savefig('Couzo_fig_scipy2014_plotting_contest.pdf')