Python MapPlot.pcolormesh方法代码示例

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

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import pcolormesh [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 pcolormesh [as 别名]
def plot():
    """Do plotting work"""
    cmap = plt.get_cmap('inferno_r')
    # cmap.set_under('black')
    # cmap.set_over('red')
    minval = (np.load('minval.npy') * units.degK).to(units.degF)
    maxval = (np.load('maxval.npy') * units.degK).to(units.degF)
    diff = maxval - minval
    lons = np.load('lons.npy')
    lats = np.load('lats.npy')
    mp = MapPlot(sector='conus',
                 title=(r"Difference between warmest 3 Oct and coldest 4 "
                        "Oct 2m Temperature"),
                 subtitle=("based on hourly NCEP Real-Time Mesoscale Analysis "
                           "(RTMA) ending midnight CDT"))
    mp.ax.text(0.5, 0.97,
               (r"Pixel Difference Range: %.1f$^\circ$F to %.1f$^\circ$F, "
                r"Domain Analysis Range: %.1f$^\circ$F to %.1f$^\circ$F"
                ) % (np.min(diff).magnitude,
                     np.max(diff).magnitude,
                     np.min(minval).magnitude,
                     np.max(maxval).magnitude),
               transform=mp.ax.transAxes, fontsize=12, ha='center',
               bbox=dict(pad=0, color='white'), zorder=50)
    mp.pcolormesh(lons, lats, diff, range(0, 61, 5),
                  cmap=cmap, clip_on=False,
                  units=r"$^\circ$F")
    mp.postprocess(filename='test.png')

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import pcolormesh [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 pcolormesh [as 别名]
def compute(valid):
    ''' Get me files '''
    prob = None
    for hr in range(-15,0):
        ts = valid + datetime.timedelta(hours=hr)
        fn = ts.strftime("hrrr.ref.%Y%m%d%H00.grib2")
        if not os.path.isfile(fn):
            continue

        grbs = pygrib.open(fn)
        gs = grbs.select(level=1000,forecastTime=(-1 * hr * 60))
        ref = generic_filter(gs[0]['values'], np.max, size=10)
        if prob is None:
            lats, lons = gs[0].latlons()
            prob = np.zeros( np.shape(ref) )
        
        prob = np.where(ref > 29, prob+1, prob)

    prob = np.ma.array(prob / 15. * 100.)
    prob.mask = np.ma.where(prob < 1, True, False)    
    
    m = MapPlot(sector='iowa',
                title='HRRR Composite Forecast 4 PM 20 May 2014 30+ dbZ Reflectivity',
                subtitle='frequency of previous 15 model runs all valid at %s, ~15km smoothed' % (valid.astimezone(pytz.timezone("America/Chicago")).strftime("%-d %b %Y %I:%M %p %Z"),))

    m.pcolormesh(lons, lats, prob, np.arange(0,101,10), units='%',
                     clip_on=False)
    m.map.drawcounties()
    m.postprocess(filename='test.ps')
    m.close()

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

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import pcolormesh [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 pcolormesh [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 pcolormesh [as 别名]
def make_plots(nc):
    ''' Generate some plots '''
    sts = compute_sts(nc)
    lats = nc.variables['lat'][:]
    lons = nc.variables['lon'][:]
    rts = (sts.astimezone(pytz.timezone("America/Chicago"))).strftime(
                                                            "%d %b %Y %H %p")
    for i, tm in enumerate(nc.variables['time'][:]):
        dt = sts + datetime.timedelta(minutes=float(tm))
        if dt.minute != 0:
            continue
        fhour = int( tm / 60.0 )
        fts = (dt.astimezone(pytz.timezone("America/Chicago"))).strftime(
                                                            "%d %b %Y %H %p")
        for pvar in PVARS:
            m = MapPlot(title='ISUMM5/Bridget Modelled %s' % (
                                                    PVARS[pvar]['title'],),
                        subtitle='Model Run: %s Forecast Valid: %s' % (rts, fts))
            vals = nc.variables[pvar][i,:,:]
            if pvar == 'bdeckt':
                vals = temperature(vals, 'K').value('F')
            m.pcolormesh(lons, lats, vals, PVARS[pvar]['levels'], units='mm')
            pqstr = "plot c %s model/frost/bridget/%02i/%s_%02i_f%03i.png bogus png" % (
                                        sts.strftime("%Y%m%d%H%M"), sts.hour,
                                        pvar, sts.hour, fhour)
            m.postprocess(pqstr=pqstr)
            m.close()

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

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

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

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import pcolormesh [as 别名]
def test_pcolormesh():
    """See if we can do pcolormesh OKish"""
    mp = MapPlot(sector='custom', north=43, east=-80, west=-96,
                 south=38, projection=reference.EPSG[2163],
                 continentalcolor='white', nocaption=True)
    lons = np.arange(-100, -80, 0.25)
    lats = np.arange(40, 50, 0.25)
    vals = np.linspace(0, 1, lats.shape[0] * lons.shape[0]
                       ).reshape([lats.shape[0], lons.shape[0]])
    lons, lats = np.meshgrid(lons, lats)
    mp.pcolormesh(lons, lats, vals, np.arange(0, 1, 0.1))
    return mp.fig

# 需要导入模块: from pyiem.plot import MapPlot [as 别名]
# 或者: from pyiem.plot.MapPlot import pcolormesh [as 别名]
def main():
    """Go Main Go"""
    nc = ncopen("/mesonet/data/iemre/1979_narr.nc")

    data = np.sum(nc.variables['apcp'][:, :, :], axis=0)

    m = MapPlot(sector='conus', axisbg='tan',
                title=(""),
                subtitle='',
                titlefontsize=16)

    t = distance(data, 'MM').value('IN')
    m.pcolormesh(nc.variables['lon'][:], nc.variables['lat'][:], t,
                 np.arange(0, 60, 2), units='F')

    m.postprocess(filename='test.png')