diff --git a/autoarray/dataset/plot/interferometer_plots.py b/autoarray/dataset/plot/interferometer_plots.py
index a05e68918..027195240 100644
--- a/autoarray/dataset/plot/interferometer_plots.py
+++ b/autoarray/dataset/plot/interferometer_plots.py
@@ -9,6 +9,44 @@
 from autoarray.structures.grids.irregular_2d import Grid2DIrregular
 
 
+def _subplot_natural_dataset(
+    dataset,
+    output_path,
+    output_filename,
+    output_format,
+    colormap,
+    use_log10,
+    title_prefix=None,
+):
+    """
+    1x2 subplot of the natural-weighted dirty image and dirty beam of an array-free
+    ``Interferometer`` (built by ``from_stream`` / ``from_sparse_terms``), which carries no
+    visibilities, uv-wavelengths or transformer, so only these real-space terms can be drawn.
+    """
+    _pf = (lambda t: f"{title_prefix.rstrip()} {t}") if title_prefix else (lambda t: t)
+
+    fig, axes = subplots(1, 2, figsize=conf_subplot_figsize(1, 2))
+
+    plot_array(
+        dataset.dirty_image_natural,
+        ax=axes[0],
+        title=_pf("Dirty Image (Natural)"),
+        colormap=colormap,
+        use_log10=use_log10,
+    )
+    plot_array(
+        dataset.dirty_beam,
+        ax=axes[1],
+        title=_pf("Dirty Beam (Natural)"),
+        colormap=colormap,
+        use_log10=use_log10,
+    )
+
+    hide_unused_axes(axes)
+    tight_layout()
+    subplot_save(fig, output_path, output_filename, output_format)
+
+
 def subplot_interferometer_dataset(
     dataset,
     output_path: Optional[str] = None,
@@ -24,6 +62,9 @@ def subplot_interferometer_dataset(
     Panels: Visibilities | UV-Wavelengths | Amplitudes vs UV-distances |
             Phases vs UV-distances | Dirty Image | Dirty S/N Map
 
+    An array-free dataset (``dataset.is_array_free``) has no visibilities, so a 1x2 subplot
+    of its ``dirty_image_natural`` and ``dirty_beam`` is written to the same filename instead.
+
     Parameters
     ----------
     dataset
@@ -39,6 +80,17 @@ def subplot_interferometer_dataset(
     use_log10
         Apply log10 normalisation to image panels.
     """
+    if dataset.is_array_free:
+        return _subplot_natural_dataset(
+            dataset,
+            output_path=output_path,
+            output_filename=output_filename,
+            output_format=output_format,
+            colormap=colormap,
+            use_log10=use_log10,
+            title_prefix=title_prefix,
+        )
+
     _pf = (lambda t: f"{title_prefix.rstrip()} {t}") if title_prefix else (lambda t: t)
 
     fig, axes = subplots(2, 3, figsize=conf_subplot_figsize(2, 3))
@@ -104,6 +156,9 @@ def subplot_interferometer_dirty_images(
     """
     1x3 subplot of dirty image, dirty noise map, and dirty S/N map.
 
+    An array-free dataset (``dataset.is_array_free``) has no visibilities, so a 1x2 subplot
+    of its ``dirty_image_natural`` and ``dirty_beam`` is written to the same filename instead.
+
     Parameters
     ----------
     dataset
@@ -119,6 +174,16 @@ def subplot_interferometer_dirty_images(
     use_log10
         Apply log10 normalisation.
     """
+    if dataset.is_array_free:
+        return _subplot_natural_dataset(
+            dataset,
+            output_path=output_path,
+            output_filename=output_filename,
+            output_format=output_format,
+            colormap=colormap,
+            use_log10=use_log10,
+        )
+
     fig, axes = subplots(1, 3, figsize=conf_subplot_figsize(1, 3))
 
     plot_array(
@@ -164,7 +229,9 @@ def fits_interferometer(
       ``uv_wavelengths_path`` to write each component to its own FITS file.
     * **Single multi-HDU file** -- pass ``file_path`` to write all components
       into one FITS file with named extensions (``data``, ``noise_map``,
-      ``uv_wavelengths``).
+      ``uv_wavelengths``). An array-free dataset (``from_stream`` /
+      ``from_sparse_terms``) has none of these, so its natural-weighted
+      ``dirty_image_natural`` and ``dirty_beam`` are written instead.
 
     Parameters
     ----------
@@ -183,8 +250,9 @@ def fits_interferometer(
         values_list = []
         ext_name_list = []
 
-        values_list.append(np.asarray(dataset.data.in_array))
-        ext_name_list.append("data")
+        if dataset.data is not None:
+            values_list.append(np.asarray(dataset.data.in_array))
+            ext_name_list.append("data")
 
         if dataset.noise_map is not None:
             values_list.append(np.asarray(dataset.noise_map.in_array))
@@ -194,13 +262,19 @@ def fits_interferometer(
             values_list.append(np.asarray(dataset.uv_wavelengths))
             ext_name_list.append("uv_wavelengths")
 
+        if dataset.is_array_free:
+            values_list.append(np.asarray(dataset.dirty_image_natural.native))
+            ext_name_list.append("dirty_image_natural")
+            values_list.append(np.asarray(dataset.dirty_beam.native))
+            ext_name_list.append("dirty_beam")
+
         hdu_list = hdu_list_for_output_from(
             values_list=values_list,
             ext_name_list=ext_name_list,
         )
         write_hdu_list(hdu_list, file_path=file_path, overwrite=overwrite)
     else:
-        if data_path is not None:
+        if dataset.data is not None and data_path is not None:
             output_to_fits(
                 values=np.asarray(dataset.data.in_array),
                 file_path=data_path, overwrite=overwrite,
diff --git a/autoarray/fit/fit_interferometer.py b/autoarray/fit/fit_interferometer.py
index b0873b8be..802d3c934 100644
--- a/autoarray/fit/fit_interferometer.py
+++ b/autoarray/fit/fit_interferometer.py
@@ -13,6 +13,54 @@
 from autoarray import type as ty
 
 
+def dirty_model_image_natural_from(dataset, image) -> Array2D:
+    """
+    Returns the naturally weighted, normalised dirty image of the model visibilities of a real-space
+    `image`, `Re(F^H W F m) / sum(w)`, computed without a transformer or any visibility-sized array.
+
+    `W~ = Re(F^H W F)` is the operator cached on the dataset's `sparse_operator` (the one the sparse
+    inversion uses for its curvature matrix), so `W~ m / sum(w)` is one FFT convolution on the real-space
+    grid. It is the model counterpart of `Interferometer.dirty_image_natural` (the natural dirty image of
+    the data, `Re(F^H W d) / sum(w)`): their difference is the natural dirty residual map. It is available
+    on both dataset types -- an array-free one built by `from_stream` / `from_sparse_terms` and an in-memory
+    one after `apply_sparse_operator()` -- and is how the visualizers draw a model on an array-free dataset.
+
+    Parameters
+    ----------
+    dataset
+        The `Interferometer` dataset, which must carry a `sparse_operator`.
+    image
+        The model image `m` on the slim masked real-space grid of the dataset's `real_space_mask`.
+    """
+    sparse_operator = getattr(dataset, "sparse_operator", None)
+
+    if sparse_operator is None:
+        raise exc.DatasetException(
+            "The natural dirty model image `W~ m / sum(w)` needs the dataset's `sparse_operator`; call "
+            "`apply_sparse_operator()` on the dataset (an array-free dataset built by from_stream / "
+            "from_sparse_terms always carries one)."
+        )
+
+    sparse_terms = getattr(dataset, "sparse_terms", None)
+
+    if sparse_terms is not None:
+        sum_weights = float(sparse_terms.sum_weights)
+    else:
+        sum_weights = float(np.sum(dataset.noise_map.array.real**-2.0))
+
+    image = np.asarray(getattr(image, "array", image), dtype=np.float64)
+
+    operated_image = sparse_operator.operated_matrix_slim_from(
+        matrix_slim=image[:, None],
+        extent_index_for_masked_pixel=dataset.real_space_mask.extent_index_for_masked_pixel,
+        xp=np,
+    )[:, 0]
+
+    return Array2D(
+        values=np.asarray(operated_image) / sum_weights, mask=dataset.real_space_mask
+    )
+
+
 class FitInterferometer(FitDataset):
     def __init__(
         self,
diff --git a/autoarray/fit/plot/fit_interferometer_plots.py b/autoarray/fit/plot/fit_interferometer_plots.py
index 2cc805b95..7472582fd 100644
--- a/autoarray/fit/plot/fit_interferometer_plots.py
+++ b/autoarray/fit/plot/fit_interferometer_plots.py
@@ -7,6 +7,78 @@
 from autoarray.plot.utils import subplots, subplot_save, symmetric_vmin_vmax, hide_unused_axes, conf_subplot_figsize, tight_layout
 
 
+def _subplot_fit_natural(
+    fit,
+    model_image,
+    output_path,
+    output_filename,
+    output_format,
+    colormap,
+    use_log10,
+    residuals_symmetric_cmap,
+):
+    """
+    Subplot of the natural-weighted dirty images of a ``FitInterferometer`` whose dataset is
+    array-free (built by ``from_stream`` / ``from_sparse_terms``), which has no visibilities,
+    transformer or visibility-space residuals.
+
+    Panels: Dirty Image (Natural) | Dirty Model Image (Natural) | Dirty Residual Map (Natural),
+    the last two only when the real-space ``model_image`` is given. The dirty model image is
+    ``W~ m / sum(w)`` (see ``dirty_model_image_natural_from``) and the dirty residual map is the
+    dirty image minus it.
+    """
+    from autoarray.fit.fit_interferometer import dirty_model_image_natural_from
+
+    dataset = fit.dataset
+    dirty_image = dataset.dirty_image_natural
+
+    if model_image is None:
+        fig, axes = subplots(1, 1, figsize=conf_subplot_figsize(1, 1))
+        axes = [axes]
+    else:
+        fig, axes = subplots(1, 3, figsize=conf_subplot_figsize(1, 3))
+
+    plot_array(
+        dirty_image,
+        ax=axes[0],
+        title="Dirty Image (Natural)",
+        colormap=colormap,
+        use_log10=use_log10,
+    )
+
+    if model_image is not None:
+        dirty_model_image = dirty_model_image_natural_from(
+            dataset=dataset, image=model_image
+        )
+        dirty_residual_map = dirty_image - dirty_model_image
+
+        if residuals_symmetric_cmap:
+            vmin_r, vmax_r = symmetric_vmin_vmax(dirty_residual_map)
+        else:
+            vmin_r = vmax_r = None
+
+        plot_array(
+            dirty_model_image,
+            ax=axes[1],
+            title="Dirty Model Image (Natural)",
+            colormap=colormap,
+            use_log10=use_log10,
+        )
+        plot_array(
+            dirty_residual_map,
+            ax=axes[2],
+            title="Dirty Residual Map (Natural)",
+            colormap=colormap,
+            use_log10=False,
+            vmin=vmin_r,
+            vmax=vmax_r,
+        )
+
+    hide_unused_axes(axes)
+    tight_layout()
+    subplot_save(fig, output_path, output_filename, output_format)
+
+
 def subplot_fit_interferometer(
     fit,
     output_path: Optional[str] = None,
@@ -15,6 +87,7 @@ def subplot_fit_interferometer(
     colormap=None,
     use_log10: bool = False,
     residuals_symmetric_cmap: bool = True,
+    model_image=None,
 ):
     """
     2×3 subplot of ``FitInterferometer`` residuals in UV-plane.
@@ -38,7 +111,24 @@ def subplot_fit_interferometer(
     residuals_symmetric_cmap
         Not used here (UV-plane residuals are scatter plots); kept for API
         consistency.
+    model_image
+        The real-space model image, used only when the fit's dataset is array-free
+        (``fit.dataset.is_array_free``): the visibility-space panels cannot be drawn, so
+        the natural-weighted dirty image, dirty model image and dirty residual map are
+        plotted to the same filename instead (the last two only if this is given).
     """
+    if fit.dataset.is_array_free:
+        return _subplot_fit_natural(
+            fit,
+            model_image=model_image,
+            output_path=output_path,
+            output_filename=output_filename,
+            output_format=output_format,
+            colormap=colormap,
+            use_log10=use_log10,
+            residuals_symmetric_cmap=residuals_symmetric_cmap,
+        )
+
     fig, axes = subplots(2, 3, figsize=conf_subplot_figsize(2, 3))
     axes = axes.flatten()
 
@@ -110,6 +200,7 @@ def subplot_fit_interferometer_dirty_images(
     colormap=None,
     use_log10: bool = False,
     residuals_symmetric_cmap: bool = True,
+    model_image=None,
 ):
     """
     2×3 subplot of ``FitInterferometer`` dirty-image components.
@@ -133,7 +224,24 @@ def subplot_fit_interferometer_dirty_images(
         Apply log10 normalisation to non-residual panels.
     residuals_symmetric_cmap
         Centre residual colour scale symmetrically around zero.
+    model_image
+        The real-space model image, used only when the fit's dataset is array-free
+        (``fit.dataset.is_array_free``): the natural-weighted dirty image, dirty model
+        image and dirty residual map are plotted to the same filename instead (the last
+        two only if this is given).
     """
+    if fit.dataset.is_array_free:
+        return _subplot_fit_natural(
+            fit,
+            model_image=model_image,
+            output_path=output_path,
+            output_filename=output_filename,
+            output_format=output_format,
+            colormap=colormap,
+            use_log10=use_log10,
+            residuals_symmetric_cmap=residuals_symmetric_cmap,
+        )
+
     fig, axes = subplots(2, 3, figsize=conf_subplot_figsize(2, 3))
     axes = axes.flatten()
 
diff --git a/autoarray/inversion/plot/inversion_plots.py b/autoarray/inversion/plot/inversion_plots.py
index 02a3f061e..22caa36ef 100644
--- a/autoarray/inversion/plot/inversion_plots.py
+++ b/autoarray/inversion/plot/inversion_plots.py
@@ -23,6 +23,18 @@
 from autoarray.inversion.plot.mapper_plots import plot_mapper
 from autoarray.structures.arrays.uniform_2d import Array2D
 
+
+def _is_array_free_interferometer(inversion) -> bool:
+    """
+    Whether ``inversion`` is an interferometer inversion of an array-free dataset (built by
+    ``Interferometer.from_stream`` / ``from_sparse_terms``), which has no transformer and so no
+    visibility-space ``mapped_reconstructed_operated_data_dict``: its reconstructed image is read
+    from ``mapped_reconstructed_data_dict`` directly. Imaging inversions have no ``transformer``
+    attribute and are unaffected.
+    """
+    return hasattr(type(inversion), "transformer") and inversion.transformer is None
+
+
 logger = logging.getLogger(__name__)
 
 
@@ -91,6 +103,9 @@ def subplot_of_mapper(
 
     # panels 1-3: reconstructed operated data (plain, log10, + mesh grid overlay)
     def _recon_array():
+        if _is_array_free_interferometer(inversion):
+            return inversion.mapped_reconstructed_data_dict[mapper]
+
         array = inversion.mapped_reconstructed_operated_data_dict[mapper]
         from autoarray.structures.visibilities import Visibilities
 
@@ -129,7 +144,7 @@ def _recon_array():
             positions=positions,
             lines=lines,
         )
-    except (AttributeError, KeyError):
+    except (AttributeError, KeyError, exc.InversionException):
         pass
 
     # panels 4-5: source reconstruction zoomed / unzoomed
@@ -401,11 +416,14 @@ def subplot_mappings(
 
     # panel 1: reconstructed operated data
     try:
-        array = inversion.mapped_reconstructed_operated_data_dict[mapper]
-        from autoarray.structures.visibilities import Visibilities
-
-        if isinstance(array, Visibilities):
+        if _is_array_free_interferometer(inversion):
             array = inversion.mapped_reconstructed_data_dict[mapper]
+        else:
+            array = inversion.mapped_reconstructed_operated_data_dict[mapper]
+            from autoarray.structures.visibilities import Visibilities
+
+            if isinstance(array, Visibilities):
+                array = inversion.mapped_reconstructed_data_dict[mapper]
         plot_array(
             array,
             ax=axes[1],
@@ -419,7 +437,7 @@ def subplot_mappings(
             region_alpha=region_alpha,
             region_labels=region_labels,
         )
-    except (AttributeError, KeyError):
+    except (AttributeError, KeyError, exc.InversionException):
         pass
 
     pixel_values = inversion.reconstruction_dict[mapper]
diff --git a/test_autoarray/dataset/plot/test_interferometer_plotters.py b/test_autoarray/dataset/plot/test_interferometer_plotters.py
index defcd0d97..b1fb937e6 100644
--- a/test_autoarray/dataset/plot/test_interferometer_plotters.py
+++ b/test_autoarray/dataset/plot/test_interferometer_plotters.py
@@ -76,3 +76,103 @@ def test__subplots_are_output(interferometer_7, plot_path, plot_patch):
     )
 
     assert str(Path(plot_path) / "dirty_images.png") in plot_patch.paths
+
+
+def _array_free_from(dataset):
+    """
+    The array-free counterpart of an in-memory `Interferometer`: the same visibilities
+    streamed through `Interferometer.from_stream` in two chunks.
+    """
+    import numpy as np
+    import autoarray as aa
+
+    uv_wavelengths = np.asarray(dataset.uv_wavelengths)
+    data = np.asarray(dataset.data.array)
+    noise_map = np.asarray(dataset.noise_map.array)
+
+    chunks = [
+        (uv_wavelengths[k0:k1], data[k0:k1], noise_map[k0:k1])
+        for k0, k1 in ((0, 3), (3, data.shape[0]))
+    ]
+
+    return aa.Interferometer.from_stream(chunks, dataset.real_space_mask)
+
+
+def test__subplots_are_output__array_free_dataset(
+    interferometer_7, plot_path, plot_patch
+):
+    pytest.importorskip("nufftax")
+
+    dataset = _array_free_from(interferometer_7)
+
+    assert dataset.is_array_free
+
+    subplot_interferometer_dataset(
+        dataset=dataset,
+        output_path=plot_path,
+        output_format="png",
+    )
+
+    assert str(Path(plot_path) / "dataset.png") in plot_patch.paths
+
+    subplot_interferometer_dirty_images(
+        dataset=dataset,
+        output_path=plot_path,
+        output_format="png",
+    )
+
+    assert str(Path(plot_path) / "dirty_images.png") in plot_patch.paths
+
+
+def test__fits_interferometer__array_free_dataset_writes_natural_terms(
+    interferometer_7, tmp_path
+):
+    pytest.importorskip("nufftax")
+
+    import numpy as np
+    from astropy.io import fits
+    from autoarray.dataset.plot.interferometer_plots import fits_interferometer
+
+    dataset = _array_free_from(interferometer_7)
+
+    file_path = tmp_path / "dataset.fits"
+
+    fits_interferometer(dataset=dataset, file_path=file_path)
+
+    with fits.open(file_path) as hdu_list:
+        ext_names = [hdu.header.get("EXTNAME") for hdu in hdu_list]
+        dirty_image = np.asarray(hdu_list["DIRTY_IMAGE_NATURAL"].data)
+        dirty_beam = np.asarray(hdu_list["DIRTY_BEAM"].data)
+
+    assert ext_names == ["DIRTY_IMAGE_NATURAL", "DIRTY_BEAM"]
+    np.testing.assert_allclose(
+        dirty_image, np.asarray(dataset.dirty_image_natural.native), rtol=1.0e-6
+    )
+    np.testing.assert_allclose(
+        dirty_beam, np.asarray(dataset.dirty_beam.native), rtol=1.0e-6
+    )
+
+    # Separate-file mode writes nothing for the absent arrays and does not raise.
+    data_path = tmp_path / "data.fits"
+
+    fits_interferometer(dataset=dataset, data_path=data_path)
+
+    assert not data_path.exists()
+
+
+def test__fits_interferometer__in_memory_dataset_unchanged(interferometer_7, tmp_path):
+    from astropy.io import fits
+    from autoarray.dataset.plot.interferometer_plots import fits_interferometer
+
+    file_path = tmp_path / "dataset.fits"
+
+    fits_interferometer(dataset=interferometer_7, file_path=file_path)
+
+    with fits.open(file_path) as hdu_list:
+        ext_names = [hdu.header.get("EXTNAME") for hdu in hdu_list]
+
+    assert [name.lower() for name in ext_names] == [
+        "data",
+        "noise_map",
+        "uv_wavelengths",
+    ]
diff --git a/test_autoarray/fit/plot/test_fit_interferometer_plotters.py b/test_autoarray/fit/plot/test_fit_interferometer_plotters.py
index 3d811050f..3e3d4991b 100644
--- a/test_autoarray/fit/plot/test_fit_interferometer_plotters.py
+++ b/test_autoarray/fit/plot/test_fit_interferometer_plotters.py
@@ -108,3 +108,74 @@ def test__fit_sub_plots(fit_interferometer_7, plot_path, plot_patch):
     )
 
     assert str(Path(plot_path) / "fit_dirty_images.png") in plot_patch.paths
+
+
+def _array_free_from(dataset):
+    """
+    The array-free counterpart of an in-memory `Interferometer`: the same visibilities
+    streamed through `Interferometer.from_stream` in two chunks.
+    """
+    import autoarray as aa
+
+    uv_wavelengths = np.asarray(dataset.uv_wavelengths)
+    data = np.asarray(dataset.data.array)
+    noise_map = np.asarray(dataset.noise_map.array)
+
+    chunks = [
+        (uv_wavelengths[k0:k1], data[k0:k1], noise_map[k0:k1])
+        for k0, k1 in ((0, 3), (3, data.shape[0]))
+    ]
+
+    return aa.Interferometer.from_stream(chunks, dataset.real_space_mask)
+
+
+def test__fit_sub_plots__array_free_dataset(
+    interferometer_7, plot_path, plot_patch, monkeypatch
+):
+    pytest.importorskip("nufftax")
+
+    import autoarray as aa
+    from autoarray.fit.plot import fit_interferometer_plots
+
+    dataset = _array_free_from(interferometer_7)
+
+    fit = aa.m.MockFitInterferometer(dataset=dataset)
+
+    model_image = aa.Array2D(
+        values=np.ones(dataset.real_space_mask.pixels_in_mask),
+        mask=dataset.real_space_mask,
+    )
+
+    titles = []
+    plot_array = fit_interferometer_plots.plot_array
+
+    def _plot_array(array, *args, title=None, **kwargs):
+        titles.append(title)
+        return plot_array(array, *args, title=title, **kwargs)
+
+    monkeypatch.setattr(fit_interferometer_plots, "plot_array", _plot_array)
+
+    aplt.subplot_fit_interferometer(
+        fit=fit,
+        output_path=plot_path,
+        output_format="png",
+        model_image=model_image,
+    )
+
+    assert str(Path(plot_path) / "fit.png") in plot_patch.paths
+    assert titles == [
+        "Dirty Image (Natural)",
+        "Dirty Model Image (Natural)",
+        "Dirty Residual Map (Natural)",
+    ]
+
+    titles.clear()
+
+    aplt.subplot_fit_interferometer_dirty_images(
+        fit=fit,
+        output_path=plot_path,
+        output_format="png",
+    )
+
+    assert str(Path(plot_path) / "fit_dirty_images.png") in plot_patch.paths
+    assert titles == ["Dirty Image (Natural)"]
diff --git a/test_autoarray/fit/test_fit_interferometer.py b/test_autoarray/fit/test_fit_interferometer.py
index d4c17f103..cb9cf0715 100644
--- a/test_autoarray/fit/test_fit_interferometer.py
+++ b/test_autoarray/fit/test_fit_interferometer.py
@@ -516,3 +516,55 @@ def fget(self):
     assert np.isfinite(fit.log_evidence)
     assert fit.figure_of_merit == fit.log_evidence
     assert touched == []
+
+
+def test__dirty_model_image_natural_from__matches_transformer_natural_dirty_image(
+    interferometer_7,
+):
+    from autoarray.fit.fit_interferometer import dirty_model_image_natural_from
+
+    dataset = interferometer_7.apply_sparse_operator()
+
+    image = aa.Array2D(
+        values=np.random.default_rng(seed=1).normal(
+            size=dataset.real_space_mask.pixels_in_mask
+        ),
+        mask=dataset.real_space_mask,
+    )
+
+    dirty_model_image = dirty_model_image_natural_from(dataset=dataset, image=image)
+
+    # The natural dirty image of the model visibilities, built exactly as
+    # `Interferometer.dirty_image_natural` builds it from the data.
+    noise_map = dataset.noise_map.array
+    visibilities = dataset.transformer.visibilities_from(image=image).array
+    weighted = aa.Visibilities(
+        visibilities=visibilities.real * noise_map.real**-2.0
+        + 1j * visibilities.imag * noise_map.imag**-2.0
+    )
+    expected = dataset.transformer.image_from(visibilities=weighted) / float(
+        np.sum(noise_map.real**-2.0)
+    )
+
+    assert isinstance(dirty_model_image, aa.Array2D)
+    assert dirty_model_image.mask is dataset.real_space_mask
+    np.testing.assert_allclose(
+        dirty_model_image.array,
+        expected.array,
+        rtol=1.0e-10,
+        atol=1.0e-10 * np.abs(expected.array).max(),
+    )
+
+
+def test__dirty_model_image_natural_from__no_sparse_operator__raises(
+    interferometer_7,
+):
+    from autoarray.fit.fit_interferometer import dirty_model_image_natural_from
+
+    image = aa.Array2D.ones(
+        shape_native=interferometer_7.real_space_mask.shape_native,
+        pixel_scales=interferometer_7.real_space_mask.pixel_scales,
+    ).apply_mask(mask=interferometer_7.real_space_mask)
+
+    with pytest.raises(aa.exc.DatasetException, match="sparse_operator"):
+        dirty_model_image_natural_from(dataset=interferometer_7, image=image)
diff --git a/test_autoarray/inversion/plot/test_inversion_plotters.py b/test_autoarray/inversion/plot/test_inversion_plotters.py
index 735ea1a2d..a88864f49 100644
--- a/test_autoarray/inversion/plot/test_inversion_plotters.py
+++ b/test_autoarray/inversion/plot/test_inversion_plotters.py
@@ -205,3 +205,105 @@ def test__subplot_mappings__pix_indexes_bypasses_the_clump_finding(
     )
 
     assert str(Path(plot_path) / "mappings_pix_indexes_0.png") in plot_patch.paths
+
+
+def _array_free_inversion():
+    """
+    A rectangular-mesh inversion of an array-free `Interferometer` (built by
+    `from_stream`), which has no transformer and no visibility-space
+    `mapped_reconstructed_operated_data_dict`.
+    """
+    import autoarray as aa
+    from autoarray.inversion.mesh.mesh.rectangular_rtu_adapt_density import (
+        overlay_grid_from,
+    )
+
+    mask = aa.Mask2D.circular(shape_native=(10, 10), pixel_scales=0.5, radius=2.0)
+
+    rng = np.random.default_rng(seed=3)
+    uv_wavelengths = rng.normal(size=(60, 2)) * 1.0e5
+    data = rng.normal(size=60) + 1j * rng.normal(size=60)
+    sigma = rng.uniform(0.5, 2.0, size=60)
+
+    dataset = aa.Interferometer.from_stream(
+        [(uv_wavelengths, data, sigma + 1j * sigma)], mask
+    )
+
+    grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=1)
+    mesh = aa.mesh.RectangularUniform(shape=(4, 4))
+    interpolator = mesh.interpolator_from(
+        source_plane_data_grid=grid,
+        source_plane_mesh_grid=aa.Grid2DIrregular(
+            overlay_grid_from(shape_native=(4, 4), grid=grid)
+        ),
+        adapt_data=None,
+    )
+    mapper = aa.Mapper(
+        interpolator=interpolator, regularization=aa.reg.Constant(coefficient=1.0)
+    )
+
+    return aa.Inversion(dataset=dataset, linear_obj_list=[mapper])
+
+
+def _record_titles(monkeypatch):
+    from autoarray.inversion.plot import inversion_plots
+
+    titles = []
+    plot_array = inversion_plots.plot_array
+
+    def _plot_array(array, *args, title=None, **kwargs):
+        titles.append(title)
+        return plot_array(array, *args, title=title, **kwargs)
+
+    monkeypatch.setattr(inversion_plots, "plot_array", _plot_array)
+
+    return titles
+
+
+def test__subplot_of_mapper__array_free_interferometer_inversion(
+    plot_path, plot_patch, monkeypatch
+):
+    pytest.importorskip("nufftax")
+
+    inversion = _array_free_inversion()
+
+    with pytest.raises(exc.InversionException):
+        inversion.mapped_reconstructed_operated_data_dict
+
+    titles = _record_titles(monkeypatch)
+
+    aplt.subplot_of_mapper(
+        inversion=inversion,
+        mapper_index=0,
+        output_path=plot_path,
+        output_format="png",
+    )
+
+    assert str(Path(plot_path) / "inversion_0.png") in plot_patch.paths
+
+    # The reconstructed-image panels are drawn from `mapped_reconstructed_data_dict`.
+    assert "Reconstructed Image" in titles
+    assert "Reconstructed Image (log10)" in titles
+    assert "Mesh Pixel Grid Overlaid" in titles
+
+
+def test__subplot_mappings__array_free_interferometer_inversion(
+    plot_path, plot_patch, monkeypatch
+):
+    pytest.importorskip("nufftax")
+
+    inversion = _array_free_inversion()
+
+    titles = _record_titles(monkeypatch)
+
+    aplt.subplot_mappings(
+        inversion=inversion,
+        pixelization_index=0,
+        pix_indexes=[[0, 1], [8]],
+        output_path=plot_path,
+        output_filename="mappings",
+        output_format="png",
+    )
+
+    assert str(Path(plot_path) / "mappings_0.png") in plot_patch.paths
+    assert "Reconstructed Image" in titles