#!/usr/bin/env python # Pollution rose plots by Zuber Farooqui, Ph.D. # zuber.farooqui [at] gmail.com # 9/17/2019 # 9/8/2020 The script is updated to include calm wind speed. User can input the # calm wind speed cut-point. Including cut-point and lesser will not be used in the plot. # import numpy as np import argparse import pandas as pd import pandas.io.sql as psql # from sqlalchemy import create_engine # getting data and arguments parser = argparse.ArgumentParser("Providing arguments for windrose plots") parser.add_argument('--ifile', dest='ifile', required=True, default='None', help='Input data file .csv format') parser.add_argument('--site', dest='site', required=True, default='site', help='CARB site number') parser.add_argument('--bdate', dest='bdate', required=True, default='bdate', help='Begin date 1-1-2018 format') parser.add_argument('--edate', dest='edate', required=True, default='edate', help='End date 12-31-2018 format') parser.add_argument('--wind-dir', dest='winddir', default='WD', help='Wind direction variable') parser.add_argument('--binwidth', dest='binwidth', type=float, default=45., help='Direction bin widths in degrees.') parser.add_argument('--wscut', dest='wscut', type=float, default=0.0, help='Wind speed cut point, excluded.') parser.add_argument('--max-pct', dest='maxpct', type=float, default=50., help='Maximum percent on plot') parser.add_argument('--fromnorth', dest='fromnorth', action='store_true', default=False, help='Bins should start from north (i.e., not bound north)') parser.add_argument('--bounds', dest='bounds', type=lambda x: np.array(eval(x)), default=np.array([10, 20, 30, 50, 100, 120, 140, np.inf]), help='Boundaries for pollutant') parser.add_argument('--outpath', dest='outpath', type=str, required=True, default=50., help='Destination to plot files') args = parser.parse_args() parser.epilog = """ Example: $ windrose_mpl.py --ifile IFILE --site 3818 --bdate bdate --edate edate --max-pct=50 --wind-dir wd --wind-speed poll --bounds="[0.5,1,2,3.5,4.5,6.5,np.inf]" --title obs --fromnorth --binwidth 22.5 --outpath outdir Description: Create windrose for site site 3818 between the dates badte and edate and on where wind-direction is name wd, and windspeed is named poll. The script will calculate calm wind speed and associated conc. at cut point. """ wd_data = pd.read_csv(args.ifile, index_col=False, engine='python') def calmws(wd_data): calm_data = wd_data.loc[wd_data['ws'] <= args.wscut] calm_per = round(calm_data.ws.size * 100 / float(wd_data.ws.size),1) return calm_per # windrose plot using matplotlib and numpy def windrose(wd_data,args): # calmness condition calculations calm_per = calmws(wd_data) # data excluding calm conditions wd_data = wd_data.loc[wd_data['ws'] > args.wscut] wd= wd_data['wd'] ws= wd_data['ws'] # from matplotlib import pyplot as plt plt.rcParams['text.usetex'] = False # Bar width width = args.binwidth halfwidth = width / 2. # Setting plot size and others fig = plt.figure(figsize=(6, 6)) ax = fig.add_axes([0.1, 0.07, 0.8, 0.8], polar=True) ax.set_theta_offset(np.radians(90)) ax.set_theta_direction(-1) # getting bounds ubs = args.bounds[1:] lbs = args.bounds[:-1] # ading colors ubcs = plt.cm.jet(np.arange(len(ubs), dtype='f') / ubs.size) color_dict = dict(zip(ubs, ubcs)) # n=masking data based on pollutant mask = np.ma.getmaskarray(ws) | np.ma.getmaskarray( wd) | (ws < args.bounds[0]) | (ws > args.bounds[-1]) ws = np.ma.masked_where(mask, ws).compressed() wd = np.ma.masked_where(mask, wd).compressed() # adding title dl=pd.to_datetime(args.bdate) dlabel=dl.year # fig.suptitle('%s' % (dlabel), size=20) if args.fromnorth: # NNE, WNE, E, SSE, SSW, WSW, WNW, NNW thetas = (np.radians((wd + halfwidth)// width * width ).astype('float64') % (2 * np.pi)).astype('float64') else: # N, NE, E, SE, S, SW, W, NW thetas = (np.radians((wd.astype('int64') + halfwidth) // width * width - halfwidth) % (2 * np.pi)).astype('float64') # uthetas = np.unique(thetas) nws = float(ws.size) for theta in uthetas: tws = ws[thetas == theta] for ubi, ub in enumerate(ubs[::-1]): ubc = color_dict[ub] pct = (tws < ub).sum() / nws * 100. ax.bar(theta, pct, width=np.radians(width), bottom=0.0, color=ubc) # labels=[] for i,l in enumerate(ubs): if (i == 0): labels.append('< %s' % l) elif (i == 'inf'): labels.append('> %s' % ubs[i-1]) else: labels.append('%s - %s' % (ubs[i-1],l)) ##to use below commented section (lines:117-124). First change line 68 and 69 to ##fig = plt.figure(figsize=(6, 7)) ##ax = fig.add_axes([0.1, 0.2, 0.8, 0.7], polar=True) #rectangles = [plt.Rectangle((0, 0), 1, 1, color=ubc) for ubc in ubcs] #fig.legend(rectangles, labels, ncol=3, # bbox_to_anchor=(.5, 0.025), loc='lower center') #ax.text(0.27, -0.09, 'PM2.5, $\mu$g/$m^3$', horizontalalignment='center', # verticalalignment='center', transform=ax.transAxes,size=12) #datelabel = '%s_%s' % (args.bdate,args.edate) #ax.text(0.95, 1.05, datelabel, horizontalalignment='center', style='italic', # verticalalignment='center', transform=ax.transAxes,size=10) # writing calm ws percentage cws_text = " Calm Winds = %s%%" % calm_per ax.text(0.15, 0.0, cws_text, horizontalalignment='center', style='italic', verticalalignment='center', transform=ax.transAxes,size=10) # ticks = np.linspace(0, args.maxpct, 6)[1:] labels = ['%s%%' % i for i in ticks] plt.yticks(ticks, labels) ax.set_rmax(args.maxpct) ax.set_clip_on(False) figpath = '%s/WRose_%s_%s_%s.png' % (args.outpath,args.site,args.bdate,args.edate) fig.savefig(figpath, transparen=True) print('Saved fig', figpath) plt.close(fig) windrose(wd_data,args)