diff --git a/pyproject.toml b/pyproject.toml
index 9ff2789c..0d87849f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -106,6 +106,10 @@ line_length = 110
 exclude = [
     "__init__.py",
 ]
+line-length = 110
+target-version = "py311"
+
+[tool.ruff.lint]
 ignore = [
     "N802",
     "N803",
@@ -124,21 +128,19 @@ ignore = [
     "D400",
     "E712",
 ]
-line-length = 110
 select = [
     "E",  # pycodestyle
     "F",  # pyflakes
     "N",  # pep8-naming
     "W",  # pycodestyle
 ]
-target-version = "py311"
 extend-select = [
     "RUF100", # Warn about unused noqa
 ]
 
-[tool.ruff.pycodestyle]
+[tool.ruff.lint.pycodestyle]
 max-doc-length = 79
 
-[tool.ruff.pydocstyle]
+[tool.ruff.lint.pydocstyle]
 convention = "numpy"
 
diff --git a/rubin_sim/maf/batches/moving_objects_batch.py b/rubin_sim/maf/batches/moving_objects_batch.py
index 5769b16f..9ff4b6fc 100644
--- a/rubin_sim/maf/batches/moving_objects_batch.py
+++ b/rubin_sim/maf/batches/moving_objects_batch.py
@@ -835,7 +835,7 @@ def _codePlot(key):
     display_dict["subgroup"] = "Completeness over time"
     display_dict["caption"] = "Completeness over time, for H values indicated in legend."
     ph.save_fig(
-        fig.number,
+        fig,
         f"{figroot}_CompletenessOverTime",
         "Combo",
         "CompletenessOverTime",
diff --git a/rubin_sim/maf/metric_bundles/metric_bundle.py b/rubin_sim/maf/metric_bundles/metric_bundle.py
index 860e0a2e..f87abe5b 100644
--- a/rubin_sim/maf/metric_bundles/metric_bundle.py
+++ b/rubin_sim/maf/metric_bundles/metric_bundle.py
@@ -790,9 +790,9 @@ def plot(self, plot_handler=None, plot_func=None, outfile_suffix=None, savefig=F
         Returns
         -------
         made_plots : `dict`
-            Dictionary of plot_type:figure number key/value pairs,
+            Dictionary of plot_type:figure key/value pairs,
             indicating what plots were created
-            and what matplotlib figure numbers were used.
+            and what matplotlib figures were used.
         """
         # Generate a plot_handler if none was set.
         if plot_handler is None:
@@ -808,10 +808,10 @@ def plot(self, plot_handler=None, plot_func=None, outfile_suffix=None, savefig=F
         plot_handler.set_plot_dicts(plot_dicts=[self.plot_dict], reset=True)
         made_plots = {}
         if plot_func is not None:
-            fignum = plot_handler.plot(plot_func, outfile_suffix=outfile_suffix)
-            made_plots[plot_func.plotType] = fignum
+            fig = plot_handler.plot(plot_func, outfile_suffix=outfile_suffix)
+            made_plots[plot_func.plotType] = fig
         else:
             for plot_func in self.plot_funcs:
-                fignum = plot_handler.plot(plot_func, outfile_suffix=outfile_suffix)
-                made_plots[plot_func.plot_type] = fignum
+                fig = plot_handler.plot(plot_func, outfile_suffix=outfile_suffix)
+                made_plots[plot_func.plot_type] = fig
         return made_plots
diff --git a/rubin_sim/maf/plots/__init__.py b/rubin_sim/maf/plots/__init__.py
index 7394531e..42dee1af 100644
--- a/rubin_sim/maf/plots/__init__.py
+++ b/rubin_sim/maf/plots/__init__.py
@@ -6,7 +6,6 @@
 from .night_pointing_plotter import *
 from .oned_plotters import *
 from .perceptual_rainbow import *
-from .plot_bundle import *
 from .plot_handler import *
 from .spatial_plotters import *
 from .special_plotters import *
diff --git a/rubin_sim/maf/plots/hg_plotters.py b/rubin_sim/maf/plots/hg_plotters.py
index cb847740..3df11d15 100644
--- a/rubin_sim/maf/plots/hg_plotters.py
+++ b/rubin_sim/maf/plots/hg_plotters.py
@@ -305,21 +305,25 @@ def _map_colors(self, values):  # pylint: disable=invalid-name, no-self-use
 
         return colors, color_mappable
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
         """Restructure the metric data to use, and build the figure.
 
         Parameters
         ----------
-        metric_value : `numpy.ndarray`
-           Metric values
-        slicer : `rubin_sim.maf.slicers.baseSlicer.BaseSlicer`
-           must have "mjd" and "duration" slice points, in units
-           of days and seconds, respectively.
-        user_plot_dict : `dict`
-           Plotting parameters
-        fignum : `int`
-           matplotlib figure number
+        metric_value : `numpy.ma.MaskedArray`
+            The metric values from the bundle.
+        slicer : `rubin_sim.maf.slicers.TwoDSlicer`
+            The slicer.
+        user_plot_dict: `dict`
+            Dictionary of plot parameters set by user
+            (overrides default values).
+        fig : `matplotlib.figure.Figure`
+            Matplotlib figure number to use. Default = None, starts new figure.
 
+        Returns
+        -------
+        fig : `matplotlib.figure.Figure`
+           Figure with the plot.
         """
         # Highest level driver for the plotter.
         # Prepares data structures and figure-wide elements
@@ -344,7 +348,8 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         self.plot_dict.update(user_plot_dict)
 
         # Generate the figure
-        fig = plt.figure(fignum, figsize=self.plot_dict["figsize"])
+        if fig is None:
+            fig = plt.figure(figsize=self.plot_dict["figsize"])
 
         # Add the plots
         color_mappable, axes = self._plot(fig, intervals)
@@ -355,7 +360,9 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
 
         # add the legend, if requested
         if self.plot_dict["legend"] and len(self.color_map) > 0:
-            self._add_figure_legend(fig, axes)
+            fig = self._add_figure_legend(fig, axes)
+
+        return fig
 
     def _add_figure_legend(self, fig, axes):
         """Creates and adds the figure legend.
@@ -370,8 +377,8 @@ def _add_figure_legend(self, fig, axes):
 
         Returns
         -------
-        figure_number : `int`
-           The matplotlib figure number
+        fig : `matplotlib.figure.Figure`
+           The matplotlib figure
         """
         # Creates and adds the figure legend.  This method follows two
         # stages: first, it adds entries for each element in the color
@@ -423,7 +430,7 @@ def _add_figure_legend(self, fig, axes):
             bbox_to_anchor=self.plot_dict["legend_bbox_to_anchor"],
         )
 
-        return fig.number
+        return fig
 
     def _add_colorbar(self, fig, color_mappable, axes):  # pylint: disable=invalid-name, no-self-use
         """Add a colorbar.
diff --git a/rubin_sim/maf/plots/hourglass_plotters.py b/rubin_sim/maf/plots/hourglass_plotters.py
index 50fbb81a..81b4bdbd 100644
--- a/rubin_sim/maf/plots/hourglass_plotters.py
+++ b/rubin_sim/maf/plots/hourglass_plotters.py
@@ -14,6 +14,7 @@ def __init__(self):
             "title": None,
             "xlabel": "Night - min(Night)",
             "ylabel": "Hours from local midnight",
+            "figsize": None,
         }
         self.filter2color = {
             "u": "purple",
@@ -24,20 +25,20 @@ def __init__(self):
             "y": "red",
         }
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
         """
         Generate the hourglass plot
         """
         if slicer.slicer_name != "HourglassSlicer":
             raise ValueError("HourglassPlot is for use with hourglass slicers")
 
-        fig = plt.figure(fignum)
-        ax = fig.add_subplot(111)
-
         plot_dict = {}
         plot_dict.update(self.default_plot_dict)
         plot_dict.update(user_plot_dict)
 
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
+        ax = fig.add_subplot(111)
         pernight = metric_value[0]["pernight"]
         perfilter = metric_value[0]["perfilter"]
 
@@ -57,9 +58,11 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
                 color=self.filter2color[key],
                 transform=ax.transAxes,
             )
-        # ax.plot(pernight['mjd'] - dmin, (pernight['twi6_rise'] - pernight['midnight']) * 24.,
+        # ax.plot(pernight['mjd'] - dmin,
+        # (pernight['twi6_rise'] - pernight['midnight']) * 24.,
         #        'blue', label=r'6$^\circ$ twilight')
-        # ax.plot(pernight['mjd'] - dmin, (pernight['twi6_set'] - pernight['midnight']) * 24.,
+        # ax.plot(pernight['mjd'] - dmin,
+        # (pernight['twi6_set'] - pernight['midnight']) * 24.,
         #        'blue')
         ax.plot(
             pernight["mjd"] - dmin,
@@ -114,4 +117,4 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
             )
             ax.axhline((pernight["twi18_set"] - pernight["midnight"]) * 24.0, color="red")
 
-        return fig.number
+        return fig
diff --git a/rubin_sim/maf/plots/mo_plotters.py b/rubin_sim/maf/plots/mo_plotters.py
index 519f6768..0e8b99b8 100644
--- a/rubin_sim/maf/plots/mo_plotters.py
+++ b/rubin_sim/maf/plots/mo_plotters.py
@@ -7,9 +7,12 @@
 
 from .plot_handler import BasePlotter
 
-# mag_sun = -27.1 # apparent r band magnitude of the sun. this sets the band for the magnitude limit.
-# see http://www.ucolick.org/~cnaw/sun.html for apparent magnitudes in other bands.
-mag_sun = -26.74  # apparent V band magnitude of the Sun (our H mags translate to V band)
+# mag_sun = -27.1
+# apparent r band magnitude of the sun.
+# this sets the band for the magnitude limit.
+# see http://www.ucolick.org/~cnaw/sun.html for apparent mags in other bands.
+mag_sun = -26.74
+# apparent V band magnitude of the Sun (our H mags translate to V band)
 km_per_au = 1.496e8
 m_per_km = 1000
 
@@ -37,7 +40,7 @@ def __init__(self):
         }
         self.min_hrange = 1.0
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
         if "linestyle" not in user_plot_dict:
             user_plot_dict["linestyle"] = "-"
         plot_dict = {}
@@ -49,11 +52,13 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
             reduce_func = np.mean
         if hvals.shape[0] == 1:
             # We have a simple set of values to plot against H.
-            # This may be due to running a summary metric, such as completeness.
+            # This may be due to running a summary metric,
+            # such as completeness.
             m_vals = metric_value[0].filled()
         elif len(hvals) == slicer.shape[1]:
             # Using cloned H distribution.
-            # Apply 'np_reduce' method directly to metric values, and plot at matching H values.
+            # Apply 'np_reduce' method directly to metric values,
+            # and plot at matching H values.
             m_vals = reduce_func(metric_value.filled(), axis=0)
         else:
             # Probably each object has its own H value.
@@ -67,7 +72,8 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
                 nbins = 30
             stepsize = hrange / float(nbins)
             bins = np.arange(min_h, min_h + hrange + stepsize / 2.0, stepsize)
-            # In each bin of H, calculate the 'np_reduce' value of the corresponding metric_values.
+            # In each bin of H, calculate the 'np_reduce' value of the
+            # corresponding metric_values.
             inds = np.digitize(hvals, bins)
             inds = inds - 1
             m_vals = np.zeros(len(bins), float)
@@ -79,7 +85,8 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
                     m_vals[i] = reduce_func(match.filled())
             hvals = bins
         # Plot the values.
-        fig = plt.figure(fignum, figsize=plot_dict["figsize"])
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
         ax = plt.gca()
         ax.plot(
             hvals,
@@ -130,7 +137,7 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         plt.xlabel(plot_dict["xlabel"])
         plt.ylabel(plot_dict["ylabel"])
         plt.tight_layout()
-        return fig.number
+        return fig
 
 
 class MetricVsOrbit(BasePlotter):
@@ -159,11 +166,12 @@ def __init__(self, xaxis="q", yaxis="e"):
             "figsize": None,
         }
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
         plot_dict = {}
         plot_dict.update(self.default_plot_dict)
         plot_dict.update(user_plot_dict)
-        fig = plt.figure(fignum, figsize=plot_dict["figsize"])
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
         xvals = slicer.slice_points["orbits"][plot_dict["xaxis"]]
         yvals = slicer.slice_points["orbits"][plot_dict["yaxis"]]
         # Set x/y bins.
@@ -249,13 +257,13 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         plt.title(plot_dict["title"])
         plt.xlabel(plot_dict["xlabel"])
         plt.ylabel(plot_dict["ylabel"])
-        return fig.number
+        return fig
 
 
 class MetricVsOrbitPoints(BasePlotter):
     """
-    Plot metric values (at a particular H value) as function of orbital parameters,
-    using points for each metric value.
+    Plot metric values (at a particular H value) as function
+    of orbital parameters, using points for each metric value.
     """
 
     def __init__(self, xaxis="q", yaxis="e"):
@@ -276,11 +284,12 @@ def __init__(self, xaxis="q", yaxis="e"):
             "figsize": None,
         }
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
         plot_dict = {}
         plot_dict.update(self.default_plot_dict)
         plot_dict.update(user_plot_dict)
-        fig = plt.figure(fignum, figsize=plot_dict["figsize"])
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
         xvals = slicer.slice_points["orbits"][plot_dict["xaxis"]]
         yvals = slicer.slice_points["orbits"][plot_dict["yaxis"]]
         # Identify the relevant metric_values for the Hvalue we want to plot.
@@ -346,4 +355,4 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         plt.title(plot_dict["title"])
         plt.xlabel(plot_dict["xlabel"])
         plt.ylabel(plot_dict["ylabel"])
-        return fig.number
+        return fig
diff --git a/rubin_sim/maf/plots/nd_plotters.py b/rubin_sim/maf/plots/nd_plotters.py
index 0b19d057..885e2434 100644
--- a/rubin_sim/maf/plots/nd_plotters.py
+++ b/rubin_sim/maf/plots/nd_plotters.py
@@ -30,28 +30,30 @@ def __init__(self):
             "clims": None,
             "cmap": perceptual_rainbow,
             "cbar_format": None,
+            "figsize": None,
         }
 
-    def __call__(self, metric_values, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_values, slicer, user_plot_dict, fig=None):
         """
         Parameters
         ----------
-        metricValue : numpy.ma.MaskedArray
-        slicer : rubin_sim.maf.slicers.NDSlicer
-        user_plot_dict: dict
-            Dictionary of plot parameters set by user (overrides default values).
-            'xaxis' and 'yaxis' values define which axes of the nd data to plot along the x/y axes.
-        fignum : int
-            Matplotlib figure number to use (default = None, starts new figure).
+        metric_value : `numpy.ma.MaskedArray`
+            The metric values from the bundle.
+        slicer : `rubin_sim.maf.slicers.TwoDSlicer`
+            The slicer.
+        user_plot_dict: `dict`
+            Dictionary of plot parameters set by user
+            (overrides default values).
+        fig : `matplotlib.figure.Figure`
+            Matplotlib figure number to use. Default = None, starts new figure.
 
         Returns
         -------
-        int
-           Matplotlib figure number used to create the plot.
+        fig : `matplotlib.figure.Figure`
+           Figure with the plot.
         """
         if slicer.slicer_name != "NDSlicer":
             raise ValueError("TwoDSubsetData plots ndSlicer metric values")
-        fig = plt.figure(fignum)
         plot_dict = {}
         plot_dict.update(self.default_plot_dict)
         plot_dict.update(user_plot_dict)
@@ -59,6 +61,9 @@ def __call__(self, metric_values, slicer, user_plot_dict, fignum=None):
             raise ValueError("xaxis and yaxis must be specified in plot_dict")
         xaxis = plot_dict["xaxis"]
         yaxis = plot_dict["yaxis"]
+
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
         # Reshape the metric data so we can isolate the values to plot
         # (just new view of data, not copy).
         newshape = []
@@ -66,7 +71,8 @@ def __call__(self, metric_values, slicer, user_plot_dict, fignum=None):
             newshape.append(len(b) - 1)
         newshape.reverse()
         md = metric_values.reshape(newshape)
-        # Sum over other dimensions. Note that masked values are not included in sum.
+        # Sum over other dimensions.
+        # Note that masked values are not included in sum.
         sumaxes = list(range(slicer.nD))
         sumaxes.remove(xaxis)
         sumaxes.remove(yaxis)
@@ -110,7 +116,7 @@ def __call__(self, metric_values, slicer, user_plot_dict, fignum=None):
         )
         cb.set_label(plot_dict["units"])
         plt.title(plot_dict["title"])
-        return fig.number
+        return fig
 
 
 class OneDSubsetData(BasePlotter):
@@ -134,33 +140,39 @@ def __init__(self):
             "alpha": 0.5,
             "cmap": perceptual_rainbow,
             "cbar_format": None,
+            "figsize": None,
         }
 
-    def plot_binned_data1_d(self, metric_values, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_values, slicer, user_plot_dict, fig=None):
         """
         Parameters
         ----------
-        metricValue : numpy.ma.MaskedArray
-        slicer : rubin_sim.maf.slicers.NDSlicer
-        user_plot_dict: dict
-            Dictionary of plot parameters set by user (overrides default values).
-            'axis' keyword identifies which axis to show in the plot (along xaxis of plot).
-        fignum : int
-            Matplotlib figure number to use (default = None, starts new figure).
+        metric_value : `numpy.ma.MaskedArray`
+            The metric values from the bundle.
+        slicer : `rubin_sim.maf.slicers.TwoDSlicer`
+            The slicer.
+        user_plot_dict: `dict`
+            Dictionary of plot parameters set by user
+            (overrides default values).
+        fig : `matplotlib.figure.Figure`
+            Matplotlib figure number to use. Default = None, starts new figure.
 
         Returns
         -------
-        int
-           Matplotlib figure number used to create the plot.
+        fig : `matplotlib.figure.Figure`
+           Figure with the plot.
         """
         if slicer.slicer_name != "NDSlicer":
             raise ValueError("TwoDSubsetData plots ndSlicer metric values")
-        fig = plt.figure(fignum)
+
         plot_dict = {}
         plot_dict.update(self.default_plot_dict)
         plot_dict.update(user_plot_dict)
         if "axis" not in plot_dict:
             raise ValueError("axis for 1-d plot must be specified in plot_dict")
+
+        if fig is None:
+            fig = plt.Figure(figsize=plot_dict["figsize"])
         # Reshape the metric data so we can isolate the values to plot
         # (just new view of data, not copy).
         axis = plot_dict["axis"]
@@ -169,7 +181,8 @@ def plot_binned_data1_d(self, metric_values, slicer, user_plot_dict, fignum=None
             newshape.append(len(b) - 1)
         newshape.reverse()
         md = metric_values.reshape(newshape)
-        # Sum over other dimensions. Note that masked values are not included in sum.
+        # Sum over other dimensions.
+        # Note that masked values are not included in sum.
         sumaxes = list(range(slicer.nD))
         sumaxes.remove(axis)
         sumaxes = tuple(sumaxes)
@@ -204,4 +217,4 @@ def plot_binned_data1_d(self, metric_values, slicer, user_plot_dict, fignum=None
         if plot_dict["histRange"] is not None:
             plt.xlim(plot_dict["histRange"])
         plt.title(plot_dict["title"])
-        return fig.number
+        return fig
diff --git a/rubin_sim/maf/plots/neo_distance_plotter.py b/rubin_sim/maf/plots/neo_distance_plotter.py
index edcba62b..7fba2e26 100644
--- a/rubin_sim/maf/plots/neo_distance_plotter.py
+++ b/rubin_sim/maf/plots/neo_distance_plotter.py
@@ -11,17 +11,21 @@
 
 class NeoDistancePlotter(BasePlotter):
     """
-    Special plotter to calculate and plot the maximum distance an H=22 NEO could be observable to,
-    in any particular particular opsim observation.
+    Special plotter to calculate and plot the maximum distance an H=22 NEO
+    could be observable to, in any particular opsim observation.
+
+    Parameters
+    ----------
+    step : `float`, optional
+        Step size to use for radial bins. Default is 0.01 AU.
+    eclip_max, eclip_min : `float`, `float`, optional
+        Range of ecliptic latitude values to include when creating the plot.
     """
 
     def __init__(self, step=0.01, eclip_max=10.0, eclip_min=-10.0):
-        """
-        eclip_min/Max:  only plot observations within X degrees of the ecliptic plane
-        step: Step size to use for radial bins. Default is 0.01 AU.
-        """
         self.plot_type = "neoxyPlotter"
         self.object_plotter = True
+
         self.default_plot_dict = {
             "title": None,
             "xlabel": "X (AU)",
@@ -31,6 +35,7 @@ def __init__(self, step=0.01, eclip_max=10.0, eclip_min=-10.0):
             "y_min": -0.25,
             "y_max": 2.5,
             "units": "Count",
+            "figsize": None,
         }
         self.filter2color = {
             "u": "purple",
@@ -45,26 +50,25 @@ def __init__(self, step=0.01, eclip_max=10.0, eclip_min=-10.0):
         self.eclip_max = np.radians(eclip_max)
         self.eclip_min = np.radians(eclip_min)
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
         """
         Parameters
         ----------
-        metric_value : numpy.ma.MaskedArray
-            Metric values calculated by rubin_sim.maf.metrics.PassMetric
-        slicer : rubin_sim.maf.slicers.UniSlicer
-        user_plot_dict: dict
-            Dictionary of plot parameters set by user (overrides default values).
-        fignum : int
-            Matplotlib figure number to use (default = None, starts new figure).
+        metric_value : `numpy.ma.MaskedArray`
+            The metric values from the bundle.
+        slicer : `rubin_sim.maf.slicers.TwoDSlicer`
+            The slicer.
+        user_plot_dict: `dict`
+            Dictionary of plot parameters set by user
+            (overrides default values).
+        fig : `matplotlib.figure.Figure`
+            Matplotlib figure number to use. Default = None, starts new figure.
 
         Returns
         -------
-        int
-           Matplotlib figure number used to create the plot.
+        fig : `matplotlib.figure.Figure`
+           Figure with the plot.
         """
-        fig = plt.figure(fignum)
-        ax = fig.add_subplot(111)
-
         in_plane = np.where(
             (metric_value[0]["eclipLat"] >= self.eclip_min) & (metric_value[0]["eclipLat"] <= self.eclip_max)
         )
@@ -73,6 +77,10 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         plot_dict.update(self.default_plot_dict)
         plot_dict.update(user_plot_dict)
 
+        if fig is None:
+            fig = plt.figure(plot_dict["figsize"])
+        ax = fig.add_subplot(111)
+
         planet_props = {"Earth": 1.0, "Venus": 0.72, "Mars": 1.52, "Mercury": 0.39}
 
         planets = []
@@ -82,7 +90,8 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         for planet in planets:
             ax.add_artist(planet)
 
-        # Let's make a 2-d histogram in polar coords, then convert and display in cartisian
+        # Let's make a 2-d histogram in polar coords,
+        # then convert and display in cartisian
 
         r_step = self.step
         rvec = np.arange(0, plot_dict["x_max"] + r_step, r_step)
@@ -105,7 +114,8 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
             theta = np.arctan2(y - 1.0, x)
             diff = np.abs(thetavec - theta)
             theta_to_use = thetavec[np.where(diff == diff.min())]
-            # This is a slow where-clause, should be possible to speed it up using
+            # This is a slow where-clause, should be possible to speed it up
+            # using
             # np.searchsorted+clever slicing or hist2d to build up the map.
             good = np.where((thetagrid == theta_to_use) & (rgrid <= dist))
             H[good] += 1
@@ -126,4 +136,4 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         ax.plot([0], [1], marker="o", color="b")
         ax.plot([0], [0], marker="o", color="y")
 
-        return fig.number
+        return fig
diff --git a/rubin_sim/maf/plots/night_pointing_plotter.py b/rubin_sim/maf/plots/night_pointing_plotter.py
index c35b0f56..cac33728 100644
--- a/rubin_sim/maf/plots/night_pointing_plotter.py
+++ b/rubin_sim/maf/plots/night_pointing_plotter.py
@@ -1,5 +1,7 @@
 __all__ = ("NightPointingPlotter",)
 
+import warnings
+
 import matplotlib.pyplot as plt
 import numpy as np
 
@@ -18,6 +20,7 @@ def __init__(self, mjd_col="observationStartMJD", alt_col="alt", az_col="az"):
             "title": None,
             "xlabel": "MJD",
             "ylabels": ["Alt", "Az"],
+            "figsize": None,
         }
         self.filter2color = {
             "u": "purple",
@@ -28,8 +31,11 @@ def __init__(self, mjd_col="observationStartMJD", alt_col="alt", az_col="az"):
             "y": "red",
         }
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
+        if fig is not None:
+            warnings.warn("Always expect to create a new plot in " "NightPointingPlotter.")
         fig, (ax1, ax2) = plt.subplots(2, sharex=True)
+
         mv = metric_value[0]
 
         u_filters = np.unique(mv["data_slice"]["filter"])
@@ -87,4 +93,4 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
                 transform=ax1.transAxes,
             )
 
-        return fig.number
+        return fig
diff --git a/rubin_sim/maf/plots/oned_plotters.py b/rubin_sim/maf/plots/oned_plotters.py
index a3b34125..ff7cf004 100644
--- a/rubin_sim/maf/plots/oned_plotters.py
+++ b/rubin_sim/maf/plots/oned_plotters.py
@@ -10,6 +10,8 @@
 
 
 class OneDBinnedData(BasePlotter):
+    """Plot data from a OneDSlicer."""
+
     def __init__(self):
         self.plot_type = "BinnedData"
         self.object_plotter = False
@@ -33,7 +35,7 @@ def __init__(self):
             "grid": True,
         }
 
-    def __call__(self, metric_values, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_values, slicer, user_plot_dict, fig=None):
         """
         Plot a set of oneD binned metric data.
         """
@@ -46,7 +48,8 @@ def __call__(self, metric_values, slicer, user_plot_dict, fignum=None):
         plot_dict = {}
         plot_dict.update(self.default_plot_dict)
         plot_dict.update(user_plot_dict)
-        fig = plt.figure(fignum, figsize=plot_dict["figsize"])
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
         # Plot the histogrammed data.
         leftedge = slicer.slice_points["bins"][:-1]
         width = np.diff(slicer.slice_points["bins"])
@@ -62,7 +65,7 @@ def __call__(self, metric_values, slicer, user_plot_dict, fignum=None):
                 color=plot_dict["color"],
             )
         else:
-            good = np.where(metric_values.mask == False)
+            good = np.where(~metric_values.mask)
             x = np.ravel(list(zip(leftedge[good], leftedge[good] + width[good])))
             y = np.ravel(list(zip(metric_values[good], metric_values[good])))
             if plot_dict["log_scale"]:
@@ -89,7 +92,8 @@ def __call__(self, metric_values, slicer, user_plot_dict, fignum=None):
             plt.ylabel(plot_dict["ylabel"], fontsize=plot_dict["fontsize"])
         if "xlabel" in plot_dict:
             plt.xlabel(plot_dict["xlabel"], fontsize=plot_dict["fontsize"])
-        # Set y limits (either from values in args, percentile_clipping or compressed data values).
+        # Set y limits (either from values in args,
+        # percentile_clipping or compressed data values).
         if plot_dict["percentile_clip"] is not None:
             y_min, y_max = percentile_clipping(
                 metric_values.compressed(), percentile=plot_dict["percentile_clip"]
@@ -115,4 +119,4 @@ def __call__(self, metric_values, slicer, user_plot_dict, fignum=None):
         if plot_dict["x_max"] is not None:
             plt.xlim(right=plot_dict["x_max"])
         plt.title(plot_dict["title"])
-        return fig.number
+        return fig
diff --git a/rubin_sim/maf/plots/plot_bundle.py b/rubin_sim/maf/plots/plot_bundle.py
deleted file mode 100644
index e4114a9d..00000000
--- a/rubin_sim/maf/plots/plot_bundle.py
+++ /dev/null
@@ -1,86 +0,0 @@
-__all__ = ("PlotBundle",)
-
-
-import matplotlib.pylab as plt
-
-from .plot_handler import PlotHandler
-
-
-class PlotBundle:
-    """
-    Object designed to help organize multiple MetricBundles that will be plotted
-    together using the PlotHandler.
-    """
-
-    def __init__(self, bundle_list=None, plot_dicts=None, plot_func=None):
-        """
-        Init object and set things if desired.
-        bundle_list: A list of bundleDict objects
-        plot_dicts: A list of dictionaries with plotting kwargs
-        plot_func: A single MAF plotting function
-        """
-        if bundle_list is None:
-            self.bundle_list = []
-        else:
-            self.bundle_list = bundle_list
-
-        if plot_dicts is None:
-            if len(self.bundle_list) > 0:
-                self.plot_dicts = [{}]
-            else:
-                self.plot_dicts = []
-        else:
-            self.plot_dicts = plot_dicts
-
-        self.plot_func = plot_func
-
-    def add_bundle(self, bundle, plot_dict=None, plot_func=None):
-        """
-        Add bundle to the object.
-        Optionally add a plot_dict and/or replace the plot_func
-        """
-        self.bundle_list.append(bundle)
-        if plot_dict is not None:
-            self.plot_dicts.append(plot_dict)
-        else:
-            self.plot_dicts.append({})
-        if plot_func is not None:
-            self.plot_func = plot_func
-
-    def increment_plot_order(self):
-        """
-        Find the maximium order number in the display dicts, and set them to +1 that
-        """
-        max_order = 0
-        for m_b in self.bundle_list:
-            if "order" in list(m_b.displayDict.keys()):
-                max_order = max([max_order, m_b.displayDict["order"]])
-
-        for m_b in self.bundle_list:
-            m_b.displayDict["order"] = max_order + 1
-
-    def percentile_legend(self):
-        """
-        Go through the bundles and change the lables if there are the correct summary stats
-        """
-        for i, mB in enumerate(self.bundle_list):
-            if mB.summary_values is not None:
-                keys = list(mB.summary_values.keys())
-                if ("25th%ile" in keys) & ("75th%ile" in keys) & ("Median" in keys):
-                    if "label" not in list(self.plot_dicts[i].keys()):
-                        self.plot_dicts[i]["label"] = ""
-                    newstr = "%0.1f/%0.1f/%0.1f " % (
-                        mB.summary_values["25th%ile"],
-                        mB.summary_values["Median"],
-                        mB.summary_values["75th%ile"],
-                    )
-                    self.plot_dicts[i]["label"] = newstr + self.plot_dicts[i]["label"]
-
-    def plot(self, out_dir="Out", results_db=None, closefigs=True):
-        ph = PlotHandler(out_dir=out_dir, results_db=results_db)
-        ph.set_metric_bundles(self.bundle_list)
-        # Auto-generate labels and things
-        ph.set_plot_dicts(plot_dicts=self.plot_dicts, plot_func=self.plot_func)
-        ph.plot(self.plot_func, plot_dicts=self.plot_dicts)
-        if closefigs:
-            plt.close("all")
diff --git a/rubin_sim/maf/plots/plot_handler.py b/rubin_sim/maf/plots/plot_handler.py
index 45cb94bb..0fb605bf 100644
--- a/rubin_sim/maf/plots/plot_handler.py
+++ b/rubin_sim/maf/plots/plot_handler.py
@@ -26,7 +26,8 @@ class BasePlotter:
 
     def __init__(self):
         self.plot_type = None
-        # This should be included in every subsequent default_plot_dict (assumed to be present).
+        # This should be included in every subsequent default_plot_dict
+        # (assumed to be present).
         self.default_plot_dict = {
             "title": None,
             "xlabel": None,
@@ -36,11 +37,51 @@ def __init__(self):
             "figsize": None,
         }
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
+        """
+        Parameters
+        ----------
+        metric_value : `numpy.ma.MaskedArray`
+            The metric values from the bundle.
+        slicer : `rubin_sim.maf.slicers.TwoDSlicer`
+            The slicer.
+        user_plot_dict: `dict`
+            Dictionary of plot parameters set by user
+            (overrides default values).
+        fig : `matplotlib.figure.Figure`
+            Matplotlib figure number to use. Default = None, starts new figure.
+
+        Returns
+        -------
+        fig : `matplotlib.figure.Figure`
+           Figure with the plot.
+        """
         pass
 
 
 class PlotHandler:
+    """Create plots from a single or series of metric bundles.
+
+    Parameters
+    ----------
+    out_dir : `str`, optional
+        Directory to save output plots.
+    results_db : `rubin_sim.maf.ResultsDb`, optional
+        ResultsDb into which to record plot location and information.
+    savefig : `bool`, optional
+        Flag for saving images to disk (versus create and return
+        to caller).
+    fig_format : `str`, optional
+        Figure format to use to save full-size output. Default PDF.
+    dpi : `int`, optional
+        DPI to save output figures to disk at. (for matplotlib figures).
+    thumbnail : `bool`, optional
+        Flag for saving thumbnails (reduced size pngs) to disk.
+    trim_whitespace : `bool`, optional
+        Flag for trimming whitespace option to matplotlib output figures.
+        Default True, usually doesn't need to be changed.
+    """
+
     def __init__(
         self,
         out_dir=".",
@@ -110,9 +151,9 @@ def set_metric_bundles(self, m_bundles):
 
     def set_plot_dicts(self, plot_dicts=None, plot_func=None, reset=False):
         """
-        Set or update (or 'reset') the plot_dict for the (possibly joint) plots.
+        Set or update the plot_dict for the (possibly joint) plots.
 
-        Resolution is:
+        Resolution is: (from lowest to higher)
         auto-generated items (colors/labels/titles)
         < anything previously set in the plot_handler
         < defaults set by the plotter
@@ -120,7 +161,8 @@ def set_plot_dicts(self, plot_dicts=None, plot_func=None, reset=False):
         < explicitly set items in the plot_dicts list passed to this method.
         """
         if reset:
-            # Have to explicitly set each dictionary to a (separate) blank dictionary.
+            # Have to explicitly set each dictionary to a (separate)
+            # blank dictionary.
             self.plot_dicts = [{} for b in self.m_bundles]
 
         if isinstance(plot_dicts, dict):
@@ -142,9 +184,10 @@ def set_plot_dicts(self, plot_dicts=None, plot_func=None, reset=False):
             tmp_plot_dict["label"] = auto_label_list[i]
             tmp_plot_dict["color"] = auto_color_list[i]
             tmp_plot_dict["cbar_format"] = auto_cbar
-            # Then update that with anything previously set in the plot_handler.
+            # Update that with anything previously set in the plot_handler.
             tmp_plot_dict.update(self.plot_dicts[i])
-            # Then override with plot_dict items set explicitly based on the plot type.
+            # Then override with plot_dict items set explicitly
+            # based on the plot type.
             if plot_func is not None:
                 tmp_plot_dict["xlabel"] = auto_xlabel
                 tmp_plot_dict["ylabel"] = auto_ylabel
@@ -154,9 +197,10 @@ def set_plot_dicts(self, plot_dicts=None, plot_func=None, reset=False):
                 for k, v in plotter_defaults.items():
                     if v is not None:
                         tmp_plot_dict[k] = v
-            # Then add/override based on the bundle plot_dict parameters if they are set.
+            # Then add/override based on the bundle plot_dict parameters
+            # if they are set.
             tmp_plot_dict.update(bundle.plot_dict)
-            # Finally, override with anything set explicitly by the user right now.
+            # Finally, override with anything set explicitly by the user.
             if plot_dicts is not None:
                 tmp_plot_dict.update(plot_dicts[i])
             # And save this new dictionary back in the class.
@@ -178,15 +222,18 @@ def _combine_metric_names(self):
         if len(self.metric_names) == 1:
             joint_name = " ".join(self.metric_names)
         else:
-            # Split each unique name into a list to see if we can merge the names.
+            # Split each unique name into a list to see
+            # if we can merge the names.
             name_lengths = [len(x.split()) for x in self.metric_names]
             name_lists = [x.split() for x in self.metric_names]
-            # If the metric names are all the same length, see if we can combine any parts.
+            # If the metric names are all the same length, see
+            # if we can combine any parts.
             if len(set(name_lengths)) == 1:
                 joint_name = []
                 for i in range(name_lengths[0]):
                     tmp = set([x[i] for x in name_lists])
-                    # Try to catch special case of filters and put them in order.
+                    # Try to catch special case of filters and
+                    # put them in order.
                     if tmp.intersection(order) == tmp:
                         filterlist = ""
                         for f in order:
@@ -197,7 +244,8 @@ def _combine_metric_names(self):
                         # Otherwise, just join and put into joint_name.
                         joint_name.append("".join(tmp))
                 joint_name = " ".join(joint_name)
-            # If the metric names are not the same length, just join everything.
+            # If the metric names are not the same length,
+            # just join everything.
             else:
                 joint_name = " ".join(self.metric_names)
         self.joint_metric_names = joint_name
@@ -225,10 +273,12 @@ def _combine_metadata(self):
         else:
             order = ["u", "g", "r", "i", "z", "y"]
             # See if there are any subcomponents we can combine,
-            # splitting on some values we expect to separate info_label clauses.
+            # splitting on some values we expect to separate
+            # info_label clauses.
             splitmetas = []
             for m in self.info_label:
-                # Try to split info_label into separate phrases (filter / proposal / constraint..).
+                # Try to split info_label into separate phrases
+                # (filter / proposal / constraint..).
                 if " and " in m:
                     m = m.split(" and ")
                 elif ", " in m:
@@ -240,10 +290,12 @@ def _combine_metadata(self):
                 # Strip white spaces from individual elements.
                 m = set([im.strip() for im in m])
                 splitmetas.append(m)
-            # Look for common elements and separate from the general info_label.
+            # Look for common elements and
+            # separate from the general info_label.
             common = set.intersection(*splitmetas)
             diff = [x.difference(common) for x in splitmetas]
-            # Now look within the 'diff' elements and see if there are any common words to split off.
+            # Now look within the 'diff' elements and
+            # see if there are any common words to split off.
             diffsplit = []
             for d in diff:
                 if len(d) > 0:
@@ -253,10 +305,12 @@ def _combine_metadata(self):
                 diffsplit.append(m)
             diffcommon = set.intersection(*diffsplit)
             diffdiff = [x.difference(diffcommon) for x in diffsplit]
-            # If the length of any of the 'differences' is 0, then we should stop and not try to subdivide.
+            # If the length of any of the 'differences' is 0,
+            # then we should stop and not try to subdivide.
             lengths = [len(x) for x in diffdiff]
             if min(lengths) == 0:
-                # Sort them in order of length (so it goes 'g', 'g dithered', etc.)
+                # Sort them in order of length
+                # (so it goes 'g', 'g dithered', etc.)
                 tmp = []
                 for d in diff:
                     tmp.append(list(d)[0])
@@ -266,7 +320,8 @@ def _combine_metadata(self):
                 diffdiff = [diff[i] for i in idx]
                 diffcommon = []
             else:
-                # diffdiff is the part where we might expect our filter values to appear;
+                # diffdiff is the part where we might expect
+                # our filter values to appear;
                 # try to put this in order.
                 diffdiff_ordered = []
                 diffdiff_end = []
@@ -304,7 +359,8 @@ def _build_title(self):
         """
         Build a plot title from the metric names, runNames and info_label.
         """
-        # Create a plot title from the unique parts of the metric/run_name/info_label.
+        # Create a plot title from the unique parts of
+        # the metric/run_name/info_label.
         plot_title = ""
         if len(self.run_names) == 1:
             plot_title += list(self.run_names)[0]
@@ -313,7 +369,8 @@ def _build_title(self):
         if len(self.metric_names) == 1:
             plot_title += ": " + list(self.metric_names)[0]
         if plot_title == "":
-            # If there were more than one of everything above, use joint info_label and metricNames.
+            # If there were more than one of everything above,
+            # use joint info_label and metricNames.
             plot_title = self.joint_metadata + " " + self.joint_metric_names
         return plot_title
 
@@ -323,7 +380,7 @@ def _build_x_ylabels(self, plot_func, len_max=25):
 
         Parameters
         ----------
-        len_max : `int` (30)
+        len_max : `int`, optional
             If the xlabel starts longer than this, add the units as a newline.
         """
         if plot_func.plot_type == "BinnedData":
@@ -372,7 +429,8 @@ def _build_x_ylabels(self, plot_func, len_max=25):
 
     def _build_legend_labels(self):
         """
-        Build a set of legend labels, using parts of the run_name/info_label/metricNames that change.
+        Build a set of legend labels,
+        using the parts of the run_name/info_label/metricNames that change.
         """
         if len(self.m_bundles) == 1:
             return [None]
@@ -449,11 +507,15 @@ def _build_cbar_format(self):
     def _build_file_root(self, outfile_suffix=None):
         """
         Build a root filename for plot outputs.
-        If there is only one metricBundle, this is equal to the metricBundle fileRoot + outfile_suffix.
-        For multiple metricBundles, this is created from the runNames, info_label and metric names.
-
-        If you do not wish to use the automatic filenames, then you could set 'savefig' to False and
-          save the file manually to disk, using the plot figure numbers returned by 'plot'.
+        If there is only one metricBundle,
+        this is equal to the metricBundle fileRoot + outfile_suffix.
+        For multiple metricBundles,
+        this is created from the runNames, info_label and metric names.
+
+        If you do not wish to use the automatic filenames,
+        then you could set 'savefig' to False and
+        save the file manually to disk,
+        using the plot figure numbers returned by 'plot'.
         """
         if len(self.m_bundles) == 1:
             outfile = self.m_bundles[0].file_root
@@ -468,7 +530,8 @@ def _build_file_root(self, outfile_suffix=None):
     def _build_display_dict(self):
         """
         Generate a display dictionary.
-        This is most useful for when there are many metricBundles being combined into a single plot.
+        This is most useful for when there are many metricBundles
+        being combined into a single plot.
         """
         if len(self.m_bundles) == 1:
             return self.m_bundles[0].display_dict
@@ -506,7 +569,8 @@ def _build_display_dict(self):
 
     def _check_plot_dicts(self):
         """
-        Check to make sure there are no conflicts in the plot_dicts that are being used in the same subplot.
+        Check to make sure there are no conflicts in the plot_dicts
+        that are being used in the same subplot.
         """
         # Check that the length is OK
         if len(self.plot_dicts) != len(self.m_bundles):
@@ -518,8 +582,10 @@ def _check_plot_dicts(self):
         # These are the keys that need to match (or be None)
         keys2_check = ["xlim", "ylim", "color_min", "color_max", "title"]
 
-        # Identify how many subplots there are. If there are more than one, just don't change anything.
-        # This assumes that if there are more than one, the plot_dicts are actually all compatible.
+        # Identify how many subplots there are.
+        # If there are more than one, just don't change anything.
+        # This assumes that if there are more than one,
+        # the plot_dicts are actually all compatible.
         subplots = set()
         for pd in self.plot_dicts:
             if "subplot" in pd:
@@ -531,7 +597,8 @@ def _check_plot_dicts(self):
             for key in keys2_check:
                 values = [pd[key] for pd in self.plot_dicts if key in pd]
                 if len(np.unique(values)) > 1:
-                    # We will reset some of the keys to the default, but for some we should do better.
+                    # We will reset some of the keys to the default,
+                    # but for some we should do better.
                     if key.endswith("Max"):
                         for pd in self.plot_dicts:
                             pd[key] = np.max(values)
@@ -548,7 +615,8 @@ def _check_plot_dicts(self):
                             + " Will reset to default value. (found values %s)" % values
                         )
                         reset_keys.append(key)
-            # Reset the most of the keys to defaults; this can generally be done safely.
+            # Reset the most of the keys to defaults;
+            # this can generally be done safely.
             for key in reset_keys:
                 for pd in self.plot_dicts:
                     pd[key] = None
@@ -562,13 +630,34 @@ def plot(
         outfile_suffix=None,
     ):
         """
-        Create plot for mBundles, using plot_func.
+        Create a plot for the active metric bundles (self.set_metric_bundles).
 
-        plot_dicts:  List of plot_dicts if one wants to use a _new_ plot_dict per MetricBundle.
+        Parameters
+        ----------
+        plot_func : `rubin_sim.plots.BasePlotter`
+            The plotter to use to make the figure.
+        plot_dicts : `list` of [`dict`], optional
+            List of plot_dicts for each metric bundle.
+            Can use these to override individual metric bundle colors, etc.
+        display_dict : `dict`, optional
+            Information to save to resultsDb to accompany the figure on the
+            show_maf pages. Generally set automatically. Includes a caption.
+        outfile_root : `str`, optional
+            Output filename. Generally set automatically, but can be
+            overriden (such as when output filenames get too long).
+        outfile_suffix : `str`, optional
+            A suffix to add to the end of the default output filename.
+            Useful when creating a series of plots, such as for a movie.
+
+        Returns
+        -------
+        fig : `matplotlib.figure.Figure`
+            The plot.
         """
         if not plot_func.object_plotter:
-            # Check that metric_values type and plotter are compatible (most are float/float, but
-            #  some plotters expect object data .. and some only do sometimes).
+            # Check that metric_values type and plotter are compatible
+            # (most are float/float, but some plotters expect object data ..
+            # and some only do sometimes).
             for m_b in self.m_bundles:
                 if m_b.metric.metric_dtype == "object":
                     metric_is_color = m_b.plot_dict.get("metric_is_color", False)
@@ -587,7 +676,7 @@ def plot(
         if len(self.m_bundles) > 1:
             plot_type = "Combo" + plot_type
         # Make plot.
-        fignum = None
+        fig = None
         for m_b, plot_dict in zip(self.m_bundles, self.plot_dicts):
             if m_b.metric_values is None:
                 # Skip this metricBundle.
@@ -598,8 +687,9 @@ def plot(
                 msg = "MetricBundle (%s) has no unmasked metric_values, skipping plots." % (m_b.file_root)
                 warnings.warn(msg)
             else:
-                fignum = plot_func(m_b.metric_values, m_b.slicer, plot_dict, fignum=fignum)
-        # Add a legend if more than one metricValue is being plotted or if legendloc is specified.
+                fig = plot_func(m_b.metric_values, m_b.slicer, plot_dict, fig=fig)
+        # Add a legend if more than one metricValue is being plotted
+        # or if legendloc is specified.
         legend_loc = None
         if "legend_loc" in self.plot_dicts[0]:
             legend_loc = self.plot_dicts[0]["legend_loc"]
@@ -609,7 +699,8 @@ def plot(
             except KeyError:
                 legend_loc = "upper right"
         if legend_loc is not None:
-            plt.figure(fignum)
+            # Activate expected figure and write legend.
+            plt.figure(fig)
             plt.legend(loc=legend_loc, fancybox=True, fontsize="smaller")
         # Add the super title if provided.
         if "suptitle" in self.plot_dicts[0]:
@@ -619,7 +710,7 @@ def plot(
             if display_dict is None:
                 display_dict = self._build_display_dict()
             self.save_fig(
-                fignum,
+                fig,
                 outfile,
                 plot_type,
                 self.joint_metric_names,
@@ -629,11 +720,11 @@ def plot(
                 self.joint_metadata,
                 display_dict,
             )
-        return fignum
+        return fig
 
     def save_fig(
         self,
-        fignum,
+        fig,
         outfile_root,
         plot_type,
         metric_name,
@@ -643,8 +734,10 @@ def save_fig(
         info_label,
         display_dict=None,
     ):
-        fig = plt.figure(fignum)
         plot_file = outfile_root + "_" + plot_type + "." + self.fig_format
+        if fig is None:
+            warnings.warn(f"Trying to save figure to {plot_file} but" "figure is None. Skipping.")
+            return
         if self.trim_whitespace:
             fig.savefig(
                 os.path.join(self.out_dir, plot_file),
@@ -661,7 +754,7 @@ def save_fig(
         # Generate a png thumbnail.
         if self.thumbnail:
             thumb_file = "thumb." + outfile_root + "_" + plot_type + ".png"
-            plt.savefig(os.path.join(self.out_dir, thumb_file), dpi=72, bbox_inches="tight")
+            fig.savefig(os.path.join(self.out_dir, thumb_file), dpi=72, bbox_inches="tight")
         # Save information about the file to results_db.
         if self.results_db:
             if display_dict is None:
diff --git a/rubin_sim/maf/plots/spatial_plotters.py b/rubin_sim/maf/plots/spatial_plotters.py
index 504b4747..f71fc67f 100644
--- a/rubin_sim/maf/plots/spatial_plotters.py
+++ b/rubin_sim/maf/plots/spatial_plotters.py
@@ -68,7 +68,7 @@ def set_color_lims(metric_value, plot_dict, key_min="color_min", key_max="color_
     # Use plotdict values if available
     color_min = plot_dict[key_min]
     color_max = plot_dict[key_max]
-    # If either is not set and we have data ..
+    # If either is not set and we have data:
     if color_min is None or color_max is None:
         if np.size(metric_value.compressed()) > 0:
             # is percentile clipping set?
@@ -102,7 +102,7 @@ def set_color_map(plot_dict):
     cmap = plot_dict["cmap"]
     if cmap is None:
         cmap = "perceptual_rainbow"
-    if type(cmap) == str:
+    if isinstance(cmap, str):
         cmap = getattr(cm, cmap)
     # Set background and masked pixel colors default healpy white and gray.
     cmap = copy.copy(cmap)
@@ -141,23 +141,24 @@ def __init__(self):
         self.ax = None
         self.im = None
 
-    def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value_in, slicer, user_plot_dict, fig=None):
         """
         Parameters
         ----------
         metric_value : `numpy.ma.MaskedArray`
-        slicer : `rubin_sim.maf.slicers.HealpixSlicer`
+            The metric values from the bundle.
+        slicer : `rubin_sim.maf.slicers.TwoDSlicer`
+            The slicer.
         user_plot_dict: `dict`
             Dictionary of plot parameters set by user
             (overrides default values).
-        fignum : `int`
-            Matplotlib figure number to use
-            (default = None, starts new figure).
+        fig : `matplotlib.figure.Figure`
+            Matplotlib figure number to use. Default = None, starts new figure.
 
         Returns
         -------
-        fignum : `int`
-           Matplotlib figure number used to create the plot.
+        fig : `matplotlib.figure.Figure`
+           Figure with the plot.
         """
         # Override the default plotting parameters with user specified values.
         plot_dict = {}
@@ -194,7 +195,8 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
             badval = slicer.badval
 
         # Generate a full-sky plot.
-        fig = plt.figure(fignum, figsize=plot_dict["figsize"])
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
         # Set up color bar limits.
         clims = set_color_lims(metric_value, plot_dict)
         cmap = set_color_map(plot_dict)
@@ -306,7 +308,7 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
         # If outputing to PDF, this fixes the colorbar white stripes
         if plot_dict["cbar_edge"]:
             cb.solids.set_edgecolor("face")
-        return fig.number
+        return fig
 
 
 class HealpixPowerSpectrum(BasePlotter):
@@ -317,7 +319,7 @@ def __init__(self):
         self.default_plot_dict.update(base_default_plot_dict)
         self.default_plot_dict.update({"maxl": None, "removeDipole": True, "linestyle": "-"})
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
         """Generate and plot the power spectrum of metric_values
         (for metrics calculated on a healpix grid).
         """
@@ -327,7 +329,8 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         plot_dict.update(self.default_plot_dict)
         plot_dict.update(user_plot_dict)
 
-        fig = plt.figure(fignum, figsize=plot_dict["figsize"])
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
         if plot_dict["subplot"] != "111":
             fig.add_subplot(plot_dict["subplot"])
         # If the mask is True everywhere (no data), just plot zeros
@@ -366,7 +369,7 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
             plt.title(plot_dict["title"])
         # Return figure number
         # (so we can reuse/add onto/save this figure if desired).
-        return fig.number
+        return fig
 
 
 class HealpixHistogram(BasePlotter):
@@ -387,7 +390,7 @@ def __init__(self):
         )
         self.base_hist = BaseHistogram()
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
         """Histogram metric_value for all healpix points."""
         if "Healpix" not in slicer.slicer_name:
             raise ValueError("HealpixHistogram is for use with healpix slicer.")
@@ -396,8 +399,8 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         plot_dict.update(user_plot_dict)
         if plot_dict["scale"] is None:
             plot_dict["scale"] = hp.nside2pixarea(slicer.nside, degrees=True) / 1000.0
-        fignum = self.base_hist(metric_value, slicer, plot_dict, fignum=fignum)
-        return fignum
+        fig = self.base_hist(metric_value, slicer, plot_dict, fig=fig)
+        return fig
 
 
 class BaseHistogram(BasePlotter):
@@ -418,11 +421,13 @@ def __init__(self):
             }
         )
 
-    def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value_in, slicer, user_plot_dict, fig=None):
         """
-        Plot a histogram of metric_values (such as would come from a spatial slicer).
+        Plot a histogram of metric_values
+        (such as would come from a spatial slicer).
         """
-        # Adjust metric values by zeropoint or norm_val, and use 'compressed' version of masked array.
+        # Adjust metric values by zeropoint or norm_val,
+        # and use 'compressed' version of masked array.
         plot_dict = {}
         plot_dict.update(self.default_plot_dict)
         plot_dict.update(user_plot_dict)
@@ -431,20 +436,25 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
         metric_value = metric_value[np.isfinite(metric_value)]
         x_min, x_max = set_color_lims(metric_value, plot_dict, key_min="x_min", key_max="x_max")
         metric_value = metric_value.compressed()
-        # Set up the bins for the histogram. User specified 'bins' overrides 'bin_size'.
-        # Note that 'bins' could be a single number or an array, simply passed to plt.histogram.
+        # Set up the bins for the histogram.
+        # User specified 'bins' overrides 'bin_size'.
+        # Note that 'bins' could be a single number or an array,
+        # simply passed to plt.histogram.
         if plot_dict["bins"] is not None:
             bins = plot_dict["bins"]
         elif plot_dict["bin_size"] is not None:
-            #  If generating a cumulative histogram, want to use full range of data (but with given bin_size).
-            #    .. but if user set histRange to be wider than full range of data, then
-            #       extend bins to cover this range, so we can make prettier plots.
+            #  If generating a cumulative histogram,
+            #  want to use full range of data (but with given bin_size).
+            #  but if user set histRange to be wider than full range of data,
+            #  then extend bins to cover this range,
+            #  so we can make prettier plots.
             if plot_dict["cumulative"]:
                 # Potentially, expand the range for the cumulative histogram.
                 bmin = np.min([metric_value.min(), x_min])
                 bmax = np.max([metric_value.max(), x_max])
                 bins = np.arange(bmin, bmax + plot_dict["bin_size"] / 2.0, plot_dict["bin_size"])
-            #  Otherwise, not cumulative so just use metric values, without potential expansion.
+            #  Otherwise, not cumulative so just use metric values,
+            #  without potential expansion.
             else:
                 bins = np.arange(
                     x_min,
@@ -459,10 +469,12 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
                     plot_dict["bin_size"],
                 )
         else:
-            # If user did not specify bins or bin_size, then we try to figure out a good number of bins.
+            # If user did not specify bins or bin_size,
+            # then we try to figure out a good number of bins.
             bins = optimal_bins(metric_value)
         # Generate plots.
-        fig = plt.figure(fignum, figsize=plot_dict["figsize"])
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
         if (
             plot_dict["subplot"] != 111
             and plot_dict["subplot"] != (1, 1, 1)
@@ -471,14 +483,16 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
             ax = fig.add_subplot(plot_dict["subplot"])
         else:
             ax = plt.gca()
-        # Check if any data falls within histRange, because otherwise histogram generation will fail.
+        # Check if any data falls within histRange,
+        # because otherwise histogram generation will fail.
         if isinstance(bins, np.ndarray):
             condition = (metric_value >= bins.min()) & (metric_value <= bins.max())
         else:
             condition = (metric_value >= x_min) & (metric_value <= x_max)
         plot_value = metric_value[condition]
         if len(plot_value) == 0:
-            # No data is within histRange/bins. So let's just make a simple histogram anyway.
+            # No data is within histRange/bins.
+            # So let's just make a simple histogram anyway.
             n, b, p = plt.hist(
                 metric_value,
                 bins=50,
@@ -489,7 +503,8 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
                 color=plot_dict["color"],
             )
         else:
-            # There is data to plot, and we've already ensured x_min/x_max/bins are more than single value.
+            # There is data to plot, and we've already ensured
+            # x_min/x_max/bins are more than single value.
             n, b, p = plt.hist(
                 metric_value,
                 bins=bins,
@@ -530,8 +545,8 @@ def mjr_formatter(y, pos):
         if plot_dict["labelsize"] is not None:
             plt.tick_params(axis="x", labelsize=plot_dict["labelsize"])
             plt.tick_params(axis="y", labelsize=plot_dict["labelsize"])
-        # Return figure number
-        return fig.number
+        # Return figure
+        return fig
 
 
 class BaseSkyMap(BasePlotter):
@@ -554,7 +569,7 @@ def __init__(self):
             }
         )
 
-    def _plot_tissot_ellipse(self, lon, lat, radius, ax=None, **kwargs):
+    def _plot_tissot_ellipse(self, lon, lat, radius, **kwargs):
         """Plot Tissot Ellipse/Tissot Indicatrix
 
         Parameters
@@ -565,7 +580,6 @@ def _plot_tissot_ellipse(self, lon, lat, radius, ax=None, **kwargs):
         latitude-like of ellipse centers (radians)
         radius : float or array_like
         radius of ellipses (radians)
-        ax : Axes object (optional)
         matplotlib axes instance on which to draw ellipses.
 
         Other Parameters
@@ -578,10 +592,8 @@ def _plot_tissot_ellipse(self, lon, lat, radius, ax=None, **kwargs):
         # Code adapted from astroML, which is BSD-licensed.
         # See http: //github.com/astroML/astroML for details.
         ellipses = []
-        if ax is None:
-            ax = plt.gca()
-        for l, b, diam in np.broadcast(lon, lat, radius * 2.0):
-            el = Ellipse((l, b), diam / np.cos(b), diam, **kwargs)
+        for lon, b, diam in np.broadcast(lon, lat, radius * 2.0):
+            el = Ellipse((lon, b), diam / np.cos(b), diam, **kwargs)
             ellipses.append(el)
         return ellipses
 
@@ -628,7 +640,7 @@ def _plot_mw_zone(
         lon = -(ra - ra_cen - np.pi) % (np.pi * 2) - np.pi
         ax.plot(lon, dec, "b.", markersize=1.8, alpha=0.4)
 
-    def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value_in, slicer, user_plot_dict, fig=None):
         """
         Plot the sky map of metric_value for a generic spatial slicer.
         """
@@ -641,7 +653,8 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
         plot_dict.update(user_plot_dict)
         metric_value = apply_zp_norm(metric_value_in, plot_dict)
 
-        fig = plt.figure(fignum, figsize=plot_dict["figsize"])
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
         # other projections available include
         # ['aitoff', 'hammer', 'lambert', 'mollweide', 'polar', 'rectilinear']
         ax = fig.add_subplot(plot_dict["subplot"], projection=plot_dict["projection"])
@@ -650,7 +663,8 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
             # Plot all data points.
             good = np.ones(len(metric_value), dtype="bool")
         else:
-            # Only plot points which are not masked. Flip numpy ma mask where 'False' == 'good'.
+            # Only plot points which are not masked.
+            # Flip numpy ma mask where 'False' == 'good'.
             good = ~metric_value.mask
 
         # Add ellipses at RA/Dec locations - but don't add colors yet.
@@ -677,7 +691,8 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
             if current:
                 ax.add_patch(current)
         else:
-            # Determine color min/max values. metricValue.compressed = non-masked points.
+            # Determine color min/max values.
+            # metricValue.compressed = non-masked points.
             clims = set_color_lims(metric_value, plot_dict)
             # Determine whether or not to use auto-log scale.
             if plot_dict["log_scale"] == "auto":
@@ -689,7 +704,8 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
                 else:
                     plot_dict["log_scale"] = False
             if plot_dict["log_scale"]:
-                # Move min/max values to things that can be marked on the colorbar.
+                # Move min/max values to things that can be marked
+                # on the colorbar.
                 # clims[0] = 10 ** (int(np.log10(clims[0])))
                 # clims[1] = 10 ** (int(np.log10(clims[1])))
                 norml = colors.LogNorm()
@@ -725,7 +741,7 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
                         format=plot_dict["cbar_format"],
                     )
                 else:
-                    # Most of the time we just want a standard horizontal colorbar
+                    # Usually we just want a standard horizontal colorbar
                     cb = plt.colorbar(
                         p,
                         shrink=0.75,
@@ -763,7 +779,7 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
                 transform=ax.transAxes,
                 fontsize=plot_dict["fontsize"],
             )
-        return fig.number
+        return fig
 
 
 class HealpixSDSSSkyMap(BasePlotter):
@@ -783,11 +799,13 @@ def __init__(self):
             }
         )
 
-    def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value_in, slicer, user_plot_dict, fig=None):
         """
-        Plot the sky map of metric_value using healpy cartview plots in thin strips.
+        Plot the sky map of metric_value using healpy cartview plots
+        in thin strips.
         raMin: Minimum RA to plot (deg)
-        raMax: Max RA to plot (deg).  Note raMin/raMax define the centers that will be plotted.
+        raMax: Max RA to plot (deg).
+        Note raMin/raMax define the centers that will be plotted.
         raLen:  Length of the plotted strips in degrees
         decMin: minimum dec value to plot
         decMax: max dec value to plot
@@ -804,8 +822,10 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
         cmap = set_color_map(plot_dict)
         racenters = np.arange(plot_dict["raMin"], plot_dict["raMax"], plot_dict["raLen"])
         nframes = racenters.size
-        fig = plt.figure(fignum)
-        # Do not specify or use plot_dict['subplot'] because this is done in each call to hp.cartview.
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
+        # Do not specify or use plot_dict['subplot']
+        # because this is done in each call to hp.cartview.
         for i, racenter in enumerate(racenters):
             if i == 0:
                 use_title = (
@@ -834,8 +854,8 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
                 fig=fig,
             )
             hp.graticule(dpar=20, dmer=20, verbose=False)
-        # Add colorbar (not using healpy default colorbar because want more tickmarks).
-        ax = fig.add_axes([0.1, 0.15, 0.8, 0.075])  # left, bottom, width, height
+        # Add colorbar (not using healpy default colorbar)
+        ax = fig.add_axes([0.1, 0.15, 0.8, 0.075])
         # Add label.
         if plot_dict["label"] is not None:
             plt.figtext(0.8, 0.9, "%s" % plot_dict["label"])
@@ -865,8 +885,7 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
         # If outputing to PDF, this fixes the colorbar white stripes
         if plot_dict["cbar_edge"]:
             cb.solids.set_edgecolor("face")
-        fig = plt.gcf()
-        return fig.number
+        return fig
 
 
 def project_lambert(longitude, latitude):
@@ -935,7 +954,7 @@ def __init__(self):
             }
         )
 
-    def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value_in, slicer, user_plot_dict, fig=None):
         if "ra" not in slicer.slice_points or "dec" not in slicer.slice_points:
             err_message = 'SpatialSlicer must contain "ra" and "dec" in slice_points metadata.'
             err_message += " SlicePoints only contains keys %s." % (slicer.slice_points.keys())
@@ -954,16 +973,19 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
             step = (clims[1] - clims[0]) / plot_dict["levels"]
             levels = np.arange(clims[0], clims[1] + step, step)
 
-        fig = plt.figure(fignum, figsize=plot_dict["figsize"])
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
         ax = fig.add_subplot(plot_dict["subplot"])
 
         x, y = project_lambert(slicer.slice_points["ra"], slicer.slice_points["dec"])
-        # Contour the plot first to remove any anti-aliasing artifacts.  Doesn't seem to work though. See:
-        # http: //stackoverflow.com/questions/15822159/aliasing-when-saving-matplotlib\
-        # -filled-contour-plot-to-pdf-or-eps
+        # Contour the plot first to remove any anti-aliasing artifacts.
+        # Doesn't seem to work though. See:
+        # http://stackoverflow.com/questions/15822159/
+        # aliasing-when-saving-matplotlib-filled-contour-plot-to-pdf-or-eps
         # tmpContour = m.contour(np.degrees(slicer.slice_points['ra']),
         #                        np.degrees(slicer.slice_points['dec']),
-        #                        metric_value.filled(np.min(clims)-1), levels, tri=True,
+        #                        metric_value.filled(np.min(clims)-1), levels,
+        #                        tri=True,
         #                        cmap=plot_dict['cmap'], ax=ax, latlon=True,
         #                        lw=1)
 
@@ -1005,4 +1027,4 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
         # If outputing to PDF, this fixes the colorbar white stripes
         if plot_dict["cbar_edge"]:
             cb.solids.set_edgecolor("face")
-        return fig.number
+        return fig
diff --git a/rubin_sim/maf/plots/special_plotters.py b/rubin_sim/maf/plots/special_plotters.py
index 18e11d6d..c9e579cf 100644
--- a/rubin_sim/maf/plots/special_plotters.py
+++ b/rubin_sim/maf/plots/special_plotters.py
@@ -31,32 +31,31 @@ def __init__(self):
             "y_max": None,
             "linewidth": 2,
             "reflinewidth": 2,
+            "figsize": None,
         }
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
         """
         Parameters
         ----------
         metric_value : `numpy.ma.MaskedArray`
-            The metric values calculated with `rubin_sim.maf.Count` and
-            a healpix slicer.
-        slicer : `rubin_sim.maf.slicers.HealpixSlicer`
+            The metric values from the bundle.
+        slicer : `rubin_sim.maf.slicers.TwoDSlicer`
+            The slicer.
         user_plot_dict: `dict`
-            Dictionary of plot parameters set by user to override defaults.
-            Note that asky and n_visits values set here and in the slicer
-            should be consistent, for plot labels and summary statistic
-            values to be consistent.
-        fignum : `int`
-            Matplotlib figure number to use. Default starts new figure.
+            Dictionary of plot parameters set by user
+            (overrides default values).
+        fig : `matplotlib.figure.Figure`
+            Matplotlib figure number to use. Default = None, starts new figure.
 
         Returns
         -------
-        fignum : `int`
-           Matplotlib figure number used to create the plot.
+        fig : `matplotlib.figure.Figure`
+           Figure with the plot.
         """
         if not hasattr(slicer, "nside"):
             raise ValueError("FOPlot to be used with healpix or healpix derived slicers.")
-        fig = plt.figure(fignum)
+
         plot_dict = {}
         plot_dict.update(self.default_plot_dict)
         plot_dict.update(user_plot_dict)
@@ -64,6 +63,9 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         if plot_dict["scale"] is None:
             plot_dict["scale"] = hp.nside2pixarea(slicer.nside, degrees=True) / 1000.0
 
+        if fig is None:
+            fig = plt.Figure(figsize=plot_dict["figsize"])
+
         # Expect metric_value to be something like number of visits
         cumulative_area = np.arange(1, metric_value.compressed().size + 1)[::-1] * plot_dict["scale"]
         plt.plot(
@@ -117,7 +119,7 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
             plt.xlim([x_min, x_max])
         if (y_min is not None) or (y_max is not None):
             plt.ylim([y_min, y_max])
-        return fig.number
+        return fig
 
 
 class SummaryHistogram(BasePlotter):
@@ -154,7 +156,7 @@ def __init__(self):
             "figsize": None,
         }
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
         """
         Parameters
         ----------
@@ -171,13 +173,13 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
             results as a step histogram (True) or as a series of values (False)
             'bins' (np.ndarray) sets the x values for the resulting plot
             and should generally match the bins used with the metric.
-        fignum : `int`
-            Matplotlib figure number to use. Default starts a new figure.
+        fig : `matplotlib.figure.Figure`
+            Matplotlib figure to use. Default starts a new figure.
 
         Returns
         -------
-        fignum: `int`
-           Matplotlib figure number used to create the plot.
+        fig: `matplotlib.figure.Figure`
+           Matplotlib figure used to create the plot.
         """
         plot_dict = {}
         plot_dict.update(self.default_plot_dict)
@@ -203,7 +205,8 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         if plot_dict["cumulative"]:
             final_hist = final_hist.cumsum()
 
-        fig = plt.figure(fignum, figsize=plot_dict["figsize"])
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
         bins = plot_dict["bins"]
         if plot_dict["histStyle"]:
             leftedge = bins[:-1]
@@ -248,4 +251,4 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         if plot_dict["xscale"] is not None:
             plt.xscale(plot_dict["xscale"])
 
-        return fig.number
+        return fig
diff --git a/rubin_sim/maf/plots/two_d_plotters.py b/rubin_sim/maf/plots/two_d_plotters.py
index 96473afd..b65124d2 100644
--- a/rubin_sim/maf/plots/two_d_plotters.py
+++ b/rubin_sim/maf/plots/two_d_plotters.py
@@ -33,24 +33,27 @@ def __init__(self):
             "aspect": "auto",
             "xextent": None,
             "origin": None,
+            "figsize": None,
         }
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
         """
         Parameters
         ----------
-        metric_value : numpy.ma.MaskedArray
-        slicer : rubin_sim.maf.slicers.BaseSpatialSlicer
-            (any spatial slicer)
-        user_plot_dict: dict
-            Dictionary of plot parameters set by user (overrides default values).
-        fignum : int
-            Matplotlib figure number to use (default = None, starts new figure).
+        metric_value : `numpy.ma.MaskedArray`
+            The metric values from the bundle.
+        slicer : `rubin_sim.maf.slicers.TwoDSlicer`
+            The slicer.
+        user_plot_dict: `dict`
+            Dictionary of plot parameters set by user
+            (overrides default values).
+        fig : `matplotlib.figure.Figure`
+            Matplotlib figure number to use. Default = None, starts new figure.
 
         Returns
         -------
-        int
-           Matplotlib figure number used to create the plot.
+        fig : `matplotlib.figure.Figure
+           Figure with the plot.
         """
         if "Healpix" in slicer.slicer_name:
             self.default_plot_dict["ylabel"] = "Healpix ID"
@@ -77,8 +80,9 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
             low_vals = np.where(metric_value.data < plot_dict["color_min"])
             metric_value.mask[low_vals] = True
 
-        figure = plt.figure(fignum)
-        ax = figure.add_subplot(111)
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
+        ax = fig.add_subplot(111)
         yextent = [0, slicer.nslice - 1]
         xextent = plot_dict["xextent"]
         extent = []
@@ -105,12 +109,31 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         # Fix white space on pdf's
         if plot_dict["cbar_edge"]:
             cb.solids.set_edgecolor("face")
-        return figure.number
+        return fig
 
 
 class VisitPairsHist(BasePlotter):
     """
-    Given an opsim2dSlicer, figure out what fraction of observations are in singles, pairs, triples, etc.
+    Given an TwoDSlicer, figure out what fraction of observations
+    are in singles, pairs, triples, etc.
+
+
+    Parameters
+    ----------
+    metric_value : `numpy.ma.MaskedArray`
+        The metric values from the bundle.
+    slicer : `rubin_sim.maf.slicers.TwoDSlicer`
+        The slicer.
+    user_plot_dict: `dict`
+        Dictionary of plot parameters set by user
+        (overrides default values).
+    fig : `matplotlib.figure.Figure`
+        Matplotlib figure number to use. Default = None, starts new figure.
+
+    Returns
+    -------
+    fig : `matplotlib.figure.Figure`
+       Figure with the plot.
     """
 
     def __init__(self):
@@ -124,24 +147,10 @@ def __init__(self):
             "color": "b",
             "xlim": [0, 20],
             "ylim": None,
+            "figsize": None,
         }
 
-    def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
-        """
-        Parameters
-        ----------
-        metric_value : numpy.ma.MaskedArray
-        slicer : rubin_sim.maf.slicers.TwoDSlicer
-        user_plot_dict: dict
-            Dictionary of plot parameters set by user (overrides default values).
-        fignum : int
-            Matplotlib figure number to use (default = None, starts new figure).
-
-        Returns
-        -------
-        int
-           Matplotlib figure number used to create the plot.
-        """
+    def __call__(self, metric_value, slicer, user_plot_dict, fig=None):
         plot_dict = {}
         plot_dict.update(self.default_plot_dict)
         # Don't clobber with None
@@ -155,8 +164,9 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         vals, bins = np.histogram(metric_value, bins)
         xvals = (bins[:-1] + bins[1:]) / 2.0
 
-        figure = plt.figure(fignum)
-        ax = figure.add_subplot(111)
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
+        ax = fig.add_subplot(111)
         ax.bar(xvals, vals * xvals, color=plot_dict["color"], label=plot_dict["label"])
         ax.set_xlabel(plot_dict["xlabel"])
         ax.set_ylabel(plot_dict["ylabel"])
@@ -164,4 +174,4 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None):
         ax.set_xlim(plot_dict["xlim"])
         ax.set_ylim(plot_dict["ylim"])
 
-        return figure.number
+        return fig
diff --git a/rubin_sim/maf/plots/xyplotter.py b/rubin_sim/maf/plots/xyplotter.py
index f5cbdd98..8837a9a6 100644
--- a/rubin_sim/maf/plots/xyplotter.py
+++ b/rubin_sim/maf/plots/xyplotter.py
@@ -12,15 +12,21 @@ class XyPlotter(BasePlotter):
     def __init__(self):
         self.object_plotter = True
         self.plot_type = "simple"
-        self.default_plot_dict = {"title": None, "xlabel": "", "ylabel": ""}
+        self.default_plot_dict = {
+            "title": None,
+            "xlabel": "",
+            "ylabel": "",
+            "figsize": None,
+        }
 
-    def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
+    def __call__(self, metric_value_in, slicer, user_plot_dict, fig=None):
         plot_dict = {}
         plot_dict.update(self.default_plot_dict)
         plot_dict.update(user_plot_dict)
         plot_dict.update(metric_value_in[0]["plot_dict"])
 
-        fig = plt.figure(fignum)
+        if fig is None:
+            fig = plt.figure(figsize=plot_dict["figsize"])
         ax = fig.add_subplot(111)
         x = metric_value_in[0]["x"]
         y = metric_value_in[0]["y"]
@@ -28,4 +34,4 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None):
         ax.set_title(plot_dict["title"])
         ax.set_xlabel(plot_dict["xlabel"])
         ax.set_ylabel(plot_dict["ylabel"])
-        return fig.number
+        return fig
diff --git a/rubin_sim/moving_objects/ooephemerides.py b/rubin_sim/moving_objects/ooephemerides.py
index 177471ec..9c61de20 100644
--- a/rubin_sim/moving_objects/ooephemerides.py
+++ b/rubin_sim/moving_objects/ooephemerides.py
@@ -51,7 +51,7 @@ class PyOrbEphemerides:
     Typical usage:
 
     >>> pyephs = PyOrbEphemerides()
-    >>> pyephs.setOrbits(orbits)
+    >>> pyephs.set_orbits(orbits)
     >>> ephs = pyephs.generateEphemerides(times, timeScale, obscode)
     """