Python MapPlot.drawcounties方法代码示例

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def doday():
    """
    Create a plot of precipitation stage4 estimates for some day
    """
    sts = mx.DateTime.DateTime(2013,5,25,12)
    ets = mx.DateTime.DateTime(2013,5,31,12)
    interval = mx.DateTime.RelativeDateTime(days=1)
    now = sts
    total = None
    while now < ets:
        fp = "/mesonet/ARCHIVE/data/%s/stage4/ST4.%s.24h.grib" % (
            now.strftime("%Y/%m/%d"), 
            now.strftime("%Y%m%d%H") )
        if os.path.isfile(fp):
            lts = now
            grbs = pygrib.open(fp)

            if total is None:
                g = grbs[1]
                total = g["values"]
                lats, lons = g.latlons()
            else:
                total += grbs[1]["values"]
            grbs.close()
        now += interval
        
    m = MapPlot(sector='iowa', title='NOAA Stage IV & Iowa ASOS Precipitation',
                subtitle='25-30 May 2013')
    m.pcolormesh(lons, lats, total / 25.4, numpy.arange(0,14.1,1), latlon=True,
                 units='inch')
    m.drawcounties()
    m.plot_values(dlons, dlats, dvals, '%.02f')
    m.postprocess(filename='test.svg')
    import iemplot
    iemplot.makefeature('test')

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def run(base, ceil, now, fn):
    """ Generate the plot """
    # Compute normal from the climate database
    sql = """SELECT station,
       sum(gddxx(%s, %s, high, low)) as gdd
       from alldata_ia WHERE year = %s and month in (5,6,7,8,9,10)
       and station != 'IA0000' and substr(station,2,1) != 'C'
       GROUP by station""" % (base, ceil, now.year)

    lats = []
    lons = []
    gdd50 = []
    ccursor.execute(sql)
    for row in ccursor:
        if row[0] not in nt.sts:
            continue
        lats.append(nt.sts[row[0]]['lat'])
        lons.append(nt.sts[row[0]]['lon'])
        gdd50.append(float(row[1]))

    m = MapPlot(title=("Iowa 1 May - %s GDD Accumulation"
                       ) % (now.strftime("%-d %B %Y"), ),
                subtitle="base %s" % (base,))
    bins = np.linspace(min(gdd50)-1, max(gdd50)+1, num=10, dtype=np.int)
    m.contourf(lons, lats, gdd50, bins)
    m.drawcounties()

    pqstr = "plot c 000000000000 summary/%s.png bogus png" % (fn,)
    m.postprocess(pqstr=pqstr)

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def main():
    """Do Something"""
    cursor = IEM.cursor()
    data = []
    cursor.execute("""SELECT ST_x(geom), ST_y(geom), tsf0, tsf1, tsf2, tsf3,
    id, rwis_subf from current c JOIN stations t on (t.iemid = c.iemid)
    WHERE c.valid > now() - '1 hour'::interval""")
    for row in cursor:
        val = cln(row[2:6])
        if val is None:
            continue
        d = dict(lat=row[1], lon=row[0], tmpf=val, id=row[6])
        if row[7] is not None and not np.isnan(row[7]):
            d['dwpf'] = row[7]
        data.append(d)

    now = datetime.datetime.now()
    m = MapPlot(axisbg='white',
                title='Iowa RWIS Average Pavement + Sub-Surface Temperature',
                subtitle=("Valid: %s (pavement in red, sub-surface in blue)"
                          "") % (now.strftime("%-d %b %Y %-I:%M %p"),))
    m.plot_station(data)
    m.drawcounties()
    pqstr = ("plot c %s rwis_sf.png rwis_sf.png png"
             "") % (datetime.datetime.utcnow().strftime("%Y%m%d%H%M"), )
    m.postprocess(view=False, pqstr=pqstr)

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def runYear(year):
    # Grab the data
    sql = """SELECT station,
        sum(case when precip >= 0.01 then 1 else 0 end) as days, max(day)
        from alldata_ia WHERE year = %s and substr(station,3,1) != 'C' 
        and station != 'IA0000' GROUP by station""" % (year,)

    lats = []
    lons = []
    vals = []
    labels = []
    ccursor.execute( sql )
    for row in ccursor:
        sid = row['station'].upper()
        if not nt.sts.has_key(sid):
            continue
        labels.append( sid[2:] )
        lats.append( nt.sts[sid]['lat'] )
        lons.append( nt.sts[sid]['lon'] )
        vals.append( row['days'] )
        maxday = row['max']

    #---------- Plot the points
    m = MapPlot(title="Days with Measurable Precipitation (%s)" % (year,),
                subtitle='Map valid January 1 - %s' % (maxday.strftime("%b %d")),
                axisbg='white')
    m.plot_values(lons, lats, vals, fmt='%.0f', labels=labels,
                  labeltextsize=8, labelcolor='tan')
    m.drawcounties()
    pqstr = "plot m %s bogus %s/summary/precip_days.png png" % (
                                        now.strftime("%Y%m%d%H%M"), year,)
    m.postprocess(pqstr=pqstr)

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
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)

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def plot(data, v):
    ''' Actually plot this data '''
    nt = NetworkTable("ISUSM")
    lats = []
    lons = []
    vals = []
    valid = None
    for sid in data.keys():
        if data[sid][v] is None:
            continue
        lats.append(nt.sts[sid]['lat'])
        lons.append(nt.sts[sid]['lon'])
        vals.append(data[sid][v])
        valid = data[sid]['valid']

    if valid is None:
        m = MapPlot(sector='iowa', axisbg='white',
                    title=('ISU Soil Moisture Network :: %s'
                           '') % (CTX[v]['title'], ),
                    figsize=(8.0, 6.4))
        m.plot_values([-95, ], [41.99, ], ['No Data Found'], '%s', textsize=30)
        m.postprocess(web=True)
        return

    m = MapPlot(sector='iowa', axisbg='white',
                title='ISU Soil Moisture Network :: %s' % (CTX[v]['title'],),
                subtitle='valid %s' % (valid.strftime("%-d %B %Y %I:%M %p"),),
                figsize=(8.0, 6.4))
    m.plot_values(lons, lats, vals, '%.1f')
    m.drawcounties()
    m.postprocess(web=True)

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def plotter(fdict):
    """ Go """
    import matplotlib
    matplotlib.use('agg')
    from pyiem.plot import MapPlot
    import matplotlib.cm as cm

    pgconn = psycopg2.connect(database='coop', host='iemdb', user='nobody')
    cursor = pgconn.cursor()

    sector = fdict.get('sector', 'IA')
    date1 = datetime.datetime.strptime(fdict.get('date1', '2015-01-01'),
                                       '%Y-%m-%d')
    date2 = datetime.datetime.strptime(fdict.get('date2', '2015-02-01'),
                                       '%Y-%m-%d')

    table = "alldata_%s" % (sector, ) if sector != 'midwest' else "alldata"
    cursor.execute("""
    WITH obs as (
        SELECT station, sday, day, precip from """ + table + """ WHERE
        day >= %s and day < %s and precip >= 0 and
        substr(station, 3, 1) != 'C' and substr(station, 3, 4) != '0000'),
    climo as (
        SELECT station, to_char(valid, 'mmdd') as sday, precip from
        climate51),
    combo as (
        SELECT o.station, o.precip - c.precip as d from obs o JOIN climo c ON
        (o.station = c.station and o.sday = c.sday)),
    deltas as (
        SELECT station, sum(d) from combo GROUP by station)

    SELECT d.station, d.sum, ST_x(t.geom), ST_y(t.geom) from deltas d
    JOIN stations t on (d.station = t.id) WHERE t.network ~* 'CLIMATE'
    """, (date1, date2))

    rows = []
    for row in cursor:
        rows.append(dict(station=row[0], delta=row[1], lon=row[2],
                         lat=row[3]))
    df = pd.DataFrame(rows)
    lons = np.array(df['lon'])
    vals = np.array(df['delta'])
    lats = np.array(df['lat'])
    sector2 = "state" if sector != 'midwest' else 'midwest'
    m = MapPlot(sector=sector2, state=sector, axisbg='white',
                title=('%s - %s Precipitation Departure [inch]'
                       ) % (date1.strftime("%d %b %Y"),
                            date2.strftime("%d %b %Y")),
                subtitle='%s vs 1950-2014 Climatology' % (date1.year,))
    rng = int(max([0 - np.min(vals), np.max(vals)]))
    cmap = cm.get_cmap('RdYlBu')
    cmap.set_bad('white')
    m.contourf(lons, lats, vals, np.linspace(0 - rng - 0.5, rng + 0.6, 10,
                                             dtype='i'),
               cmap=cmap, units='inch')
    m.plot_values(lons, lats, vals, fmt='%.2f')
    if sector == 'iowa':
        m.drawcounties()

    return m.fig, df

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def main():
    """Go Main"""
    pgconn = get_dbconn('postgis')
    df = read_postgis("""
    select geom, issue from sbw where wfo = 'PUB' and phenomena = 'TO'
    and significance = 'W' and status = 'NEW' and issue > '2007-10-01'
    and issue < '2019-01-01'
    """, pgconn, geom_col='geom', crs={'init': 'epsg:4326', 'no_defs': True})

    bounds = df['geom'].total_bounds
    # bounds = [-102.90293903,   40.08745967,  -97.75622311,   43.35172981]
    bbuf = 0.25
    mp = MapPlot(
        sector='custom', west=bounds[0] - bbuf,
        south=bounds[1] - bbuf,
        east=bounds[2] + bbuf, north=bounds[3] + bbuf,
        continentalcolor='white',  # '#b3242c',
        title='NWS Pueblo Issued Tornado Warnings [2008-2018]',
        subtitle='%s warnings plotted' % (len(df.index), ))
    crs_new = ccrs.Mercator()
    crs = ccrs.PlateCarree()
    new_geometries = [crs_new.project_geometry(ii, src_crs=crs)
                      for ii in df['geom'].values]
    # mp.draw_cwas()
    mp.ax.add_geometries(new_geometries, crs=crs_new, lw=0.5,
                         edgecolor='red', facecolor='None', alpha=1,
                         zorder=5)
    mp.drawcounties()
    mp.postprocess(filename='test.png')

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def main():
    """Go!"""
    title = 'NOAA MRMS Q3: RADAR + Guage Corrected Rainfall Estimates + NWS Storm Reports'
    mp = MapPlot(sector='custom',
                 north=42.3, east=-93.0, south=41.65, west=-94.1,
                 axisbg='white',
                 titlefontsize=14,
                 title=title,
                 subtitle='Valid: 14 June 2018')

    shp = shapefile.Reader('cities.shp')
    for record in shp.shapeRecords():
        geo = shape(record.shape)
        mp.ax.add_geometries([geo], ccrs.PlateCarree(), zorder=Z_OVERLAY2,
                             facecolor='None', edgecolor='k', lw=2)

    grbs = pygrib.open('MRMS_GaugeCorr_QPE_24H_00.00_20180614-200000.grib2')
    grb = grbs.message(1)
    pcpn = distance(grb['values'], 'MM').value('IN')
    lats, lons = grb.latlons()
    lons -= 360.
    clevs = [0.01, 0.1, 0.3, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 8, 10]
    cmap = nwsprecip()
    cmap.set_over('k')

    mp.pcolormesh(lons, lats, pcpn, clevs, cmap=cmap, latlon=True,
                  units='inch')
    lons, lats, vals, labels = get_data()
    mp.drawcounties()
    mp.plot_values(lons, lats, vals, "%s", labels=labels,
                   labelbuffer=1, labelcolor='white')

    mp.drawcities(labelbuffer=5, minarea=0.2)
    mp.postprocess(filename='test.png')

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def main():
    """Go Main"""
    grbs = pygrib.open('ds.snow.bin')
    # skip 1-off first field
    total = None
    lats = lons = None
    for grb in grbs[1:]:
        if lats is None:
            lats, lons = grb.latlons()
            total = grb['values']
            continue
        total += grb['values']
    # TODO tz-hack here
    analtime = grb.analDate - datetime.timedelta(hours=5)

    mp = MapPlot(
        sector='custom', west=-100, east=-92, north=45, south=41,
        axisbg='tan',
        title=("NWS Forecasted Accumulated Snowfall "
               "thru 7 PM 12 April 2019"),
        subtitle='NDFD Forecast Issued %s' % (
            analtime.strftime("%-I %p %-d %B %Y"), )
    )
    cmap = nwssnow()
    cmap.set_bad('tan')
    mp.pcolormesh(
        lons, lats, total * 39.3701,
        [0.01, 1, 2, 3, 4, 6, 8, 12, 18, 24, 30, 36],
        cmap=cmap,
        units='inch')

    mp.drawcounties()
    mp.drawcities()
    mp.postprocess(filename='test.png')
    mp.close()

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def run(year):
    cursor.execute("""
    WITH obs as (
    SELECT station, sum(precip) from alldata_ia where year = %s GROUP by station
    ), climate as (
    SELECT station, sum(precip) from climate51 GROUP by station
    )
    SELECT o.station, o.sum - c.sum as diff from obs o JOIN climate c 
    on (c.station = o.station) ORDER by diff ASC
    """, (year,))
    lats = []
    lons = []
    vals = []
    for row in cursor:
        if not nt.sts.has_key(row[0]) or row[0] in rejs or row[0][2] == 'C':
            continue
        print row
        lats.append( nt.sts[row[0]]['lat'])
        lons.append( nt.sts[row[0]]['lon'])
        vals.append( row[1] )
    
    m = MapPlot(title='%s Precipitation Departure' % (year,))
    cmap = cm.get_cmap('BrBG')
    #cmap.set_over('blue')
    #cmap.set_under('red')
    m.contourf(lons, lats, vals, np.arange(-24,24.1,2), cmap=cmap, units='inch')
    #m.plot_values(lons, lats, vals, '%.02f')
    m.drawcounties()
    m.postprocess(filename='%s.png' % (year,))

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def doday(ts, realtime):
    """
    Create a plot of precipitation stage4 estimates for some day

    We should total files from 1 AM to midnight local time
    """
    sts = ts.replace(hour=1)
    ets = sts + datetime.timedelta(hours=24)
    interval = datetime.timedelta(hours=1)
    now = sts
    total = None
    lts = None
    while now < ets:
        gmt = now.astimezone(pytz.timezone("UTC"))
        fn = gmt.strftime(("/mesonet/ARCHIVE/data/%Y/%m/%d/"
                           +"stage4/ST4.%Y%m%d%H.01h.grib"))
        if os.path.isfile(fn):
            lts = now
            grbs = pygrib.open(fn)

            if total is None:
                g = grbs[1]
                total = g["values"]
                lats, lons = g.latlons()
            else:
                total += grbs[1]["values"]
            grbs.close()
        now += interval

    if lts is None and ts.hour > 1:
        print 'stage4_today_total.py found no data!'
    if lts is None:
        return
    lts = lts - datetime.timedelta(minutes=1)
    subtitle = "Total between 12:00 AM and %s" % (lts.strftime("%I:%M %p %Z"),)
    routes = 'ac'
    if not realtime:
        routes = 'a'
    for sector in ['iowa', 'midwest', 'conus']:
        pqstr = "plot %s %s00 %s_stage4_1d.png %s_stage4_1d.png png" % (routes,
                ts.strftime("%Y%m%d%H"), sector, sector )
        
        m = MapPlot(sector=sector,
                    title="%s NCEP Stage IV Today's Precipitation" % (
                                                    ts.strftime("%-d %b %Y"),),
                    subtitle=subtitle)
            
        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
    
        m.pcolormesh(lons, lats, total / 24.5, clevs, units='inch')
    
        #map.drawstates(zorder=2)
        if sector == 'iowa':
            m.drawcounties()
        m.postprocess(pqstr=pqstr)
        m.close()

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def plotter(fdict):
    """ Go """
    import matplotlib
    matplotlib.use('agg')
    pgconn = psycopg2.connect(database='coop', host='iemdb', user='nobody')

    state = fdict.get('state', 'IA')[:2]
    varname = fdict.get('var', 'total_precip')
    sector = fdict.get('sector', 'state')
    opt = fdict.get('opt', 'both')
    over = fdict.get('over', 'monthly')
    month = int(fdict.get('month', datetime.date.today().month))

    df = read_sql("""
    WITH data as (
        SELECT station, extract(month from valid) as month,
        sum(precip) as total_precip, avg(high) as avg_high,
        avg(low) as avg_low, avg((high+low)/2.) as avg_temp
        from ncdc_climate81 GROUP by station, month)

    SELECT station, ST_X(geom) as lon, ST_Y(geom) as lat, month,
    total_precip, avg_high, avg_low, avg_temp from data d JOIN stations t
    ON (d.station = t.id) WHERE t.network = 'NCDC81' and
    t.state in ('IA', 'ND', 'SD', 'NE', 'KS', 'MO', 'IL', 'WI', 'MN', 'MI',
    'IN', 'OH', 'KY')
    """, pgconn, index_col=['station', 'month'])

    if over == 'monthly':
        title = "%s %s" % (calendar.month_name[month], PDICT3[varname])
        df.reset_index(inplace=True)
        df2 = df[df['month'] == month]
    else:
        title = "Annual %s" % (PDICT3[varname], )
        if varname == 'total_precip':
            df2 = df.sum(axis=0, level='station')
        else:
            df2 = df.mean(axis=0, level='station')
        df2['lat'] = df['lat'].mean(axis=0, level='station')
        df2['lon'] = df['lon'].mean(axis=0, level='station')
    m = MapPlot(sector=sector, state=state, axisbg='white',
                title=('NCEI 1981-2010 Climatology of %s'
                       ) % (title,),
                subtitle=('based on National Centers for '
                          'Environmental Information (NCEI) 1981-2010'
                          ' Climatology'))
    levels = np.linspace(df2[varname].min(), df2[varname].max(), 10)
    levels = [round(x, PRECISION[varname]) for x in levels]
    if opt in ['both', 'contour']:
        m.contourf(df2['lon'].values, df2['lat'].values,
                   df2[varname].values, levels, units=UNITS[varname])
    if sector == 'state':
        m.drawcounties()
    if opt in ['both', 'values']:
        m.plot_values(df2['lon'].values, df2['lat'].values,
                      df2[varname].values,
                      fmt='%%.%if' % (PRECISION[varname],))

    return m.fig, df

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def main():
    """GO!"""
    pgconn = get_dbconn('idep')
    cursor = pgconn.cursor()

    scenario = int(sys.argv[1])
    mp = MapPlot(sector='iowa', axisbg='white', nologo=True,
                 subtitle='1 Jan 2014 thru 31 Dec 2014',
                 caption='Daily Erosion Project',
                 title=('Harvest Index 0.8 Change in 2014 Soil Delivery '
                        'from Baseline'))

    cursor.execute("""
    with baseline as (
        SELECT huc_12, sum(avg_delivery) * 4.463 as loss from results_by_huc12
        where valid between '2014-01-01' and '2015-01-01' and
        scenario = 0 GROUP by huc_12),
    scenario as (
        SELECT huc_12, sum(avg_delivery) * 4.463 as loss from results_by_huc12
        where valid between '2014-01-01' and '2015-01-01' and
        scenario = %s GROUP by huc_12),
    agg as (
        SELECT b.huc_12, b.loss as baseline_loss, s.loss as scenario_loss from
        baseline b LEFT JOIN scenario s on (b.huc_12 = s.huc_12))

     SELECT ST_Transform(simple_geom, 4326),
     (scenario_loss  - baseline_loss) / 1.0 as val, i.huc_12
     from huc12 i JOIN agg d on (d.huc_12 = i.huc_12)
     WHERE i.states ~* 'IA' ORDER by val DESC

    """, (scenario, ))

    # bins = np.arange(0, 101, 10)
    bins = [-5, -2, -1, -0.5, 0, 0.5, 1, 2, 5]
    cmap = plt.get_cmap("BrBG_r")
    cmap.set_under('purple')
    cmap.set_over('black')
    norm = mpcolors.BoundaryNorm(bins, cmap.N)

    for row in cursor:
        # print "%s,%s" % (row[2], row[1])
        polygon = loads(row[0].decode('hex'))
        arr = np.asarray(polygon.exterior)
        points = mp.ax.projection.transform_points(ccrs.Geodetic(),
                                                   arr[:, 0], arr[:, 1])
        val = float(row[1])
        # We have very small negative numbers that should just be near a
        # positive zero
        if val < 0 and val > -0.1:
            val = 0.001
        color = cmap(norm([val, ]))[0]
        poly = Polygon(points[:, :2], fc=color, ec='k', zorder=2, lw=0.1)
        mp.ax.add_patch(poly)

    mp.draw_colorbar(bins, cmap, norm, units='T/a/yr')

    mp.drawcounties()
    mp.postprocess(filename='test.png')

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import drawcounties [as 别名]
def do(ts, hours):
    """
    Create a plot of precipitation stage4 estimates for some day
    """
    ts = ts.replace(minute=0)
    sts = ts - datetime.timedelta(hours=hours)
    ets = ts 
    interval = datetime.timedelta(hours=1)
    now = sts
    total = None
    lts = None
    while now < ets:
        fn = "/mesonet/ARCHIVE/data/%s/stage4/ST4.%s.01h.grib" % (
            now.strftime("%Y/%m/%d"), 
            now.strftime("%Y%m%d%H") )

        if os.path.isfile(fn):
            lts = now
            grbs = pygrib.open(fn)

            if total is None:
                g = grbs[1]
                total = g["values"]
                lats, lons = g.latlons()
            else:
                total += grbs[1]["values"]
            grbs.close()
        now += interval
    
    if lts is None and ts.hour > 1:
        print 'Missing StageIV data!'
    if lts is None:
        return
    
    cmap = cm.get_cmap("jet")
    cmap.set_under('white')
    cmap.set_over('black')
    clevs = [0.01,0.1,0.25,0.5,1,2,3,5,8,9.9]
    localtime = (ts - datetime.timedelta(minutes=1)).astimezone(
                                        pytz.timezone("America/Chicago"))

    for sector in ['iowa', 'midwest', 'conus']:
        m = MapPlot(sector=sector,
                    title='NCEP Stage IV %s Hour Precipitation' % (hours,),
                    subtitle='Total up to %s' % (
                                    localtime.strftime("%d %B %Y %I %p %Z"),))
        m.pcolormesh(lons, lats, total / 24.5, clevs, units='inch')
        pqstr = "plot %s %s00 %s_stage4_%sh.png %s_stage4_%sh_%s.png png" % (
                                'ac', ts.strftime("%Y%m%d%H"), sector, hours,
                                sector, hours, ts.strftime("%H"))
        if sector == 'iowa':
            m.drawcounties()
        m.postprocess(pqstr=pqstr)
        m.close()