Python Basemap.contour方法代码示例

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import contour [as 别名]
lons, lats = P.meshgrid(lons, lats)

# Example to show how to make multi-panel plots.

# 2-panel plot, oriented vertically, colorbar on bottom.

rcParams["figure.subplot.hspace"] = 0.4  # more height between subplots
rcParams["figure.subplot.wspace"] = 0.5  # more width between subplots

# create new figure
fig = P.figure()
# panel 1
mnh = Basemap(lon_0=-105, boundinglat=20.0, resolution="c", area_thresh=10000.0, projection="nplaea")
xnh, ynh = mnh(lons, lats)
ax = fig.add_subplot(211)
CS = mnh.contour(xnh, ynh, hgt, 15, linewidths=0.5, colors="k")
CS = mnh.contourf(xnh, ynh, hgt, 15, cmap=P.cm.Spectral)
# colorbar on bottom.
l, b, w, h = ax.get_position()
cax = P.axes([l, b - 0.05, w, 0.025])  # setup colorbar axes
P.colorbar(cax=cax, orientation="horizontal", ticks=CS.levels[0::4])  # draw colorbar
P.axes(ax)  # make the original axes current again
mnh.drawcoastlines(linewidth=0.5)
delat = 30.0
circles = P.arange(0.0, 90.0, delat).tolist() + P.arange(-delat, -90, -delat).tolist()
mnh.drawparallels(circles, labels=[1, 0, 0, 0])
delon = 45.0
meridians = P.arange(0, 360, delon)
mnh.drawmeridians(meridians, labels=[1, 0, 0, 1])
P.title("NH 500 hPa Height (cm.Spectral)")

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import contour [as 别名]
# read in data on lat/lon grid.
hgt = load('500hgtdata.gz')
lons = load('500hgtlons.gz')
lats = load('500hgtlats.gz')
# shift data so lons go from -180 to 180 instead of 0 to 360.
hgt,lons = shiftgrid(180.,hgt,lons,start=False)
lons, lats = meshgrid(lons, lats)

# create new figure
fig=figure()
# setup of sinusoidal basemap
m = Basemap(resolution='c',projection='sinu',lon_0=0)
ax = fig.add_axes([0.1,0.1,0.7,0.7])
# make a filled contour plot.
x, y = m(lons, lats)
CS = m.contour(x,y,hgt,15,linewidths=0.5,colors='k')
CS = m.contourf(x,y,hgt,15,cmap=cm.jet)
l,b,w,h=ax.get_position()
cax = axes([l+w+0.075, b, 0.05, h]) # setup colorbar axes
colorbar(drawedges=True, cax=cax) # draw colorbar
axes(ax)  # make the original axes current again
# draw coastlines and political boundaries.
m.drawcoastlines()
m.drawmapboundary()
m.fillcontinents()
# draw parallels and meridians.
parallels = arange(-60.,90,30.)
m.drawparallels(parallels,labels=[1,0,0,0])
meridians = arange(-360.,360.,30.)
m.drawmeridians(meridians)
title('Sinusoidal Filled Contour Demo')

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import contour [as 别名]
meridians = N.arange(-120., -60.0, 10.0)
fields = (cress1, cress2, cress3, cress_187, barnes_analyses[0],
  barnes_analyses[1], barnes_analyses[2], barnes_3pass)
names = ('Cressman 1st Pass', 'Cressman 2nd Pass', 'Cressman 3rd Pass',
  'Cressman R=1.87d', r'$\rm{Barnes} \gamma=1.0$',r'$\rm{Barnes} \gamma=0.4$',
  r'$\rm{Barnes} \gamma=0.2$', r'$\rm{Barnes 3-pass} \gamma=1.0$')
filenames = ('cress1', 'cress2', 'cress3', 'cress187', 'barnes1', 'barnes2',
  'barnes3', 'barnes3pass')
for field,name,filename in zip(fields, names, filenames):
  f = P.figure()
  bm_ps.drawstates()
  bm_ps.drawcountries()
  bm_ps.drawcoastlines()
  bm_ps.drawparallels(parallels, labels = [1,1,0,0])
  bm_ps.drawmeridians(meridians, labels = [0,0,1,1])
  cp = bm_ps.contour(x_bm, y_bm, field, contours, cmap=colormap)
  P.clabel(cp, fontsize=fontsize, fmt='%.1f')
  f.text(0.5,0.95,name,horizontalalignment='center',fontsize=16)
  if save_work:
    P.savefig(filename + '.png', dpi = 300)
  else:
    P.show()

filename = 'cress12diff'
f = P.figure()
bm_ps.drawstates()
bm_ps.drawcountries()
bm_ps.drawcoastlines()
bm_ps.drawparallels(parallels, labels = [1,1,0,0])
bm_ps.drawmeridians(meridians, labels = [0,0,1,1])
cp = bm_ps.contour(x_bm, y_bm, cress1 - cress2, cmap=colormap)

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import contour [as 别名]
# draw the edge of the map projection region (the projection limb)
map.drawmapboundary()
# draw lat/lon grid lines every 30 degrees.
map.drawmeridians(p.arange(0,360,30))
map.drawparallels(p.arange(-90,90,30))
# lat/lon coordinates of five cities.
lats=[40.02,32.73,38.55,48.25,17.29]
lons=[-105.16,-117.16,-77.00,-114.21,-88.10]
cities=['Boulder, CO','San Diego, CA',
        'Washington, DC','Whitefish, MT','Belize City, Belize']
# compute the native map projection coordinates for cities.
x,y = map(lons,lats)
# plot filled circles at the locations of the cities.
map.plot(x,y,'bo')
# plot the names of those five cities.
for name,xpt,ypt in zip(cities,x,y):
    p.text(xpt+50000,ypt+50000,name,fontsize=9)
# make up some data on a regular lat/lon grid.
nlats = 73; nlons = 145; delta = 2.*p.pi/(nlons-1)
lats = (0.5*p.pi-delta*p.indices((nlats,nlons))[0,:,:])
lons = (delta*p.indices((nlats,nlons))[1,:,:])
wave = 0.75*(p.sin(2.*lats)**8*p.cos(4.*lons))
mean = 0.5*p.cos(2.*lats)*((p.sin(2.*lats))**2 + 2.)
# compute native map projection coordinates of lat/lon grid.
x, y = map(lons*180./p.pi, lats*180./p.pi)
# contour data over the map.
cs = map.contour(x,y,wave+mean,15,linewidths=1.5)
p.show()
#p.savefig('wiki_example.ps')

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import contour [as 别名]
fig = p.figure()
ax = fig.add_axes([0.1,0.1,0.7,0.7])
# set desired contour levels.
clevs = p.arange(960,1061,5)
# compute native x,y coordinates of grid.
x, y = m(lons, lats)
# define parallels and meridians to draw.
parallels = p.arange(-80.,90,20.)
meridians = p.arange(0.,360.,20.)
# number of repeated frames at beginning and end is n1.
nframe = 0; n1 = 10
l,b,w,h=ax.get_position()
# loop over times, make contour plots, draw coastlines, 
# parallels, meridians and title.
for nt,date in enumerate(datelabels[1:]):
    CS = m.contour(x,y,slp[nt,:,:],clevs,linewidths=0.5,colors='k',animated=True)
    CS = m.contourf(x,y,slp[nt,:,:],clevs,cmap=p.cm.RdBu_r,animated=True)
    # plot wind vectors on lat/lon grid.
    # rotate wind vectors to map projection coordinates.
    #urot,vrot = m.rotate_vector(u[nt,:,:],v[nt,:,:],lons,lats)
    # plot wind vectors over map.
    #Q = m.quiver(x,y,urot,vrot,scale=500) 
    # plot wind vectors on projection grid (looks better).
    # first, shift grid so it goes from -180 to 180 (instead of 0 to 360
    # in longitude).  Otherwise, interpolation is messed up.
    ugrid,newlons = shiftgrid(180.,u[nt,:,:],longitudes,start=False)
    vgrid,newlons = shiftgrid(180.,v[nt,:,:],longitudes,start=False)
    # transform vectors to projection grid.
    urot,vrot,xx,yy = m.transform_vector(ugrid,vgrid,newlons,latitudes,51,51,returnxy=True,masked=True)
    # plot wind vectors over map.
    Q = m.quiver(xx,yy,urot,vrot,scale=500)

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import contour [as 别名]
# Transform the grid to basemap space for plotting
cm_per_m = 100.0
x_bm,y_bm = to_basemap_xy(bm_ps, x_grid / cm_per_m, y_grid / cm_per_m,
  -110.0, 90.0, 1.0/15000000.0)

#Specify contour intervals, as well as what parallels and meridians to show
contours = N.arange(5300., 6000., 60.0)
parallels = N.arange(25., 60.0, 10.0)
meridians = N.arange(-120., -60.0, 10.0)

#These draw the various maps (self explanatory)
bm_ps.drawstates()
bm_ps.drawcountries()
bm_ps.drawcoastlines()

#Draw grid lines, labels specify where to label the grid lines
#[left, right, top, bottom]
bm_ps.drawparallels(parallels, labels = [1,1,0,0])
bm_ps.drawmeridians(meridians, labels = [0,0,1,1])

#Remove the cmap argument below to get the normal, colored colortable
cp = bm_ps.contour(x_bm, y_bm, data, contours, cmap=M.cm.get_cmap('gist_gray'))

#Label contours
P.clabel(cp)
#Normally just P.title('My plot'), but the longitude labels get overlapped
P.gcf().text(0.5,0.95,'My Plot',horizontalalignment='center',fontsize=16)

#P.savefig('plot.png', dpi = 300) # This saves to a file
P.show()

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import contour [as 别名]
contours = N.arange(5300., 6000., 60.0)
parallels = N.arange(25., 60.0, 10.0)
meridians = N.arange(-120., -60.0, 10.0)
fields = (heights_uniform, heights_cress, heights_uniforbm_ps, heights_cress2)
names = ('Uniform Weights', 'Cressman Weights',\
  'Uniform weights with wind correction',\
  'Cressman weights with wind correction')
filenames = ('uniform1','cress1','uniform2','cress2')
for field,name,filename in zip(fields, names, filenames):
  f = P.figure()
  bm_ps.drawstates()
  bm_ps.drawcountries()
  bm_ps.drawcoastlines()
  bm_ps.drawparallels(parallels, labels = [1,1,0,0])
  bm_ps.drawmeridians(meridians, labels = [0,0,1,1])
  cp = bm_ps.contour(x_bm, y_bm, field, contours, cmap=colormap)
  P.clabel(cp, fontsize=fontsize, fmt='%.1f')
  f.text(0.5,0.95,name,horizontalalignment='center',fontsize=16)
  if save_work:
    P.savefig(filename + '.png', dpi = 300)
  else:
    P.show()
  height_diff = list()
  if save_work:
    sys.stdout = open(filename + '.txt', 'w')
  print name
  print "\nAnalysis Array:"
  for n in N.arange(y.size):
    print "\t%d" % (n + 1),
  print
  for (num, row) in enumerate(field.T):

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import contour [as 别名]
def plot_ncdf_output(id, datafn, **kwds):
    """Plot model field id from the data in netCDF file datafn.

    Positional Input Parameter:
    * id:  Name of the id of the field to plot.  String.

    * datafn:  Filename containing the output data to plot.  String.

    Input keyword parameter descriptions are found in the docstring
    for Qtcm methods ploti, plotm, and other methods that call this
    private method.  In general, those methods take the keyword
    parameters they receive and pass it along unchanged as keyword
    parameters to this function.  In that sense, this function is
    seldom used as a stand-alone function, but rather is usually
    used coupled with a Qtcm instance.

    The data fields read in from the netCDF output file are dimensioned
    (time, lat, lon).  This is different than how the data is stored
    in the compiled QTCM model fields (lon, lat, time), and at the
    Python level (lon, lat).  The reason this is the case is that
    f2py automatically makes the arrays passed between the Python
    and Fortran levels match.

    For a lat vs. lon plot, the contour plot is superimposed onto
    a cylindrical projection map of the Earth with continents drawn
    and labeled meridians and parallels.  The title also includes
    the model time, and x- and y-axis labels are not drawn.

    All numerical data used for plotting come from the netCDF output
    file for consistency (e.g., the dimensions of u1).  Currently
    this method only works for 3-D data arrays (two in space, one
    in time).
    """
    #- Accomodate other ids.  The id that this routine will use
    #  (i.e., iduse corresponds to the name in the netCDF output file)
    #  is called iduse.  Set this to id, except for the case where
    #  some aliases of ids are entered in which iduse is the alias:

    if id == 'Qc': iduse = 'Prec'
    elif id == 'FLWut': iduse = 'OLR'
    elif id == 'STYPE': iduse = 'stype'
    else: iduse = id


    #- Set defined keyword defaults.  All are set to None except for
    #  nlatlon which gets an integer:

    plotkwds_ids = ['lat', 'lon', 'time', 'fn', 'levels', 'title',
                    'xlabel', 'ylabel', 
                    'filled', 'nlatlon', 'tmppreview']

    plotkwds = {}
    for ikey in plotkwds_ids:
        if kwds.has_key(ikey):
            plotkwds[ikey] = copy.copy(kwds[ikey])
        else:
            plotkwds[ikey] = None

    if not kwds.has_key('nlatlon'):
        plotkwds['nlatlon'] = 8


    #- Get data and dimensions of iduse to plot:

    fileobj = S.NetCDFFile(datafn, mode='r')
    data = N.array(fileobj.variables[iduse].getValue())
    data_name = fileobj.variables[iduse].long_name
    data_units = fileobj.variables[iduse].units

    dim = {}
    dimname = {}
    dimunits = {}

    dim['lat'] = N.array(fileobj.variables['lat'].getValue())
    dimname['lat'] = fileobj.variables['lat'].long_name
    dimunits['lat'] = fileobj.variables['lat'].units

    dim['lon'] = N.array(fileobj.variables['lon'].getValue())
    dimname['lon'] = fileobj.variables['lon'].long_name
    dimunits['lon'] = fileobj.variables['lon'].units
    
    dim['time'] = N.array(fileobj.variables['time'].getValue())
    dimname['time'] = fileobj.variables['time'].long_name
    dimunits['time'] = fileobj.variables['time'].units

    fileobj.close()


    #- Alter data long name to remove any units.  The definition
    #  of units as the substring within the [] is the same as used in
    #  defVar in output.F90 of the compiled QTCM model.  Remove 
    #  underscores and extra whitespace in data_name and data_units, 
    #  replacing with a single whitespace character between words:

    idx1 = data_name.find('[')
    idx2 = data_name.find(']')
    if idx1 != -1 and idx2 != -1:
        data_name = data_name[:idx1] + data_name[idx2+1:]
    data_name = data_name.strip()

#.........这里部分代码省略.........

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import contour [as 别名]
			'green': ((0.,)*3,(.5,1.,1.),(1.,0.,0.)),
			'blue': ((0.,1.,1.),(.5,1.,1.),(1.,0.,0.))}
		return matplotlib.colors.LinearSegmentedColormap('bwr',cdict,256)

	figure(figsize=(8,8))
	subplot(ncol,nrow,nslice)
	for i in xrange(nslice):
		lon = composites.getLongitude()
		lat = composites.getLatitude()
		xx,yy = meshgrid(lon,lat)
		m = Basemap(resolution='l',lat_0=N.average(lat),lon_0=N.average(lon),
			llcrnrlon=min(lon),llcrnrlat=min(lat),
			urcrnrlon=max(lon),urcrnrlat=max(lat))
		levels = vcs.mkscale(genutil.minmax(composites))
		m.contourf(lon,lat,composites[i],levels=levels,cmap=bwr)
		clabel(m.contour(lon,lat,composites[i],levels=levels))
		title("Phase %i/%i" (i+1,6))
		m.drawcoastlines()
		m.drawcountries()
		m.fillcontinents(color='coral')
		m.drawparallels(vcs.mkscale(genutil.minmax(lat)),labels=[1,0,0,0])
		m.drawmeridians(vcs.mkscale(genutil.minmax(lon)),labels=[0,0,0,1])
	figtext(.5,1.,"\n%s [%s]" % ("El Nino phase composites\nSea surface temperature", composites.units))
	savefig(sys.argv[0].replace(".py",".png"))
	show()
	

# Fall back to vcs because we have it!
except:
	# Plot 1 phase over two, then a time series
	# TODO: we must do something nicer!!

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import contour [as 别名]
plots = ['contour','pcolor']
#plots = ['contour','imshow']

for np,plot in enumerate(plots):

    fig.add_subplot(1,2,np+1)
    ax = gca()

    # plot data.
    print plot+' plot ...'
    if plot == 'pcolor':
        m.pcolor(x,y,hgt,shading='flat')
    elif plot == 'imshow':
        im = m.imshow(hgt)
    elif plot == 'contour':
        levels, colls = m.contour(x,y,hgt,15,linewidths=0.5,colors='k')
        levels, colls = m.contourf(x,y,hgt,15,cmap=cm.jet,colors=None)

    # set size of plot to match aspect ratio of map.
    l,b,w,h = ax.get_position()
    b = 0.5 - 0.5*w*m.aspect; h = w*m.aspect
    ax.set_position([l,b,w,h])
    
    # draw map.
    m.drawcoastlines()
	
    # draw parallels
    delat = 30.
    delon = 90.
    circles = arange(10.,90.+delat,delat).tolist()
    m.drawparallels(circles,labels=[0,0,1,1])

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import contour [as 别名]
x_bm,y_bm = ps2.to_basemap_xy(m2, x_grid / cm_per_m, y_grid / cm_per_m)
#Radius factor used in all of the analyses
R = 4.32

heights_cress = oban.analyze_grid(data.heights, x_grid, y_grid, data.x, data.y,
  oban.cressman_weights, R)
contours = N.arange(5300., 6000., 60.0)
parallels = N.arange(25., 60.0, 10.0)
meridians = N.arange(-120., -60.0, 10.0)

##m2.drawstates()
##m2.drawcountries()
##m2.drawcoastlines()
##m2.drawparallels(parallels, labels = [1,1,0,0])
##m2.drawmeridians(meridians, labels = [0,0,1,1])
ob_x, ob_y = ps2.to_basemap_xy(m2, data.x / cm_per_m, data.y / cm_per_m)
#m2.plot(ob_x, ob_y, 'bx')
#m2.plot(x_bm, y_bm, 'g.')
for name in M.cm.cmapnames:
  f = P.figure()
  m2.drawstates()
  m2.drawcountries()
  m2.drawcoastlines()
  m2.drawparallels(parallels, labels = [1,1,0,0])
  m2.drawmeridians(meridians, labels = [0,0,1,1])
  c = m2.contour(x_bm, y_bm, heights_cress, contours, cmap=M.cm.get_cmap(name))
  P.title(name)
  #P.clabel(c)
  P.show()

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import contour [as 别名]
lat_150=[]
for xy in verts:
	lon_150.append(xy[0])
	lat_150.append(xy[1])

close(f)
dd_box=[array(lon_150).min(), array(lon_150).max(), array(lat_150).min(), array(lat_150).max()]
fat=0.05
box=[dd_box[0]-fat, dd_box[1]+fat, dd_box[2]-fat, dd_box[3]+fat]
f=figure()
mapobj= Basemap(projection='merc',lat_0=(box[2]+box[3])/2.0, lon_0=(box[0]+box[1])/2.0,llcrnrlat=box[2], llcrnrlon=box[0], urcrnrlat=box[3] , urcrnrlon=box[1], resolution='l',area_thresh=1., lat_ts=(box[2]+box[3])/2.0)
longr, latgr=meshgrid(lons,lats)
xx, yy = mapobj(longr, latgr)
mapobj.drawmapboundary()
mapobj.readshapefile('/flurry/home/scollis/shapes/cstntcd_r','coast',drawbounds=True,linewidth=0.5,color='k',antialiased=1,ax=None)
mapobj.contour(xx,yy,angs, levels=[150,30], colors=['r'])
mapobj.drawmeridians(array([130.2, 130.4, 130.6, 130.8,131.0,131.2, 131.4]), labels=[1,0,0,1])
mapobj.drawparallels(array([--12.8, -12.6, -12.4, -12.2, -12.0, -11.8, -11.6, -11.4]), labels=[1,0,0,1])
dd_box_xx, dd_box_yy=mapobj(array([dd_box[0], dd_box[1]]),array([dd_box[2], dd_box[3]]))
xy=dd_box_xx[0], dd_box_yy[0]
width, height= dd_box_xx[1]-dd_box_xx[0], dd_box_yy[1]-dd_box_yy[0]
my_patch=matplotlib.patches.Rectangle(xy, width, height, edgecolor='blue', facecolor='white')
ax=gca()
ax.add_patch(my_patch)
radar_xx, radar_yy=mapobj(array([ber_loc[1], gp_loc[1]]), array([ber_loc[0], gp_loc[0]]))
mapobj.plot(radar_xx, radar_yy, 'bo')
ax.text(radar_xx[0]+1000.0, radar_yy[0]-3000.0, 'Berrimah')
ax.text(radar_xx[1]+1000.0, radar_yy[1]-3000.0, 'Gunn Point')
savefig('/flurry/home/scollis/results/area_vis.png')
close(f)