Menu
▾
▴
SF.net SVN: matplotlib:[6962] trunk/toolkits/basemap
|
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.
|
