def do_month(year, month, routes):
    """ Generate a MRMS plot for the month!"""

    sts = datetime.datetime(year,month,1)
    ets = sts + datetime.timedelta(days=35)
    ets = ets.replace(day=1)

    today = datetime.datetime.now()
    if ets > today:
        ets = today

    idx0 = iemre.daily_offset(sts)
    idx1 = iemre.daily_offset(ets)

    nc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_mrms_daily.nc" % (year,),
                          'r')

    lats = nc.variables['lat'][:]
    lons = nc.variables['lon'][:]
    p01d = np.sum(nc.variables['p01d'][idx0:idx1,:,:],0) / 24.5
    nc.close()

    m = MapPlot(sector='iowa', title='MRMS %s - %s Total Precipitation' % (
            sts.strftime("%-d %b"), 
            (ets - datetime.timedelta(days=1)).strftime("%-d %b %Y")),
            subtitle='Data from NOAA MRMS Project')
    x,y = np.meshgrid(lons, lats)
    bins = [0.01, 0.1, 0.5, 1, 1.5, 2, 3, 4, 5, 6, 7, 8, 12, 16, 20]
    m.pcolormesh(x, y, p01d, bins, units='inches')
    m.drawcounties()
    currentfn = "summary/iowa_mrms_q3_month.png"
    archivefn = sts.strftime("%Y/%m/summary/iowa_mrms_q3_month.png")
    pqstr = "plot %s %s00 %s %s png" % (
                routes, sts.strftime("%Y%m%d%H"), currentfn, archivefn)
    m.postprocess(pqstr=pqstr)

def main():
    """Go Main"""
    total = None
    years = 0.
    for yr in range(1981, 2018):
        print(yr)
        ncfn = "/mesonet/data/prism/%s_daily.nc" % (yr, )
        nc = netCDF4.Dataset(ncfn)
        if total is None:
            lons = nc.variables['lon'][:]
            lats = nc.variables['lat'][:]
            total = np.zeros(nc.variables['tmax'].shape[1:], np.float)
        days = np.zeros(nc.variables['tmax'].shape[1:], np.float)
        sidx = daily_offset(datetime.date(yr, 1, 1))
        eidx = daily_offset(datetime.date(yr, 7, 4))
        for idx in range(sidx, eidx):
            days += np.where(nc.variables['tmax'][idx, :, :] > THRESHOLD,
                             1, 0)
        nc.close()
        years += 1.
        total += days

    val = days - (total / years)
    print(np.max(val))
    print(np.min(val))
    mp = MapPlot(sector='conus',
                 title=("OSU PRISM 2017 Days with High >= 90$^\circ$F "
                        "Departure"),
                 subtitle=("2017 thru 4 July against 1981-2016 "
                           "Year to Date Average"))
    mp.contourf(lons, lats, val, np.arange(-25, 26, 5),
                units='days', cmap=plt.get_cmap('seismic'))
    mp.postprocess(filename='test.png')

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 test_daily_offset():
    """ Compute the offsets """
    ts = utc(2013, 1, 1, 0, 0)
    offset = iemre.daily_offset(ts)
    assert offset == 0

    ts = datetime.date(2013, 2, 1)
    offset = iemre.daily_offset(ts)
    assert offset == 31

    ts = utc(2013, 1, 5, 12, 0)
    offset = iemre.daily_offset(ts)
    assert offset == 4

def do_precip12(nc, ts):
    """Compute the 24 Hour precip at 12 UTC, we do some more tricks though"""
    offset = iemre.daily_offset(ts)
    ets = datetime.datetime(ts.year, ts.month, ts.day, 12)
    ets = ets.replace(tzinfo=pytz.timezone("UTC"))
    sts = ets - datetime.timedelta(hours=24)
    offset1 = iemre.hourly_offset(sts)
    offset2 = iemre.hourly_offset(ets)
    if ts.month == 1 and ts.day == 1:
        print(("p01d_12z  for %s [idx:%s] %s(%s)->%s(%s) SPECIAL"
               ) % (ts, offset, sts.strftime("%Y%m%d%H"), offset1,
                    ets.strftime("%Y%m%d%H"), offset2))
        ncfn = "/mesonet/data/iemre/%s_mw_hourly.nc" % (ets.year,)
        if not os.path.isfile(ncfn):
            print("Missing %s" % (ncfn,))
            return
        hnc = netCDF4.Dataset(ncfn)
        phour = np.sum(hnc.variables['p01m'][:offset2, :, :], 0)
        hnc.close()
        hnc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_hourly.nc" % (
            sts.year,))
        phour += np.sum(hnc.variables['p01m'][offset1:, :, :], 0)
        hnc.close()
    else:
        print(("p01d_12z  for %s [idx:%s] %s(%s)->%s(%s)"
               ) % (ts, offset, sts.strftime("%Y%m%d%H"), offset1,
                    ets.strftime("%Y%m%d%H"), offset2))
        ncfn = "/mesonet/data/iemre/%s_mw_hourly.nc" % (ts.year,)
        if not os.path.isfile(ncfn):
            print("Missing %s" % (ncfn,))
            return
        hnc = netCDF4.Dataset(ncfn)
        phour = np.sum(hnc.variables['p01m'][offset1:offset2, :, :], 0)
        hnc.close()
    nc.variables['p01d_12z'][offset] = phour

def do_hrrr(ts):
    """Convert the hourly HRRR data to IEMRE grid"""
    total = None
    xaxis = None
    yaxis = None
    for hr in range(5, 23):  # Only need 5 AM to 10 PM for solar
        utc = ts.replace(hour=hr).astimezone(pytz.timezone("UTC"))
        fn = utc.strftime(("/mesonet/ARCHIVE/data/%Y/%m/%d/model/hrrr/%H/"
                           "hrrr.t%Hz.3kmf00.grib2"))
        if not os.path.isfile(fn):
            # print 'HRRR file %s missing' % (fn,)
            continue
        grbs = pygrib.open(fn)
        try:
            if utc >= SWITCH_DATE:
                grb = grbs.select(name='Downward short-wave radiation flux')
            else:
                grb = grbs.select(parameterNumber=192)
        except ValueError:
            print 'coop/hrrr_solarrad.py %s had no solar rad' % (fn,)
            continue
        if len(grb) == 0:
            print 'Could not find SWDOWN in HRR %s' % (fn,)
            continue
        g = grb[0]
        if total is None:
            total = g.values
            lat1 = g['latitudeOfFirstGridPointInDegrees']
            lon1 = g['longitudeOfFirstGridPointInDegrees']
            llcrnrx, llcrnry = LCC(lon1, lat1)
            nx = g['Nx']
            ny = g['Ny']
            dx = g['DxInMetres']
            dy = g['DyInMetres']
            xaxis = llcrnrx + dx * np.arange(nx)
            yaxis = llcrnry + dy * np.arange(ny)
        else:
            total += g.values

    if total is None:
        print 'coop/hrrr_solarrad.py found no HRRR data for %s' % (
                                                    ts.strftime("%d %b %Y"), )
        return

    # We wanna store as W m-2, so we just average out the data by hour
    total = total / 24.0

    nc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_daily.nc" % (ts.year,),
                         'a')
    offset = iemre.daily_offset(ts)
    data = nc.variables['rsds'][offset, :, :]
    for i, lon in enumerate(iemre.XAXIS):
        for j, lat in enumerate(iemre.YAXIS):
            (x, y) = LCC(lon, lat)
            i2 = np.digitize([x], xaxis)[0]
            j2 = np.digitize([y], yaxis)[0]
            data[j, i] = total[j2, i2]

    nc.variables['rsds'][offset] = data
    nc.close()

def plot_daily(date, interval, plotvar, mc, mckey):
    """ Generate the plot, please """
    opts = PLOT_OPS[plotvar]
    offset = iemre.daily_offset(date)
    
    nc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_daily.nc" % (date.year,))
    data = nc.variables[opts['ncvar_daily']][offset] / 25.4 # inches
    lons = nc.variables['lon'][:]
    lats = nc.variables['lat'][:]
    extra = lons[-1] + (lons[-1] - lons[-2])
    lons = np.concatenate([lons, [extra,]])

    extra = lats[-1] + (lats[-1] - lats[-2])
    lats = np.concatenate([lats, [extra,]])
    x,y = np.meshgrid(lons, lats)

    nc.close()

    p = plot.MapPlot(sector='midwest',
                 title='%s IEM Reanalysis %s [%s]' % (date.strftime("%-d %b %Y"),
                                                 opts['title'], opts['units'])
                 )
    p.pcolormesh(x, y, data, opts['clevs'],
                 units=opts['units'])
    p.postprocess(web=True, memcache=mc, memcachekey=mckey, memcacheexpire=0)

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 run(ts):
    ''' Actually do the work, please '''
    nc = netCDF4.Dataset('/mesonet/data/iemre/%s_mw_mrms_daily.nc' % (
                                                            ts.year,),
                         'a')
    offset = iemre.daily_offset(ts)
    ncprecip = nc.variables['p01d']
    ts += datetime.timedelta(hours=24)
    gmtts = ts.astimezone(pytz.timezone("UTC"))

    fn = gmtts.strftime(("/mesonet/ARCHIVE/data/%Y/%m/%d/GIS/q2/"
                        "p24h_%Y%m%d%H00.png"))
    img = Image.open(fn)
    data = np.asarray(img)
    # data is 3500,7000 , starting at upper L
    data = np.flipud(data)
    # Anything over 254 is bad
    res = np.where(data > 254, 0, data)
    res = np.where(np.logical_and(data >= 0, data < 100), data * 0.25, res)
    res = np.where(np.logical_and(data >= 100, data < 180),
                   25. + ((data - 100) * 1.25), res)
    res = np.where(np.logical_and(data >= 180, data < 255),
                   125. + ((data - 180) * 5.), res)

    y1 = (iemre.NORTH - mrms.SOUTH) * 100.0
    y0 = (iemre.SOUTH - mrms.SOUTH) * 100.0
    x0 = (iemre.WEST - mrms.WEST) * 100.0
    x1 = (iemre.EAST - mrms.WEST) * 100.0
    ncprecip[offset, :, :] = res[y0:y1, x0:x1]
    nc.close()

def do_coop(ts):
    """Use COOP solar radiation data"""
    pgconn = psycopg2.connect(database='coop', host='iemdb', user='nobody')
    cursor = pgconn.cursor()

    cursor.execute("""SELECT ST_x(geom), ST_y(geom),
        coalesce(narr_srad, merra_srad) from alldata a JOIN stations t
        ON (a.station = t.id) WHERE
        day = %s and t.network ~* 'CLIMATE' and substr(id, 3, 1) != 'C'
        and substr(id, 3, 4) != '0000'
    """, (ts.strftime("%Y-%m-%d"), ))
    lons = []
    lats = []
    vals = []
    for row in cursor:
        if row[2] is None or row[2] < 0:
            continue
        lons.append(row[0])
        lats.append(row[1])
        vals.append(row[2])

    nn = NearestNDInterpolator((np.array(lons), np.array(lats)),
                               np.array(vals))
    xi, yi = np.meshgrid(iemre.XAXIS, iemre.YAXIS)

    nc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_daily.nc" % (ts.year,),
                         'a')
    offset = iemre.daily_offset(ts)
    # Data above is MJ / d / m-2, we want W / m-2
    nc.variables['rsds'][offset, :, :] = nn(xi, yi) * 1000000. / 86400.
    nc.close()

def do_precip(nc, ts):
    """Compute the 6 UTC to 6 UTC precip

    We need to be careful here as the timestamp sent to this app is today,
    we are actually creating the analysis for yesterday
    """
    sts = datetime.datetime(ts.year, ts.month, ts.day, 6)
    sts = sts.replace(tzinfo=pytz.timezone("UTC"))
    ets = sts + datetime.timedelta(hours=24)
    offset = iemre.daily_offset(ts)
    offset1 = iemre.hourly_offset(sts)
    offset2 = iemre.hourly_offset(ets)
    if ts.month == 12 and ts.day == 31:
        print(("p01d      for %s [idx:%s] %s(%s)->%s(%s) SPECIAL"
               ) % (ts, offset, sts.strftime("%Y%m%d%H"), offset1,
                    ets.strftime("%Y%m%d%H"), offset2))
        hnc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_hourly.nc" % (
            ets.year,))
        phour = np.sum(hnc.variables['p01m'][:offset2, :, :], 0)
        hnc.close()
        hnc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_hourly.nc" % (
            sts.year,))
        phour += np.sum(hnc.variables['p01m'][offset1:, :, :], 0)
        hnc.close()
    else:
        print(("p01d      for %s [idx:%s] %s(%s)->%s(%s)"
               ) % (ts, offset, sts.strftime("%Y%m%d%H"), offset1,
                    ets.strftime("%Y%m%d%H"), offset2))
        hnc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_hourly.nc" % (
            sts.year,))
        phour = np.sum(hnc.variables['p01m'][offset1:offset2, :, :], 0)
        hnc.close()
    nc.variables['p01d'][offset] = phour

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 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 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 do_var(varname):
    """
    Run our estimator for a given variable
    """
    currentnc = None
    sql = """select day, station from alldata_%s WHERE %s is null
        and day >= '1893-01-01' ORDER by day ASC""" % (state.lower(), varname)
    ccursor.execute(sql)
    for row in ccursor:
        day = row[0]
        station = row[1]
        if station not in nt.sts:
            continue

        sql = """
            SELECT station, %s from alldata_%s WHERE %s is not NULL
            and station in %s and day = '%s'
            """ % (varname, state, varname, tuple(friends[station]), day)
        ccursor2.execute(sql)
        weight = []
        value = []
        for row2 in ccursor2:
            idx = friends[station].index(row2[0])
            weight.append(weights[station][idx])
            value.append(row2[1])

        if len(weight) < 3:
            # Nearest neighbors failed, so lets look at our grided analysis
            # and sample from it
            if currentnc is None or currentnc.title.find(str(day.year)) == -1:
                currentnc = netCDF4.Dataset(("/mesonet/data/iemre/"
                                             "%s_mw_daily.nc") % (day.year,))
            tidx = iemre.daily_offset(datetime.datetime(day.year, day.month,
                                                        day.day))
            iidx, jidx = iemre.find_ij(nt.sts[station]['lon'],
                                       nt.sts[station]['lat'])
            iemreval = currentnc.variables[vnameconv[varname]][tidx, jidx,
                                                               iidx]
            if varname in ('high', 'low'):
                interp = temperature(iemreval, 'K').value('F')
            else:
                interp = distance(iemreval, 'MM').value('IN')
            print '--> Neighbor failure, %s %s %s' % (station, day, varname)
        else:
            mass = sum(weight)
            interp = np.sum(np.array(weight) * np.array(value) / mass)

        dataformat = '%.2f'
        if varname in ['high', 'low']:
            dataformat = '%.0f'
        print(('Set station: %s day: %s varname: %s value: %s'
               ) % (station, day, varname, dataformat % (interp,)))
        sql = """
            UPDATE alldata_%s SET estimated = true, %s = %s WHERE
            station = '%s' and day = '%s'
            """ % (state.lower(), varname,
                   dataformat % (interp,), station, day)
        sql = sql.replace(' nan ', ' null ')
        ccursor2.execute(sql)

def do_process(valid, fn):
    """Process this file, please """
    # shape of data is (1, 621, 1405)
    data = rasterio.open(fn).read()
    varname = fn.split("_")[1]
    nc = netCDF4.Dataset("/mesonet/data/prism/%s_daily.nc" % (valid.year,),
                         'a')
    idx = daily_offset(valid)
    nc.variables[varname][idx] = np.flipud(data[0])
    nc.close()

def plotter(fdict):
    """ Go """
    import matplotlib
    matplotlib.use('agg')
    from pyiem.plot import MapPlot
    ptype = fdict.get('ptype', 'c')
    date = datetime.datetime.strptime(fdict.get('date', '2015-01-01'),
                                      '%Y-%m-%d')
    varname = fdict.get('var', 'rsds')

    idx0 = iemre.daily_offset(date)
    nc = netCDF4.Dataset(("/mesonet/data/iemre/%s_mw_daily.nc"
                          ) % (date.year, ), 'r')
    lats = nc.variables['lat'][:]
    lons = nc.variables['lon'][:]
    if varname == 'rsds':
        # Value is in W m**-2, we want MJ
        data = nc.variables[varname][idx0, :, :] * 86400. / 1000000.
        units = 'MJ d-1'
        clevs = np.arange(0, 37, 3.)
        clevs[0] = 0.01
        clevstride = 1
    elif varname in ['p01d', 'p01d_12z']:
        # Value is in W m**-2, we want MJ
        data = nc.variables[varname][idx0, :, :] / 25.4
        units = 'inch'
        clevs = np.arange(0, 0.25, 0.05)
        clevs = np.append(clevs, np.arange(0.25, 3., 0.25))
        clevs = np.append(clevs, np.arange(3., 10.0, 1))
        clevs[0] = 0.01
        clevstride = 1
    elif varname in ['high_tmpk', 'low_tmpk', 'high_tmpk_12z', 'low_tmpk_12z']:
        # Value is in W m**-2, we want MJ
        data = temperature(nc.variables[varname][idx0, :, :], 'K').value('F')
        units = 'F'
        clevs = np.arange(-30, 120, 2)
        clevstride = 5
    nc.close()

    title = date.strftime("%-d %B %Y")
    m = MapPlot(sector='midwest', axisbg='white', nocaption=True,
                title='IEM Reanalysis of %s for %s' % (PDICT.get(varname),
                                                       title),
                subtitle='Data derived from various NOAA datasets'
                )
    if np.ma.is_masked(np.max(data)):
        return 'Data Unavailable'
    x, y = np.meshgrid(lons, lats)
    if ptype == 'c':
        m.contourf(x, y, data, clevs, clevstride=clevstride, units=units)
    else:
        m.pcolormesh(x, y, data, clevs, clevstride=clevstride, units=units)

    return m.fig

def do_precip(nc, ts):
    """Compute the precip totals based on the hourly analysis totals"""
    offset = iemre.daily_offset(ts)
    ets = ts.replace(hour=12, tzinfo=pytz.timezone("UTC"))
    sts = ets - datetime.timedelta(hours=24)
    offset1 = iemre.hourly_offset(sts)
    offset2 = iemre.hourly_offset(ets)
    hnc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_hourly.nc" % (ts.year,))
    phour = np.sum(hnc.variables['p01m'][offset1:offset2, :, :], 0)
    nc.variables['p01d'][offset] = phour
    hnc.close()

def run(ts):
    """ Actually do the work, please """
    nc = netCDF4.Dataset(('/mesonet/data/iemre/%s_mw_mrms_daily.nc'
                          '') % (ts.year,), 'a')
    offset = iemre.daily_offset(ts)
    ncprecip = nc.variables['p01d']

    # We want this mrms variable to replicate the netcdf file, so the
    # origin is the southwestern corner
    ts += datetime.timedelta(hours=24)
    gmtts = ts.astimezone(pytz.timezone("UTC"))

    gribfn = gmtts.strftime(("/mnt/a4/data/%Y/%m/%d/mrms/ncep/"
                             "RadarOnly_QPE_24H/"
                             "RadarOnly_QPE_24H_00.00_%Y%m%d-%H%M00.grib2.gz"))
    if not os.path.isfile(gribfn):
        print("merge_mrms_q3.py MISSING %s" % (gribfn,))
        return

    fp = gzip.GzipFile(gribfn, 'rb')
    (_, tmpfn) = tempfile.mkstemp()
    tmpfp = open(tmpfn, 'wb')
    tmpfp.write(fp.read())
    tmpfp.close()
    grbs = pygrib.open(tmpfn)
    grb = grbs[1]
    lats, _ = grb.latlons()
    os.unlink(tmpfn)

    val = grb['values']
    # Anything less than zero, we set to zero
    val = np.where(val < 0, 0, val)

    # CAREFUL HERE!  The MRMS grid is North to South
    # set top (smallest y)
    y0 = int((lats[0, 0] - iemre.NORTH) * 100.0)
    y1 = int((lats[0, 0] - iemre.SOUTH) * 100.0)
    x0 = int((iemre.WEST - mrms.WEST) * 100.0)
    x1 = int((iemre.EAST - mrms.WEST) * 100.0)
    # print 'y0:%s y1:%s x0:%s x1:%s' % (y0, y1, x0, x1)
    ncprecip[offset, :, :] = np.flipud(val[y0:y1, x0:x1])
    # m = MapPlot(sector='midwest')
    # x, y = np.meshgrid(nc.variables['lon'][:], nc.variables['lat'][:])
    # m.pcolormesh(x, y, ncprecip[offset,:,:], range(10), latlon=True)
    # m.postprocess(filename='test.png')
    # (fig, ax) = plt.subplots()
    # ax.imshow(mrms)
    # fig.savefig('test.png')
    # (fig, ax) = plt.subplots()
    # ax.imshow(mrms[y0:y1,x0:x1])
    # fig.savefig('test2.png')
    nc.close()

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")