Python datatypes.speed函数代码示例

def plotter(fdict):
    """ Go """
    import matplotlib
    matplotlib.use('agg')
    import matplotlib.pyplot as plt
    ASOS = psycopg2.connect(database='asos', host='iemdb', user='nobody')
    cursor = ASOS.cursor(cursor_factory=psycopg2.extras.DictCursor)

    station = fdict.get('station', 'AMW')
    units = fdict.get('units', 'mph')
    network = fdict.get('network', 'IA_ASOS')

    nt = NetworkTable(network)

    cursor.execute("""
    SELECT extract(doy from valid), sknt * 0.514, drct from alldata
    where station = %s and sknt >= 0 and drct >= 0
    """, (station, ))

    uwnd = np.zeros((366,), 'f')
    vwnd = np.zeros((366,), 'f')
    cnt = np.zeros((366,), 'f')
    for row in cursor:
        u, v = uv(row[1], row[2])
        uwnd[int(row[0]) - 1] += u
        vwnd[int(row[0]) - 1] += v
        cnt[int(row[0]) - 1] += 1

    u = speed(uwnd / cnt, 'MPS').value(units.upper())
    v = speed(vwnd / cnt, 'mps').value(units.upper())

    df = pd.DataFrame(dict(u=pd.Series(u),
                           v=pd.Series(v),
                           day_of_year=pd.Series(np.arange(1, 366))))

    (fig, ax) = plt.subplots(1, 1)

    ax.plot(np.arange(1, 366), smooth(u[:-1], 14, 'hamming'), color='r',
            label='u, West(+) : East(-) component')
    ax.plot(np.arange(1, 366), smooth(v[:-1], 14, 'hamming'), color='b',
            label='v, South(+) : North(-) component')
    ax.set_xticks([1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 365])
    ax.set_xticklabels(calendar.month_abbr[1:])
    ax.legend(ncol=2, fontsize=11, loc=(0., -0.15))
    ax.grid(True)
    ax.set_xlim(0, 366)
    ax.set_title(("[%s] %s Daily Average Component Wind Speed\n"
                  "[%s-%s] 14 day smooth filter applied, %.0f obs found"
                  "") % (station, nt.sts[station]['name'],
                         nt.sts[station]['archive_begin'].year,
                         datetime.datetime.now().year, np.sum(cnt)))
    ax.set_ylabel("Average Wind Speed %s" % (PDICT.get(units), ))

    box = ax.get_position()
    ax.set_position([box.x0, box.y0 + box.height * 0.1, box.width,
                     box.height * 0.9])

    return fig, df

def main():
    """Go Main"""
    pgconn = get_dbconn('asos')
    df = read_sql("""
    SELECT valid - '1 hour'::interval as valid,
    drct, sknt, gust_sknt, pres1, tmpf, dwpf
    from t2018_1minute
    where station = %s and valid >= '2018-06-14 08:30' and
    valid <= '2018-06-14 10:15' ORDER by valid ASC
    """, pgconn, params=('PHP', ), index_col='valid')
    xticks = []
    xticklabels = []
    for valid in df.index.values:
        if pd.to_datetime(valid).minute % 15 == 0:
            xticks.append(valid)
            ts = pd.to_datetime(valid) - datetime.timedelta(hours=5)
            xticklabels.append(ts.strftime("%-H:%M\n%p"))

    fig = plt.figure(figsize=(8, 9))
    ax = fig.add_axes([0.1, 0.55, 0.75, 0.35])
    ax.plot(df.index.values, df['tmpf'], label='Air Temp')
    ax.plot(df.index.values, df['dwpf'], label='Dew Point')
    ax.legend()
    ax.grid(True)
    ax.set_ylabel("Temperature $^\circ$F")
    ax.set_xticks(xticks)
    ax.set_xticklabels(xticklabels)
    ax.set_title(("Philip, SD (KPHP) ASOS 1 Minute Interval Data for 14 Jun 2018\n"
                  "Heat Burst Event, data missing in NCEI files 8:02 to 8:10 AM"))

    ax = fig.add_axes([0.1, 0.08, 0.75, 0.35])

    ax.bar(df.index.values, speed(df['gust_sknt'], 'KT').value('MPH'),
           width=1/1440., color='red')
    ax.bar(df.index.values, speed(df['sknt'], 'KT').value('MPH'),
           width=1/1440., color='tan')
    ax.set_ylabel("Wind Speed (tan) & Gust (red) [mph]")
    ax.grid(True, zorder=5)
    ax.set_ylim(0, 60)

    ax2 = ax.twinx()
    ax2.plot(df.index.values, pressure(df['pres1'], 'IN').value('MB'),
             color='g', lw=2)
    ax2.set_ylabel("Air Pressure [hPa]", color='green')
    ax2.set_xticks(xticks)
    ax2.set_xticklabels(xticklabels)
    ax.set_xlabel("14 June 2018 MDT")
    ax2.set_ylim(923, 926)
    ax2.set_yticks(np.arange(923, 926.1, 0.5))
    # ax2.set_zorder(ax.get_zorder()-1)
    # ax2.set_ylim(0, 360)
    # ax2.set_yticks(range(0, 361, 45))
    # ax2.set_yticklabels(['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW', 'N'])

    fig.savefig('test.png')

def read_excel(siteid, fn):
    df = pd.read_excel(fn, skiprows=[1, ])
    newcols = {}
    for k in df.columns:
        newcols[k] = XREF.get(k, k)
    df.rename(columns=newcols, inplace=True)
    df['valid'] = df['valid'] + datetime.timedelta(hours=TZREF[siteid])
    # do some conversions
    print("ALERT: doing windspeed unit conv")
    df['windspeed_mps'] = speed(df['windspeed_mps'].values, 'KMH').value('MPS')
    print("ALERT: doing windgustunit conv")
    df['windgust_mps'] = speed(df['windgust_mps'].values, 'KMH').value('MPS')
    return df

def test_uv():
    """ Test calculation of uv wind components """
    speed = datatypes.speed([10, ], 'KT')
    mydir = datatypes.direction([0, ], 'DEG')
    u, v = meteorology.uv(speed, mydir)
    assert u.value("KT") == 0.
    assert v.value("KT") == -10.

    speed = datatypes.speed([10, 20, 15], 'KT')
    mydir = datatypes.direction([90, 180, 135], 'DEG')
    u, v = meteorology.uv(speed, mydir)
    assert u.value("KT")[0] == -10
    assert v.value("KT")[1] == 20.
    assert abs(v.value("KT")[2] - 10.6) < 0.1

def rabbit_tracks(row):
    """Generate a rabbit track for this attr"""
    res = ""
    if row['sknt'] is None or row['sknt'] <= 5 or row['drct'] is None:
        return res
    # 5 carrots at six minutes to get 30 minutes?
    lat0 = row['lat']
    lon0 = row['lon']
    drct = row['drct']
    sknt = row['sknt']
    x0, y0 = P3857(lon0, lat0)
    smps = speed(sknt, 'KTS').value('MPS')
    angle = dir2ccwrot(drct)
    rotation = (drct + 180) % 360
    rad = math.radians(angle)
    x = x0 + math.cos(rad) * smps * SECONDS
    y = y0 + math.sin(rad) * smps * SECONDS
    # Draw white line out 30 minutes
    lons, lats = P3857(x, y, inverse=True)
    res += ("Line: 1, 0, \"Cell [%s]\"\n"
            "%.4f, %.4f\n"
            "%.4f, %.4f\n"
            "END:\n") % (row['storm_id'],
                         lat0, lon0, lats[-1], lons[-1])
    for i in range(3):
        res += ("Icon: %.4f,%.4f,%.0f,1,10,\"+%.0f min\"\n"
                ) % (lats[i], lons[i], rotation, (i+1)*15)
    return res

def plotter(fdict):
    """ Go """
    import matplotlib
    matplotlib.use('agg')
    import matplotlib.pyplot as plt
    import matplotlib.patheffects as PathEffects
    pgconn = psycopg2.connect(database='iem', host='iemdb', user='nobody')
    cursor = pgconn.cursor(cursor_factory=psycopg2.extras.DictCursor)

    station = fdict.get('zstation', 'AMW')
    network = fdict.get('network', 'IA_ASOS')
    units = fdict.get('units', 'MPH').upper()
    if units not in PDICT:
        units = 'MPH'
    year = int(fdict.get('year', datetime.datetime.now().year))
    month = int(fdict.get('month', datetime.datetime.now().month))
    sts = datetime.date(year, month, 1)
    ets = (sts + datetime.timedelta(days=35)).replace(day=1)
    nt = NetworkTable(network)

    cursor.execute("""
      SELECT day, avg_sknt, vector_avg_drct from summary s JOIN stations t
      ON (t.iemid = s.iemid) WHERE t.id = %s and t.network = %s and
      s.day >= %s and s.day < %s ORDER by day ASC
    """, (station, network, sts, ets))
    days = []
    drct = []
    sknt = []
    for row in cursor:
        if row[1] is None:
            continue
        days.append(row[0].day)
        drct.append(row[2])
        sknt.append(row[1])
    if len(sknt) == 0:
        return "ERROR: No Data Found"
    df = pd.DataFrame(dict(day=pd.Series(days),
                           drct=pd.Series(drct),
                           sknt=pd.Series(sknt)))
    sknt = speed(np.array(sknt), 'KT').value(units)
    (fig, ax) = plt.subplots(1, 1)
    ax.bar(np.array(days)-0.4, sknt, ec='green', fc='green')
    pos = max([min(sknt) / 2.0, 0.5])
    for d, _, r in zip(days, sknt, drct):
        draw_line(plt, d, max(sknt)+0.5, (270. - r) / 180. * np.pi)
        txt = ax.text(d, pos, drct2text(r), ha='center', rotation=90,
                      color='white', va='center')
        txt.set_path_effects([PathEffects.withStroke(linewidth=2,
                                                     foreground="k")])
    ax.grid(True, zorder=11)
    ax.set_title(("%s [%s]\n%s Daily Average Wind Speed and Direction"
                  ) % (nt.sts[station]['name'], station,
                       sts.strftime("%b %Y")))
    ax.set_xlim(0.5, max(days)+0.5)
    ax.set_ylim(top=max(sknt)+2)

    ax.set_ylabel("Average Wind Speed [%s]" % (PDICT.get(units),))

    return fig, df

def uv(speed, direction):
    """
    Compute the u and v components of the wind
    @param wind speed in whatever units
    @param dir wind direction with zero as north
    @return u and v components
    """
    if (not isinstance(speed, dt.speed) or
            not isinstance(direction, dt.direction)):
        raise InvalidArguments(("uv() needs speed and direction "
                                "objects as args"))
    # Get radian units
    rad = direction.value("RAD")
    if rad is None or speed.value() is None:
        return None, None
    u = (0 - speed.value()) * np.sin(rad)
    v = (0 - speed.value()) * np.cos(rad)
    return (dt.speed(u, speed.get_units()), dt.speed(v, speed.get_units()))

def plotter(fdict):
    """ Go """
    import matplotlib
    matplotlib.use('agg')
    from pyiem.plot import calendar_plot
    pgconn = psycopg2.connect(database='iem', host='iemdb', user='nobody')
    cursor = pgconn.cursor(cursor_factory=psycopg2.extras.DictCursor)

    ctx = get_autoplot_context(fdict, get_description())
    station = ctx['station']
    varname = ctx['var']
    network = ctx['network']
    sdate = ctx['sdate']
    edate = ctx['edate']

    nt = NetworkTable(network)

    # Get Climatology
    cdf = read_sql("""SELECT to_char(valid, 'mmdd') as sday, high, low,
    precip from ncdc_climate81 WHERE station = %s
    """, psycopg2.connect(database='coop', host='iemdb', user='nobody'),
                   params=(nt.sts[station]['ncdc81'],), index_col='sday')

    cursor.execute("""
    SELECT day, max_tmpf, min_tmpf, max_dwpf, min_dwpf,
    pday, coalesce(avg_sknt, 0) as avg_sknt from summary s JOIN stations t
    on (t.iemid = s.iemid) WHERE s.day >= %s and s.day <= %s and
    t.id = %s and t.network = %s ORDER by day ASC
    """, (sdate, edate, station, network))
    rows = []
    data = {}
    for row in cursor:
        hd = row['max_tmpf'] - cdf.at[row[0].strftime("%m%d"), 'high']
        ld = row['min_tmpf'] - cdf.at[row[0].strftime("%m%d"), 'low']
        rows.append(dict(day=row['day'], max_tmpf=row['max_tmpf'],
                         avg_smph=speed(row['avg_sknt'], 'KT').value('MPH'),
                         min_dwpf=row['min_dwpf'], max_dwpf=row['max_dwpf'],
                         high_departure=hd, low_departure=ld,
                         min_tmpf=row['min_tmpf'], pday=row['pday']))
        data[row[0]] = {'val': safe(rows[-1], varname)}
        if varname == 'high_departure':
            data[row[0]]['color'] = 'b' if hd < 0 else 'r'
        elif varname == 'low_departure':
            data[row[0]]['color'] = 'b' if ld < 0 else 'r'
    df = pd.DataFrame(rows)

    title = ('[%s] %s Daily %s\n%s thru %s'
             ) % (station, nt.sts[station]['name'],
                  PDICT.get(varname), sdate.strftime("%-d %b %Y"),
                  edate.strftime("%-d %b %Y"))

    fig = calendar_plot(sdate, edate, data,
                        title=title)
    return fig, df

def googlesheet(siteid, sheetkey):
    """Harvest a google sheet, please"""
    rows = []
    config = util.get_config()
    sheets = util.get_sheetsclient(config, "td")
    f = sheets.spreadsheets().get(spreadsheetId=sheetkey, includeGridData=True)
    j = util.exponential_backoff(f.execute)
    for sheet in j['sheets']:
        # sheet_title = sheet['properties']['title']
        for griddata in sheet['data']:
            for row, rowdata in enumerate(griddata['rowData']):
                if 'values' not in rowdata:  # empty sheet
                    continue
                if row == 1:  # skip units
                    continue
                if row == 0:
                    header = []
                    for col, celldata in enumerate(rowdata['values']):
                        header.append(celldata['formattedValue'])
                    continue
                data = {}
                for col, celldata in enumerate(rowdata['values']):
                    data[header[col]] = fmt(celldata.get('formattedValue'))
                rows.append(data)
    df = pd.DataFrame(rows)
    print("googlesheet has columns: %s" % (repr(df.columns.values),))
    newcols = {}
    for k in df.columns:
        newcols[k] = XREF.get(k, k)
    df.rename(columns=newcols, inplace=True)
    df['valid'] = pd.to_datetime(df['valid'], errors='raise',
                                 format='%m/%d/%y %H:%M')
    df['valid'] = df['valid'] + datetime.timedelta(hours=TZREF[siteid])

    # do some conversions
    print("ALERT: doing windspeed unit conv")
    df['windspeed_mps'] = speed(df['windspeed_mps'].values, 'KMH').value('MPS')
    print("ALERT: doing windgustunit conv")
    df['windgust_mps'] = speed(df['windgust_mps'].values, 'KMH').value('MPS')
    return df

def make_rwis(i, j, initts, oldncout):
    """ Generate spinup file """
    i = i - IOFFSET
    j = j - JOFFSET

    o = open('rwis.xml', 'w')
    o.write("""<?xml version="1.0"?>
<observation>
 <header>
  <filetype>rwis-observation</filetype>
  <version>1.0</version>
  <road-station>oaa</road-station>
  </header>
  <measure-list>""")
    if oldncout is None:
        fake_rwis(o, initts)
        return

    ts0 = find_initts(oldncout)
    # at Air Temp in C
    tmpc = dt.temperature(oldncout.variables['tmpk'][:, i, j], 'K').value('C')
    # td Dew point in C
    dwpc = dt.temperature(oldncout.variables['dwpk'][:, i, j], 'K').value('C')
    # pi presence of precipitation 0: No -- 1: Yes
    # ws wind speed in km / hr
    ws = dt.speed(oldncout.variables['wmps'][:, i, j], 'MPS').value('KMH')
    # sc condition code  1=DryCond 2=Wet 3=Ice 4=MixWaterSnow
    #                    5=dew 6=Meltsnow 7=Frost 8=Ice
    # Was set to 33 for SSI ?
    icond = oldncout.variables['icond'][:, i, j]
    # st road surface temp
    bridgec = dt.temperature(
        oldncout.variables['bdeckt'][:, i, j], 'K').value('C')
    # sst sub surface temp
    subsfc = dt.temperature(
        oldncout.variables['subsfct'][:, i, j], 'K').value('C')
    t1 = initts + datetime.timedelta(hours=12)
    for tstep in range(4, len(oldncout.dimensions['time']), 4):
        ts = ts0 + datetime.timedelta(
                                minutes=int(oldncout.variables['time'][tstep]))
        if ts > t1:
            break
        o.write("""<measure><observation-time>%s</observation-time>
<at>%.2f</at><td>%.2f</td><pi>0</pi><ws>%.2f</ws><sc>%s</sc><st>%.2f</st>
<sst>%.2f</sst></measure>
      """ % (ts.strftime("%Y-%m-%dT%H:%MZ"), tmpc[tstep], dwpc[tstep],
             ws[tstep], icond[tstep], bridgec[tstep], subsfc[tstep]))

    o.write("</measure-list></observation>")
    o.close()

def plotter(fdict):
    """ Go """
    import matplotlib
    matplotlib.use('agg')
    import matplotlib.pyplot as plt
    pgconn = psycopg2.connect(database='asos', host='iemdb', user='nobody')
    ctx = get_autoplot_context(fdict, get_description())

    station = ctx['zstation']
    network = ctx['network']
    units = ctx['units']

    nt = NetworkTable(network)

    df = read_sql("""
        select date_trunc('hour', valid) as ts, avg(sknt) as sknt,
        max(drct) as drct from alldata
        WHERE station = %s and sknt is not null and drct is not null
        GROUP by ts
        """, pgconn, params=(station, ), parse_dates=('ts',),
                  index_col=None)
    sknt = speed(df['sknt'].values, 'KT')
    drct = direction(df['drct'].values, 'DEG')
    df['u'], df['v'] = [x.value('MPS') for x in meteorology.uv(sknt, drct)]
    df['month'] = df['ts'].dt.month
    grp = df[['month', 'u', 'v', 'sknt']].groupby('month').mean()
    grp['u_%s' % (units,)] = speed(grp['u'].values, 'KT').value(units.upper())
    grp['v_%s' % (units,)] = speed(grp['u'].values, 'KT').value(units.upper())
    grp['sped_%s' % (units,)] = speed(grp['sknt'].values,
                                      'KT').value(units.upper())
    drct = meteorology.drct(speed(grp['u'].values, 'KT'),
                            speed(grp['v'].values, 'KT'))
    grp['drct'] = drct.value('DEG')
    maxval = grp['sped_%s' % (units,)].max()
    (fig, ax) = plt.subplots(1, 1)
    ax.barh(grp.index.values, grp['sped_%s' % (units,)].values,
            align='center')
    ax.set_xlabel("Average Wind Speed [%s]" % (UNITS[units],))
    ax.set_yticks(grp.index.values)
    ax.set_yticklabels(calendar.month_abbr[1:])
    ax.grid(True)
    ax.set_xlim(0, maxval * 1.2)
    for mon, row in grp.iterrows():
        ax.text(maxval * 1.1, mon, drct2text(row['drct']), ha='center',
                va='center', bbox=dict(color='white'))
        ax.text(row['sped_%s' % (units,)] * 0.98, mon,
                "%.1f" % (row['sped_%s' % (units,)],), ha='right',
                va='center', bbox=dict(color='white',
                                       boxstyle='square,pad=0.03',))
    ax.set_ylim(12.5, 0.5)
    ax.set_title(("[%s] %s [%s-%s]\nMonthly Average Wind Speed and"
                  " Vector Average Direction"
                  ) % (station, nt.sts[station]['name'],
                       df['ts'].min().year,
                       df['ts'].max().year))

    return fig, grp

def computeOthers(d):
    r = {}
    # Need something to compute other values needed for output
    for sid in d.keys():
        ob = d[sid]
        ob["ticks"] = calendar.timegm(ob['utc_valid'].timetuple())
        if ob['sknt'] is not None:
            ob["sped"] = ob["sknt"] * 1.17
        if ob.get('tmpf') is not None and ob.get('dwpf') is not None:
            tmpf = temperature(ob['tmpf'], 'F')
            dwpf = temperature(ob['dwpf'], 'F')
            ob["relh"] = meteorology.relh(tmpf, dwpf).value('%')
        else:
            ob['relh'] = None
        if ob['relh'] == 'M':
            ob['relh'] = None

        if (ob.get('tmpf') is not None and ob.get('dwpf') is not None and
                ob.get('sped') is not None):
            tmpf = temperature(ob['tmpf'], 'F')
            dwpf = temperature(ob['dwpf'], 'F')
            sknt = speed(ob['sped'], 'MPH')
            ob["feel"] = meteorology.feelslike(tmpf, dwpf, sknt).value("F")
        else:
            ob['feel'] = None
        if ob['feel'] == 'M':
            ob['feel'] = None

        ob["altiTend"] = 'S'
        ob["drctTxt"] = util.drct2text(ob["drct"])
        if ob["max_drct"] is None:
            ob["max_drct"] = 0
        ob["max_drctTxt"] = util.drct2text(ob["max_drct"])
        ob["20gu"] = 0
        if ob['gust'] is not None:
            ob["gmph"] = ob["gust"] * 1.17
        if ob['max_gust'] is not None:
            ob["max_sped"] = ob["max_gust"] * 1.17
        else:
            ob['max_sped'] = 0
        ob['pday'] = 0 if ob['pday'] is None else ob['pday']
        ob['pmonth'] = 0 if ob['pmonth'] is None else ob['pmonth']
        ob["gtim"] = "0000"
        ob["gtim2"] = "12:00 AM"
        if ob["max_gust_ts"] is not None and ob["max_gust_ts"] != "null":
            ob["gtim"] = ob["max_gust_ts"].strftime("%H%M")
            ob["gtim2"] = ob["max_gust_ts"].strftime("%-I:%M %p")
        r[sid] = ob
    return r

def grid_day(nc, ts):
    """
    """
    offset = iemre.daily_offset(ts)
    print(('cal hi/lo for %s [idx:%s]') % (ts, offset))
    sql = """
       SELECT ST_x(s.geom) as lon, ST_y(s.geom) as lat, s.state,
       s.name, s.id as station,
       (CASE WHEN pday >= 0 then pday else null end) as precipdata,
       (CASE WHEN max_tmpf > -50 and max_tmpf < 130
           then max_tmpf else null end) as highdata,
       (CASE WHEN min_tmpf > -50 and min_tmpf < 95
           then min_tmpf else null end) as lowdata,
       (CASE WHEN max_dwpf > -50 and max_dwpf < 130
           then max_dwpf else null end) as highdwpf,
       (CASE WHEN min_dwpf > -50 and min_dwpf < 95
           then min_dwpf else null end) as lowdwpf,
        (CASE WHEN avg_sknt >= 0 and avg_sknt < 100
         then avg_sknt else null end) as avgsknt
       from summary_%s c, stations s WHERE day = '%s' and
       s.network in ('IA_ASOS', 'MN_ASOS', 'WI_ASOS', 'IL_ASOS', 'MO_ASOS',
        'KS_ASOS', 'NE_ASOS', 'SD_ASOS', 'ND_ASOS', 'KY_ASOS', 'MI_ASOS',
        'OH_ASOS', 'AWOS') and c.iemid = s.iemid
        """ % (ts.year, ts.strftime("%Y-%m-%d"))
    df = read_sql(sql, pgconn)

    if len(df.index) > 4:
        res = generic_gridder(df, 'highdata')
        nc.variables['high_tmpk'][offset] = datatypes.temperature(
                                                res, 'F').value('K')
        res = generic_gridder(df, 'lowdata')
        nc.variables['low_tmpk'][offset] = datatypes.temperature(
                                            res, 'F').value('K')
        hres = generic_gridder(df, 'highdwpf')
        lres = generic_gridder(df, 'lowdwpf')
        nc.variables['avg_dwpk'][offset] = datatypes.temperature(
                                            (hres + lres) / 2., 'F').value('K')
        res = generic_gridder(df, 'avgsknt')
        res = np.where(res < 0, 0, res)
        nc.variables['wind_speed'][offset] = datatypes.speed(
                                            res, 'KT').value('MPS')
    else:
        print "%s has %02i entries, FAIL" % (ts.strftime("%Y-%m-%d"),
                                             cursor.rowcount)

def test_drct():
    """Conversion of u and v to direction"""
    r = meteorology.drct(
        datatypes.speed(np.array([10, 20]), 'KT'),
        datatypes.speed(np.array([10, 20]), 'KT')
    ).value("DEG")
    assert r[0] == 225
    r = meteorology.drct(
        datatypes.speed(-10, 'KT'),
        datatypes.speed(10, 'KT')
    ).value("DEG")
    assert r == 135
    r = meteorology.drct(
        datatypes.speed(-10, 'KT'),
        datatypes.speed(-10, 'KT')
    ).value("DEG")
    assert r == 45
    r = meteorology.drct(
        datatypes.speed(10, 'KT'),
        datatypes.speed(-10, 'KT')
    ).value("DEG")
    assert r == 315

def do_windalerts(obs):
    """Iterate through the obs and do wind alerts where appropriate"""
    for sid in obs:
        # Problem sites with lightning issues
        if sid in ['RBFI4', 'RTMI4', 'RWII4', 'RCAI4', 'RDYI4',
                   'RDNI4', 'RCDI4', 'RCII4', 'RCLI4']:
            continue
        ob = obs[sid]
        # screening
        if ob.get('gust', 0) < 40:
            continue
        if np.isnan(ob['gust']):
            continue
        smph = speed(ob['gust'], 'KT').value('MPH')
        if smph < 50:
            continue
        if smph > 100:
            print(('process_rwis did not relay gust %.1f MPH from %s'
                  '') % (smph, sid))
            continue
        # Use a hacky tmp file to denote a wind alert that was sent
        fn = "/tmp/iarwis.%s.%s" % (sid, ob['valid'].strftime("%Y%m%d%H%M"))
        if os.path.isfile(fn):
            continue
        o = open(fn, 'w')
        o.write(" ")
        o.close()
        lts = ob['valid'].astimezone(pytz.timezone("America/Chicago"))
        stname = NT.sts[sid]['name']
        msg = ("At %s, a wind gust of %.1f mph (%.1f kts) was recorded "
               "at the %s (%s) Iowa RWIS station"
               "") % (lts.strftime("%I:%M %p %d %b %Y"), smph, ob['gust'],
                      stname, sid)
        mt = MIMEText(msg)
        mt['From'] = '[email protected]'
        # mt['To'] = '[email protected]'
        mt['To'] = '[email protected]'
        mt['Subject'] = 'Iowa RWIS Wind Gust %.0f mph %s' % (smph, stname)
        s = smtplib.SMTP('mailhub.iastate.edu')
        s.sendmail(mt['From'], [mt['To']], mt.as_string())
        s.quit()