`earthpy.plot.plot_bands` method ignores `ax` keyword argument if raster has multiple bands

[ahasha]

Describe the bug earthpy.plot.plot_bands method ignores ax keyword argument if raster has multiple bands

To Reproduce Steps to reproduce the behavior:

If raster is an array with raster.ndim > 2 and raster.shape[0] > 1, then the ax keyword argument has no effect and the plot_bands method produces a new axes array. The user has no access to the figure object to customize the plot.

import earthpy.plot as ep
import matplotlib.pyplot as plt

fig, axs = plt.subplots(3, 1)
ep.plot_bands(raster, ax=axs)

will produce an output that looks like

Expected behavior

Ideally, if an axis grid consistent with the number of bands is passed in the ax keyword argument, the method would use these axes objects as expected. If this is not feasible, an error message indicating that the ax keyword is not supported when multi-band arrays are provided is needed to avoid confusion.

**What Operating System Are you Running? **

• Mac

[ahasha]

I have a branch in my fork with a potential solution to this issue: https://github.com/ahasha/earthpy/tree/issue_796

Please let me know if a PR would be welcome. Thanks!

[nkorinek]

Hey @ahasha! Seeing this now, the solution on your branch looks awesome. A PR would be very welcome, might ask for some changes before I merge it but go ahead and get it started!