From f988c1a56d72fe5b21f147fea37c9b8c8dd83e5a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?V=C3=ADctor=20Molina=20Garc=C3=ADa?= Date: Wed, 19 Oct 2022 19:16:46 +0200 Subject: [PATCH 1/7] Fix URLs in examples using broken addcyclic --- examples/fcstmaps.py | 12 +++++------- examples/fcstmaps_axesgrid.py | 13 +++++-------- examples/plothighsandlows.py | 12 ++++-------- 3 files changed, 14 insertions(+), 23 deletions(-) diff --git a/examples/fcstmaps.py b/examples/fcstmaps.py index dcfe6802d..e907c299d 100644 --- a/examples/fcstmaps.py +++ b/examples/fcstmaps.py @@ -3,27 +3,25 @@ from __future__ import unicode_literals # this example reads today's numerical weather forecasts # from the NOAA OpenDAP servers and makes a multi-panel plot. +import datetime as dt import numpy as np import matplotlib.pyplot as plt import sys import numpy.ma as ma -import datetime from mpl_toolkits.basemap import Basemap, addcyclic from netCDF4 import Dataset as NetCDFFile, num2date # today's date is default. if len(sys.argv) > 1: - YYYYMMDD = sys.argv[1] + date = dt.datetime.strptime(sys.argv[1], "%Y%m%d") else: - YYYYMMDD = datetime.datetime.today().strftime('%Y%m%d') + date = dt.datetime.today() # set OpenDAP server URL. try: - URLbase="http://nomads.ncep.noaa.gov:9090/dods/gfs/gfs" - URL=URLbase+YYYYMMDD+'/gfs_00z' - print(URL) - data = NetCDFFile(URL) + urlbase = "http://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs%Y%m%d/gfs_0p25_00z" + data = NetCDFFile(date.strftime(urlbase)) except: msg = """ opendap server not providing the requested data. diff --git a/examples/fcstmaps_axesgrid.py b/examples/fcstmaps_axesgrid.py index d4d3d44df..9c7b673af 100644 --- a/examples/fcstmaps_axesgrid.py +++ b/examples/fcstmaps_axesgrid.py @@ -4,11 +4,11 @@ # this example reads today's numerical weather forecasts # from the NOAA OpenDAP servers and makes a multi-panel plot. # This version demonstrates the use of the AxesGrid toolkit. +import datetime as dt import numpy as np import matplotlib.pyplot as plt import sys import numpy.ma as ma -import datetime from mpl_toolkits.basemap import Basemap, addcyclic from mpl_toolkits.axes_grid1 import AxesGrid from netCDF4 import Dataset as NetCDFFile, num2date @@ -16,23 +16,20 @@ # today's date is default. if len(sys.argv) > 1: - YYYYMMDD = sys.argv[1] + date = dt.datetime.strptime(sys.argv[1], "%Y%m%d") else: - YYYYMMDD = datetime.datetime.today().strftime('%Y%m%d') + date = dt.datetime.today() # set OpenDAP server URL. try: - URLbase="http://nomads.ncep.noaa.gov:9090/dods/gfs/gfs" - URL=URLbase+YYYYMMDD+'/gfs_00z' - print(URL) - data = NetCDFFile(URL) + urlbase = "http://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs%Y%m%d/gfs_0p25_00z" + data = NetCDFFile(date.strftime(urlbase)) except: msg = """ opendap server not providing the requested data. Try another date by providing YYYYMMDD on command line.""" raise IOError(msg) - # read lats,lons,times. print(data.variables.keys()) diff --git a/examples/plothighsandlows.py b/examples/plothighsandlows.py index 9159e13b0..a0ea1c34e 100644 --- a/examples/plothighsandlows.py +++ b/examples/plothighsandlows.py @@ -4,9 +4,9 @@ plot H's and L's on a sea-level pressure map (uses scipy.ndimage.filters and netcdf4-python) """ +import datetime as dt import numpy as np import matplotlib.pyplot as plt -from datetime import datetime from mpl_toolkits.basemap import Basemap, addcyclic from scipy.ndimage.filters import minimum_filter, maximum_filter from netCDF4 import Dataset @@ -22,15 +22,11 @@ def extrema(mat,mode='wrap',window=10): return np.nonzero(mat == mn), np.nonzero(mat == mx) # plot 00 UTC today. -date = datetime.now().strftime('%Y%m%d')+'00' +urlbase = "http://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs%Y%m%d/gfs_0p25_00z" +date = dt.datetime.now() # open OpenDAP dataset. -#data=Dataset("http://nomads.ncep.noaa.gov:9090/dods/gfs/gfs/%s/gfs_%sz_anl" %\ -# (date[0:8],date[8:10])) -data=Dataset("http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd%s/gfs_hd_%sz"%\ - (date[0:8],date[8:10])) - - +data = Dataset(date.strftime(urlbase)) # read lats,lons. lats = data.variables['lat'][:] From 3313c180ac7bf1cffe7fa14cc4e23cce299d3365 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?V=C3=ADctor=20Molina=20Garc=C3=ADa?= Date: Wed, 19 Oct 2022 19:26:35 +0200 Subject: [PATCH 2/7] Fix addcyclic array indexing with list instead of tuple --- packages/basemap/src/mpl_toolkits/basemap/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/basemap/src/mpl_toolkits/basemap/__init__.py b/packages/basemap/src/mpl_toolkits/basemap/__init__.py index bdacd1c96..653813f34 100644 --- a/packages/basemap/src/mpl_toolkits/basemap/__init__.py +++ b/packages/basemap/src/mpl_toolkits/basemap/__init__.py @@ -5127,7 +5127,7 @@ def _addcyclic(a): except IndexError: raise ValueError('The specified axis does not correspond to an ' 'array dimension.') - return npsel.concatenate((a,a[slicer]),axis=axis) + return npsel.concatenate((a,a[tuple(slicer)]),axis=axis) def _addcyclic_lon(a): """addcyclic function for a single longitude array""" # select the right numpy functions From 4824629169fcdb4eb74791649c38e09dcb1bbcfb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?V=C3=ADctor=20Molina=20Garc=C3=ADa?= Date: Wed, 19 Oct 2022 19:35:29 +0200 Subject: [PATCH 3/7] Add bugfix for issue 555 to CHANGELOG --- CHANGELOG.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8a6e9b9e7..1f4415f45 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,12 @@ https://keepachangelog.com/en/1.0.0/ https://semver.org/spec/v2.0.0.html +## [1.3.5] + +### Fixed +- Fix broken array slicing inside `addcyclic` (PR [#559], solves issue + [#555], thanks to @fragkoul). + ## [1.3.4] - 2022-08-10 ### Changed @@ -924,6 +930,10 @@ https://semver.org/spec/v2.0.0.html - Fix glitches in drawing of parallels and meridians. +[#559]: +https://github.com/matplotlib/basemap/pull/559 +[#555]: +https://github.com/matplotlib/basemap/issues/555 [#548]: https://github.com/matplotlib/basemap/pull/548 [#547]: From 62e94a2a1b77597976ec8cefb6ad7b16144f3a36 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?V=C3=ADctor=20Molina=20Garc=C3=ADa?= Date: Wed, 19 Oct 2022 20:09:47 +0200 Subject: [PATCH 4/7] Cleanse examples/fcstmaps.py --- examples/fcstmaps.py | 197 ++++++++++++++++++++++++------------------- 1 file changed, 110 insertions(+), 87 deletions(-) diff --git a/examples/fcstmaps.py b/examples/fcstmaps.py index e907c299d..c942c339b 100644 --- a/examples/fcstmaps.py +++ b/examples/fcstmaps.py @@ -1,90 +1,113 @@ -from __future__ import (absolute_import, division, print_function) +"""Make a multi-panel plot from numerical weather forecast in NOAA OPeNDAP.""" +from __future__ import print_function -from __future__ import unicode_literals -# this example reads today's numerical weather forecasts -# from the NOAA OpenDAP servers and makes a multi-panel plot. -import datetime as dt +import netCDF4 import numpy as np import matplotlib.pyplot as plt -import sys -import numpy.ma as ma -from mpl_toolkits.basemap import Basemap, addcyclic -from netCDF4 import Dataset as NetCDFFile, num2date - - -# today's date is default. -if len(sys.argv) > 1: - date = dt.datetime.strptime(sys.argv[1], "%Y%m%d") -else: - date = dt.datetime.today() - -# set OpenDAP server URL. -try: - urlbase = "http://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs%Y%m%d/gfs_0p25_00z" - data = NetCDFFile(date.strftime(urlbase)) -except: - msg = """ -opendap server not providing the requested data. -Try another date by providing YYYYMMDD on command line.""" - raise IOError(msg) - - -# read lats,lons,times. - -print(data.variables.keys()) -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 = num2date(times,units=fcsttimes.units,calendar='standard') -# make a list of YYYYMMDDHH strings. -verifdates = [fdate.strftime('%Y%m%d%H') for fdate in fdates] -# convert times to forecast hours. -fcsthrs = [] -for fdate in fdates: - fdiff = fdate-fdates[0] - fcsthrs.append(fdiff.days*24. + fdiff.seconds/3600.) -print(fcsthrs) -print(verifdates) -lats = latitudes[:] -nlats = len(lats) -lons1 = longitudes[:] -nlons = len(lons1) - -# unpack 2-meter temp forecast data. - -t2mvar = data.variables['tmp2m'] -t2m = np.zeros((ntimes,nlats,nlons+1),np.float32) -# create Basemap instance for Orthographic projection. -m = Basemap(lon_0=-90,lat_0=60,projection='ortho') -# add wrap-around point in longitude. -for nt in range(ntimes): - t2m[nt,:,:], lons = addcyclic(t2mvar[nt,:,:], lons1) -# convert to celsius. -t2m = t2m-273.15 -# contour levels -clevs = np.arange(-30,30.1,2.) -lons, lats = np.meshgrid(lons, lats) -x, y = m(lons, lats) -# create figure. -fig=plt.figure(figsize=(6,8)) -# make subplots. -for nt,fcsthr in enumerate(fcsthrs): - ax = fig.add_subplot(321+nt) - cs = m.contourf(x,y,t2m[nt,:,:],clevs,cmap=plt.cm.jet,extend='both') - m.drawcoastlines(linewidth=0.5) - m.drawcountries() - m.drawparallels(np.arange(-80,81,20)) - m.drawmeridians(np.arange(0,360,20)) - # panel title - plt.title('%d-h forecast valid '%fcsthr+verifdates[nt],fontsize=9) -# figure title -plt.figtext(0.5,0.95, - "2-m temp (\N{DEGREE SIGN}C) forecasts from %s"%verifdates[0], - horizontalalignment='center',fontsize=14) -# a single colorbar. -cax = plt.axes([0.1, 0.05, 0.8, 0.025]) -plt.colorbar(cax=cax, orientation='horizontal') -plt.show() +from mpl_toolkits.basemap import Basemap +from mpl_toolkits.basemap import addcyclic + + +def main(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") + if verbose: + print("Data variables:") + print(sorted(data.variables)) + except OSError as err: + 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") for fdate in fdates] + if verbose: + print("Forecast datetime strings:") + print(verifdates) + + # Convert times to forecast hours. + fcsthrs = [] + for fdate in fdates: + fdiff = fdate - fdates[0] + fcsthrs.append(fdiff.days * 24. + fdiff.seconds / 3600.) + if verbose: + 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) + for nt in range(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. + for nt, fcsthr in enumerate(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__": + + import sys + import datetime as dt + + # Parse input date (default: today). + if len(sys.argv) > 1: + dateobj = dt.datetime.strptime(sys.argv[1], "%Y%m%d") + else: + dateobj = dt.datetime.today() + + sys.exit(main(dateobj, verbose=True)) From 34d83cede083ddf86fd763f712e3e48523ebc613 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?V=C3=ADctor=20Molina=20Garc=C3=ADa?= Date: Wed, 19 Oct 2022 20:32:29 +0200 Subject: [PATCH 5/7] Cleanse examples/fcstmaps_axesgrid.py --- examples/fcstmaps_axesgrid.py | 221 +++++++++++++++++++--------------- 1 file changed, 123 insertions(+), 98 deletions(-) diff --git a/examples/fcstmaps_axesgrid.py b/examples/fcstmaps_axesgrid.py index 9c7b673af..2748a1bac 100644 --- a/examples/fcstmaps_axesgrid.py +++ b/examples/fcstmaps_axesgrid.py @@ -1,102 +1,127 @@ -from __future__ import (absolute_import, division, print_function) +"""Make a multi-panel plot from numerical weather forecast in NOAA OPeNDAP. -from __future__ import unicode_literals -# this example reads today's numerical weather forecasts -# from the NOAA OpenDAP servers and makes a multi-panel plot. -# This version demonstrates the use of the AxesGrid toolkit. -import datetime as dt +This version demonstrates the use of the AxesGrid toolkit. +""" +from __future__ import print_function + +import netCDF4 import numpy as np import matplotlib.pyplot as plt -import sys -import numpy.ma as ma -from mpl_toolkits.basemap import Basemap, addcyclic +from mpl_toolkits.basemap import Basemap +from mpl_toolkits.basemap import addcyclic from mpl_toolkits.axes_grid1 import AxesGrid -from netCDF4 import Dataset as NetCDFFile, num2date - - -# today's date is default. -if len(sys.argv) > 1: - date = dt.datetime.strptime(sys.argv[1], "%Y%m%d") -else: - date = dt.datetime.today() - -# set OpenDAP server URL. -try: - urlbase = "http://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs%Y%m%d/gfs_0p25_00z" - data = NetCDFFile(date.strftime(urlbase)) -except: - msg = """ -opendap server not providing the requested data. -Try another date by providing YYYYMMDD on command line.""" - raise IOError(msg) - -# read lats,lons,times. - -print(data.variables.keys()) -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 = num2date(times,units=fcsttimes.units,calendar='standard') -# make a list of YYYYMMDDHH strings. -verifdates = [fdate.strftime('%Y%m%d%H') for fdate in fdates] -# convert times to forecast hours. -fcsthrs = [] -for fdate in fdates: - fdiff = fdate-fdates[0] - fcsthrs.append(fdiff.days*24. + fdiff.seconds/3600.) -print(fcsthrs) -print(verifdates) -lats = latitudes[:] -nlats = len(lats) -lons1 = longitudes[:] -nlons = len(lons1) - -# unpack 2-meter temp forecast data. - -t2mvar = data.variables['tmp2m'] - -# create figure, set up AxesGrid. -fig=plt.figure(figsize=(6,8)) -grid = AxesGrid(fig, [0.05,0.01,0.9,0.9], - nrows_ncols=(3, 2), - axes_pad=0.25, - cbar_mode='single', - cbar_pad=0.3, - cbar_size=0.1, - cbar_location='top', - share_all=True, - ) - -# create Basemap instance for Orthographic projection. -m = Basemap(lon_0=-90,lat_0=60,projection='ortho') -# add wrap-around point in longitude. -t2m = np.zeros((ntimes,nlats,nlons+1),np.float32) -for nt in range(ntimes): - t2m[nt,:,:], lons = addcyclic(t2mvar[nt,:,:], lons1) -# convert to celsius. -t2m = t2m-273.15 -# contour levels -clevs = np.arange(-30,30.1,2.) -lons, lats = np.meshgrid(lons, lats) -x, y = m(lons, lats) -# make subplots. -for nt,fcsthr in enumerate(fcsthrs): - ax = grid[nt] - m.ax = ax - cs = m.contourf(x,y,t2m[nt,:,:],clevs,cmap=plt.cm.jet,extend='both') - m.drawcoastlines(linewidth=0.5) - m.drawcountries() - m.drawparallels(np.arange(-80,81,20)) - m.drawmeridians(np.arange(0,360,20)) - # panel title - ax.set_title('%d-h forecast valid '%fcsthr+verifdates[nt],fontsize=9) -# figure title -plt.figtext(0.5,0.95, - "2-m temp (\N{DEGREE SIGN}C) forecasts from %s"%verifdates[0], - horizontalalignment='center',fontsize=14) -# a single colorbar. -cbar = fig.colorbar(cs, cax=grid.cbar_axes[0], orientation='horizontal') -plt.show() + + +def main(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") + if verbose: + print("Data variables:") + print(sorted(data.variables)) + except OSError as err: + 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") for fdate in fdates] + if verbose: + print("Forecast datetime strings:") + print(verifdates) + + # Convert times to forecast hours. + fcsthrs = [] + for fdate in fdates: + fdiff = fdate - fdates[0] + fcsthrs.append(fdiff.days * 24. + fdiff.seconds / 3600.) + if verbose: + 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) + for nt in range(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 and AxesGrid instance. + fig = plt.figure(figsize=(6, 8)) + grid = AxesGrid( + fig, + [0.05, 0.01, 0.9, 0.9], + nrows_ncols=(3, 2), + axes_pad=0.5, + cbar_mode="single", + cbar_pad=0.75, + cbar_size=0.1, + cbar_location="top", + share_all=True) + + # Make subplots. + for nt, fcsthr in enumerate(fcsthrs): + bmap.ax = grid[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. + bmap.ax.set_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 = grid.cbar_axes[0] + fig.colorbar(cs, cax=cax, orientation="horizontal") + plt.show() + + +if __name__ == "__main__": + + import sys + import datetime as dt + + # Parse input date (default: today). + if len(sys.argv) > 1: + dateobj = dt.datetime.strptime(sys.argv[1], "%Y%m%d") + else: + dateobj = dt.datetime.today() + + sys.exit(main(dateobj, verbose=True)) From 2cc0d9f388cc9a8322c3c42847e93a567a5eb119 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?V=C3=ADctor=20Molina=20Garc=C3=ADa?= Date: Wed, 19 Oct 2022 20:50:21 +0200 Subject: [PATCH 6/7] Cleanse examples/plothighsandlows.py --- examples/plothighsandlows.py | 192 +++++++++++++++++++---------------- 1 file changed, 106 insertions(+), 86 deletions(-) diff --git a/examples/plothighsandlows.py b/examples/plothighsandlows.py index a0ea1c34e..edacc6875 100644 --- a/examples/plothighsandlows.py +++ b/examples/plothighsandlows.py @@ -1,90 +1,110 @@ -from __future__ import (absolute_import, division, print_function) +"""Plot H's and L's on a sea-level pressure map.""" +from __future__ import print_function -""" -plot H's and L's on a sea-level pressure map -(uses scipy.ndimage.filters and netcdf4-python) -""" import datetime as dt +import netCDF4 import numpy as np import matplotlib.pyplot as plt -from mpl_toolkits.basemap import Basemap, addcyclic -from scipy.ndimage.filters import minimum_filter, maximum_filter -from netCDF4 import Dataset - -def extrema(mat,mode='wrap',window=10): - """find the indices of local extrema (min and max) - in the input array.""" - 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 - # (mat == mn) true if pixel is equal to the local in - # Return the indices of the maxima, minima - return np.nonzero(mat == mn), np.nonzero(mat == mx) - -# plot 00 UTC today. -urlbase = "http://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs%Y%m%d/gfs_0p25_00z" -date = dt.datetime.now() - -# open OpenDAP dataset. -data = Dataset(date.strftime(urlbase)) - -# 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=50) -# create Basemap instance. -m =\ -Basemap(llcrnrlon=0,llcrnrlat=-80,urcrnrlon=360,urcrnrlat=80,projection='mill') -# 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=(8,4.5)) -ax = fig.add_axes([0.05,0.05,0.9,0.85]) -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),labels=[1,1,0,0]) -m.drawmeridians(np.arange(0,360,60),labels=[0,0,0,1]) -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. -yoffset = 0.022*(m.ymax-m.ymin) -dmin = yoffset -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', - ha='center',va='center',color='b') - plt.text(x,y-yoffset,repr(int(p)),fontsize=9, - ha='center',va='top',color='b', - bbox = dict(boxstyle="square",ec='None',fc=(1,1,1,0.5))) - 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', - ha='center',va='center',color='r') - plt.text(x,y-yoffset,repr(int(p)),fontsize=9, - ha='center',va='top',color='r', - bbox = dict(boxstyle="square",ec='None',fc=(1,1,1,0.5))) - xyplotted.append((x,y)) -plt.title('Mean Sea-Level Pressure (with Highs and Lows) %s' % date) -plt.show() +from mpl_toolkits.basemap import Basemap +from mpl_toolkits.basemap import addcyclic +from scipy.ndimage import minimum_filter +from scipy.ndimage import maximum_filter + + +def extrema(mat, mode="wrap", window=10): + """Find the indices of local extrema (min and max) in the input array.""" + + minimum = minimum_filter(mat, size=window, mode=mode) + maximum = maximum_filter(mat, size=window, mode=mode) + + # Return the indices of the maxima, minima. + # (mat == maximum) true if pixel is equal to the local max. + # (mat == minimum) true if pixel is equal to the local in. + return np.nonzero(mat == minimum), np.nonzero(mat == maximum) + + +def main(): + """Main function.""" + + # Plot 00 UTC today. + url = "http://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs%Y%m%d/gfs_0p25_00z" + date = dt.datetime.now() + + # Open OPeNDAP dataset. + data = netCDF4.Dataset(date.strftime(url)) + + # Read lats and lons. + lats = data.variables["lat"][:] + lons1 = data.variables["lon"][:] + + # Read prmsl and convert to hPa (mbar). + 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=50) + + # Create Basemap instance. + bmap = Basemap(projection="mill", + llcrnrlon=0, llcrnrlat=-80, + urcrnrlon=360, urcrnrlat=80) + + # Add wrap-around point in longitude. + prmsl, lons = addcyclic(prmsl, lons1) + + # Define contour levels. + clevs = np.arange(900, 1100., 5.) + + # Find x, y of map projection grid. + lons, lats = np.meshgrid(lons, lats) + x, y = bmap(lons, lats) + + # Create figure. + fig = plt.figure(figsize=(8, 4.5)) + fig.add_axes([0.05, 0.05, 0.9, 0.85]) + bmap.contour(x, y, prmsl, clevs, colors="k", linewidths=1.0) + bmap.drawcoastlines(linewidth=1.25) + bmap.fillcontinents(color="0.8") + bmap.drawparallels(np.arange(-80, 81, 20), labels=[1, 1, 0, 0]) + bmap.drawmeridians(np.arange(0, 360, 60), labels=[0, 0, 0, 1]) + xlows, xhighs = x[local_min], x[local_max] + ylows, yhighs = y[local_min], y[local_max] + lowvals, highvals = prmsl[local_min], prmsl[local_max] + + # Plot lows as blue L's, with min pressure value underneath. + # Do not plot if there is already a L or H within dmin meters. + xyplotted = [] + yoffset = 0.022 * (bmap.ymax - bmap.ymin) + dmin = yoffset + for x, y, p in zip(xlows, ylows, lowvals): + if bmap.xmin < x < bmap.xmax and bmap.ymin < y < bmap.ymax: + dist = [np.sqrt((x - x0)**2 + (y - y0)**2) for x0, y0 in xyplotted] + if not dist or min(dist) > dmin: + bbox = dict(boxstyle="square", ec="None", fc=(1, 1, 1, 0.5)) + plt.text(x, y, "L", fontsize=14, fontweight="bold", + ha="center", va="center", color="b") + plt.text(x, y - yoffset, repr(int(p)), fontsize=9, + ha="center", va="top", color="b", bbox=bbox) + 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 bmap.xmin < x < bmap.xmax and bmap.ymin < y < bmap.ymax: + dist = [np.sqrt((x - x0)**2 + (y - y0)**2) for x0, y0 in xyplotted] + if not dist or min(dist) > dmin: + bbox = dict(boxstyle="square", ec="None", fc=(1, 1, 1, 0.5)) + plt.text(x, y, "H", fontsize=14, fontweight="bold", + ha="center", va="center", color="r") + plt.text(x, y - yoffset, repr(int(p)), fontsize=9, + ha="center", va="top", color="r", bbox=bbox) + xyplotted.append((x, y)) + + # Set plot title and show. + plt.title("Mean Sea-Level Pressure (with Highs and Lows) %s" % date) + plt.show() + + +if __name__ == "__main__": + + import sys + sys.exit(main()) From 2e2544c3915a01852504e4f2ccfb3399bb7fc823 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?V=C3=ADctor=20Molina=20Garc=C3=ADa?= Date: Wed, 19 Oct 2022 20:55:32 +0200 Subject: [PATCH 7/7] Simplify main call in some of the examples --- examples/fcstmaps.py | 3 +-- examples/fcstmaps_axesgrid.py | 3 +-- examples/plothighsandlows.py | 4 +--- 3 files changed, 3 insertions(+), 7 deletions(-) diff --git a/examples/fcstmaps.py b/examples/fcstmaps.py index c942c339b..12796fcae 100644 --- a/examples/fcstmaps.py +++ b/examples/fcstmaps.py @@ -109,5 +109,4 @@ def main(date, verbose=True): dateobj = dt.datetime.strptime(sys.argv[1], "%Y%m%d") else: dateobj = dt.datetime.today() - - sys.exit(main(dateobj, verbose=True)) + main(dateobj, verbose=True) diff --git a/examples/fcstmaps_axesgrid.py b/examples/fcstmaps_axesgrid.py index 2748a1bac..7e2f226cf 100644 --- a/examples/fcstmaps_axesgrid.py +++ b/examples/fcstmaps_axesgrid.py @@ -123,5 +123,4 @@ def main(date, verbose=True): dateobj = dt.datetime.strptime(sys.argv[1], "%Y%m%d") else: dateobj = dt.datetime.today() - - sys.exit(main(dateobj, verbose=True)) + main(dateobj, verbose=True) diff --git a/examples/plothighsandlows.py b/examples/plothighsandlows.py index edacc6875..f86fddda6 100644 --- a/examples/plothighsandlows.py +++ b/examples/plothighsandlows.py @@ -105,6 +105,4 @@ def main(): if __name__ == "__main__": - - import sys - sys.exit(main()) + main()