fcstmaps.py

"""Make a multi-panel plot from numerical weather forecast in NOAA OPeNDAP."""
from __future__ importprint_function

importnetCDF4
importnumpyasnp
importmatplotlib.pyplotasplt
frommpl_toolkits.basemapimportBasemap
frommpl_toolkits.basemapimportaddcyclic


defmain(date, verbose=True):
"""Main function."""

# Open dataset from OPeNDAP URL.
url="http://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs%Y%m%d/gfs_0p25_00z"
try:
data=netCDF4.Dataset(date.strftime(url), "r")
ifverbose:
print("Data variables:")
print(sorted(data.variables))
exceptOSErroraserr:
err.args= (err.args[0], "date not found in OPeNDAP server")
raise

# Read lats, lons, and times.
latitudes=data.variables["lat"]
longitudes=data.variables["lon"]
fcsttimes=data.variables["time"]
times=fcsttimes[0:6] # First 6 forecast times.
ntimes=len(times)

# Convert times for datetime instances.
fdates=netCDF4.num2date(
times, units=fcsttimes.units, calendar="standard")

# Make a list of YYYYMMDDHH strings.
verifdates= [fdate.strftime("%Y%m%d%H") forfdateinfdates]
ifverbose:
print("Forecast datetime strings:")
print(verifdates)

# Convert times to forecast hours.
fcsthrs= []
forfdateinfdates:
fdiff=fdate-fdates[0]
fcsthrs.append(fdiff.days*24.+fdiff.seconds/3600.)
ifverbose:
print("Forecast datetime hours:")
print(fcsthrs)

# Unpack 2-meter temp forecast data.
lats=latitudes[:]
nlats=len(lats)
lons1=longitudes[:]
nlons=len(lons1)
t2mvar=data.variables["tmp2m"]

# Create Basemap instance for orthographic projection.
bmap=Basemap(lon_0=-90, lat_0=60, projection="ortho")

# Add wrap-around point in longitude.
t2m=np.zeros((ntimes, nlats, nlons+1), np.float32)
forntinrange(ntimes):
t2m[nt, :, :], lons=addcyclic(t2mvar[nt, :, :], lons1)

# Convert to Celsius.
t2m=t2m-273.15

# Define contour levels.
clevs=np.arange(-30, 30.1, 2.0)
lons, lats=np.meshgrid(lons, lats)
x, y=bmap(lons, lats)

# Create figure.
fig=plt.figure(figsize=(6, 8))

# Make subplots.
fornt, fcsthrinenumerate(fcsthrs):
fig.add_subplot(321+nt)
cs=bmap.contourf(x, y, t2m[nt, :, :], clevs,
cmap=plt.cm.jet, extend="both")
bmap.drawcoastlines(linewidth=0.5)
bmap.drawcountries()
bmap.drawparallels(np.arange(-80, 81, 20))
bmap.drawmeridians(np.arange(0, 360, 20))
# Set panel title.
plt.title(
"%d-h forecast valid "%fcsthr+verifdates[nt], fontsize=9)

# Set figure title.
plt.figtext(
0.5, 0.95,
"2-m temp (\N{DEGREE SIGN}C) forecasts from %s"%verifdates[0],
horizontalalignment="center", fontsize=14)

# Draw a single colorbar.
cax=plt.axes([0.1, 0.05, 0.8, 0.025])
plt.colorbar(cs, cax=cax, orientation="horizontal")
plt.show()


if__name__=="__main__":

importsys
importdatetimeasdt

# Parse input date (default: today).
iflen(sys.argv) >1:
dateobj=dt.datetime.strptime(sys.argv[1], "%Y%m%d")
else:
dateobj=dt.datetime.today()
main(dateobj, verbose=True)