|
1 | | -from __future__ import (absolute_import, division, print_function) |
| 1 | +"""Plot H's and L's on a sea-level pressure map.""" |
| 2 | +from __future__ import print_function |
2 | 3 |
|
3 | | -""" |
4 | | -plot H's and L's on a sea-level pressure map |
5 | | -(uses scipy.ndimage.filters and netcdf4-python) |
6 | | -""" |
| 4 | +import datetime as dt |
| 5 | +import netCDF4 |
7 | 6 | import numpy as np |
8 | 7 | import matplotlib.pyplot as plt |
9 | | -from datetime import datetime |
10 | | -from mpl_toolkits.basemap import Basemap, addcyclic |
11 | | -from scipy.ndimage.filters import minimum_filter, maximum_filter |
12 | | -from netCDF4 import Dataset |
13 | | - |
14 | | -def extrema(mat,mode='wrap',window=10): |
15 | | - """find the indices of local extrema (min and max) |
16 | | - in the input array.""" |
17 | | - mn = minimum_filter(mat, size=window, mode=mode) |
18 | | - mx = maximum_filter(mat, size=window, mode=mode) |
19 | | - # (mat == mx) true if pixel is equal to the local max |
20 | | - # (mat == mn) true if pixel is equal to the local in |
21 | | - # Return the indices of the maxima, minima |
22 | | - return np.nonzero(mat == mn), np.nonzero(mat == mx) |
23 | | - |
24 | | -# plot 00 UTC today. |
25 | | -date = datetime.now().strftime('%Y%m%d')+'00' |
26 | | - |
27 | | -# open OpenDAP dataset. |
28 | | -#data=Dataset("http://nomads.ncep.noaa.gov:9090/dods/gfs/gfs/%s/gfs_%sz_anl" %\ |
29 | | -# (date[0:8],date[8:10])) |
30 | | -data=Dataset("http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd%s/gfs_hd_%sz"%\ |
31 | | - (date[0:8],date[8:10])) |
32 | | - |
33 | | - |
34 | | - |
35 | | -# read lats,lons. |
36 | | -lats = data.variables['lat'][:] |
37 | | -lons1 = data.variables['lon'][:] |
38 | | -nlats = len(lats) |
39 | | -nlons = len(lons1) |
40 | | -# read prmsl, convert to hPa (mb). |
41 | | -prmsl = 0.01*data.variables['prmslmsl'][0] |
42 | | -# the window parameter controls the number of highs and lows detected. |
43 | | -# (higher value, fewer highs and lows) |
44 | | -local_min, local_max = extrema(prmsl, mode='wrap', window=50) |
45 | | -# create Basemap instance. |
46 | | -m =\ |
47 | | -Basemap(llcrnrlon=0,llcrnrlat=-80,urcrnrlon=360,urcrnrlat=80,projection='mill') |
48 | | -# add wrap-around point in longitude. |
49 | | -prmsl, lons = addcyclic(prmsl, lons1) |
50 | | -# contour levels |
51 | | -clevs = np.arange(900,1100.,5.) |
52 | | -# find x,y of map projection grid. |
53 | | -lons, lats = np.meshgrid(lons, lats) |
54 | | -x, y = m(lons, lats) |
55 | | -# create figure. |
56 | | -fig=plt.figure(figsize=(8,4.5)) |
57 | | -ax = fig.add_axes([0.05,0.05,0.9,0.85]) |
58 | | -cs = m.contour(x,y,prmsl,clevs,colors='k',linewidths=1.) |
59 | | -m.drawcoastlines(linewidth=1.25) |
60 | | -m.fillcontinents(color='0.8') |
61 | | -m.drawparallels(np.arange(-80,81,20),labels=[1,1,0,0]) |
62 | | -m.drawmeridians(np.arange(0,360,60),labels=[0,0,0,1]) |
63 | | -xlows = x[local_min]; xhighs = x[local_max] |
64 | | -ylows = y[local_min]; yhighs = y[local_max] |
65 | | -lowvals = prmsl[local_min]; highvals = prmsl[local_max] |
66 | | -# plot lows as blue L's, with min pressure value underneath. |
67 | | -xyplotted = [] |
68 | | -# don't plot if there is already a L or H within dmin meters. |
69 | | -yoffset = 0.022*(m.ymax-m.ymin) |
70 | | -dmin = yoffset |
71 | | -for x,y,p in zip(xlows, ylows, lowvals): |
72 | | - if x < m.xmax and x > m.xmin and y < m.ymax and y > m.ymin: |
73 | | - dist = [np.sqrt((x-x0)**2+(y-y0)**2) for x0,y0 in xyplotted] |
74 | | - if not dist or min(dist) > dmin: |
75 | | - plt.text(x,y,'L',fontsize=14,fontweight='bold', |
76 | | - ha='center',va='center',color='b') |
77 | | - plt.text(x,y-yoffset,repr(int(p)),fontsize=9, |
78 | | - ha='center',va='top',color='b', |
79 | | - bbox = dict(boxstyle="square",ec='None',fc=(1,1,1,0.5))) |
80 | | - xyplotted.append((x,y)) |
81 | | -# plot highs as red H's, with max pressure value underneath. |
82 | | -xyplotted = [] |
83 | | -for x,y,p in zip(xhighs, yhighs, highvals): |
84 | | - if x < m.xmax and x > m.xmin and y < m.ymax and y > m.ymin: |
85 | | - dist = [np.sqrt((x-x0)**2+(y-y0)**2) for x0,y0 in xyplotted] |
86 | | - if not dist or min(dist) > dmin: |
87 | | - plt.text(x,y,'H',fontsize=14,fontweight='bold', |
88 | | - ha='center',va='center',color='r') |
89 | | - plt.text(x,y-yoffset,repr(int(p)),fontsize=9, |
90 | | - ha='center',va='top',color='r', |
91 | | - bbox = dict(boxstyle="square",ec='None',fc=(1,1,1,0.5))) |
92 | | - xyplotted.append((x,y)) |
93 | | -plt.title('Mean Sea-Level Pressure (with Highs and Lows) %s' % date) |
94 | | -plt.show() |
| 8 | +from mpl_toolkits.basemap import Basemap |
| 9 | +from mpl_toolkits.basemap import addcyclic |
| 10 | +from scipy.ndimage import minimum_filter |
| 11 | +from scipy.ndimage import maximum_filter |
| 12 | + |
| 13 | + |
| 14 | +def extrema(mat, mode="wrap", window=10): |
| 15 | + """Find the indices of local extrema (min and max) in the input array.""" |
| 16 | + |
| 17 | + minimum = minimum_filter(mat, size=window, mode=mode) |
| 18 | + maximum = maximum_filter(mat, size=window, mode=mode) |
| 19 | + |
| 20 | + # Return the indices of the maxima, minima. |
| 21 | + # (mat == maximum) true if pixel is equal to the local max. |
| 22 | + # (mat == minimum) true if pixel is equal to the local in. |
| 23 | + return np.nonzero(mat == minimum), np.nonzero(mat == maximum) |
| 24 | + |
| 25 | + |
| 26 | +def main(): |
| 27 | + """Main function.""" |
| 28 | + |
| 29 | + # Plot 00 UTC today. |
| 30 | + url = "http://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs%Y%m%d/gfs_0p25_00z" |
| 31 | + date = dt.datetime.now() |
| 32 | + |
| 33 | + # Open OPeNDAP dataset. |
| 34 | + data = netCDF4.Dataset(date.strftime(url)) |
| 35 | + |
| 36 | + # Read lats and lons. |
| 37 | + lats = data.variables["lat"][:] |
| 38 | + lons1 = data.variables["lon"][:] |
| 39 | + |
| 40 | + # Read prmsl and convert to hPa (mbar). |
| 41 | + prmsl = 0.01 * data.variables["prmslmsl"][0] |
| 42 | + |
| 43 | + # The window parameter controls the number of highs and lows detected |
| 44 | + # (higher value, fewer highs and lows). |
| 45 | + local_min, local_max = extrema(prmsl, mode="wrap", window=50) |
| 46 | + |
| 47 | + # Create Basemap instance. |
| 48 | + bmap = Basemap(projection="mill", |
| 49 | + llcrnrlon=0, llcrnrlat=-80, |
| 50 | + urcrnrlon=360, urcrnrlat=80) |
| 51 | + |
| 52 | + # Add wrap-around point in longitude. |
| 53 | + prmsl, lons = addcyclic(prmsl, lons1) |
| 54 | + |
| 55 | + # Define contour levels. |
| 56 | + clevs = np.arange(900, 1100., 5.) |
| 57 | + |
| 58 | + # Find x, y of map projection grid. |
| 59 | + lons, lats = np.meshgrid(lons, lats) |
| 60 | + x, y = bmap(lons, lats) |
| 61 | + |
| 62 | + # Create figure. |
| 63 | + fig = plt.figure(figsize=(8, 4.5)) |
| 64 | + fig.add_axes([0.05, 0.05, 0.9, 0.85]) |
| 65 | + bmap.contour(x, y, prmsl, clevs, colors="k", linewidths=1.0) |
| 66 | + bmap.drawcoastlines(linewidth=1.25) |
| 67 | + bmap.fillcontinents(color="0.8") |
| 68 | + bmap.drawparallels(np.arange(-80, 81, 20), labels=[1, 1, 0, 0]) |
| 69 | + bmap.drawmeridians(np.arange(0, 360, 60), labels=[0, 0, 0, 1]) |
| 70 | + xlows, xhighs = x[local_min], x[local_max] |
| 71 | + ylows, yhighs = y[local_min], y[local_max] |
| 72 | + lowvals, highvals = prmsl[local_min], prmsl[local_max] |
| 73 | + |
| 74 | + # Plot lows as blue L's, with min pressure value underneath. |
| 75 | + # Do not plot if there is already a L or H within dmin meters. |
| 76 | + xyplotted = [] |
| 77 | + yoffset = 0.022 * (bmap.ymax - bmap.ymin) |
| 78 | + dmin = yoffset |
| 79 | + for x, y, p in zip(xlows, ylows, lowvals): |
| 80 | + if bmap.xmin < x < bmap.xmax and bmap.ymin < y < bmap.ymax: |
| 81 | + dist = [np.sqrt((x - x0)**2 + (y - y0)**2) for x0, y0 in xyplotted] |
| 82 | + if not dist or min(dist) > dmin: |
| 83 | + bbox = dict(boxstyle="square", ec="None", fc=(1, 1, 1, 0.5)) |
| 84 | + plt.text(x, y, "L", fontsize=14, fontweight="bold", |
| 85 | + ha="center", va="center", color="b") |
| 86 | + plt.text(x, y - yoffset, repr(int(p)), fontsize=9, |
| 87 | + ha="center", va="top", color="b", bbox=bbox) |
| 88 | + xyplotted.append((x, y)) |
| 89 | + # Plot highs as red H's, with max pressure value underneath. |
| 90 | + xyplotted = [] |
| 91 | + for x, y, p in zip(xhighs, yhighs, highvals): |
| 92 | + if bmap.xmin < x < bmap.xmax and bmap.ymin < y < bmap.ymax: |
| 93 | + dist = [np.sqrt((x - x0)**2 + (y - y0)**2) for x0, y0 in xyplotted] |
| 94 | + if not dist or min(dist) > dmin: |
| 95 | + bbox = dict(boxstyle="square", ec="None", fc=(1, 1, 1, 0.5)) |
| 96 | + plt.text(x, y, "H", fontsize=14, fontweight="bold", |
| 97 | + ha="center", va="center", color="r") |
| 98 | + plt.text(x, y - yoffset, repr(int(p)), fontsize=9, |
| 99 | + ha="center", va="top", color="r", bbox=bbox) |
| 100 | + xyplotted.append((x, y)) |
| 101 | + |
| 102 | + # Set plot title and show. |
| 103 | + plt.title("Mean Sea-Level Pressure (with Highs and Lows) %s" % date) |
| 104 | + plt.show() |
| 105 | + |
| 106 | + |
| 107 | +if __name__ == "__main__": |
| 108 | + main() |
0 commit comments