diff --git a/CHANGES.rst b/CHANGES.rst
index 9a60dc27c8..c4cf256166 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -21,6 +21,10 @@ Mosviz
 Specviz
 ^^^^^^^
 
+- Spectrum that has incompatible flux unit with what is already loaded
+  will no longer be loaded as ghost spectrum. It will now be rejected
+  with an error message on the snackbar. [#2485]
+
 Specviz2d
 ^^^^^^^^^
 
diff --git a/jdaviz/app.py b/jdaviz/app.py
index ba44417cf3..7b30de92d0 100644
--- a/jdaviz/app.py
+++ b/jdaviz/app.py
@@ -1930,9 +1930,14 @@ def set_data_visibility(self, viewer_reference, data_label, visible=True, replac
                     f"No data item found with label '{data_label}'. Label must be one "
                     "of:\n\t" + "\n\t".join(dc_labels))
 
-            [data] = [x for x in self.data_collection if x.label == data_label]
+            data = self.data_collection[data_label]
 
             viewer.add_data(data, percentile=95, color=viewer.color_cycler())
+
+            # Specviz removes the data from collection in viewer.py if flux unit incompatible.
+            if data_label not in self.data_collection:
+                return
+
             viewer.set_plot_axes()
 
             add_data_message = AddDataMessage(data, viewer,
diff --git a/jdaviz/configs/specviz/plugins/viewers.py b/jdaviz/configs/specviz/plugins/viewers.py
index 153b12a907..f41fb50aa6 100644
--- a/jdaviz/configs/specviz/plugins/viewers.py
+++ b/jdaviz/configs/specviz/plugins/viewers.py
@@ -4,17 +4,19 @@
 from astropy import table
 from astropy import units as u
 from astropy.nddata import StdDevUncertainty, VarianceUncertainty, InverseVariance
+from echo import delay_callback
+from glue.config import data_translator
 from glue.core import BaseData
-from glue_astronomy.spectral_coordinates import SpectralCoordinates
+from glue.core.exceptions import IncompatibleAttribute
+from glue.core.units import UnitConverter
 from glue.core.subset import Subset
 from glue.core.subset_group import GroupedSubset
-from glue.config import data_translator
+from glue_astronomy.spectral_coordinates import SpectralCoordinates
 from glue_jupyter.bqplot.profile import BqplotProfileView
-from glue.core.exceptions import IncompatibleAttribute
 from matplotlib.colors import cnames
 from specutils import Spectrum1D
 
-from jdaviz.core.events import SpectralMarksChangedMessage, LineIdentifyMessage
+from jdaviz.core.events import SpectralMarksChangedMessage, LineIdentifyMessage, SnackbarMessage
 from jdaviz.core.registries import viewer_registry
 from jdaviz.core.marks import SpectralLine, LineUncertainties, ScatterMask, OffscreenLinesMarks
 from jdaviz.core.linelists import load_preset_linelist, get_available_linelists
@@ -24,6 +26,7 @@
 
 __all__ = ['SpecvizProfileView']
 
+uc = UnitConverter()
 
 uncertainty_str_to_cls_mapping = {
     "std": StdDevUncertainty,
@@ -366,6 +369,19 @@ def add_data(self, data, color=None, alpha=None, **layer_state):
         if len(self.layers) == 0:
             reset_plot_axes = True
         else:
+            # Check if the new data flux unit is actually compatible since flux not linked.
+            try:
+                uc.to_unit(data, data.find_component_id("flux"), [1, 1],
+                           u.Unit(self.state.y_display_unit))  # Error if incompatible
+            except Exception as err:
+                # Raising exception here introduces a dirty state that messes up next load_data
+                # but not raising exception also causes weird behavior unless we remove the data
+                # completely.
+                self.session.hub.broadcast(SnackbarMessage(
+                    f"Failed to load {data.label}, so removed it: {repr(err)}",
+                    sender=self, color='error'))
+                self.jdaviz_app.data_collection.remove(data)
+                return False
             reset_plot_axes = False
 
         # The base class handles the plotting of the main
@@ -375,8 +391,9 @@ def add_data(self, data, color=None, alpha=None, **layer_state):
         if reset_plot_axes:
             x_units = data.get_component(self.state.x_att.label).units
             y_units = data.get_component("flux").units
-            self.state.x_display_unit = x_units if len(x_units) else None
-            self.state.y_display_unit = y_units if len(y_units) else None
+            with delay_callback(self.state, "x_display_unit", "y_display_unit"):
+                self.state.x_display_unit = x_units if len(x_units) else None
+                self.state.y_display_unit = y_units if len(y_units) else None
             self.set_plot_axes()
 
         self._plot_uncertainties()
@@ -541,22 +558,23 @@ def set_plot_axes(self):
         else:
             spectral_axis_unit_type = str(x_unit.physical_type).title()
 
-        self.figure.axes[0].label = f"{spectral_axis_unit_type} [{self.state.x_display_unit}]"
-        self.figure.axes[1].label = f"{flux_unit_type} [{self.state.y_display_unit}]"
-
-        # Make it so axis labels are not covering tick numbers.
-        self.figure.fig_margin["left"] = 95
-        self.figure.fig_margin["bottom"] = 60
-        self.figure.send_state('fig_margin')  # Force update
-        self.figure.axes[0].label_offset = "40"
-        self.figure.axes[1].label_offset = "-70"
-        # NOTE: with tick_style changed below, the default responsive ticks in bqplot result
-        # in overlapping tick labels.  For now we'll hardcode at 8, but this could be removed
-        # (default to None) if/when bqplot auto ticks react to styling options.
-        self.figure.axes[1].num_ticks = 8
-
-        # Set Y-axis to scientific notation
-        self.figure.axes[1].tick_format = '0.1e'
-
-        for i in (0, 1):
-            self.figure.axes[i].tick_style = {'font-size': 15, 'font-weight': 600}
+        with self.figure.hold_sync():
+            self.figure.axes[0].label = f"{spectral_axis_unit_type} [{self.state.x_display_unit}]"
+            self.figure.axes[1].label = f"{flux_unit_type} [{self.state.y_display_unit}]"
+
+            # Make it so axis labels are not covering tick numbers.
+            self.figure.fig_margin["left"] = 95
+            self.figure.fig_margin["bottom"] = 60
+            self.figure.send_state('fig_margin')  # Force update
+            self.figure.axes[0].label_offset = "40"
+            self.figure.axes[1].label_offset = "-70"
+            # NOTE: with tick_style changed below, the default responsive ticks in bqplot result
+            # in overlapping tick labels.  For now we'll hardcode at 8, but this could be removed
+            # (default to None) if/when bqplot auto ticks react to styling options.
+            self.figure.axes[1].num_ticks = 8
+
+            # Set Y-axis to scientific notation
+            self.figure.axes[1].tick_format = '0.1e'
+
+            for i in (0, 1):
+                self.figure.axes[i].tick_style = {'font-size': 15, 'font-weight': 600}
diff --git a/jdaviz/configs/specviz/tests/test_helper.py b/jdaviz/configs/specviz/tests/test_helper.py
index 68ae8870c6..8ee076a0a2 100644
--- a/jdaviz/configs/specviz/tests/test_helper.py
+++ b/jdaviz/configs/specviz/tests/test_helper.py
@@ -443,3 +443,18 @@ def test_spectra_partial_overlap(specviz_helper):
                                          'Wave 7.02222e+03 Angstrom (2 pix)',
                                          'Flux 6.00000e+01 nJy')
     assert label_mouseover.icon == 'b'
+
+
+def test_spectra_incompatible_flux(specviz_helper):
+    """https://github.com/spacetelescope/jdaviz/issues/2459"""
+    wav = [1.1, 1.2, 1.3] * u.um
+    sp1 = Spectrum1D(flux=[1, 1.1, 1] * (u.MJy / u.sr), spectral_axis=wav)
+    sp2 = Spectrum1D(flux=[1, 1, 1.1] * (u.MJy), spectral_axis=wav)
+    flux3 = ([1, 1.1, 1] * u.MJy).to(u.erg / u.s / u.cm / u.cm / u.AA, u.spectral_density(wav))
+    sp3 = Spectrum1D(flux=flux3, spectral_axis=wav)
+
+    specviz_helper.load_data(sp2, data_label="2")  # OK
+    specviz_helper.load_data(sp1, data_label="1")  # Not OK
+    specviz_helper.load_data(sp3, data_label="3")  # OK
+
+    assert specviz_helper.app.data_collection.labels == ["2", "3"]