Python Basemap.barbs方法代码示例

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
def test_plot_all_map_vels():
    def gen_tracks():
        for fn in glob('tracks/*.json'):
            with open(fn, 'rb') as fp:
                try:
                    yield RKJSON(fp)
                except BadInputException:
                    pass

    tracks = list(gen_tracks())

    print "Considering {} tracks".format(len(tracks))

    lllat, urlat, lllon, urlon = Grid.calculate_bounds(tracks)

    m = Basemap(projection='cyl',
                llcrnrlon=lllon, llcrnrlat=lllat,
                urcrnrlon=urlon, urcrnrlat=urlat,
                resolution='h')
    m.drawcoastlines()

    for t in tracks:
        vels = t.calculate_vels()
        x, y = m(vels.lon, vels.lat)
        m.barbs(x, y, vels.u, vels.v)

    plt.annotate("Docklands", (144.948, -37.815))
    plt.annotate("Brunswick", (144.960, -37.767))
    plt.annotate("Brunswick East", (144.979, -37.769))
    plt.annotate("Richmond", (144.999, -37.819))
    plt.annotate("Fitzroy", (144.978, -37.800))

    plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
def test_plot_map_vels(xml):
    track = GPX(xml)
    vels = track.calculate_vels(smooth_vels=True)

    lllat, lllon = min(track.lat) - 0.01, min(track.lon) - 0.01
    urlat, urlon = max(track.lat) + 0.01, max(track.lon) + 0.01

    # find the centre point
    lat_0 = (urlat - lllat) / 2 + lllat
    lon_0 = (urlon - lllon) / 2 + lllon

    # FIXME: rsphere required because my Proj is screwy
    m = Basemap(projection='cyl',
                llcrnrlon=lllon, llcrnrlat=lllat,
                urcrnrlon=urlon, urcrnrlat=urlat,
                lat_0=lat_0, lon_0=lon_0,
                resolution='h')
                # rsphere=(6378137.00, 6356752.3142))

    x, y = m(vels.lon, vels.lat)

    m.drawcoastlines()
    m.drawrivers()
    m.barbs(x, y, vels.u, vels.v) #, vels.anom, cmap=plt.get_cmap('RdBu_r'))
    plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
def test_plot_map_vels(json):
    track = RKJSON(json)
    vels = track.calculate_vels()

    lllat, lllon = min(track.lat) - 0.01, min(track.lon) - 0.01
    urlat, urlon = max(track.lat) + 0.01, max(track.lon) + 0.01

    print vels.bearing

    # find the centre point
    lat_0 = (urlat - lllat) / 2 + lllat
    lon_0 = (urlon - lllon) / 2 + lllon

    # FIXME: rsphere required because my Proj is screwy
    m = Basemap(projection='cyl',
                llcrnrlon=lllon, llcrnrlat=lllat,
                urcrnrlon=urlon, urcrnrlat=urlat,
                lat_0=lat_0, lon_0=lon_0,
                resolution='h')

    x, y = m(vels.lon, vels.lat)
    plt.annotate("Start", (x[0], y[0]))

    m.barbs(x, y, vels.u, vels.v) #, vels.anom, cmap=plt.get_cmap('RdBu_r'))
    plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
def plot_gini_grid_barbs(gini_grid, box = [-110, -70, 20, 52], resolution = 'l',
                   parallels = None,
                   meridians = None,
                   vmin = None, vmax = None,
                   fld = 'IR', title = None,
                   degrade = 5, u='u', v='v'):
    m = Basemap(llcrnrlon = box[0] ,llcrnrlat = box[2] , urcrnrlon = box[1],
                   urcrnrlat = box[3] , projection = 'mill', area_thresh =1000 ,
                   resolution = resolution)
    
    x, y = m(gini_grid.fields['lon']['data'], gini_grid.fields['lat']['data'])
    # create figure.
    if parallels == None:
        parallels = np.linspace(10,50, 9)
    if meridians == None:
        meridians = np.linspace(-110, -80,7)
    m.drawparallels(parallels, labels=[1,1,0,0])
    m.drawmeridians(meridians,labels=[0,0,0,1]) 
    pc = m.pcolormesh(x, y , gini_grid.fields[fld]['data'][0,:], cmap=plt.get_cmap('gray'),
                      vmin = vmin, vmax = vmax)
    qq = m.barbs(x[::degrade,::degrade], y[::degrade,::degrade], 
                  gini_grid.fields[u]['data'][0,::degrade,::degrade], 
                  gini_grid.fields[v]['data'][0,::degrade,::degrade])

    plt.title(title)
    m.drawcoastlines(linewidth=1.25)
    m.drawstates()
    plt.colorbar(mappable=pc, label = 'Counts')

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
def sfc_plot(starttime, endtime, variables, variablest, locations, 
             met, xi, yi, xmin, xmax, ymin, ymax):
    ''' Script for plotting the mesonet data with wind barbs over a 
    county map in a given time interval 
    '''
    interval = int((endtime - starttime).total_seconds()/300)
    z_max = np.max(met[variablest[1]])                                                                  
    z_min = np.min(met[variablest[1]])                                                                  
    levels = np.arange(z_min, z_max+0.1, 0.1)
    shapefile = 'UScounties/UScounties'
    if not os.path.exists('%s' %(variables)):
        os.makedirs('%s' %(variables))
    for i in range(interval):
        time_selection = starttime + dt.timedelta(minutes=5*i)
        zi = interpolate.griddata((met.ix[time_selection]['Lon'], 
                                   met.ix[time_selection]['Lat']), 
                                   met.ix[time_selection][variablest[1]],
                                   (xi, yi), method='linear')
        maps = Basemap(llcrnrlon=xmin, llcrnrlat=ymin, 
                       urcrnrlon=xmax, urcrnrlat=ymax, projection='cyl')
        maps.readshapefile(shapefile, name='counties')
        if (variables == 'dew_point'):
            maps.contourf(xi, yi, zi, levels, cmap=plt.cm.gist_earth_r)
        if (variables == 'temperature'):
            maps.contourf(xi, yi, zi, levels, cmap=plt.cm.jet)
        if variables == 'rainfall':
            maps.contourf(xi, yi, zi, levels, cmap=plt.cm.YlGn)
        if ((variables == 'pressure') or (variables == 'wind_speed') or 
            (variables == 'gust_speed')):
            maps.contourf(xi, yi, zi, levels, cmap=plt.cm.gist_earth)
        c = plt.colorbar()  
        c.set_label(variablest[0])  
        maps.scatter(met.ix[time_selection]['Lon'], 
                     met.ix[time_selection]['Lat'], marker='o', c='b', s=5)
        maps.barbs(met.ix[time_selection]['Lon'], 
                   met.ix[time_selection]['Lat'], 
                   met.ix[time_selection]['u'].values, 
                   met.ix[time_selection]['v'].values)
        plt.title(variablest[1]) 
        filename = '%s_%s.png' % (variables, 
                                  time_selection.strftime('%Y%m%d_%H%M'))
        plt.tight_layout()
        plt.savefig(variables + '/' + filename, dpi=150)
        plt.clf()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
def CO2nWind(file,outputdir):
  """
  Outputs plots of CO2 in the lowest level and 10m winds
  """
  import matplotlib
  matplotlib.use('Agg')
  import numpy as np
  import sys
  from netCDF4 import Dataset
  import matplotlib.pyplot as plt
  from mpl_toolkits.basemap import Basemap
  # read in file and vars
  f = Dataset(file,'r')
  wrfeta = f.variables['ZNU'][0][:]
  times = f.variables['Times'][:]
  wrflats = f.variables['XLAT'][0][:]
  wrflons = f.variables['XLONG'][0][:]
  var = f.variables['CO2_ANT'][:,0,:,:]
  u = f.variables['U'][:,0,:,:]
  v = f.variables['V'][:,0,:,:]
  # destagger u/v
  u = (u[:,:,:-1] + u[:,:,1:])/2.
  v = (v[:,:-1,:] + v[:,1:,:])/2.

  # four corners of domain
  wrflat_s = wrflats[0,len(wrflats)-1]
  wrflat_n = wrflats[len(wrflats)-1,len(wrflons[0])-1]
  wrflon_w = wrflons[0,0]
  wrflon_e = wrflons[len(wrflats)-1,len(wrflons[0])-1]

  z = 0 # assuming lowest level of model

  # set up map
  map = Basemap(projection='merc',llcrnrlon=wrflon_w,urcrnrlon=wrflon_e,llcrnrlat=wrflat_s,urcrnrlat=wrflat_n,resolution='i')
  map.drawstates()
  map.drawcounties()
  map.drawcoastlines()
  x,y = map(wrflons,wrflats)

  # loop through times
  for t in range(len(times)):
    timestr = ''.join(times[t,:])
    map.drawstates(color='gray',linewidth=1)
    map.drawcounties(color='white')
    map.drawcoastlines(color='gray',linewidth=1)
    plt1 = map.pcolormesh(x,y,var[t,:,:],vmin=380,vmax=450)
    #plt1 = map.pcolormesh(x,y,var[t,:,:],vmin=np.amin(var),vmax=np.amax(var))
    winds = map.barbs(x[::20,::20],y[::20,::20],u[t,::20,::20]*1.94,v[t,::20,::20]*1.94,length=6,color='white') # *1.94 to convert m/s to knots (barb convention)
    colorbar = map.colorbar(plt1,"right", size="5%",pad="2%")
    colorbar.set_label(f.variables['CO2_ANT'].description+' '+f.variables['CO2_ANT'].units)
    plt.title('WRF output valid: '+timestr)
    plt.savefig(outputdir+'/%03d_' % (t) +timestr+'_CO2_wind.png')
    plt.clf()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
def plot_gridded_ensemble(gfsx, contour_units=None, max_level=None,
        barb_units='knot', cols=2, save_path=None, save_fmt='svg'):
    """
    Plots the average windspeed deviation of ensemble forecast members over
    the published GSF forecast for each grid point and forecast time.
    The result is (for each forecast time) a color filled contour plot
    indicating the ensemble deviation superimposed over the usual wind barbs
    with the published GFS forecast vectors.

    Parameters:
    -----------
    fcst_gfsx: xray.Dataset
        Dataset with the GFS forecast for windspeed (uwnd and vwnd) and
        ensemble forecast deviation (other forecast variables will be ignored
        if present).
    contour_units: str
        The units in which to plot the ensemble spread indicator (heatmap). If
        specified, the data sets current unit must be convertible into the
        desired unit by slocum.lib.units. If not specified the exsiting units in
        the data set will be used.
    max_level: float
        Heatmap/contour plot maximum level. Minimum level is assumed to be 0.
        Values greater tha max_level will be mapped onto the max_level color.
        If not specified, the maximum value found in the data set will be used
        as the upper end of the color bar.
    barb_units: str
        Units in which to plot the GFS wind speed/direction forecast data (wind
        bars). If specified, the data sets current unit must be convertible
        into the desired unit by slocum.lib.units.
    cols: int
        Number of columns for subplot layout.
    save_path: str
        If specified the plot will be saved into this directory with file
        name ``ens_<bounding box>_<t0>.svg`` where *bounding box* is
        specified as *ll_lat, ll_lon - ur_lat, ur_lon*.  If not specified or
        ``None`` the plot will be displayed interactively.
    save_fmt: str
        Format under which to save the image (only relevant if *save_path* is
        specified). Can be any image file extension that plt.savefig() will
        recognize as a valid format.
    """
    # adjust units as requested:
    for v in (conv.UWND, conv.VWND):
        units.convert_units(gfsx[v], barb_units)
    if contour_units:
        units.convert_units(gfsx[conv.ENS_SPREAD_WS], contour_units)
    if not max_level:
        max_level = gfsx[conv.ENS_SPREAD_WS].max()

    if isinstance(gfsx, np.datetime64): # True is gfsx has not been packed
        f_times = gfsx[conv.TIME].values
    else:                               # time variable has int offsets
        f_times = xray.conventions.decode_cf_datetime(
                gfsx['time'], gfsx['time'].attrs['units'])

    lats = gfsx[conv.LAT].values
    lons = gfsx[conv.LON].values

    # Basemap.transform_vector requires lats and lons each to be in ascending
    # order:
    lat_inds = range(len(lats))
    lon_inds = range(len(lons))
    if NautAngle(lats[0]).is_north_of(NautAngle(lats[-1])):
        lat_inds.reverse()
        lats = lats[lat_inds]
    if NautAngle(lons[0]).is_east_of(NautAngle(lons[-1])):
        lon_inds.reverse()
        lons = lons[lon_inds]

    # determine layout:
    fig_width_inches = 10
    plot_aspect_ratio = (abs(lons[-1] - lons[0]) /
                         float(abs(lats[-1] - lats[0])))
    rows = int(np.ceil(gfsx.dimensions[conv.TIME] / float(cols)))
    fig_height_inches = (fig_width_inches / (float(cols) * plot_aspect_ratio)
                         * rows + 2)
    fig = plt.figure(figsize=(fig_width_inches, fig_height_inches))
    fig.suptitle("GFS wind forecast in %s and "
                 "ensemble wind speed deviation" % barb_units,
                 fontsize=12)
    grid = AxesGrid(fig, [0.01, 0.01, 0.95, 0.93],
                    nrows_ncols=(rows, cols),
                    axes_pad=0.8,
                    cbar_mode='single',
                    cbar_pad=0.0,
                    cbar_size=0.2,
                    cbar_location='bottom',
                    share_all=True,)

    # size for lat/lon labels, timestamp:
    label_fontsize = 'medium' if cols <= 2 else 'small'
    # format string for colorbar labels:
    decimals = max(0, int(2 - np.floor(np.log10(max_level))))
    cb_label_fmt = '%.' + '%d' % decimals + 'f'

    # heatmap color scaling and levels
    spread_levels = np.linspace(0., max_level, 50)

    m = Basemap(projection='merc', llcrnrlon=lons[0], llcrnrlat=lats[0],
            urcrnrlon=lons[-1], urcrnrlat=lats[-1], resolution='l')
#.........这里部分代码省略.........

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
lons, lats = np.meshgrid(lons, lats)

v10 = np.ones((lons.shape)) * 15
u10 = np.zeros((lons.shape))

u10_rot, v10_rot, x, y = map.rotate_vector(u10, v10, lons, lats, returnxy=True)

ax = fig.add_subplot(121)
ax.set_title('Without rotation')

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#cc9955', lake_color='aqua', zorder = 0)
map.drawcoastlines(color = '0.15')


map.barbs(x, y, u10, v10, 
    pivot='middle', barbcolor='#333333')

ax = fig.add_subplot(122)
ax.set_title('Rotated vectors')

map.drawmapboundary(fill_color='9999FF')
map.fillcontinents(color='#ddaa66', lake_color='9999FF', zorder = 0)
map.drawcoastlines(color = '0.15')

map.barbs(x, y, u10_rot, v10_rot, 
    pivot='middle', barbcolor='#ff7777')

plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
norm = mpl.colors.Normalize(vmin=-35, vmax=80)
ax.text(0.5, .97, title, transform=plt.gca().transAxes,fontsize=16, horizontalalignment='center', bbox=dict(facecolor='white', alpha=0.9, boxstyle='round'))                                                     
#cax = m.pcolormesh(x, y, data, cmap=cmap, norm=norm)
cax = m.imshow(data, extent = (x.min(), x.max(), y.min(), y.max()), cmap=cmap, norm=norm, origin="upper")
m.drawcoastlines(color='#FFFFFF')
m.drawstates(color='#FFFFFF')
m.drawcountries(color='#FFFFFF')
m.drawcounties(color='#FFFFFF', linewidth=0.09)
cbar = m.colorbar(cax)

x_grid, y_grid = m(cape_grid_lon, cape_grid_lat)
c = m.contour(x_grid, y_grid, CAPE, np.arange(500,6500,500), cmap='spring_r')
plt.clabel(c, fmt='%4.0f')
stride = 10
#m.barbs(x_grid[::stride, ::stride], y_grid[::stride, ::stride], u_shear[::stride,::stride], v_shear[::stride,::stride], mag[::stride,::stride], cmap='RdPu') 
m.barbs(grid_lon, grid_lat, u_shear.T, v_shear.T, mag.T, clim=[30,70], latlon=True, cmap='RdPu') 

plt.tight_layout()
# Make the newest graphic.
out_filename = path_to_imgs + 'nexrad_rap_composite_' + today.strftime("%Y%m%d_%H%M.png")
plt.savefig(out_filename, bbox_inches='tight')

# Delete the oldest graphic.
old_filename = path_to_imgs + 'nexrad_rap_composite_' + end.strftime("%Y%m%d_%H%M.png")
os.system('/usr/bin/rm ' + old_filename)

print "Making looper"
path_to_looper = '/data/soundings/blumberg/programs/wx_maps/looper/make_looper.py'
os.system('/data/soundings/anaconda/bin/python ' + path_to_looper + ' ' + path_to_imgs + ' nexrad_rap')

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
qk = plt.quiverkey(Q, 0.1, 0.1, 20, '20 m/s', labelpos='W')
# draw coastlines, parallels, meridians.
m.drawcoastlines(linewidth=1.5)
m.drawparallels(parallels)
m.drawmeridians(meridians)
# add colorbar
cb = m.colorbar(CS2,"bottom", size="5%", pad="2%")
cb.set_label('hPa')
# set plot title
ax.set_title('SLP and Wind Vectors '+str(date))
fig1.savefig('plotwindvec1.png')

# create 2nd figure, add axes
fig2 = plt.figure(figsize=(8,10))
ax = fig2.add_axes([0.1,0.1,0.8,0.8])
# plot SLP contours
CS1 = m.contour(x,y,slp,clevs,linewidths=0.5,colors='k',animated=True)
CS2 = m.contourf(x,y,slp,clevs,cmap=plt.cm.RdBu_r,animated=True)
# plot wind barbs over map.
barbs = m.barbs(xx,yy,uproj,vproj,length=5,barbcolor='k',flagcolor='r',linewidth=0.5)
# draw coastlines, parallels, meridians.
m.drawcoastlines(linewidth=1.5)
m.drawparallels(parallels)
m.drawmeridians(meridians)
# add colorbar
cb = m.colorbar(CS2,"bottom", size="5%", pad="2%")
cb.set_label('hPa')
# set plot title.
ax.set_title('SLP and Wind Barbs '+str(date))
fig2.savefig('plotwindvec2.png')

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
              lat_0=0, lon_0=0)

lons = np.linspace(-180, 180, 8)
lats = np.linspace(-90, 90, 8)

v10 = np.ones((lons.shape)) * 15
u10 = np.zeros((lons.shape))

u10, v10 = np.meshgrid(u10, v10)

u10_rot, v10_rot, x_rot, y_rot = map.transform_vector(u10, v10, 
                                                      lons, lats, 
                                                      15, 15, 
                                                      returnxy=True)

map.drawmapboundary(fill_color='#9999FF')
map.fillcontinents(color='#ddaa66', lake_color='#9999FF', zorder = 0)
map.drawcoastlines(color = '0.15')

#Drawing the original points
lons, lats = np.meshgrid(lons, lats)
x, y = map(lons, lats)
map.barbs(x, y, u10, v10, 
    pivot='middle', barbcolor='#333333')

#Drawing the rotated & interpolated points
map.barbs(x_rot, y_rot, u10_rot, v10_rot, 
    pivot='middle', barbcolor='#ff5555')

plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
from osgeo import gdal
import numpy as np


map = Basemap(llcrnrlon=-93.7, llcrnrlat=28., urcrnrlon=-66.1, urcrnrlat=39.5,
              projection='lcc', lat_1=30., lat_2=60., lat_0=34.83158, lon_0=-98.)

ds = gdal.Open("../sample_files/wrf.tiff")
lons = ds.GetRasterBand(4).ReadAsArray()
lats = ds.GetRasterBand(5).ReadAsArray()
u10 = ds.GetRasterBand(1).ReadAsArray()
v10 = ds.GetRasterBand(2).ReadAsArray()

x, y = map(lons, lats)

yy = np.arange(0, y.shape[0], 4)
xx = np.arange(0, x.shape[1], 4)

points = np.meshgrid(yy, xx)

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#cc9955', lake_color='aqua', zorder = 0)
map.drawcoastlines(color = '0.15')

map.barbs(x[points], y[points], u10[points], v10[points], 
    pivot='middle', barbcolor='#333333')


plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
lon = np.squeeze(file.variables['XLONG'][:,:]);
lat = np.squeeze(file.variables['XLAT'][:,:]);

U,V = np.multiply(U,1.94384449),np.multiply(V,1.94384449); 
XA,YA = m(lon,lat); 

os.chdir('/home/disk/manta14/awsmith/umcm/umcm_runs/katrina_awo_gfs10/output/')
file = Dataset('umwmout_2005-08-29_12:00:00.nc','r'); 
SWH = np.squeeze(file.variables['swh'][:,:]);
lonw = np.squeeze(file.variables['lon'][:,:]); 
latw = np.squeeze(file.variables['lat'][:,:]); 

XW,YW = m(lonw,latw); 

m.contourf(XW,YW,SWH,np.arange(0,31,1),cmap='jet')
m.barbs(XA[0:-1:8,0:-1:8],YA[0:-1:8,0:-1:8],U[0:-1:8,0:-1:8],V[0:-1:8,0:-1:8],barbcolor='k',flagcolor='r',linewidth=0.5)
plt.title('Hurricane Katrina (2005) SWH (m) & Surface Winds (kt), August 29 1200 UTC',fontsize=8)
## Add land mask over 
m.fillcontinents(color='w')
m.drawcountries()
m.drawstates()

## Add ocean mask over
plt.figure(num=5,figsize=(6,6))
m = Basemap(llcrnrlon=-95.,llcrnrlat=20.,urcrnrlon=-80.,urcrnrlat=35.,projection='lcc',lat_1=20.,lat_2=40.,lon_0=-60.,resolution ='l',area_thresh=1000.)
m.drawcoastlines()
m.drawcountries()
m.drawparallels(np.arange(20,36,2),labels=[1,0,0,0])
m.drawmeridians(np.arange(-95,-78,2),labels=[0,0,0,1])
import mpl_toolkits.basemap
REFD_OM = mpl_toolkits.basemap.maskoceans(lon,lat,REFD,'True','l');

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
map = Basemap(llcrnrlon=-82., llcrnrlat=28., urcrnrlon=-79., urcrnrlat=29.5,
              projection='lcc', lat_1=30., lat_2=60., lat_0=34.83158, lon_0=-98.,
              resolution='i')

ds = gdal.Open("../sample_files/wrf.tiff")
lons = ds.GetRasterBand(4).ReadAsArray()
lats = ds.GetRasterBand(5).ReadAsArray()
u10 = ds.GetRasterBand(1).ReadAsArray()
v10 = ds.GetRasterBand(2).ReadAsArray()

x, y = map(lons, lats)

x2 = np.linspace(x[0][0],x[0][-1],x.shape[1]*2)
y2 = np.linspace(y[0][0],y[-1][0],y.shape[0]*2)

x2, y2 = np.meshgrid(x2, y2)

u10_2 = interp(u10,  x[0], np.flipud(y[:, 0]), x2, np.flipud(y2),order=1)
v10_2 = interp(v10,  x[0], np.flipud(y[:, 0]), x2, np.flipud(y2),order=1)

map.drawmapboundary(fill_color='#9999FF')
map.fillcontinents(color='#cc9955', lake_color='#9999FF', zorder = 0)
map.drawcoastlines(color = '0.15')

map.barbs(x, y, u10, v10, 
    pivot='middle', barbcolor='#555555')

map.barbs(x2, y2, u10_2, v10_2, 
    pivot='middle', barbcolor='#FF3333')

plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import barbs [as 别名]
    llcrnrlon=bot_left_lon,llcrnrlat=bot_left_lat,\
    urcrnrlon=top_right_lon,urcrnrlat=top_right_lat,)

plt.figure(figsize=[16*1.5,10*1.5])
m.drawcoastlines()
m.drawcountries()
m.drawstates()

x,y = m(lon,lat)    

# Now plot wind and geopotential height
m.pcolormesh(x,y,speed,cmap="gist_stern_r")
cb = plt.colorbar(shrink=.7)
cb.set_label("Wind Speed (m/s)")

interval = 50
m.barbs(x[::50,::50],y[::50,::50],U[::50,::50],V[::50,::50])  

# Just becuase we can, add 500 mb Geopotential Height Contours on top.
CS = m.contour(x,y,Z,colors='black',levels=np.arange(0,6000,60),linewidths=3)
plt.clabel(CS, inline=1, fontsize=12,fmt='%1.0f')      


validDate = grb.validDate

plt.title('HRRR Analysis Valid '+validDate.strftime('%Y-%b-%d %H:%M:%S')+'\n'+
          str(grb.level)+' mb '+Zname+'and Wind')

plt.savefig('HRRR_anl_'+str(mb)+'.png',dpi=300,bbox_inches='tight')
plt.show()