Repository: OSOceanAcoustics/echopype (language: python)

Issue:
Retire the `visualize` subpackage?
I am wondering if we should retire the subpackage `echopype.visualize` (`epviz`), and keep the EK500 colormap import. [Current docs here](https://echopype.readthedocs.io/en/stable/viz.html).

There are a few main reasons:
- The current epviz functionality tries to guess automatically what users want to do if they miss some input arguments, but the handling is brittle if the data are not in one of the predicted form. Given the flexibility of EK80 configuration, the "guess" often breaks.
- It seems that instead of having a specialized subpackage, we can just do a better job in the docs to show users how to plot data using native xarray plotting functions. This way users will likely be able to explore it in a more flexible way.
- Now we have the [`Echoshader`](https://github.com/OSOceanAcoustics/echoshader) package that provides interactive visualization functionalities via the accessor mechanism that is more powerful for the users, so we can direct users there.

@valentina-s @lsetiawan: Any thoughts? I know you spent quite a bit of time developing epviz and writing the docs too.

### Tasks
- [x] deprecation warning at v0.8.4
- [x] remove subpackage at v0.9.0
  - [x] make sure to retain EK500 colormap via import
- [x] update docs
  - [x] add notice for removal (but retaining the page until a later version)
  - [x] add notice for how to import EK500 colormap



Discussion / hints:
Now that all the plotting functionalities are mostly within Echoshader, I agree that this `visualize` subpackage can be retired. If there's a way to port the EK500 colormap that would be really great!

> It seems that instead of having a specialized subpackage, we can just do a better job in the docs to show users how to plot data using native xarray plotting functions.

I really like this idea of showing users how to use xarray plotting functions to build their own static plots!
Ok, I'll put in a deprecating warning for v0.8.4 and redirect people to Echoshader, and remove it at v0.9.0 (since its current state is buggy and will be very involved to fix).

For docs, it'll take a bit of time to add in examples, but the first thing we can do is to link to xarray dataset/dataarray intros and its plotting documentation, so that people have a good place to find references.

Fix patch:
diff --git a/echopype/colormap/__init__.py b/echopype/colormap/__init__.py
new file mode 100644
index 000000000..03be43b2f
--- /dev/null
+++ b/echopype/colormap/__init__.py
@@ -0,0 +1,9 @@
+"""
+Colormaps for plotting echogram.
+
+**NOTE: To use this subpackage. `Matplotlib` and `cmocean` package must be installed.**
+"""
+
+from . import cm
+
+__all__ = ["cm"]
diff --git a/echopype/visualize/cm.py b/echopype/colormap/cm.py
similarity index 63%
rename from echopype/visualize/cm.py
rename to echopype/colormap/cm.py
index 15199169b..618e45747 100644
--- a/echopype/visualize/cm.py
+++ b/echopype/colormap/cm.py
@@ -1,9 +1,9 @@
-import numpy as np
 import matplotlib as mpl
+import numpy as np
 
 __cmap_colors = {
-    'ek500': {
-        'rgb': (
+    "ek500": {
+        "rgb": (
             np.array(
                 [
                     [159, 159, 159],  # light grey
@@ -21,8 +21,8 @@
             )
             / 255
         ),
-        'under': '1',  # white
-        'over': np.array([120, 60, 40]) / 255,  # dark brown
+        "under": "1",  # white
+        "over": np.array([120, 60, 40]) / 255,  # dark brown
     }
 }
 
@@ -39,46 +39,44 @@ def _create_cmap(rgb, under=None, over=None):
 
 
 cmap_d = {}
-cmapnames = ['ek500']
+cmapnames = ["ek500"]
 
 # add colormaps and reversed to dictionary
 for cmapname in cmapnames:
     colors_d = __cmap_colors[cmapname]
-    rgb = colors_d['rgb']
+    rgb = colors_d["rgb"]
     cmap_d[cmapname] = _create_cmap(
-        rgb, under=colors_d.get('under', None), over=colors_d.get('over', None)
+        rgb, under=colors_d.get("under", None), over=colors_d.get("over", None)
     )
     cmap_d[cmapname].name = cmapname
-    cmap_d[cmapname + '_r'] = _create_cmap(
+    cmap_d[cmapname + "_r"] = _create_cmap(
         rgb[::-1, :],
-        under=colors_d.get('over', None),
-        over=colors_d.get('under', None),
+        under=colors_d.get("over", None),
+        over=colors_d.get("under", None),
     )
-    cmap_d[cmapname + '_r'].name = cmapname + '_r'
+    cmap_d[cmapname + "_r"].name = cmapname + "_r"
 
     # Register the cmap with matplotlib
     rgb_with_alpha = np.zeros((rgb.shape[0], 4))
     rgb_with_alpha[:, :3] = rgb
     rgb_with_alpha[:, 3] = 1.0  # set alpha channel to 1
-    reg_map = mpl.colors.ListedColormap(
-        rgb_with_alpha, 'ep.' + cmapname, rgb.shape[0]
-    )
-    if 'under' in colors_d:
-        reg_map.set_under(colors_d['under'])
+    reg_map = mpl.colors.ListedColormap(rgb_with_alpha, "ep." + cmapname, rgb.shape[0])
+    if "under" in colors_d:
+        reg_map.set_under(colors_d["under"])
 
-    if 'over' in colors_d:
-        reg_map.set_over(colors_d['over'])
+    if "over" in colors_d:
+        reg_map.set_over(colors_d["over"])
     mpl.colormaps.register(cmap=reg_map)
 
     # Register the reversed map
     reg_map_r = mpl.colors.ListedColormap(
-        rgb_with_alpha[::-1, :], 'ep.' + cmapname + '_r', rgb.shape[0]
+        rgb_with_alpha[::-1, :], "ep." + cmapname + "_r", rgb.shape[0]
     )
-    if 'under' in colors_d:
-        reg_map_r.set_over(colors_d['under'])
+    if "under" in colors_d:
+        reg_map_r.set_over(colors_d["under"])
 
-    if 'over' in colors_d:
-        reg_map_r.set_under(colors_d['over'])
+    if "over" in colors_d:
+        reg_map_r.set_under(colors_d["over"])
     mpl.colormaps.register(cmap=reg_map_r)
 
 # make colormaps available to call
diff --git a/echopype/visualize/__init__.py b/echopype/visualize/__init__.py
deleted file mode 100644
index 9cb6855d7..000000000
--- a/echopype/visualize/__init__.py
+++ /dev/null
@@ -1,21 +0,0 @@
-"""
-Visualization module to quickly plot raw, Sv, and MVBS dataset.
-
-**NOTE: To use this subpackage. `Matplotlib` and `cmocean` package must be installed.**
-"""
-import warnings
-from .api import create_echogram
-from . import cm
-
-__all__ = ["create_echogram", "cm"]
-
-warnings.warn(
-    (
-        "`echopype.visualization` is deprecated and will be removed in an upcoming release. \n"
-        "Use Echoshader for visualization functions instead. \n"
-        "Repository: https://github.com/OSOceanAcoustics/echoshader \n"
-        "Echogram plotting: https://echoshader.readthedocs.io/en/latest/version_0.1.0/echogram_examples.html"  # noqa
-    ),
-    DeprecationWarning,
-    2,
-)
diff --git a/echopype/visualize/api.py b/echopype/visualize/api.py
deleted file mode 100644
index 9135f4d4c..000000000
--- a/echopype/visualize/api.py
+++ /dev/null
@@ -1,264 +0,0 @@
-from typing import Optional, Union, List, Type
-
-import xarray as xr
-
-from .plot import _plot_echogram, FacetGrid, QuadMesh
-from ..echodata import EchoData
-from ..calibrate.calibrate_ek import CalibrateEK60, CalibrateEK80
-from ..calibrate.calibrate_azfp import CalibrateAZFP
-from ..utils.log import _init_logger
-
-logger = _init_logger(__name__)
-
-
-def create_echogram(
-    data: Union[EchoData, xr.Dataset],
-    channel: Union[str, List[str], None] = None,
-    frequency: Union[str, List[str], None] = None,
-    get_range: Optional[bool] = None,
-    range_kwargs: dict = {},
-    vertical_offset: Union[int, float, xr.DataArray, bool, None] = None,
-    **kwargs,
-) -> List[Union[FacetGrid, QuadMesh]]:
-    """Create an Echogram from an EchoData object or Sv and MVBS Dataset.
-
-    Parameters
-    ----------
-    data : EchoData or xr.Dataset
-        Echodata or Xarray Dataset to be plotted
-    channel : str or list of str, optional
-        The channel to be plotted.
-        Otherwise all channels will be plotted.
-    frequency : int, float, or list of float or ints, optional
-        The frequency to be plotted.
-        If not specified, all frequency will be plotted.
-    get_range : bool, optional
-        Flag as to whether range (``echo_range``) should be computed or not,
-        by default it will just plot `range_sample`` as the yaxis.
-
-        Note that for data that is "Sv" xarray dataset, `get_range` defaults
-        to `True`.
-    range_kwargs : dict
-        Keyword arguments dictionary for computing range (``echo_range``).
-        Keys are `env_params`, `waveform_mode`, and `encode_mode`.
-    vertical_offset : int, float, xr.DataArray, or bool, optional
-        Water level data array for platform water level correction.
-        Note that auto addition of water level can be performed
-        when data is an EchoData object by setting this argument
-        to `True`. Currently because the water level information
-        is not available as part of the Sv dataset, a warning is issued
-        when `vertical_offset=True` in this case and no correction is
-        performed. This behavior will change in the future when the
-        default content of Sv dataset is updated to include this information.
-    **kwargs: optional
-        Additional keyword arguments for xarray plot pcolormesh.
-
-    Notes
-    -----
-    The EK80 echosounder can be configured to transmit
-    either broadband (``waveform_mode="BB"``)
-    or narrowband (``waveform_mode="CW"``) signals.
-    When transmitting in broadband mode, the returned echoes are
-    encoded as complex samples (``encode_mode="complex"``).
-    When transmitting in narrowband mode, the returned echoes can be encoded
-    either as complex samples (``encode_mode="complex"``)
-    or as power/angle combinations (``encode_mode="power"``) in a format
-    similar to those recorded by EK60 echosounders.
-
-    """
-    range_attrs = {
-        'long_name': 'Range',
-        'units': 'm',
-    }
-
-    if channel and frequency:
-        logger.warning(
-            "Both channel and frequency are specified. Channel filtering will be used."
-        )
-
-    if isinstance(channel, list) and len(channel) == 1:
-        channel = channel[0]
-    elif isinstance(frequency, list) and len(frequency) == 1:
-        frequency = frequency[0]
-
-    if isinstance(data, EchoData):
-        if data.sonar_model.lower() == 'ad2cp':
-            raise ValueError(
-                "Visualization for AD2CP sonar model is currently unsupported."
-            )
-        yaxis = 'range_sample'
-        variable = 'backscatter_r'
-        ds = data["Sonar/Beam_group1"]
-        if 'ping_time' in ds:
-            _check_ping_time(ds.ping_time)
-        if get_range is True:
-            yaxis = 'echo_range'
-
-            if data.sonar_model.lower() == 'azfp':
-                if 'azfp_cal_type' not in range_kwargs:
-                    range_kwargs['azfp_cal_type'] = 'Sv'
-                if 'env_params' not in range_kwargs:
-                    raise ValueError(
-                        "Please provide env_params in range_kwargs!"
-                    )
-            elif data.sonar_model.lower() == 'ek60':
-                if 'waveform_mode' not in range_kwargs:
-                    range_kwargs['waveform_mode'] = 'CW'
-                elif range_kwargs['waveform_mode'] != 'CW':
-                    raise ValueError(
-                        f"waveform_mode {range_kwargs['waveform_mode']} is invalid. EK60 waveform_mode must be 'CW'."  # noqa
-                    )
-
-                if 'encode_mode' not in range_kwargs:
-                    range_kwargs['encode_mode'] = 'power'
-                elif range_kwargs['encode_mode'] != 'power':
-                    raise ValueError(
-                        f"encode_mode {range_kwargs['encode_mode']} is invalid. EK60 encode_mode must be 'power'."  # noqa
-                    )
-            elif data.sonar_model.lower() == 'ek80':
-                if not all(
-                    True if mode in range_kwargs else False
-                    for mode in ['waveform_mode', 'encode_mode']
-                ):
-                    raise ValueError(
-                        "Please provide waveform_mode and encode_mode in range_kwargs for EK80 sonar model."  # noqa
-                    )
-                waveform_mode = range_kwargs['waveform_mode']
-                encode_mode = range_kwargs['encode_mode']
-
-                if waveform_mode not in ("BB", "CW"):
-                    raise ValueError(
-                        f"waveform_mode {waveform_mode} is invalid. EK80 waveform_mode must be 'BB' or 'CW'."  # noqa
-                    )
-                elif encode_mode not in ("complex", "power"):
-                    raise ValueError(
-                        f"encode_mode {encode_mode} is invalid. EK80 waveform_mode must be 'complex' or 'power'."  # noqa
-                    )
-                elif waveform_mode == "BB" and encode_mode == "power":
-                    raise ValueError(
-                        "Data from broadband ('BB') transmission must be recorded as complex samples"  # noqa
-                    )
-
-            # Compute range via calibration objects
-            if data.sonar_model == "AZFP":
-                cal_obj = CalibrateAZFP(
-                    echodata=data,
-                    env_params=range_kwargs.get("env_params", {}),
-                    cal_params=None,
-                    waveform_mode=None,
-                    encode_mode=None,
-                )
-                if range_kwargs["azfp_cal_type"] is None:
-                    raise ValueError("azfp_cal_type must be specified when sonar_model is AZFP")
-                cal_obj.compute_echo_range(cal_type=range_kwargs["azfp_cal_type"])
-            elif data.sonar_model in ("EK60", "EK80", "ES70", "ES80", "EA640"):
-                if data.sonar_model in ["EK60", "ES70"]:
-                    cal_obj = CalibrateEK60(
-                        echodata=data,
-                        env_params=range_kwargs.get("env_params", {}),
-                        cal_params=None,
-                        ecs_file=None,
-                    )
-                else:
-                    cal_obj = CalibrateEK80(
-                        echodata=data,
-                        env_params=range_kwargs.get("env_params", {}),
-                        cal_params=None,
-                        ecs_file=None,
-                        waveform_mode=range_kwargs.get("waveform_mode", "CW"),
-                        encode_mode=range_kwargs.get("encode_mode", "power"),
-                    )
-            range_in_meter = cal_obj.range_meter
-
-            range_in_meter.attrs = range_attrs
-            if vertical_offset is not None:
-                range_in_meter = _add_vertical_offset(
-                    range_in_meter=range_in_meter,
-                    vertical_offset=vertical_offset,
-                    data_type=EchoData,
-                    platform_data=data["Platform"],
-                )
-            ds = ds.assign_coords({'echo_range': range_in_meter})
-            ds.echo_range.attrs = range_attrs
-
-    elif isinstance(data, xr.Dataset):
-        if 'ping_time' in data:
-            _check_ping_time(data.ping_time)
-        variable = 'Sv'
-        ds = data
-        yaxis = 'echo_range'
-        if 'echo_range' not in data.dims and get_range is False:
-            # Range in dims indicates that data is MVBS.
-            yaxis = 'range_sample'
-
-        # If depth is available in ds, use it.
-        ds = ds.set_coords('echo_range')
-        if vertical_offset is not None:
-            ds['echo_range'] = _add_vertical_offset(
-                range_in_meter=ds.echo_range,
-                vertical_offset=vertical_offset,
-                data_type=xr.Dataset,
-            )
-        ds.echo_range.attrs = range_attrs
-    else:
-        raise ValueError(f"Unsupported data type: {type(data)}")
-
-    return _plot_echogram(
-        ds,
-        xaxis='ping_time',
-        yaxis=yaxis,
-        variable=variable,
-        channel=channel,
-        frequency=frequency,
-        **kwargs,
-    )
-
-
-def _check_ping_time(ping_time):
-    if ping_time.shape[0] < 2:
-        raise ValueError("Ping time must have a length that is greater or equal to 2")
-
-
-def _add_vertical_offset(
-    range_in_meter: xr.DataArray,
-    vertical_offset: Union[int, float, xr.DataArray, bool],
-    data_type: Union[Type[xr.Dataset], Type[EchoData]],
-    platform_data: Optional[xr.Dataset] = None,
-) -> xr.DataArray:
-    # Below, we rename time2 to ping_time because range_in_meter is in ping_time
-    if isinstance(vertical_offset, bool):
-        if vertical_offset is True:
-            if data_type == xr.Dataset:
-                logger.warning(
-                    "Boolean type found for vertical_offset. Ignored since data is an xarray dataset."
-                )
-                return range_in_meter
-            elif data_type == EchoData:
-                if (
-                    isinstance(platform_data, xr.Dataset)
-                    and 'vertical_offset' in platform_data
-                ):
-                    return range_in_meter + platform_data.vertical_offset.rename({'time2': 'ping_time'})
-                else:
-                    logger.warning(
-                        "Boolean type found for vertical_offset. Please provide platform data with vertical_offset in it or provide a separate vertical_offset data."  # noqa
-                    )
-                    return range_in_meter
-        logger.warning(f"vertical_offset value of {vertical_offset} is ignored.")
-        return range_in_meter
-    if isinstance(vertical_offset, xr.DataArray):
-        check_dims = list(range_in_meter.dims)
-        check_dims.remove('channel')
-        if 'time2' in vertical_offset:
-            vertical_offset = vertical_offset.rename({'time2': 'ping_time'})
-
-        if not any(
-            True if d in vertical_offset.dims else False for d in check_dims
-        ):
-            raise ValueError(
-                f"vertical_offset must have any of these dimensions: {', '.join(check_dims)}"
-            )
-        # Adds vertical_offset to range if it exists
-        return range_in_meter + vertical_offset
-    elif isinstance(vertical_offset, (int, float)):
-        return range_in_meter + vertical_offset
diff --git a/echopype/visualize/plot.py b/echopype/visualize/plot.py
deleted file mode 100644
index b2a3a9b23..000000000
--- a/echopype/visualize/plot.py
+++ /dev/null
@@ -1,267 +0,0 @@
-import matplotlib.pyplot as plt
-import matplotlib.cm
-import xarray as xr
-import numpy as np
-from xarray.plot.facetgrid import FacetGrid
-from matplotlib.collections import QuadMesh
-from typing import Optional, Union, List
-from .cm import cmap_d
-from ..utils.log import _init_logger
-
-logger = _init_logger(__name__)
-
-
-def _format_axis_label(axis_variable):
-    return axis_variable.replace('_', " ").title()
-
-
-def _set_label(
-    fg: Union[FacetGrid, QuadMesh, None] = None,
-    filter_var: str = 'channel',
-    filter_val: Union[str, int, float, None] = None,
-    col: Optional[str] = None,
-):
-    props = {'boxstyle': 'square', 'facecolor': 'white', 'alpha': 0.7}
-    if isinstance(fg, FacetGrid):
-        text_pos = [0.02, 0.06]
-        fontsize = 14
-        if col == 'beam':
-            if isinstance(filter_val, list) or filter_val is None:
-                for rl in fg.row_labels:
-                    if rl is not None:
-                        rl.set_text('')
-
-                for idx, cl in enumerate(fg.col_labels):
-                    if cl is not None:
-                        cl.set_text(f'Beam {fg.col_names[idx]}')
-
-            text_pos = [0.04, 0.06]
-            fontsize = 10
-
-        for idx, ax in enumerate(fg.axes.flat):
-            name_dicts = fg.name_dicts.flatten()
-            if filter_var in name_dicts[idx]:
-                chan = name_dicts[idx][filter_var]
-                if col == filter_var:
-                    ax.set_title('')
-            else:
-                chan = filter_val
-                ax.set_title(f'Beam {fg.col_names[idx]}')
-            axtext = chan
-            if filter_var == 'frequency':
-                axtext = f"{int(chan / 1000)} kHz"
-            ax.text(
-                *text_pos,
-                axtext,
-                transform=ax.transAxes,
-                fontsize=fontsize,
-                verticalalignment='bottom',
-                bbox=props,
-            )
-    else:
-        if filter_val is None:
-            raise ValueError(
-                f'{filter_var.title()} value is missing for single echogram plotting.'
-            )
-        axtext = filter_val
-        if filter_var == 'frequency':
-            axtext = f"{int(axtext / 1000)} kHz"
-        ax = fg.axes
-        ax.text(
-            0.02,
-            0.04,
-            axtext,
-            transform=ax.transAxes,
-            fontsize=13,
-            verticalalignment='bottom',
-            bbox=props,
-        )
-        plt.title('')
-        plt.tight_layout()
-
-
-def _set_plot_defaults(kwargs):
-    plot_defaults = {
-        'cmap': 'jet',
-        'figsize': (15, 10),
-        'robust': False,
-        'yincrease': False,
-        'col_wrap': 1,
-    }
-
-    # Set plot defaults if not passed in kwargs
-    for k, v in plot_defaults.items():
-        if k not in kwargs:
-            kwargs[k] = v
-        elif k == 'cmap' and k in kwargs:
-            cmap = kwargs[k]
-            try:
-                if cmap in cmap_d:
-                    cmap = f'ep.{cmap}'
-                    kwargs[k] = cmap
-                matplotlib.cm.get_cmap(cmap)
-            except:
-                import cmocean
-
-                if cmap.startswith('cmo'):
-                    _, cmap = cmap.split('.')
-
-                if cmap in cmocean.cm.cmap_d:
-                    kwargs[k] = f'cmo.{cmap}'
-                else:
-                    raise ValueError(f"{cmap} is not a valid colormap.")
-
-    # Remove extra plotting attributes that should be set
-    # by echopype devs
-    exclude_attrs = ['x', 'y', 'col', 'row']
-    for attr in exclude_attrs:
-        if attr in kwargs:
-            logger.warning(f"{attr} in kwargs. Removing.")
-            kwargs.pop(attr)
-
-    return kwargs
-
-
-def _plot_echogram(
-    ds: xr.Dataset,
-    channel: Union[str, List[str], None] = None,
-    frequency: Union[str, List[str], None] = None,
-    variable: str = 'backscatter_r',
-    xaxis: str = 'ping_time',
-    yaxis: str = 'echo_range',
-    **kwargs,
-) -> Union[FacetGrid, QuadMesh]:
-    kwargs = _set_plot_defaults(kwargs)
-
-    row = None
-    col = None
-    filter_var = 'channel'
-    filter_val = None
-    # perform frequency filtering
-    if channel is not None:
-        if 'channel' not in ds.dims:
-            raise ValueError("Channel filtering is not available because channel is not a dimension for your dataset!")
-        ds = ds.sel(channel=channel)
-        filter_val = channel
-    elif frequency is not None:
-        duplicates = False
-        if 'channel' in ds.dims:
-            if len(np.unique(ds.frequency_nominal)) != len(ds.frequency_nominal):
-                duplicates = True
-            ds = ds.where(ds.frequency_nominal.isin(frequency), drop=True)
-        else:
-            if len(np.unique(ds.frequency)) != len(ds.frequency):
-                duplicates = True
-            ds = ds.sel(frequency=frequency)
-
-        if duplicates:
-            raise ValueError("Duplicate frequency found, please use channel for filtering.")
-        filter_val = frequency
-        filter_var = 'frequency'
-
-    if 'backscatter_i' in ds.variables:
-        col = 'beam'
-        kwargs.setdefault('figsize', (15, 5))
-        kwargs.update(
-            {
-                'col_wrap': None,
-            }
-        )
-        filtered_ds = np.abs(ds.backscatter_r + 1j * ds.backscatter_i)
-    else:
-        filtered_ds = ds[variable]
-        if 'beam' in filtered_ds.dims:
-            filtered_ds = filtered_ds.isel(beam=0).drop('beam')
-
-    if 'channel' in filtered_ds.dims and frequency is not None:
-        filtered_ds = filtered_ds.assign_coords({'frequency': ds.frequency_nominal})
-        filtered_ds = filtered_ds.swap_dims({'channel': 'frequency'})
-        if filtered_ds.frequency.size == 1:
-            filtered_ds = filtered_ds.isel(frequency=0)
-
-    # figure out channel/frequency size
-    # to determine plotting method
-    if filtered_ds[filter_var].size > 1:
-        if col is None:
-            col = filter_var
-        else:
-            row = filter_var
-
-    filtered_ds[xaxis].attrs = {
-        'long_name': filtered_ds[xaxis].attrs.get(
-            'long_name', _format_axis_label(xaxis)
-        ),
-        'units': filtered_ds[xaxis].attrs.get('units', ''),
-    }
-    filtered_ds[yaxis].attrs = {
-        'long_name': filtered_ds[yaxis].attrs.get(
-            'long_name', _format_axis_label(yaxis)
-        ),
-        'units': filtered_ds[yaxis].attrs.get('units', ''),
-    }
-
-    plots = []
-    if not filtered_ds[filter_var].shape:
-        if (
-            np.any(filtered_ds.isnull()).values == np.array(True)
-            and 'echo_range' in filtered_ds.coords
-            and 'range_sample' in filtered_ds.dims
-            and variable in ['backscatter_r', 'Sv']
-        ):
-            # Handle the nans for echodata and Sv
-            filtered_ds = filtered_ds.sel(
-                ping_time=filtered_ds.echo_range.dropna(dim='ping_time', how='all').ping_time
-            )
-            filtered_ds = filtered_ds.sel(
-                range_sample=filtered_ds.echo_range.dropna(dim='range_sample').range_sample
-            )
-        plot = filtered_ds.plot.pcolormesh(
-            x=xaxis,
-            y=yaxis,
-            col=col,
-            row=row,
-            **kwargs,
-        )
-        _set_label(plot, filter_var=filter_var, filter_val=filter_val, col=col)
-        plots.append(plot)
-    else:
-        # Scale plots
-        num_chan = len(filtered_ds[filter_var])
-        chan_scaling = (-0.06, -0.16)
-        figsize_scale = tuple(
-            [1 + (scale * num_chan) for scale in chan_scaling]
-        )
-        new_size = tuple(
-            [
-                size * figsize_scale[idx]
-                for idx, size in enumerate(kwargs.get('figsize'))
-            ]
-        )
-        kwargs.update({'figsize': new_size})
-
-        for f in filtered_ds[filter_var]:
-            d = filtered_ds[filtered_ds[filter_var] == f.values]
-            if (
-                np.any(d.isnull()).values == np.array(True)
-                and 'echo_range' in d.coords
-                and 'range_sample' in d.dims
-                and variable in ['backscatter_r', 'Sv']
-            ):
-                # Handle the nans for echodata and Sv
-                d = d.sel(
-                    ping_time=d.echo_range.dropna(dim='ping_time', how='all').ping_time
-                )
-                d = d.sel(
-                    range_sample=d.echo_range.dropna(dim='range_sample').range_sample
-                )
-
-            plot = d.plot.pcolormesh(
-                x=xaxis,
-                y=yaxis,
-                col=col,
-                row=row,
-                **kwargs,
-            )
-            _set_label(plot, filter_var=filter_var, filter_val=filter_val, col=col)
-            plots.append(plot)
-    return plots


Test patch:
diff --git a/echopype/tests/colormap/test_colormap.py b/echopype/tests/colormap/test_colormap.py
new file mode 100644
index 000000000..c991230ac
--- /dev/null
+++ b/echopype/tests/colormap/test_colormap.py
@@ -0,0 +1,18 @@
+import echopype.colormap
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+import numpy as np
+
+import pytest
+
+
+@pytest.mark.parametrize(
+    "cmap_name",
+    [
+        ("ep.ek500"),
+    ],
+    ids=["ek500"],
+)
+def test_colormap(cmap_name):
+    fig = plt.imshow(np.random.rand(10,10), cmap=cmap_name)
+    assert isinstance(fig, mpl.image.AxesImage) is True
diff --git a/echopype/tests/visualize/test_plot.py b/echopype/tests/visualize/test_plot.py
deleted file mode 100644
index 5c5dc6290..000000000
--- a/echopype/tests/visualize/test_plot.py
+++ /dev/null
@@ -1,327 +0,0 @@
-import echopype
-import echopype.visualize
-from echopype.testing import TEST_DATA_FOLDER
-from echopype.calibrate.calibrate_ek import CalibrateEK60
-
-import pytest
-from xarray.plot.facetgrid import FacetGrid
-from matplotlib.collections import QuadMesh
-import xarray as xr
-import numpy as np
-
-ek60_path = TEST_DATA_FOLDER / "ek60"
-ek80_path = TEST_DATA_FOLDER / "ek80_new"
-azfp_path = TEST_DATA_FOLDER / "azfp"
-ad2cp_path = TEST_DATA_FOLDER / "ad2cp"
-
-param_args = ("filepath", "sonar_model", "azfp_xml_path", "range_kwargs")
-param_testdata = [
-    (
-        ek60_path / "ncei-wcsd" / "Summer2017-D20170719-T211347.raw",
-        "EK60",
-        None,
-        {},
-    ),
-    (
-        ek60_path / "DY1002_EK60-D20100318-T023008_rep_freq.raw",
-        "EK60",
-        None,
-        {},
-    ),
-    (
-        ek80_path / "echopype-test-D20211004-T235930.raw",
-        "EK80",
-        None,
-        {'waveform_mode': 'BB', 'encode_mode': 'complex'},
-    ),
-    (
-        ek80_path / "D20211004-T233354.raw",
-        "EK80",
-        None,
-        {'waveform_mode': 'CW', 'encode_mode': 'power'},
-    ),
-    (
-        ek80_path / "D20211004-T233115.raw",
-        "EK80",
-        None,
-        {'waveform_mode': 'CW', 'encode_mode': 'complex'},
-    ),
-    (
-        azfp_path / "17082117.01A",
-        "AZFP",
-        azfp_path / "17041823.XML",
-        {},
-    ),  # Will always need env variables
-    pytest.param(
-        ad2cp_path / "raw" / "090" / "rawtest.090.00001.ad2cp",
-        "AD2CP",
-        None,
-        {},
-        marks=pytest.mark.xfail(
-            run=False,
-            reason="Not supported at the moment",
-        ),
-    ),
-]
-
-
-@pytest.mark.parametrize(param_args, param_testdata)
-def test_plot_multi(
-    filepath,
-    sonar_model,
-    azfp_xml_path,
-    range_kwargs,
-):
-    # TODO: Need to figure out how to compare the actual rendered plots
-    ed = echopype.open_raw(filepath, sonar_model, azfp_xml_path)
-    plots = echopype.visualize.create_echogram(ed)
-    assert isinstance(plots, list) is True
-    assert all(isinstance(plot, FacetGrid) for plot in plots) is True
-
-
-@pytest.mark.parametrize(param_args, param_testdata)
-def test_plot_single(
-    filepath,
-    sonar_model,
-    azfp_xml_path,
-    range_kwargs,
-):
-    # TODO: Need to figure out how to compare the actual rendered plots
-    ed = echopype.open_raw(filepath, sonar_model, azfp_xml_path)
-    plots = echopype.visualize.create_echogram(
-        ed, channel=ed["Sonar/Beam_group1"].channel[0].values
-    )
-    assert isinstance(plots, list) is True
-    if (
-        sonar_model.lower() == 'ek80'
-        and range_kwargs['encode_mode'] == 'complex'
-    ):
-        assert all(isinstance(plot, FacetGrid) for plot in plots) is True
-    else:
-        assert all(isinstance(plot, QuadMesh) for plot in plots) is True
-
-
-@pytest.mark.parametrize(param_args, param_testdata)
-def test_plot_multi_get_range(
-    filepath,
-    sonar_model,
-    azfp_xml_path,
-    range_kwargs,
-):
-    # TODO: Need to figure out how to compare the actual rendered plots
-    ed = echopype.open_raw(filepath, sonar_model, azfp_xml_path)
-    if ed.sonar_model.lower() == 'azfp':
-        avg_temperature = ed["Environment"]['temperature'].values.mean()
-        env_params = {
-            'temperature': avg_temperature,
-            'salinity': 27.9,
-            'pressure': 59,
-        }
-        range_kwargs['env_params'] = env_params
-    plots = echopype.visualize.create_echogram(
-        ed, get_range=True, range_kwargs=range_kwargs
-    )
-    assert isinstance(plots, list) is True
-    assert all(isinstance(plot, FacetGrid) for plot in plots) is True
-
-    # Beam shape check
-    if (
-        sonar_model.lower() == 'ek80'
-        and range_kwargs['encode_mode'] == 'complex'
-    ):
-        assert plots[0].axes.shape[-1] > 1
-    else:
-        assert plots[0].axes.shape[-1] == 1
-
-    # Channel shape check
-    assert ed["Sonar/Beam_group1"].channel.shape[0] == len(plots)
-
-
-@pytest.mark.parametrize(param_args, param_testdata)
-def test_plot_Sv(
-    filepath,
-    sonar_model,
-    azfp_xml_path,
-    range_kwargs,
-):
-    # TODO: Need to figure out how to compare the actual rendered plots
-    ed = echopype.open_raw(filepath, sonar_model, azfp_xml_path)
-    if ed.sonar_model.lower() == 'azfp':
-        avg_temperature = ed["Environment"]['temperature'].values.mean()
-        env_params = {
-            'temperature': avg_temperature,
-            'salinity': 27.9,
-            'pressure': 59,
-        }
-        range_kwargs['env_params'] = env_params
-        if 'azfp_cal_type' in range_kwargs:
-            range_kwargs.pop('azfp_cal_type')
-    Sv = echopype.calibrate.compute_Sv(ed, **range_kwargs)
-    plots = echopype.visualize.create_echogram(Sv)
-    assert isinstance(plots, list) is True
-    assert all(isinstance(plot, FacetGrid) for plot in plots) is True
-
-
-@pytest.mark.parametrize(param_args, param_testdata)
-def test_plot_mvbs(
-    filepath,
-    sonar_model,
-    azfp_xml_path,
-    range_kwargs,
-):
-
-    # TODO: Need to figure out how to compare the actual rendered plots
-    ed = echopype.open_raw(filepath, sonar_model, azfp_xml_path)
-    if ed.sonar_model.lower() == 'azfp':
-        avg_temperature = ed["Environment"]['temperature'].values.mean()
-        env_params = {
-            'temperature': avg_temperature,
-            'salinity': 27.9,
-            'pressure': 59,
-        }
-        range_kwargs['env_params'] = env_params
-        if 'azfp_cal_type' in range_kwargs:
-            range_kwargs.pop('azfp_cal_type')
-    Sv = echopype.calibrate.compute_Sv(ed, **range_kwargs)
-    mvbs = echopype.commongrid.compute_MVBS(Sv, ping_time_bin='10S')
-
-    plots = []
-    try:
-        plots = echopype.visualize.create_echogram(mvbs)
-    except Exception as e:
-        assert isinstance(e, ValueError)
-        assert str(e) == "Ping time must have a length that is greater or equal to 2"  # noqa
-
-    if len(plots) > 0:
-        assert all(isinstance(plot, FacetGrid) for plot in plots) is True
-
-
-@pytest.mark.parametrize(
-    ("vertical_offset", "expect_warning"),
-    [
-        (True, False),
-        ([True], True),
-        (False, True),
-        (xr.DataArray(np.array(50.0)), False),
-        ([10, 30.5], False),
-        (10, False),
-        (30.5, False),
-    ],
-)
-def test_vertical_offset_echodata(vertical_offset, expect_warning, caplog):
-    from echopype.echodata import EchoData
-    from echopype.visualize.api import _add_vertical_offset
-    echopype.verbose()
-
-    filepath = ek60_path / "ncei-wcsd" / "Summer2017-D20170719-T211347.raw"
-    sonar_model = "EK60"
-    range_kwargs = {}
-
-    echodata = echopype.open_raw(
-        sonar_model=sonar_model, raw_file=filepath, xml_path=None
-    )
-
-    cal_obj = CalibrateEK60(
-        echodata=echodata,
-        env_params=range_kwargs.get("env_params", {}),
-        cal_params=None,
-        ecs_file=None,
-    )
-    range_in_meter = cal_obj.range_meter
-
-    single_array = range_in_meter.sel(channel='GPT  18 kHz 009072058c8d 1-1 ES18-11').isel(ping_time=0).values
-    no_input_vertical_offset = False
-
-    if isinstance(vertical_offset, list):
-        vertical_offset = vertical_offset[0]
-        echodata["Platform"] = echodata["Platform"].drop_vars('vertical_offset')
-        no_input_vertical_offset = True
-
-    if isinstance(vertical_offset, xr.DataArray):
-        original_array = single_array + vertical_offset.values
-    elif isinstance(vertical_offset, bool) and vertical_offset is True:
-        if not no_input_vertical_offset:
-            original_array = (
-                single_array
-                + echodata["Platform"].vertical_offset.isel(time2=0).values
-            )
-        else:
-            original_array = single_array
-    elif vertical_offset is not False and isinstance(vertical_offset, (int, float)):
-        original_array = single_array + vertical_offset
-    else:
-        original_array = single_array
-
-    results = None
-    try:
-        results = _add_vertical_offset(
-            range_in_meter=range_in_meter,
-            vertical_offset=vertical_offset,
-            data_type=EchoData,
-            platform_data=echodata["Platform"],
-        )
-        if expect_warning:
-            assert 'WARNING' in caplog.text
-    except Exception as e:
-        assert isinstance(e, ValueError)
-        assert str(e) == 'vertical_offset must have any of these dimensions: ping_time, range_sample'  # noqa
-
-    if isinstance(results, xr.DataArray):
-        final_array = results.sel(channel='GPT  18 kHz 009072058c8d 1-1 ES18-11').isel(ping_time=0).values
-        print(f"original_array = {original_array}")
-        print(f"results = {results}")
-        assert np.array_equal(original_array, final_array)
-
-
-@pytest.mark.parametrize(
-    ("vertical_offset", "expect_warning"),
-    [
-        (True, True),
-        (False, True),
-        (xr.DataArray(np.array(50.0)), False),
-        (10, False),
-        (30.5, False),
-    ],
-)
-def test_vertical_offset_Sv_dataset(vertical_offset, expect_warning, caplog):
-    from echopype.visualize.api import _add_vertical_offset
-    echopype.verbose()
-
-    filepath = ek60_path / "ncei-wcsd" / "Summer2017-D20170719-T211347.raw"
-    sonar_model = "EK60"
-    range_kwargs = {}
-
-    echodata = echopype.open_raw(
-        sonar_model=sonar_model, raw_file=filepath, xml_path=None
-    )
-    Sv = echopype.calibrate.compute_Sv(echodata, **range_kwargs)
-    ds = Sv.set_coords('echo_range')
-    range_in_meter = ds.echo_range
-
-    single_array = range_in_meter.sel(channel='GPT  18 kHz 009072058c8d 1-1 ES18-11').isel(ping_time=0).values
-
-    if isinstance(vertical_offset, xr.DataArray):
-        original_array = single_array + vertical_offset.values
-    elif not isinstance(vertical_offset, bool) and isinstance(vertical_offset, (int, float)):
-        original_array = single_array + vertical_offset
-    else:
-        original_array = single_array
-
-    results = None
-    try:
-        results = _add_vertical_offset(
-            range_in_meter=range_in_meter,
-            vertical_offset=vertical_offset,
-            data_type=xr.Dataset,
-        )
-
-        if expect_warning:
-            assert 'WARNING' in caplog.text
-    except Exception as e:
-        assert isinstance(e, ValueError)
-        assert str(e) == 'vertical_offset must have any of these dimensions: ping_time, range_sample'  # noqa
-
-    if isinstance(results, xr.DataArray):
-        final_array = results.sel(channel='GPT  18 kHz 009072058c8d 1-1 ES18-11').isel(ping_time=0).values
-
-        assert np.array_equal(original_array, final_array)