Python plot.MapPlot类代码示例

def plotter(fdict):
    """ Go """
    import matplotlib
    matplotlib.use('agg')
    from pyiem.plot import MapPlot
    utc = datetime.datetime.utcnow()
    bins = [0, 1, 14, 31, 91, 182, 273, 365, 730, 1460, 2920, 3800]
    pgconn = psycopg2.connect(database='postgis', host='iemdb', user='nobody')
    cursor = pgconn.cursor(cursor_factory=psycopg2.extras.DictCursor)

    phenomena = fdict.get('phenomena', 'TO')
    significance = fdict.get('significance', 'W')

    cursor.execute("""
     select wfo,  extract(days from ('TODAY'::date - max(issue))) as m
     from warnings where significance = %s and phenomena = %s
     GROUP by wfo ORDER by m ASC
    """, (significance, phenomena))
    data = {}
    rows = []
    for row in cursor:
        wfo = row[0] if row[0] != 'JSJ' else 'SJU'
        rows.append(dict(wfo=wfo, days=row[1]))
        data[wfo] = max([row[1], 0])
    df = pd.DataFrame(rows)

    m = MapPlot(sector='nws', axisbg='white', nocaption=True,
                title='Days since Last %s %s by NWS Office' % (
                        vtec._phenDict.get(phenomena, phenomena),
                        vtec._sigDict.get(significance, significance)),
                subtitle='Valid %s' % (utc.strftime("%d %b %Y %H%M UTC"),))
    m.fill_cwas(data, bins=bins, ilabel=True, units='Days',
                lblformat='%.0f')

    return m.fig, df

def draw_map():
    """make maps, not war."""
    m = MapPlot(sector='conus',
                title='4 March 2019 :: DEP Precip Points')
    update_grid(m.ax)
    m.postprocess(filename='/tmp/map_clipoints.png')
    m.close()

def main():
    """Map some CLI data"""
    pgconn = get_dbconn('iem')

    df = read_sql("""
    WITH data as (
        SELECT station, snow_jul1 - snow_jul1_normal as s
        from cli_data where valid = '2019-02-18' and snow_jul1 > 0
        and snow_jul1_normal > 0)

    select station, st_x(geom) as lon, st_y(geom) as lat, c.s as val from
    data c JOIN stations s on (s.id = c.station)
    WHERE s.network = 'NWSCLI'
    """, pgconn, index_col=None)
    df['color'] = '#ff0000'
    df.loc[df['val'] > 0, 'color'] = '#0000ff'

    mp = MapPlot(sector='midwest', axisbg='white',
                 title=("2018-2019 Snowfall Total Departure "
                        "from Average [inches]"),
                 subtitle='18 Feb 2019 Based on NWS CLI Reporting Sites')
    mp.plot_values(
        df['lon'].values, df['lat'].values,
        df['val'].values, fmt='%.1f', textsize=12, color=df['color'].values,
        labelbuffer=1)
    mp.postprocess(filename='test.png')

def test_drawugcs2():
    """3 filled zones"""
    mp = MapPlot(sector='iowa', title='Zones, 3 filled in Iowa, label',
                 subtitle='test_drawugcs2',
                 nocaption=True)
    mp.fill_ugcs({"IAZ001": 10, "IAZ003": 20, "IAZ005": 30}, ilabel=True)
    return mp.fig

def test_fillstates():
    """Can we fill states"""
    data = {'AK': 10, 'HI': 30, 'IA': 40, 'NY': 80}
    mp = MapPlot(sector='nws', title='Fill AK, HI, IA, NY States',
                 subtitle='test_fillstates', nocaption=True)
    mp.fill_states(data, ilabel=True)
    return mp.fig

def runYear(year):
    # Grab the data
    sql = """SELECT station, sum(precip) as total, max(day)
           from alldata_ia WHERE year = %s and
           station != 'IA0000' and
           substr(station,3,1) != 'C' and
           precip is not null GROUP by station""" % (year,)

    lats = []
    lons = []
    vals = []
    labels = []
    ccursor.execute( sql )
    for row in ccursor:
        sid = row['station']
        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['total'] )
        maxday = row['max']

    m = MapPlot(title="Total Precipitation [inch] (%s)" % (year,),
                subtitle='1 January - %s' % (maxday.strftime("%d %B"),),
                axisbg='white')
    m.plot_values(lons, lats, vals, labels=labels, fmt='%.2f',
                  labeltextsize=8, labelcolor='tan')
    pqstr = "plot m %s bogus %s/summary/total_precip.png png" % (
                                        now.strftime("%Y%m%d%H%M"), year,)
    m.postprocess(pqstr=pqstr)

def plot():
    """Do plotting work"""
    cmap1 = plt.get_cmap('inferno_r')
    colors = list(cmap1(np.arange(10) / 10.))
    cmap2 = plt.get_cmap('Pastel1')
    colors.extend(list(cmap2(np.arange(2) / 2.)))
    cmap = ListedColormap(colors)
    
    cmap.set_under('tan')
    cmap.set_over('white')
    minval = np.load('minval.npy')
    maxval = np.load('maxval.npy')
    diff = maxval - minval
    lons = np.load('lons.npy')
    lats = np.load('lats.npy')
    mp = MapPlot(sector='midwest', statebordercolor='white',
                 title=(r"Diff between coldest wind chill and warmest "
                        "air temp 29 Jan - 3 Feb 2019"),
                 subtitle=("based on hourly NCEP Real-Time Mesoscale Analysis "
                           "(RTMA) ending midnight CST"))

    levels = list(range(0, 101, 10))
    levels.extend([105, 110])
    mp.pcolormesh(lons, lats, diff, levels,
                  cmap=cmap, clip_on=False,
                  units=r"$^\circ$F", spacing='proportional')
    mp.postprocess(filename='test.png')

def test_contourf():
    ''' Test the contourf plot with labels specified '''
    mp = MapPlot(sector='iowa', nocaption=True)
    mp.contourf(np.arange(-94, -89), np.arange(40, 45),
                np.arange(5), np.arange(5),
                clevlabels=['a', 'b', 'c', 'd', 'e'])
    return mp.fig

def test_michigan():
    """See what we do with Michigan"""
    mp = MapPlot(sector='state', state='MI', nocaption=True)
    mp.contourf(np.arange(-84, -75), np.arange(36, 45), np.arange(9),
                np.arange(9),
                clevlabels=['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i'])
    return mp.fig

def main():
    """Go Main"""
    df = get_database_data()
    print(df)
    vals = {}
    labels = {}
    for wfo, row in df.iterrows():
        if wfo == 'JSJ':
            wfo = 'SJU'
        vals[wfo] = row['percent']
        labels[wfo] = '%.0f%%' % (row['percent'], )
        #if row['count'] == 0:
        #    labels[wfo] = '-'

    bins = np.arange(0, 101, 10)    
    #bins = [1, 25, 50, 75, 100, 125, 150, 200, 300]
    #bins = [-50, -25, -10, -5, 0, 5, 10, 25, 50]
    # bins[0] = 1
    #clevlabels = ['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW', 'N']
    cmap = plt.get_cmap('PuOr')
    mp = MapPlot(sector='nws', continentalcolor='white', figsize=(12., 9.),
                 title=("2018 Percentage of Time with 1+ Flood Warning Active"),
                 subtitle=('1 January - 30 September 2018, based on IEM archives'))
    mp.fill_cwas(vals, bins=bins, lblformat='%s', labels=labels,
                 cmap=cmap, ilabel=True,  # clevlabels=clevlabels,
                 units='percent')
    
    mp.postprocess(filename='test.png')

def main():
    """Go MAin"""
    df = pd.read_csv('flood_emergencies.csv')
    df2 = df[['source', 'eventid', 'phenomena', 'significance', 'year']
             ].drop_duplicates()
    gdf = df2.groupby('source').count()
    vals = {}
    labels = {}
    for wfo, row in gdf.iterrows():
        if wfo == 'TJSJ':
            wfo = 'SJU'
        else:
            wfo = wfo[1:]
        vals[wfo] = int(row['eventid'])
        labels[wfo] = "%s" % (row['eventid'], )

    bins = list(range(0, 31, 3))
    bins[0] = 1.
    cmap = plt.get_cmap('plasma_r')
    cmap.set_over('black')
    cmap.set_under('white')
    mp = MapPlot(sector='nws', continentalcolor='white', figsize=(12., 9.),
                 title=("2003-2018 Flash Flood Emergency Events"),
                 subtitle=('based on unofficial IEM archives, searching '
                           '"FFS", "FLW", "FFS".'))
    mp.fill_cwas(vals, bins=bins, lblformat='%s', labels=labels,
                 cmap=cmap, ilabel=True,  # clevlabels=month_abbr[1:],
                 units='count')
    mp.postprocess(filename='test.png')

def main():
    """Go Main"""
    pgconn = get_dbconn('postgis')
    df = read_postgis("""
    select geom, issue from bot_warnings where wfo = 'PUB'
    """, 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',
                 title='Bot Issued Tornado Warnings [2008-2018] for PUB',
                 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,
                         edgecolor='r', facecolor='None', alpha=1., lw=0.5,
                         zorder=10)
    mp.postprocess(filename='test.png')

def run(ts, routes):
    """ Run for a given UTC timestamp """
    fn = ts.strftime(("/mesonet/ARCHIVE/data/%Y/%m/%d/model/rtma/%H/"
                      "rtma.t%Hz.awp2p5f000.grib2"))
    if not os.path.isfile(fn):
        print 'wind_power.py missing', fn
        return

    grb = pygrib.open(fn)
    try:
        u = grb.select(name='10 metre U wind component')[0]
        v = grb.select(name='10 metre V wind component')[0]
    except:
        print('Missing u/v wind for wind_power.py\nFN: %s' % (fn,))
        return
    mag = (u['values']**2 + v['values']**2)**.5

    mag = (mag * 1.35)**3 * 0.002641
    # 0.002641

    lats, lons = u.latlons()
    lts = ts.astimezone(pytz.timezone("America/Chicago"))
    pqstr = ("plot %s %s00 midwest/rtma_wind_power.png "
             "midwest/rtma_wind_power_%s00.png png"
             ) % (routes, ts.strftime("%Y%m%d%H"), ts.strftime("%H"))
    m = MapPlot(sector='midwest',
                title=(r'Wind Power Potential :: '
                       '(speed_mps_10m * 1.35)$^3$ * 0.002641'),
                subtitle=('valid: %s based on NOAA Realtime '
                          'Mesoscale Analysis'
                          ) % (lts.strftime("%d %b %Y %I %p")))
    m.pcolormesh(lons, lats, mag, numpy.array(levels), units='MW')

    m.postprocess(pqstr=pqstr)

def run(ts, routes):
    """ Run for a given UTC timestamp """
    fn = ts.strftime("/mesonet/ARCHIVE/data/%Y/%m/%d/model/rtma/%H/rtma.t%Hz.awp2p5f000.grib2")
    if not os.path.isfile(fn):
        print "wind_power.py missing", fn
        return

    grb = pygrib.open(fn)
    u = grb.select(name="10 metre U wind component")[0]
    v = grb.select(name="10 metre V wind component")[0]
    mag = (u["values"] ** 2 + v["values"] ** 2) ** 0.5

    mag = (mag * 1.35) ** 3 * 0.002641
    # 0.002641

    lats, lons = u.latlons()
    lts = ts.astimezone(pytz.timezone("America/Chicago"))
    pqstr = "plot %s %s00 midwest/rtma_wind_power.png midwest/rtma_wind_power_%s00.png png" % (
        routes,
        ts.strftime("%Y%m%d%H"),
        ts.strftime("%H"),
    )
    m = MapPlot(
        sector="midwest",
        title=r"Wind Power Potential :: (speed_mps_10m * 1.35)$^3$ * 0.002641",
        subtitle="valid: %s based on NOAA Realtime Mesoscale Analysis" % (lts.strftime("%d %b %Y %I %p")),
    )
    m.pcolormesh(lons, lats, mag, numpy.array(levels), units="MW")

    m.postprocess(pqstr=pqstr)

def test_drawiowawfo():
    """Iowa Contour Plot"""
    mp = MapPlot(sector='iowawfo', title='Iowa Contour plot', nocaption=True)
    mp.contourf(np.arange(-94, -85), np.arange(36, 45), np.arange(9),
                np.arange(9),
                clevlabels=['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i'])
    return mp.fig