From: <js...@us...> - 2009-03-07 21:52:38
|
Revision: 6962 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=6962&view=rev Author: jswhit Date: 2009-03-07 21:52:21 +0000 (Sat, 07 Mar 2009) Log Message: ----------- add new example showing how to plot H's and L's on a pressure map (uses scipy) Modified Paths: -------------- trunk/toolkits/basemap/MANIFEST.in trunk/toolkits/basemap/examples/README trunk/toolkits/basemap/examples/run_all.py Added Paths: ----------- trunk/toolkits/basemap/examples/plothighsandlows.py Modified: trunk/toolkits/basemap/MANIFEST.in =================================================================== --- trunk/toolkits/basemap/MANIFEST.in 2009-03-07 19:32:39 UTC (rev 6961) +++ trunk/toolkits/basemap/MANIFEST.in 2009-03-07 21:52:21 UTC (rev 6962) @@ -11,6 +11,7 @@ include setup.cfg include setupegg.py include src/* +include examples/plothighsandlows.py include examples/save_background.py include examples/embedding_map_in_wx.py include examples/cubed_sphere.py Modified: trunk/toolkits/basemap/examples/README =================================================================== --- trunk/toolkits/basemap/examples/README 2009-03-07 19:32:39 UTC (rev 6961) +++ trunk/toolkits/basemap/examples/README 2009-03-07 21:52:21 UTC (rev 6962) @@ -125,3 +125,6 @@ figure (without having to redraw coastlines). maskoceans.py shows how to mask 'wet' areas on a plot. + +plothighsandlows.py shows to plot H's and L's at the local extrema of a pressure map +(requires scipy). Added: trunk/toolkits/basemap/examples/plothighsandlows.py =================================================================== --- trunk/toolkits/basemap/examples/plothighsandlows.py (rev 0) +++ trunk/toolkits/basemap/examples/plothighsandlows.py 2009-03-07 21:52:21 UTC (rev 6962) @@ -0,0 +1,89 @@ +""" +plot H's and L's on a sea-level pressure map +""" +import numpy as np +import matplotlib.pyplot as plt +import sys +from mpl_toolkits.basemap import Basemap, NetCDFFile, addcyclic +from scipy.ndimage.filters import minimum_filter, maximum_filter + +def extrema(mat,mode='wrap',window=10): + mn = minimum_filter(mat, size=window, mode=mode) + mx = maximum_filter(mat, size=window, mode=mode) + # (mat == mx) true if pixel is equal to the local max + # The next computation suppresses responses where + # the function is flat. + local_maxima = ((mat == mx) & (mat != mn)) + local_minima = ((mat == mn) & (mat != mx)) + # Get the indices of the maxima, minima + return np.nonzero(local_minima), np.nonzero(local_maxima) + +if len(sys.argv) < 2: + print 'enter date to plot (YYYYMMDDHH) on command line' + raise SystemExit + +# get date from command line. +YYYYMMDDHH = sys.argv[1] + +# open OpenDAP dataset. +data=NetCDFFile("http://nomad1.ncep.noaa.gov:9090/dods/gdas/rotating/gdas"+YYYYMMDDHH) + +# read lats,lons. +lats = data.variables['lat'][:] +lons1 = data.variables['lon'][:] +nlats = len(lats) +nlons = len(lons1) +# read prmsl, convert to hPa (mb). +prmsl = 0.01*data.variables['prmslmsl'][0] +# the window parameter controls the number of highs and lows detected. +# (higher value, fewer highs and lows) +local_min, local_max = extrema(prmsl, mode='wrap', window=25) +# create Basemap instance. +m = Basemap(llcrnrlon=0,llcrnrlat=-65,urcrnrlon=360,urcrnrlat=65,projection='merc') +# add wrap-around point in longitude. +prmsl, lons = addcyclic(prmsl, lons1) +# contour levels +clevs = np.arange(900,1100.,5.) +# find x,y of map projection grid. +lons, lats = np.meshgrid(lons, lats) +x, y = m(lons, lats) +# create figure. +fig=plt.figure(figsize=(12,6)) +cs = m.contour(x,y,prmsl,clevs,colors='k',linewidths=1.) +m.drawcoastlines(linewidth=1.25) +m.fillcontinents(color='0.8') +m.drawparallels(np.arange(-80,81,20)) +m.drawmeridians(np.arange(0,360,60)) +xlows = x[local_min]; xhighs = x[local_max] +ylows = y[local_min]; yhighs = y[local_max] +lowvals = prmsl[local_min]; highvals = prmsl[local_max] +# plot lows as blue L's, with min pressure value underneath. +xyplotted = [] +# don't plot if there is already a L or H within dmin meters. +dmin = 500000 +for x,y,p in zip(xlows, ylows, lowvals): + if x < m.xmax and x > m.xmin and y < m.ymax and y > m.ymin: + dist = [np.sqrt((x-x0)**2+(y-y0)**2) for x0,y0 in xyplotted] + if not dist or min(dist) > dmin: + plt.text(x,y,'L',fontsize=14,fontweight='bold', + horizontalalignment='center', + verticalalignment='center',color='blue') + plt.text(x,y-400000,repr(int(p)),fontsize=9, + horizontalalignment='center', + verticalalignment='top',color='blue') + xyplotted.append((x,y)) +# plot highs as red H's, with max pressure value underneath. +xyplotted = [] +for x,y,p in zip(xhighs, yhighs, highvals): + if x < m.xmax and x > m.xmin and y < m.ymax and y > m.ymin: + dist = [np.sqrt((x-x0)**2+(y-y0)**2) for x0,y0 in xyplotted] + if not dist or min(dist) > dmin: + plt.text(x,y,'H',fontsize=14,fontweight='bold', + horizontalalignment='center', + verticalalignment='center',color='red') + plt.text(x,y-400000,repr(int(p)),fontsize=9, + horizontalalignment='center', + verticalalignment='top',color='red') + xyplotted.append((x,y)) +plt.title('Mean Sea-Level Pressure (with Highs and Lows) %s' % YYYYMMDDHH) +plt.show() Modified: trunk/toolkits/basemap/examples/run_all.py =================================================================== --- trunk/toolkits/basemap/examples/run_all.py 2009-03-07 19:32:39 UTC (rev 6961) +++ trunk/toolkits/basemap/examples/run_all.py 2009-03-07 21:52:21 UTC (rev 6962) @@ -6,7 +6,8 @@ test_files.remove('pnganim.py') test_files.remove('geos_demo_2.py') test_files.remove('plotsst.py') -test_files.remove('embedding_map_in_wx.py') +test_files.remove('embedding_map_in_wx.py') # requires wx +test_files.remove('plothighsandlows.py') # requires scipy print test_files py_path = os.environ.get('PYTHONPATH') if py_path is None: This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |