def main():
    """Go Main Go"""
    pgconn = get_dbconn('scan')
    for station in ['S2004', 'S2196', 'S2002', 'S2072', 'S2068',
                    'S2031', 'S2001', 'S2047']:
        df = read_sql("""
        select extract(year from valid + '2 months'::interval) as wy,
        tmpf, dwpf from alldata where station = %s and tmpf is not null
        and dwpf is not null
        """, pgconn, params=(station, ), index_col=None)
        df['mixingratio'] = meteorology.mixing_ratio(
            temperature(df['dwpf'].values, 'F')).value('KG/KG')
        df['vapor_pressure'] = mcalc.vapor_pressure(
            1000. * units.mbar,
            df['mixingratio'].values * units('kg/kg')).to(units('kPa'))
        df['saturation_mixingratio'] = (
            meteorology.mixing_ratio(
                temperature(df['tmpf'].values, 'F')).value('KG/KG'))
        df['saturation_vapor_pressure'] = mcalc.vapor_pressure(
            1000. * units.mbar,
            df['saturation_mixingratio'].values * units('kg/kg')).to(units('kPa'))
        df['vpd'] = df['saturation_vapor_pressure'] - df['vapor_pressure']
        means = df.groupby('wy').mean()
        counts = df.groupby('wy').count()
        for yr, row in means.iterrows():
            print(("%s,%s,%.0f,%.3f"
                   ) % (yr, station, counts.at[yr, 'vpd'], row['vpd']))

def grid_day(nc, ts):
    """
    """
    offset = iemre.daily_offset(ts)
    icursor.execute("""
       SELECT ST_x(s.geom) as lon, ST_y(s.geom) as lat, 
       (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 
       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")))

    if icursor.rowcount > 4:
        res = generic_gridder(icursor, 'highdata')
        nc.variables['high_tmpk'][offset] = datatypes.temperature(res, 'F').value('K')
        icursor.scroll(0, mode='absolute')
        res = generic_gridder(icursor, 'lowdata')
        nc.variables['low_tmpk'][offset] = datatypes.temperature(res, 'F').value('K')
        icursor.scroll(0, mode='absolute')
        #res = generic_gridder(icursor, 'precipdata')
        #nc.variables['p01d'][offset] = res * 25.4
    else:
        print "%s has %02i entries, FAIL" % (ts.strftime("%Y-%m-%d"), 
            icursor.rowcount)

def get(station):
    cursor.execute(
        """
      SELECT valid, sknt, tmpf, dwpf from alldata where station = %s 
      and tmpf is not null and dwpf is not null 
      and valid > '1971-01-01' ORDER by valid ASC 
    """,
        (station,),
    )

    hits = {}
    running = False
    startr = None
    for row in cursor:
        relh = met.relh(temperature(row[2], "F"), temperature(row[3], "F")).value("%")
        if relh > 25 or row[1] < (25.0 / 1.15):
            if running:
                delta = (row[0] - startr).seconds
                if delta >= 60 * 60 * 1:
                    # print station, delta, row
                    hits[row[0].strftime("%Y%m%d")] = 1
            running = False
        else:
            running = True
            startr = row[0]

    return len(hits.keys())

def replace_forecast(df, location):
    """Replace dataframe data with forecast for this location"""
    pgconn = psycopg2.connect(database='coop', host='iemdb', user='nobody')
    cursor = pgconn.cursor()
    today = datetime.date.today()
    nextjan1 = datetime.date(today.year + 1, 1, 1)
    coop = XREF[location]['climodat']
    years = [int(y) for y in np.arange(df.index.values.min().year,
                                       df.index.values.max().year + 1)]
    cursor.execute("""
        SELECT day, high, low, precip from alldata_forecast WHERE
        modelid = (SELECT id from forecast_inventory WHERE model = 'NDFD'
        ORDER by modelts DESC LIMIT 1) and station = %s and day >= %s
    """, (coop, today))
    rcols = ['maxt', 'mint', 'rain']
    for row in cursor:
        valid = row[0]
        maxc = temperature(row[1], 'F').value('C')
        minc = temperature(row[2], 'F').value('C')
        rain = distance(row[3], 'IN').value('MM')
        for year in years:
            df.loc[valid.replace(year=year), rcols] = (maxc, minc, rain)

    # Need to get radiation from CFS
    cursor.execute("""
        SELECT day, srad from alldata_forecast WHERE
        modelid = (SELECT id from forecast_inventory WHERE model = 'CFS'
        ORDER by modelts DESC LIMIT 1) and station = %s and day >= %s
        and day < %s
    """, (coop, today, nextjan1))
    for row in cursor:
        valid = row[0]
        for year in years:
            df.loc[valid.replace(year=year), 'radn'] = row[1]

def load(dirname, location, sdate):
    """ Read a file please """
    data = []
    idx = []
    mindoy = int(sdate.strftime("%j"))
    for line in open("%s/%s.met" % (dirname, location)):
        line = line.strip()
        if not line.startswith('19') and not line.startswith('20'):
            continue
        tokens = line.split()
        if int(tokens[1]) < mindoy:
            continue
        data.append(tokens)
        ts = (datetime.date(int(tokens[0]), 1, 1) +
              datetime.timedelta(days=int(tokens[1])-1))
        idx.append(ts)
    if len(data[0]) < 10:
        cols = ['year', 'doy', 'radn', 'maxt', 'mint', 'rain']
    else:
        cols = ['year', 'doy', 'radn', 'maxt', 'mint',
                'rain', 'gdd', 'st4', 'st12', 'st24',
                'st50', 'sm12', 'sm24', 'sm50']
    df = pd.DataFrame(data, index=idx,
                      columns=cols)
    for col in cols:
        df[col] = pd.to_numeric(df[col], errors='coerce')
    if len(data[0]) < 10:
        df['gdd'] = gdd(temperature(df['maxt'].values, 'C'),
                        temperature(df['mint'].values, 'C'))
    df['gddcum'] = df.groupby(['year'])['gdd'].apply(lambda x: x.cumsum())
    df['raincum'] = distance(
        df.groupby(['year'])['rain'].apply(lambda x: x.cumsum()),
        'MM').value('IN')
    return df

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
       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')
    else:
        print "%s has %02i entries, FAIL" % (ts.strftime("%Y-%m-%d"),
                                             cursor.rowcount)

def figure(val, qcval):
    if qcval > 1000:
        return None
    if np.ma.is_masked(val) or np.ma.is_masked(qcval):
        return None
    return temperature(val + qcval,
                       'K').value('F') - temperature(val, 'K').value('F')

def estimate_hilo(ts):
    """Estimate the High and Low Temperature based on gridded data"""
    idx = iemre.daily_offset(ts)
    nc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_daily.nc" % (ts.year, ),
                         'r')
    highgrid12 = temperature(nc.variables['high_tmpk_12z'][idx, :, :],
                             'K').value('F')
    lowgrid12 = temperature(nc.variables['low_tmpk_12z'][idx, :, :],
                            'K').value('F')
    highgrid00 = temperature(nc.variables['high_tmpk'][idx, :, :],
                             'K').value('F')
    lowgrid00 = temperature(nc.variables['low_tmpk'][idx, :, :],
                            'K').value('F')
    nc.close()

    for sid in nt.sts.keys():
        if nt.sts[sid]['temp24_hour'] in [0, 22, 23]:
            val = highgrid00[nt.sts[sid]['gridj'], nt.sts[sid]['gridi']]
        else:
            val = highgrid12[nt.sts[sid]['gridj'], nt.sts[sid]['gridi']]
        if val > -80 and val < 140:
            nt.sts[sid]['high'] = "%.0f" % (val, )

        if nt.sts[sid]['temp24_hour'] in [0, 22, 23]:
            val = lowgrid00[nt.sts[sid]['gridj'], nt.sts[sid]['gridi']]
        else:
            val = lowgrid12[nt.sts[sid]['gridj'], nt.sts[sid]['gridi']]
        if val > -80 and val < 140:
            nt.sts[sid]['low'] = "%.0f" % (val, )

def grid_day(nc, ts):
    """
    I proctor the gridding of data on an hourly basis
    @param ts Timestamp of the analysis, we'll consider a 20 minute window
    """
    offset = iemre.daily_offset(ts)
    if ts.day == 29 and ts.month == 2:
        ts = datetime.datetime(2000, 3, 1)

    sql = """SELECT * from ncdc_climate71 WHERE valid = '%s' and
             substr(station,3,4) != '0000' and substr(station,3,1) != 'C'
             """ % (ts.strftime("%Y-%m-%d"), )
    cursor.execute(sql)
    if cursor.rowcount > 4:
        res = generic_gridder(cursor, 'high')
        if res is not None:
            nc.variables['high_tmpk'][offset] = datatypes.temperature(res, 'F').value('K')
        cursor.scroll(0, mode='absolute')
        res = generic_gridder(cursor, 'low')
        if res is not None:
            nc.variables['low_tmpk'][offset] = datatypes.temperature(res, 'F').value('K')
        cursor.scroll(0, mode='absolute')
        res = generic_gridder(cursor, 'precip')
        if res is not None:
            nc.variables['p01d'][offset] = res * 25.4
    else:
        print "%s has %02i entries, FAIL" % (ts.strftime("%Y-%m-%d"),
                                             cursor.rowcount)

def main():
    """Do Something Fun!"""
    form = cgi.FormContent()
    ts = datetime.datetime.strptime(form["date"][0], "%Y-%m-%d")
    lat = float(form["lat"][0])
    lon = float(form["lon"][0])
    fmt = form["format"][0]
    if fmt != 'json':
        sys.stdout.write("Content-type: text/plain\n\n")
        sys.stdout.write("ERROR: Service only emits json at this time")
        return

    i, j = iemre.find_ij(lon, lat)
    offset = iemre.daily_offset(ts)

    res = {'data': [], }

    fn = "/mesonet/data/iemre/%s_mw_daily.nc" % (ts.year,)

    sys.stdout.write('Content-type: application/json\n\n')
    if not os.path.isfile(fn):
        sys.stdout.write(json.dumps(res))
        sys.exit()

    if i is None or j is None:
        sys.stdout.write(json.dumps({'error': 'Coordinates outside of domain'}
                                    ))
        return

    nc = netCDF4.Dataset(fn, 'r')

    c2000 = ts.replace(year=2000)
    coffset = iemre.daily_offset(c2000)

    cnc = netCDF4.Dataset("/mesonet/data/iemre/mw_dailyc.nc", 'r')

    res['data'].append({
        'daily_high_f': myrounder(
           datatypes.temperature(
                nc.variables['high_tmpk'][offset, j, i], 'K').value('F'), 1),
        'climate_daily_high_f': myrounder(
           datatypes.temperature(
                cnc.variables['high_tmpk'][coffset, j, i], 'K').value("F"), 1),
        'daily_low_f': myrounder(
           datatypes.temperature(
                nc.variables['low_tmpk'][offset, j, i], 'K').value("F"), 1),
        'climate_daily_low_f': myrounder(
           datatypes.temperature(
                cnc.variables['low_tmpk'][coffset, j, i], 'K').value("F"), 1),
        'daily_precip_in': myrounder(
           nc.variables['p01d'][offset, j, i] / 25.4, 2),
        'climate_daily_precip_in': myrounder(
           cnc.variables['p01d'][coffset, j, i] / 25.4, 2),
      })
    nc.close()
    cnc.close()

    sys.stdout.write(json.dumps(res))

def do_salus(ctx):
    """ Generate SALUS
    StationID, Year, DOY, SRAD, Tmax, Tmin, Rain, DewP, Wind, Par, dbnum
    CTRL, 1981, 1, 5.62203, 2.79032, -3.53361, 5.43766, NaN, NaN, NaN, 2
    CTRL, 1981, 2, 3.1898, 1.59032, -6.83361, 1.38607, NaN, NaN, NaN, 3
    """
    if len(ctx['stations']) > 1:
        ssw(("ERROR: SALUS output is only "
             "permitted for one station at a time."))
        return

    dbconn = get_database()
    cursor = dbconn.cursor(cursor_factory=psycopg2.extras.DictCursor)

    scenario_year = 2030
    asts = datetime.date(2030, 1, 1)
    if ctx['scenario'] == 'yes':
        # Tricky!
        scenario_year = ctx['scenario_year']
        today = datetime.date.today()
        asts = datetime.date(scenario_year, today.month, today.day)

    table = get_tablename(ctx['stations'])
    station = ctx['stations'][0]
    thisyear = datetime.datetime.now().year
    cursor.execute("""
    WITH scenario as (
        SELECT
 ('""" + str(thisyear) + """-'||month||'-'||extract(day from day))::date as day,
        high, low, precip, station,
        coalesce(narr_srad, merra_srad, hrrr_srad) as srad
        from """ + table + """ WHERE station = %s and
        day >= %s and year = %s
    ), obs as (
        SELECT day,
        high, low, precip,  station,
        coalesce(narr_srad, merra_srad, hrrr_srad) as srad
        from """ + table + """ WHERE station = %s and
        day >= %s and day <= %s ORDER by day ASC
    ), total as (
        SELECT *, extract(doy from day) as doy from obs
        UNION SELECT * from scenario
    )

    SELECT * from total ORDER by day ASC
    """, (station, asts, scenario_year, station, ctx['sts'], ctx['ets']))
    ssw(("StationID, Year, DOY, SRAD, Tmax, Tmin, Rain, DewP, "
         "Wind, Par, dbnum\n"))
    for i, row in enumerate(cursor):
        srad = -99 if row['srad'] is None else row['srad']
        ssw(("%s, %s, %s, %.4f, %.2f, %.2f, %.2f, , , , %s\n"
             ) % (
                  station[:4], row["day"].year,
                  int(row["doy"]), srad,
                  temperature(row["high"], 'F').value('C'),
                  temperature(row["low"], 'F').value('C'),
                  row["precip"] * 25.4, i + 2))

def test_gdd_with_nans():
    """Can we properly deal with nan's and not emit warnings?"""
    highs = np.ma.array([70, 80, np.nan, 90],
                        mask=[False, False, True, False])
    lows = highs - 10
    r = meteorology.gdd(datatypes.temperature(highs, 'F'),
                        datatypes.temperature(lows, 'F'),
                        50, 86)
    assert np.ma.is_masked(r[2])

def dbsave(ts, data):
    """Save the data! """
    pgconn = psycopg2.connect(database="coop", host="iemdb")
    cursor = pgconn.cursor()
    # Check to see if we already have data for this date
    cursor.execute(
        """SELECT id from forecast_inventory
      WHERE model = 'CFS' and modelts = %s""",
        (ts,),
    )
    if cursor.rowcount > 0:
        modelid = cursor.fetchone()[0]
        cursor.execute(
            """DELETE from alldata_forecast where
        modelid = %s""",
            (modelid,),
        )
        if cursor.rowcount > 0:
            print("Removed %s previous entries" % (cursor.rowcount,))
    else:
        cursor.execute(
            """INSERT into forecast_inventory(model, modelts)
        VALUES ('CFS', %s) RETURNING id""",
            (ts,),
        )
        modelid = cursor.fetchone()[0]

    for date in data["fx"].keys():
        d = data["fx"][date]
        if d["high"] is None or d["low"] is None or d["precip"] is None or d["srad"] is None:
            print("Missing data for date: %s" % (date,))
            del (data["fx"][date])

    for sid in nt.sts.keys():
        # Skip virtual stations
        if sid[2:] == "0000" or sid[2] == "C":
            continue
        # Careful here, lon is 0-360 for this file
        i = np.digitize([nt.sts[sid]["lon"] + 360], data["x"])[0]
        j = np.digitize([nt.sts[sid]["lat"]], data["y"])[0]
        for date in data["fx"]:
            d = data["fx"][date]
            high = bnds(temperature(d["high"][j, i], "K").value("F"), -70, 140)
            low = bnds(temperature(d["low"][j, i], "K").value("F"), -90, 120)
            precip = bnds(round(float(d["precip"][j, i] / 25.4), 2), 0, 30)
            srad = bnds(d["srad"][j, i] / 1000000.0, 0, 50)
            if high is None or low is None or precip is None or srad is None:
                continue
            cursor.execute(
                """INSERT into alldata_forecast(modelid,
            station, day, high, low, precip, srad)
            VALUES (%s, %s, %s, %s, %s, %s, %s)
            """,
                (modelid, sid, date, high, low, precip, srad),
            )
    cursor.close()
    pgconn.commit()

def load(dirname, location, sdate):
    """ Read a file please """
    data = []
    idx = []
    for line in open("%s/%s.met" % (dirname, location)):
        line = line.strip()
        if not line.startswith("19") and not line.startswith("20"):
            continue
        tokens = line.split()
        data.append(tokens)
        ts = datetime.date(int(tokens[0]), 1, 1) + datetime.timedelta(days=int(tokens[1]) - 1)
        idx.append(ts)
    if len(data[0]) < 10:
        cols = ["year", "doy", "radn", "maxt", "mint", "rain"]
    else:
        cols = [
            "year",
            "doy",
            "radn",
            "maxt",
            "mint",
            "rain",
            "gdd",
            "st4",
            "st12",
            "st24",
            "st50",
            "sm12",
            "sm24",
            "sm50",
        ]
    df = pd.DataFrame(data, index=idx, columns=cols)
    for col in cols:
        df[col] = pd.to_numeric(df[col], errors="coerce")
    if len(data[0]) < 10:
        df["gdd"] = gdd(temperature(df["maxt"].values, "C"), temperature(df["mint"].values, "C"))
    bins = []
    today = datetime.date.today()
    for valid, _ in df.iterrows():
        if valid >= today:
            bins.append(0)
            continue
        if sdate == "nov1" and valid.month >= 11:
            bins.append(valid.year + 1)
            continue
        if valid.month < today.month:
            bins.append(valid.year)
            continue
        if valid.month == today.month and valid.day < today.day:
            bins.append(valid.year)
            continue
        bins.append(0)
    df["bin"] = bins
    df["rain"] = distance(df["rain"].values, "MM").value("IN")
    df["avgt"] = temperature((df["maxt"] + df["mint"]) / 2.0, "C").value("F")
    return df

def test_heatindex():
    ''' Test our heat index calculations '''
    t = datatypes.temperature(80.0, 'F')
    td = datatypes.temperature(70.0, 'F')
    hdx = meteorology.heatindex(t, td)
    assert abs(hdx.value("F") - 83.93) < 0.01

    t = datatypes.temperature(30.0, 'F')
    hdx = meteorology.heatindex(t, td)
    assert abs(hdx.value("F") - 30.00) < 0.01

def process(model, lon, lat):
    """ Generate a plot for this given combination """
    (fig, ax) = plt.subplots(1,1)
    modelts = get_latest_time(model)
    if modelts is None:
        ax.text(0.5, 0.5, "No Data Found to Plot!", ha='center')
        sys.stdout.write("Content-Type: image/png\n\n")
        fig.savefig( sys.stdout, format="png")
        return
    nc = netCDF4.Dataset(
            modelts.strftime("/mesonet/share/frost/"+model+"/%Y%m%d%H%M_iaoutput.nc"),'r')
    times = get_times(nc)
    i, j = get_ij(lon, lat, nc)
    
    ax.plot(times, temperature(nc.variables['bdeckt'][:,i,j],'K').value('F'),
            color='k', label='Bridge Deck Temp' if model == 'bridget' else 'Pavement')
    ax.plot(times, temperature(nc.variables['tmpk'][:,i,j], 'K').value("F"),
            color='r', label='Air Temp')
    ax.plot(times, temperature(nc.variables['dwpk'][:,i,j], 'K').value("F"),
            color='g', label='Dew Point')
    #ax.set_ylim(-30,150)
    ax.set_title(("ISUMM5 %s Timeseries\n"
                 +"i: %s j:%s lon: %.2f lat: %.2f Model Run: %s") % (model,
                    i, j, nc.variables['lon'][i,j], nc.variables['lat'][i,j],
                    modelts.astimezone(pytz.timezone("America/Chicago")).strftime(
                            "%-d %b %Y %-I:%M %p")))
    
    ax.xaxis.set_major_locator(
                               mdates.DayLocator(interval=1,
                                        tz=pytz.timezone("America/Chicago"))
                               )
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%d %b\n%Y',
                                        tz=pytz.timezone("America/Chicago")))
    ax.axhline(32, linestyle='-.')
    ax.grid(True)
    ax.set_ylabel("Temperature $^\circ$F")

    (ymin, ymax) = ax.get_ylim()

    for i2, ifrost in enumerate(nc.variables['ifrost'][:-1,i,j]):
        ax.barh(ymax-1, 1.0/24.0/4.0, left=times[i2], 
                fc=get_ifrost_color(ifrost), ec='none')
    for i2, icond in enumerate(nc.variables['icond'][:-1,i,j]):
        ax.barh(ymax-2, 1.0/24.0/4.0, left=times[i2], 
                fc=get_icond_color(model, icond), ec='none')

    # Shrink current axis's height by 10% on the bottom
    box = ax.get_position()
    ax.set_position([box.x0, box.y0 + box.height * 0.1,
                 box.width, box.height * 0.9])
    ax.legend(loc='upper center',
              bbox_to_anchor=(0.5, -0.08),fancybox=True, shadow=True, ncol=3)
    
    sys.stdout.write("Content-Type: image/png\n\n")
    fig.savefig( sys.stdout, format="png")

def daily_process(nwsli, maxts):
    """ Process the daily file """
    # print '-------------- DAILY PROCESS ----------------'
    fn = "%s%s" % (BASE, STATIONS[nwsli]['daily'].split("/")[1])
    if not os.path.isfile(fn):
        return 0
    lines = open(fn).readlines()
    if len(lines) < 5:
        return 0
    # Read header....
    headers = []
    for col in lines[1].strip().replace('"', '').split(","):
        headers.append(VARCONV.get(col.lower(), col.lower()))
    # Read data
    processed = 0
    for i in range(len(lines)-1, 3, -1):
        tokens = lines[i].strip().replace('"','').split(",")
        if len(tokens) != len(headers):
            continue
        valid = datetime.datetime.strptime(tokens[ headers.index('timestamp')][:10],
                                           '%Y-%m-%d')
        valid = valid.date() - datetime.timedelta(days=1)
        if valid < maxts:
            break
        if valid == maxts: # Reprocess
            icursor.execute("""DELETE from sm_daily WHERE valid = '%s' and
            station = '%s' """ % (valid.strftime("%Y-%m-%d"), nwsli))
        # We are ready for dbinserting!
        dbcols = "station,valid," + ",".join(headers[2:])
        dbvals = "'%s','%s'," % (nwsli, valid.strftime("%Y-%m-%d"))
        for v in tokens[2:]:
            dbvals += "%s," % (formatter(v),)
        sql = "INSERT into sm_daily (%s) values (%s)" % (dbcols, dbvals[:-1])
        icursor.execute(sql)

        # Need a timezone
        valid = datetime.datetime(valid.year, valid.month, valid.day, 12, 0)
        valid = valid.replace(tzinfo=pytz.timezone("America/Chicago"))
        ob = Observation(nwsli, 'ISUSM', valid)
        ob.data['max_tmpf'] = temperature(
                    float(tokens[headers.index('tair_c_max')]), 'C').value('F')
        ob.data['min_tmpf'] = temperature(
                    float(tokens[headers.index('tair_c_min')]), 'C').value('F')
        ob.data['pday'] = round(
            float(tokens[headers.index('rain_mm_tot')]) / 24.5, 2)
        ob.data['et_inch'] = float(tokens[headers.index('dailyet')]) / 24.5
        ob.data['srad_mj'] = float(tokens[headers.index('slrmj_tot')])
        ob.data['max_sknt'] = float(tokens[headers.index('ws_mps_max')]) * 1.94
        ob.save(accesstxn)
        #    print 'soilm_ingest.py station: %s ts: %s daily updated no data?' % (
        #                                nwsli, valid.strftime("%Y-%m-%d"))
        processed += 1
    return processed

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 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 load_iemre():
    """Use IEM Reanalysis for non-precip data

    24km product is smoothed down to the 0.01 degree grid
    """
    printt("load_iemre() called")
    xaxis = np.arange(MYWEST, MYEAST, 0.01)
    yaxis = np.arange(MYSOUTH, MYNORTH, 0.01)
    xi, yi = np.meshgrid(xaxis, yaxis)

    fn = iemre.get_daily_ncname(VALID.year)
    if not os.path.isfile(fn):
        printt("Missing %s for load_solar, aborting" % (fn,))
        sys.exit()
    nc = netCDF4.Dataset(fn, 'r')
    offset = iemre.daily_offset(VALID)
    lats = nc.variables['lat'][:]
    lons = nc.variables['lon'][:]
    lons, lats = np.meshgrid(lons, lats)

    # Storage is W m-2, we want langleys per day
    data = nc.variables['rsds'][offset, :, :] * 86400. / 1000000. * 23.9
    # Default to a value of 300 when this data is missing, for some reason
    nn = NearestNDInterpolator((np.ravel(lons), np.ravel(lats)),
                               np.ravel(data))
    SOLAR[:] = iemre_bounds_check('rsds', nn(xi, yi), 0, 1000)

    data = temperature(nc.variables['high_tmpk'][offset, :, :], 'K').value('C')
    nn = NearestNDInterpolator((np.ravel(lons), np.ravel(lats)),
                               np.ravel(data))
    HIGH_TEMP[:] = iemre_bounds_check('high_tmpk', nn(xi, yi), -60, 60)

    data = temperature(nc.variables['low_tmpk'][offset, :, :], 'K').value('C')
    nn = NearestNDInterpolator((np.ravel(lons), np.ravel(lats)),
                               np.ravel(data))
    LOW_TEMP[:] = iemre_bounds_check('low_tmpk', nn(xi, yi), -60, 60)

    data = temperature(nc.variables['avg_dwpk'][offset, :, :], 'K').value('C')
    nn = NearestNDInterpolator((np.ravel(lons), np.ravel(lats)),
                               np.ravel(data))
    DEWPOINT[:] = iemre_bounds_check('avg_dwpk', nn(xi, yi), -60, 60)

    data = nc.variables['wind_speed'][offset, :, :]
    nn = NearestNDInterpolator((np.ravel(lons), np.ravel(lats)),
                               np.ravel(data))
    WIND[:] = iemre_bounds_check('wind_speed', nn(xi, yi), 0, 30)

    nc.close()
    printt("load_iemre() finished")