Plot on curvelinear grid

Preface

If you are working with radar station data, it is almost ever only available as polar data. This means you have a 2D-array, one dimension holding the azimuth (PPI) or elevation (RHI) angle values and the other holding the range values.

In \(\omega radlib\) it is assumed that the first dimension is over the azimuth/elevation angles, while the second dimension is over the range bins.

Create Curvelinear Grid

The creation process of the curvelinear grid is bundled in the helper function wradlib.vis.create_cg(). I will not dwell too much on that, just this far wradlib.vis.create_cg() uses a derived Axes implementation.

wradlib.vis.create_cg() takes scan type (‘PPI’ or ‘RHI’) as argument, figure object and grid definition are optional. The grid creation process generates three axes objects and set some reasonable starting values for labeling.

The returned objects are cgax, caax and paax.

  • cgax: matplotlib toolkit axisartist Axes object, Curvelinear Axes which holds the angle-range-grid

  • caax: matplotlib Axes object (twin to cgax), Cartesian Axes (x-y-grid) for plotting cartesian data

  • paax: matplotlib Axes object (parasite to cgax), The parasite axes object for plotting polar data

A typical invocation of wradlib.vis.create_cg() for a PPI is:

# create curvelinear axes
cgax, caax, paax = create_cg('PPI', fig, subplot)

For plotting actual polar data two functions exist, depending on whether your data holds a PPI wradlib.vis.plot_ppi() or an RHI (wradlib.vis.plot_rhi()).

Note

  1. Other than most plotting functions you cannot give an axes object as an argument. All necessary axes objects are created on the fly. You may give an figure object and/or an subplot specification as parameter. For further information on howto plot multiple cg plots in one figure, have a look at the special section Plotting on Grids.

  2. When using the refrac keyword with wradlib.vis.plot_rhi() the data is plotted to the cartesian axis caax.

Plotting on Curvelinear Grids

Plot CG PPI

wradlib.vis.plot_ppi() with keyword cg=True is used in this section.

Simple CG PPI

First we will look into plotting a PPI. We start with importing the necessary modules:

[1]:
import wradlib as wrl
import matplotlib.pyplot as pl
import warnings

warnings.filterwarnings("ignore")
try:
    get_ipython().magic("matplotlib inline")
except:
    pl.ion()
import numpy as np
/home/runner/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Next, we will load a polar scan from the WRADLIB_DATA folder and prepare it:

[2]:
# load a polar scan
filename = wrl.util.get_wradlib_data_file("misc/polar_dBZ_tur.gz")
data = np.loadtxt(filename)

# create range and azimuth arrays accordingly
r = np.arange(0, data.shape[1], dtype=np.float)
r += (r[1] - r[0]) / 2.0
r *= 1000.0
az = np.arange(0, data.shape[0], dtype=np.float)
az += (az[1] - az[0]) / 2.0

# mask data array for better presentation
mask_ind = np.where(data <= np.nanmin(data))
data[mask_ind] = np.nan
ma = np.ma.array(data, mask=np.isnan(data))

For this simple example, we do not need the returned axes. The plotting routine would be invoked like this:

[3]:
fig = pl.figure(figsize=(10, 8))
ax, pm = wrl.vis.plot_ppi(ma, fig=fig, proj="cg")
caax = ax.parasites[0]
paax = ax.parasites[1]
ax.parasites[1].set_aspect("equal")
t = pl.title("Simple CG PPI", y=1.05)
../../_images/notebooks_visualisation_wradlib_plot_curvelinear_grids_17_0.png

Decorated CG PPI

Now we will make use of some of the capabilities of this curvelinear axes.

You see, that for labeling x- and y-axis the cartesian axis is used. The azimuth label is set via :func:text. Also a colorbar is easily added. The plotting routine would be invoked like this, adding range and azimuth arrays:

[4]:
fig = pl.figure(figsize=(10, 8))
cgax, pm = wrl.vis.plot_ppi(ma, r=r, az=az, fig=fig, proj="cg")
caax = cgax.parasites[0]
paax = cgax.parasites[1]

pl.title("Decorated CG PPI", y=1.05)
cbar = pl.colorbar(pm, pad=0.075, ax=paax)
caax.set_xlabel("x_range [km]")
caax.set_ylabel("y_range [km]")
pl.text(1.0, 1.05, "azimuth", transform=caax.transAxes, va="bottom", ha="right")
cbar.set_label("reflectivity [dBZ]")
../../_images/notebooks_visualisation_wradlib_plot_curvelinear_grids_20_0.png

And, we will use cg keyword to set the starting value for the curvelinear grid. This is because data at the center of the image is obscured by the gridlines. We also adapt the radial_spacing to better align the two grids.

[5]:
cg = {"radial_spacing": 14.0, "latmin": 10e3}
fig = pl.figure(figsize=(10, 8))
cgax, pm = wrl.vis.plot_ppi(ma, r=r, az=az, fig=fig, proj="cg")
caax = cgax.parasites[0]
paax = cgax.parasites[1]

t = pl.title("Decorated CG PPI", y=1.05)
cbar = pl.gcf().colorbar(pm, pad=0.075, ax=paax)
caax.set_xlabel("x_range [km]")
caax.set_ylabel("y_range [km]")
pl.text(1.0, 1.05, "azimuth", transform=caax.transAxes, va="bottom", ha="right")
cbar.set_label("reflectivity [dBZ]")
../../_images/notebooks_visualisation_wradlib_plot_curvelinear_grids_22_0.png

Sector CG PPI

What if I want to plot only an interesting sector of the whole PPI? Not as easy, one might think. Note, that we can use infer_intervals = True here to get nice grid cell alignment. We also can generate a so called floating axis using the cgax now. Here we go:

[6]:
cg = {"angular_spacing": 20.0}
fig = pl.figure(figsize=(10, 8))
cgax, pm = wrl.vis.plot_ppi(
    ma[200:250, 40:80],
    r=r[40:80],
    az=az[200:250],
    fig=fig,
    proj=cg,
    rf=1e3,
    infer_intervals=True,
)
caax = cgax.parasites[0]
paax = cgax.parasites[1]

t = pl.title("Decorated Sector CG PPI", y=1.05)
cbar = pl.gcf().colorbar(pm, pad=0.075, ax=paax)
caax.set_xlabel("x_range [km]")
caax.set_ylabel("y_range [km]")
pl.text(1.0, 1.05, "azimuth", transform=caax.transAxes, va="bottom", ha="right")
cbar.set_label("reflectivity [dBZ]")

# add floating axis
cgax.axis["lat"] = cgax.new_floating_axis(0, 240)
cgax.axis["lat"].set_ticklabel_direction("-")
cgax.axis["lat"].label.set_text("range [km]")
cgax.axis["lat"].label.set_rotation(180)
cgax.axis["lat"].label.set_pad(10)
../../_images/notebooks_visualisation_wradlib_plot_curvelinear_grids_25_0.png

Special Markers

One more good thing about curvelinear axes is that you can plot polar as well as cartesian data. However, you have to be careful, where to plot. Polar data has to be plottet to the parasite axis (paax). Cartesian data can be plottet to caax, although you can also plot cartesian data to the main cgax.

Anyway, it is easy to overlay your polar data, with other station data (e.g. gauges). Taking the former sector example, we can plot some additional stations:

[7]:
fig = pl.figure(figsize=(10, 8))
cg = {"angular_spacing": 20.0}
cgax, pm = wrl.vis.plot_ppi(
    ma[200:250, 40:80],
    r[40:80],
    az[200:250],
    fig=fig,
    proj=cg,
    rf=1000.0,
    infer_intervals=True,
)
caax = cgax.parasites[0]
paax = cgax.parasites[1]
t = pl.title("Decorated Sector CG PPI", y=1.05)
cbar = pl.gcf().colorbar(pm, pad=0.075, ax=paax)
caax.set_xlabel("x_range [km]")
caax.set_ylabel("y_range [km]")
pl.text(1.0, 1.05, "azimuth", transform=caax.transAxes, va="bottom", ha="right")
cbar.set_label("reflectivity [dBZ]")
cgax.axis["lat"] = cgax.new_floating_axis(0, 240)
cgax.axis["lat"].set_ticklabel_direction("-")
cgax.axis["lat"].label.set_text("range [km]")
cgax.axis["lat"].label.set_rotation(180)
cgax.axis["lat"].label.set_pad(10)
# plot on cartesian axis
caax.plot(-60, -60, "ro", label="caax")
caax.plot(-50, -70, "ro")
# plot on polar axis
paax.plot(220, 88, "bo", label="paax")
# plot on cg axis (same as on cartesian axis)
cgax.plot(-60, -70, "go", label="cgax")
# legend on main cg axis
cgax.legend()
[7]:
<matplotlib.legend.Legend at 0x7fee7d0d1150>
../../_images/notebooks_visualisation_wradlib_plot_curvelinear_grids_28_1.png

Special Specials

But there is more to know, when using the curvelinear grids! As an example, you can get access to the underlying cgax and grid_helper to change azimuth and range resolution as well as tick labels:

[8]:
from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter

# cg = {'lon_cycle': 360.}
cg = {"angular_spacing": 20.0}
fig = pl.figure(figsize=(10, 8))
cgax, pm = wrl.vis.plot_ppi(
    ma[200:250, 40:80],
    r[40:80],
    az[200:250],
    rf=1e3,
    fig=fig,
    proj=cg,
    infer_intervals=True,
)
caax = cgax.parasites[0]
paax = cgax.parasites[1]

t = pl.title("Decorated Sector CG PPI", y=1.05)
t.set_y(1.05)
cbar = pl.gcf().colorbar(pm, pad=0.075, ax=paax)
caax.set_xlabel("x_range [km]")
caax.set_ylabel("y_range [km]")
pl.text(1.0, 1.05, "azimuth", transform=caax.transAxes, va="bottom", ha="right")
cbar.set_label("reflectivity [dBZ]")
gh = cgax.get_grid_helper()
# set azimuth resolution to 15deg
locs = [i for i in np.arange(0.0, 360.0, 5.0)]
gh.grid_finder.grid_locator1 = FixedLocator(locs)
gh.grid_finder.tick_formatter1 = DictFormatter(
    dict([(i, r"${0:.0f}^\circ$".format(i)) for i in locs])
)
gh.grid_finder.grid_locator2._nbins = 20
gh.grid_finder.grid_locator2._steps = [1, 1.5, 2, 2.5, 5, 10]
cgax.axis["lat"] = cgax.new_floating_axis(0, 240)
cgax.axis["lat"].set_ticklabel_direction("-")
cgax.axis["lat"].label.set_text("range [km]")
cgax.axis["lat"].label.set_rotation(180)
cgax.axis["lat"].label.set_pad(10)
../../_images/notebooks_visualisation_wradlib_plot_curvelinear_grids_31_0.png

The use of FixedLocator and DictFormatter should be clear. The use of _nbins and _steps is a bit of head-twisting. With _steps you can set the possible divisions of the range. In connection with the _nbins the range grid is created depending on maximum range. In the above situation with _nbins set to 10 we get an range grid resolution of 25 (divider 2.5). When setting steps to 20 we get a resolution of 15 (divider 1.5). Choosing 30 lead to resolution of 10 (divider 1/10). So it may be good to play around a bit, for wanted results.

As you might have noticed the cartesian grid remained the same and the azimuth labels are bit overplottet. But matplotlib would be not matplotlib if there would be no solution. First we take care of the labeling. We push the title a bit higher to get space and toggle the caax labels to right and top:

t = plt.title('Very Special Sector CG PPI', y=1.1)
caax.toggle_axisline()

Then we toggle “left” and “right” and “top” and “bottom” axis behaviour. We also have to put the colorbar a bit to the side and alter the location of the azimuth label. And, not to forgot to adapt the ticklabels of the cartesian axes. With little effort we got a better (IMHO) representation.

[9]:
fig = pl.figure(figsize=(12, 10), constrained_layout=True)
cg = {"angular_spacing": 20.0}
cgax, pm = wrl.vis.plot_ppi(
    ma[200:251, 40:81],
    r[40:81],
    az[200:251],
    rf=1e3,
    fig=fig,
    proj=cg,
    infer_intervals=True,
)
caax = cgax.parasites[0]
paax = cgax.parasites[1]

t = pl.title("Very Special Sector CG PPI", y=1.1)
cbar = pl.gcf().colorbar(pm, pad=0.1, ax=paax)
pl.text(0.5, 1.05, "x_range [km]", transform=caax.transAxes, va="bottom", ha="center")
pl.text(
    1.1,
    0.5,
    "y_range [km]",
    transform=caax.transAxes,
    va="bottom",
    ha="center",
    rotation="vertical",
)
caax.set_xlabel("x_range [km]")
caax.set_ylabel("y_range [km]")
caax.toggle_axisline()


# make ticklabels of right and top axis visible
caax.axis["top", "right"].set_visible(True)
caax.axis["top", "right"].major_ticklabels.set_visible(True)
caax.grid(True)

from matplotlib.ticker import MaxNLocator

caax.xaxis.set_major_locator(MaxNLocator(15))
caax.yaxis.set_major_locator(MaxNLocator(15))

# make ticklabels of left and bottom axis visible,
cgax.axis["left"].major_ticklabels.set_visible(True)
cgax.axis["bottom"].major_ticklabels.set_visible(True)
cgax.axis["left"].get_helper().nth_coord_ticks = 0
cgax.axis["bottom"].get_helper().nth_coord_ticks = 0
# and also set tickmarklength to zero for better presentation
cgax.axis["right"].major_ticks.set_ticksize(0)
cgax.axis["top"].major_ticks.set_ticksize(0)
# make ticklabels of right and top axis unvisible,
cgax.axis["right"].major_ticklabels.set_visible(False)
cgax.axis["top"].major_ticklabels.set_visible(False)
# and also set tickmarklength to zero for better presentation
cgax.axis["right"].major_ticks.set_ticksize(0)
cgax.axis["top"].major_ticks.set_ticksize(0)
pl.text(0.5, -0.065, "azimuth", transform=caax.transAxes, va="bottom", ha="center")
pl.text(
    -0.1,
    0.5,
    "azimuth",
    transform=caax.transAxes,
    va="bottom",
    ha="center",
    rotation="vertical",
)
cbar.set_label("reflectivity [dBZ]")

gh = cgax.get_grid_helper()
# set azimuth resolution to 5deg
locs = [i for i in np.arange(0.0, 360.0, 5.0)]
gh.grid_finder.grid_locator1 = FixedLocator(locs)
gh.grid_finder.tick_formatter1 = DictFormatter(
    dict([(i, r"${0:.0f}^\circ$".format(i)) for i in locs])
)
gh.grid_finder.grid_locator2._nbins = 30
gh.grid_finder.grid_locator2._steps = [1, 1.5, 2, 2.5, 5, 10]
cgax.axis["lat"] = cgax.new_floating_axis(0, 240)
cgax.axis["lat"].set_ticklabel_direction("-")
cgax.axis["lat"].label.set_text("range [km]")
cgax.axis["lat"].label.set_rotation(180)
cgax.axis["lat"].label.set_pad(10)
Error in callback <function _draw_all_if_interactive at 0x7fee8965b250> (for post_execute):
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/artist.py:1378, in _get_tightbbox_for_layout_only(obj, *args, **kwargs)
   1377 try:
-> 1378     return obj.get_tightbbox(*args, **{**kwargs, "for_layout_only": True})
   1379 except TypeError:

TypeError: HostAxesBase.get_tightbbox() got an unexpected keyword argument 'for_layout_only'

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/pyplot.py:119, in _draw_all_if_interactive()
    117 def _draw_all_if_interactive():
    118     if matplotlib.is_interactive():
--> 119         draw_all()

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/_pylab_helpers.py:132, in Gcf.draw_all(cls, force)
    130 for manager in cls.get_all_fig_managers():
    131     if force or manager.canvas.figure.stale:
--> 132         manager.canvas.draw_idle()

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/backend_bases.py:2054, in FigureCanvasBase.draw_idle(self, *args, **kwargs)
   2052 if not self._is_idle_drawing:
   2053     with self._idle_draw_cntx():
-> 2054         self.draw(*args, **kwargs)

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/backends/backend_agg.py:405, in FigureCanvasAgg.draw(self)
    401 # Acquire a lock on the shared font cache.
    402 with RendererAgg.lock, \
    403      (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
    404       else nullcontext()):
--> 405     self.figure.draw(self.renderer)
    406     # A GUI class may be need to update a window using this draw, so
    407     # don't forget to call the superclass.
    408     super().draw()

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/artist.py:74, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     72 @wraps(draw)
     73 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 74     result = draw(artist, renderer, *args, **kwargs)
     75     if renderer._rasterizing:
     76         renderer.stop_rasterizing()

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/artist.py:51, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     48     if artist.get_agg_filter() is not None:
     49         renderer.start_filter()
---> 51     return draw(artist, renderer)
     52 finally:
     53     if artist.get_agg_filter() is not None:

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/figure.py:3065, in Figure.draw(self, renderer)
   3063 if self.axes and self.get_layout_engine() is not None:
   3064     try:
-> 3065         self.get_layout_engine().execute(self)
   3066     except ValueError:
   3067         pass

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/layout_engine.py:255, in ConstrainedLayoutEngine.execute(self, fig)
    252 w_pad = self._params['w_pad'] / width
    253 h_pad = self._params['h_pad'] / height
--> 255 return do_constrained_layout(fig, w_pad=w_pad, h_pad=h_pad,
    256                              wspace=self._params['wspace'],
    257                              hspace=self._params['hspace'],
    258                              rect=self._params['rect'],
    259                              compress=self._compress)

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/_constrained_layout.py:119, in do_constrained_layout(fig, h_pad, w_pad, hspace, wspace, rect, compress)
    109     return
    111 for _ in range(2):
    112     # do the algorithm twice.  This has to be done because decorations
    113     # change size after the first re-position (i.e. x/yticklabels get
   (...)
    117     # make margins for all the axes and subfigures in the
    118     # figure.  Add margins for colorbars...
--> 119     make_layout_margins(layoutgrids, fig, renderer, h_pad=h_pad,
    120                         w_pad=w_pad, hspace=hspace, wspace=wspace)
    121     make_margin_suptitles(layoutgrids, fig, renderer, h_pad=h_pad,
    122                           w_pad=w_pad)
    124     # if a layout is such that a columns (or rows) margin has no
    125     # constraints, we need to make all such instances in the grid
    126     # match in margin size.

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/_constrained_layout.py:371, in make_layout_margins(layoutgrids, fig, renderer, w_pad, h_pad, hspace, wspace)
    367     return
    369 margin = get_margin_from_padding(ax, w_pad=w_pad, h_pad=h_pad,
    370                                  hspace=hspace, wspace=wspace)
--> 371 pos, bbox = get_pos_and_bbox(ax, renderer)
    372 # the margin is the distance between the bounding box of the axes
    373 # and its position (plus the padding from above)
    374 margin['left'] += pos.x0 - bbox.x0

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/_constrained_layout.py:600, in get_pos_and_bbox(ax, renderer)
    598 # pos is in panel co-ords, but we need in figure for the layout
    599 pos = pos.transformed(fig.transSubfigure - fig.transFigure)
--> 600 tightbbox = martist._get_tightbbox_for_layout_only(ax, renderer)
    601 if tightbbox is None:
    602     bbox = pos

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/artist.py:1380, in _get_tightbbox_for_layout_only(obj, *args, **kwargs)
   1378     return obj.get_tightbbox(*args, **{**kwargs, "for_layout_only": True})
   1379 except TypeError:
-> 1380     return obj.get_tightbbox(*args, **kwargs)

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/mpl_toolkits/axes_grid1/parasite_axes.py:223, in HostAxesBase.get_tightbbox(self, renderer, call_axes_locator, bbox_extra_artists)
    218 def get_tightbbox(self, renderer=None, call_axes_locator=True,
    219                   bbox_extra_artists=None):
    220     bbs = [
    221         *[ax.get_tightbbox(renderer, call_axes_locator=call_axes_locator)
    222           for ax in self.parasites],
--> 223         super().get_tightbbox(renderer,
    224                               call_axes_locator=call_axes_locator,
    225                               bbox_extra_artists=bbox_extra_artists)]
    226     return Bbox.union([b for b in bbs if b.width != 0 or b.height != 0])

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/axes/_base.py:4451, in _AxesBase.get_tightbbox(self, renderer, call_axes_locator, bbox_extra_artists, for_layout_only)
   4448     bbox_artists = self.get_default_bbox_extra_artists()
   4450 for a in bbox_artists:
-> 4451     bbox = a.get_tightbbox(renderer)
   4452     if (bbox is not None
   4453             and 0 < bbox.width < np.inf
   4454             and 0 < bbox.height < np.inf):
   4455         bb.append(bbox)

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/artist.py:337, in Artist.get_tightbbox(self, renderer)
    322 def get_tightbbox(self, renderer=None):
    323     """
    324     Like `.Artist.get_window_extent`, but includes any clipping.
    325
   (...)
    335         The enclosing bounding box (in figure pixel coordinates).
    336     """
--> 337     bbox = self.get_window_extent(renderer)
    338     if self.get_clip_on():
    339         clip_box = self.get_clip_box()

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/collections.py:306, in Collection.get_window_extent(self, renderer)
    303 def get_window_extent(self, renderer=None):
    304     # TODO: check to ensure that this does not fail for
    305     # cases other than scatter plot legend
--> 306     return self.get_datalim(transforms.IdentityTransform())

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/collections.py:262, in Collection.get_datalim(self, transData)
    259 offsets = self.get_offsets()
    261 paths = self.get_paths()
--> 262 if not len(paths):
    263     # No paths to transform
    264     return transforms.Bbox.null()
    266 if not transform.is_affine:

TypeError: object of type 'NoneType' has no len()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/artist.py:1378, in _get_tightbbox_for_layout_only(obj, *args, **kwargs)
   1377 try:
-> 1378     return obj.get_tightbbox(*args, **{**kwargs, "for_layout_only": True})
   1379 except TypeError:

TypeError: HostAxesBase.get_tightbbox() got an unexpected keyword argument 'for_layout_only'

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/IPython/core/formatters.py:338, in BaseFormatter.__call__(self, obj)
    336     pass
    337 else:
--> 338     return printer(obj)
    339 # Finally look for special method names
    340 method = get_real_method(obj, self.print_method)

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/IPython/core/pylabtools.py:152, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    149     from matplotlib.backend_bases import FigureCanvasBase
    150     FigureCanvasBase(fig)
--> 152 fig.canvas.print_figure(bytes_io, **kw)
    153 data = bytes_io.getvalue()
    154 if fmt == 'svg':

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/backend_bases.py:2314, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2308     renderer = _get_renderer(
   2309         self.figure,
   2310         functools.partial(
   2311             print_method, orientation=orientation)
   2312     )
   2313     with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2314         self.figure.draw(renderer)
   2316 if bbox_inches:
   2317     if bbox_inches == "tight":

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/artist.py:74, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     72 @wraps(draw)
     73 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 74     result = draw(artist, renderer, *args, **kwargs)
     75     if renderer._rasterizing:
     76         renderer.stop_rasterizing()

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/artist.py:51, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     48     if artist.get_agg_filter() is not None:
     49         renderer.start_filter()
---> 51     return draw(artist, renderer)
     52 finally:
     53     if artist.get_agg_filter() is not None:

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/figure.py:3065, in Figure.draw(self, renderer)
   3063 if self.axes and self.get_layout_engine() is not None:
   3064     try:
-> 3065         self.get_layout_engine().execute(self)
   3066     except ValueError:
   3067         pass

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/layout_engine.py:255, in ConstrainedLayoutEngine.execute(self, fig)
    252 w_pad = self._params['w_pad'] / width
    253 h_pad = self._params['h_pad'] / height
--> 255 return do_constrained_layout(fig, w_pad=w_pad, h_pad=h_pad,
    256                              wspace=self._params['wspace'],
    257                              hspace=self._params['hspace'],
    258                              rect=self._params['rect'],
    259                              compress=self._compress)

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/_constrained_layout.py:119, in do_constrained_layout(fig, h_pad, w_pad, hspace, wspace, rect, compress)
    109     return
    111 for _ in range(2):
    112     # do the algorithm twice.  This has to be done because decorations
    113     # change size after the first re-position (i.e. x/yticklabels get
   (...)
    117     # make margins for all the axes and subfigures in the
    118     # figure.  Add margins for colorbars...
--> 119     make_layout_margins(layoutgrids, fig, renderer, h_pad=h_pad,
    120                         w_pad=w_pad, hspace=hspace, wspace=wspace)
    121     make_margin_suptitles(layoutgrids, fig, renderer, h_pad=h_pad,
    122                           w_pad=w_pad)
    124     # if a layout is such that a columns (or rows) margin has no
    125     # constraints, we need to make all such instances in the grid
    126     # match in margin size.

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/_constrained_layout.py:371, in make_layout_margins(layoutgrids, fig, renderer, w_pad, h_pad, hspace, wspace)
    367     return
    369 margin = get_margin_from_padding(ax, w_pad=w_pad, h_pad=h_pad,
    370                                  hspace=hspace, wspace=wspace)
--> 371 pos, bbox = get_pos_and_bbox(ax, renderer)
    372 # the margin is the distance between the bounding box of the axes
    373 # and its position (plus the padding from above)
    374 margin['left'] += pos.x0 - bbox.x0

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/_constrained_layout.py:600, in get_pos_and_bbox(ax, renderer)
    598 # pos is in panel co-ords, but we need in figure for the layout
    599 pos = pos.transformed(fig.transSubfigure - fig.transFigure)
--> 600 tightbbox = martist._get_tightbbox_for_layout_only(ax, renderer)
    601 if tightbbox is None:
    602     bbox = pos

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/artist.py:1380, in _get_tightbbox_for_layout_only(obj, *args, **kwargs)
   1378     return obj.get_tightbbox(*args, **{**kwargs, "for_layout_only": True})
   1379 except TypeError:
-> 1380     return obj.get_tightbbox(*args, **kwargs)

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/mpl_toolkits/axes_grid1/parasite_axes.py:223, in HostAxesBase.get_tightbbox(self, renderer, call_axes_locator, bbox_extra_artists)
    218 def get_tightbbox(self, renderer=None, call_axes_locator=True,
    219                   bbox_extra_artists=None):
    220     bbs = [
    221         *[ax.get_tightbbox(renderer, call_axes_locator=call_axes_locator)
    222           for ax in self.parasites],
--> 223         super().get_tightbbox(renderer,
    224                               call_axes_locator=call_axes_locator,
    225                               bbox_extra_artists=bbox_extra_artists)]
    226     return Bbox.union([b for b in bbs if b.width != 0 or b.height != 0])

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/axes/_base.py:4451, in _AxesBase.get_tightbbox(self, renderer, call_axes_locator, bbox_extra_artists, for_layout_only)
   4448     bbox_artists = self.get_default_bbox_extra_artists()
   4450 for a in bbox_artists:
-> 4451     bbox = a.get_tightbbox(renderer)
   4452     if (bbox is not None
   4453             and 0 < bbox.width < np.inf
   4454             and 0 < bbox.height < np.inf):
   4455         bb.append(bbox)

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/artist.py:337, in Artist.get_tightbbox(self, renderer)
    322 def get_tightbbox(self, renderer=None):
    323     """
    324     Like `.Artist.get_window_extent`, but includes any clipping.
    325
   (...)
    335         The enclosing bounding box (in figure pixel coordinates).
    336     """
--> 337     bbox = self.get_window_extent(renderer)
    338     if self.get_clip_on():
    339         clip_box = self.get_clip_box()

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/collections.py:306, in Collection.get_window_extent(self, renderer)
    303 def get_window_extent(self, renderer=None):
    304     # TODO: check to ensure that this does not fail for
    305     # cases other than scatter plot legend
--> 306     return self.get_datalim(transforms.IdentityTransform())

File ~/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/matplotlib/collections.py:262, in Collection.get_datalim(self, transData)
    259 offsets = self.get_offsets()
    261 paths = self.get_paths()
--> 262 if not len(paths):
    263     # No paths to transform
    264     return transforms.Bbox.null()
    266 if not transform.is_affine:

TypeError: object of type 'NoneType' has no len()
<Figure size 1200x1000 with 2 Axes>

Plot CG RHI

wradlib.vis.plot_rhi() is used in this section. An CG RHI plot is a little different compared to an CG PPI plot. I covers only one quadrant and the data is plottet counterclockwise from “east” (3 o’clock) to “north” (12 o’clock).

Everything else is much the same and you can do whatever you want as shown in the section Plot CG PPI.

So just a quick example of an cg rhi plot with some decorations. Note, the grid_locator1 for the theta angles is overwritten and now the grid is much finer.

[10]:
from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter

# reading in GAMIC hdf5 file
filename = wrl.util.get_wradlib_data_file("hdf5/2014-06-09--185000.rhi.mvol")
data1, metadata = wrl.io.read_gamic_hdf5(filename)
data1 = data1["SCAN0"]["ZH"]["data"]
r = metadata["SCAN0"]["r"]
th = metadata["SCAN0"]["el"]
az = metadata["SCAN0"]["az"]
site = (
    metadata["VOL"]["Longitude"],
    metadata["VOL"]["Latitude"],
    metadata["VOL"]["Height"],
)
# mask data array for better presentation
mask_ind = np.where(data1 <= np.nanmin(data1))
data1[mask_ind] = np.nan
ma1 = np.ma.array(data1, mask=np.isnan(data1))
[11]:
fig = pl.figure(figsize=(10, 8))
cgax, pm = wrl.vis.plot_rhi(ma1, r=r, th=th, rf=1e3, fig=fig, ax=111, proj="cg")
caax = cgax.parasites[0]
paax = cgax.parasites[1]

t = pl.title("Decorated CG RHI", y=1.05)
cgax.set_ylim(0, 14)
cbar = pl.gcf().colorbar(pm, pad=0.05, ax=paax)
cbar.set_label("reflectivity [dBZ]")
caax.set_xlabel("x_range [km]")
caax.set_ylabel("y_range [km]")
pl.text(1.0, 1.05, "azimuth", transform=caax.transAxes, va="bottom", ha="right")
gh = cgax.get_grid_helper()

# set theta to some nice values
locs = [
    0.0,
    1.0,
    2.0,
    3.0,
    4.0,
    5.0,
    6.0,
    7.0,
    8.0,
    9.0,
    10.0,
    11.0,
    12.0,
    13.0,
    14.0,
    15.0,
    16.0,
    17.0,
    18.0,
    20.0,
    22.0,
    25.0,
    30.0,
    35.0,
    40.0,
    50.0,
    60.0,
    70.0,
    80.0,
    90.0,
]
gh.grid_finder.grid_locator1 = FixedLocator(locs)
gh.grid_finder.tick_formatter1 = DictFormatter(
    dict([(i, r"${0:.0f}^\circ$".format(i)) for i in locs])
)
../../_images/notebooks_visualisation_wradlib_plot_curvelinear_grids_38_0.png

Plotting on Grids

There are serveral possibilities to plot multiple cg plots in one figure. Since both plotting routines are equipped with the same mechanisms it is concentrated mostly on RHI plots.

Note

Using the :func:tight_layout and :func:subplots_adjust functions most alignment problems can be avoided.

The Built-In Method

Using the matplotlib grid definition for the parameter subplot, we can easily plot two or more plots in one figure on a regular grid:

[12]:
subplots = [221, 222, 223, 224]
fig = pl.figure(figsize=(10, 8))
fig.subplots_adjust(wspace=0.2, hspace=0.35)
for sp in subplots:
    cgax, pm = wrl.vis.plot_rhi(ma1, r, th, rf=1e3, ax=sp, proj="cg")
    caax = cgax.parasites[0]
    paax = cgax.parasites[1]
    t = pl.title("CG RHI #%(sp)d" % locals(), y=1.1)
    cgax.set_ylim(0, 15)
    cbar = pl.gcf().colorbar(pm, pad=0.125, ax=paax)
    caax.set_xlabel("range [km]")
    caax.set_ylabel("height [km]")
    gh = cgax.get_grid_helper()
    # set theta to some nice values
    locs = [0.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, 60.0, 90.0]
    gh.grid_finder.grid_locator1 = FixedLocator(locs)
    gh.grid_finder.tick_formatter1 = DictFormatter(
        dict([(i, r"${0:.0f}^\circ$".format(i)) for i in locs])
    )
    cbar.set_label("reflectivity [dBZ]")
../../_images/notebooks_visualisation_wradlib_plot_curvelinear_grids_45_0.png

The GridSpec Method

Here the abilities of Matplotlib GridSpec are used. Now we can also plot on irregular grids. Just create your grid and take the GridSpec object as an input to the parameter ax as follows (some padding has to be adjusted to get a nice plot):

[13]:
import matplotlib.gridspec as gridspec

gs = gridspec.GridSpec(3, 3, hspace=0.75, wspace=0.4)
subplots = [gs[0, :], gs[1, :-1], gs[1:, -1], gs[-1, 0], gs[-1, -2]]
cbarpad = [0.05, 0.075, 0.2, 0.2, 0.2]
labelpad = [1.25, 1.25, 1.1, 1.25, 1.25]
fig = pl.figure(figsize=(10, 8))
for i, sp in enumerate(subplots):
    cgax, pm = wrl.vis.plot_rhi(ma1, r, th, rf=1e3, ax=sp, proj="cg")
    caax = cgax.parasites[0]
    paax = cgax.parasites[1]
    t = pl.title("CG RHI #%(i)d" % locals(), y=labelpad[i])
    cgax.set_ylim(0, 15)
    cbar = fig.colorbar(pm, pad=cbarpad[i], ax=paax)
    caax.set_xlabel("range [km]")
    caax.set_ylabel("height [km]")
    gh = cgax.get_grid_helper()
    # set theta to some nice values
    locs = [0.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, 60.0, 90.0]
    gh.grid_finder.grid_locator1 = FixedLocator(locs)
    gh.grid_finder.tick_formatter1 = DictFormatter(
        dict([(i, r"${0:.0f}^\circ$".format(i)) for i in locs])
    )
    cbar.set_label("reflectivity [dBZ]")
../../_images/notebooks_visualisation_wradlib_plot_curvelinear_grids_48_0.png

The AxesDivider Method

Here the capabilities of Matplotlib AxesGrid1 are used.

We make a PPI now, it matches much better. Just plot your PPI data and create an axes divider:

from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.ticker import NullFormatter, FuncFormatter, MaxNLocator
divider = make_axes_locatable(cgax)

Now you can easily append more axes to plot some other things, eg a maximum intensity projection:

axMipX = divider.append_axes("top", size=1.2, pad=0.1, sharex=cgax))
axMipY = divider.append_axes("right", size=1.2, pad=0.1, sharey=cgax))

OK, we have to create the mip data, we use the wradlib.georef.polar.maximum_intensity_projection():

[14]:
# angle of *cut* through ppi and scan elev.
angle = 0.0
elev = 0.0

filename = wrl.util.get_wradlib_data_file("misc/polar_dBZ_tur.gz")
data2 = np.loadtxt(filename)
# we need to have meter here for the georef function inside mip
d1 = np.arange(data2.shape[1], dtype=np.float) * 1000
d2 = np.arange(data2.shape[0], dtype=np.float)
data2 = np.roll(data2, (d2 >= angle).nonzero()[0][0], axis=0)

# calculate max intensity proj
xs, ys, mip1 = wrl.util.maximum_intensity_projection(
    data2, r=d1, az=d2, angle=angle, elev=elev
)
xs, ys, mip2 = wrl.util.maximum_intensity_projection(
    data2, r=d1, az=d2, angle=90 + angle, elev=elev
)

We also need a new formatter:

[15]:
def mip_formatter(x, pos):
    x = x / 1000.0
    fmt_str = "{:g}".format(x)
    if np.abs(x) > 0 and np.abs(x) < 1:
        return fmt_str.replace("0", "", 1)
    else:
        return fmt_str

Warning

As of matplotlib version 3.3.0 AxesDivider Method is broken for plots with twin/parasite axes. The following plot is done with normal axes.

[16]:
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.ticker import NullFormatter, FuncFormatter, MaxNLocator

fig = pl.figure(figsize=(10, 8))
# normal cg plot
cgax, pm = wrl.vis.plot_ppi(data2, r=d1, az=d2, fig=fig)  # , proj={'latmin': 10000.,
#       'radial_spacing': 12})
caax = cgax  # .parasites[0]
paax = cgax  # .parasites[1]

cgax.grid(True)

cgax.set_xlim(-np.max(d1), np.max(d1))
cgax.set_ylim(-np.max(d1), np.max(d1))
caax.xaxis.set_major_formatter(FuncFormatter(mip_formatter))
caax.yaxis.set_major_formatter(FuncFormatter(mip_formatter))
caax.set_xlabel("x_range [km]")
caax.set_ylabel("y_range [km]")

# axes divider section
divider = make_axes_locatable(cgax)
axMipX = divider.append_axes("top", size=1.2, pad=0.1, sharex=cgax)
axMipY = divider.append_axes("right", size=1.2, pad=0.1, sharey=cgax)

# special handling for labels etc.
# cgax.axis["right"].major_ticklabels.set_visible(False)
# cgax.axis["top"].major_ticklabels.set_visible(False)
axMipX.xaxis.set_major_formatter(NullFormatter())
axMipX.yaxis.set_major_formatter(FuncFormatter(mip_formatter))
axMipX.yaxis.set_major_locator(MaxNLocator(5))
axMipY.yaxis.set_major_formatter(NullFormatter())
axMipY.xaxis.set_major_formatter(FuncFormatter(mip_formatter))
axMipY.xaxis.set_major_locator(MaxNLocator(5))

# plot max intensity proj
ma = np.ma.array(mip1, mask=np.isnan(mip1))
axMipX.pcolormesh(xs, ys, ma)
ma = np.ma.array(mip2, mask=np.isnan(mip2))
axMipY.pcolormesh(ys.T, xs.T, ma.T)

# set labels, limits etc
axMipX.set_xlim(-np.max(d1), np.max(d1))
er = 6370000.0
axMipX.set_ylim(0, wrl.georef.bin_altitude(d1[-2], elev, 0, re=er))
axMipY.set_xlim(0, wrl.georef.bin_altitude(d1[-2], elev, 0, re=er))
axMipY.set_ylim(-np.max(d1), np.max(d1))
axMipX.set_ylabel("height [km]")
axMipY.set_xlabel("height [km]")
axMipX.grid(True)
axMipY.grid(True)
t = pl.gcf().suptitle("AxesDivider MIP Example")
t.set_y(0.925)
../../_images/notebooks_visualisation_wradlib_plot_curvelinear_grids_56_0.png