pmarshwx · June 24, 2013 13:59 · Jun 24, 2013
diff --git a/M1CP.txt b/M1CP.txt
diff --git a/test.py b/test.py
@@ -0,0 +1,103 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.basemap import Basemap
+
+### At least for the time being, the sfcoa_lonlats.npz file
+### can be downloaded here: 
+### http://www.pmarshwx.com/tmp/sfcoa_lonlats.npz
+
+
+def parse(lines, mask=-9999.):
+    '''
+    A simple routine to parse the output of GEMPAK's GDLIST output.
+
+    Parameters
+    ----------
+    lines : string, sequence
+        The output from GDLIST. If a string, will parse on the '\n' character.
+        If a list, will assume the file has already been split.
+    mask : scalar
+        The value of the mask. If no mask, use None or False.
+
+    Returns
+    -------
+    An 2D numpy array
+
+    '''
+    begin = False
+    lonsize = None; latsize = None
+    vals = []
+    if isinstance(lines, str):
+        lines = lines.split('\n')
+    for line in lines:
+        if not lonsize:
+            if 'GRID SIZE' not in line: continue
+            parts = line.split()
+            parts = [p.strip() for p in parts]
+            lonsize = int(parts[-2])
+            latsize = int(parts[-1])
+        if 'ROW'+str(latsize) in line or (begin and 'ROW' in line):
+            begin = True
+            val_tmp = line.split('ROW')[1].split()[1:]
+        elif begin:
+            val_tmp = line.split()
+        else:
+            continue
+        for val in val_tmp:
+            vals.append(float(val))
+    vals = np.array(vals)
+    if mask:
+        vals = np.ma.masked_where(vals == mask, vals)
+    vals = vals.reshape(latsize, -1)
+    return vals
+
+
+f = np.load('sfcoa_lonlats.npz')
+lons = f['lons']
+lats = f['lats']
+f.close()
+
+lines = open('M1CP.txt').read().split('\n')
+mlcape = parse(lines)
+mlcape = mlcape.filled(0)
+
+
+projection = 'lcc'
+resolution = 'l'
+llcrnrlon = -121.71
+llcrnrlat = 22.77
+urcrnrlon = -65.6
+urcrnrlat = 50.15
+lon_0 = -95
+lat_0 = 35.
+lat_1 = 25.
+lat_2 = 25.
+
+m = Basemap(projection=projection, resolution=resolution, llcrnrlon=llcrnrlon,
+            llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat,
+            lon_0=lon_0, lat_0=lat_0, lat_1=lat_1, lat_2=lat_2)
+
+m.drawcoastlines()
+m.drawcountries()
+m.drawstates()
+m.contourf(lons, lats, mlcape, levels=range(0, 5001, 500), latlon=True)
+plt.show()
+
+m.drawcoastlines()
+m.drawcountries()
+m.drawstates()
+x, y = m(lons, lats)
+m.contourf(x, y, mlcape, levels=range(0, 5001, 500))
+plt.show()
+
+m.drawcoastlines()
+m.drawcountries()
+m.drawstates()
+m.pcolormesh(lons, lats, mlcape, latlon=True)
+plt.show()
+
+m.drawcoastlines()
+m.drawcountries()
+m.drawstates()
+plt.contourf(x, y, mlcape, levels=range(0, 5001, 500))
+plt.show()