diff --git a/jdaviz/app.py b/jdaviz/app.py
index 4376fbf055..3f9a9d4e0e 100644
--- a/jdaviz/app.py
+++ b/jdaviz/app.py
@@ -118,6 +118,9 @@ def to_unit(self, data, cid, values, original_units, target_units):
                 else:
                     spectral_values = spec.spectral_axis
 
+                # By default include spectral density equivalency
+                eqv = u.spectral_density(spectral_values)
+
                 # Ensure a spectrum passed through Spectral Extraction plugin
                 if '_pixel_scale_factor' in spec.meta:
 
@@ -128,21 +131,17 @@ def to_unit(self, data, cid, values, original_units, target_units):
                     if (u.sr in u.Unit(original_units).bases) and \
                            (u.sr not in u.Unit(target_units).bases):
                         # Surface Brightness -> Flux
-                        eqv = [(u.MJy / u.sr,
-                                u.MJy,
-                                lambda x: (x * spec.meta['_pixel_scale_factor']),
-                                lambda x: x)]
+                        eqv += [(u.MJy / u.sr,
+                                 u.MJy,
+                                 lambda x: (x * spec.meta['_pixel_scale_factor']),
+                                 lambda x: x)]
                     elif (u.sr not in u.Unit(original_units).bases) and \
                             (u.sr in u.Unit(target_units).bases):
                         # Flux -> Surface Brightness
-                        eqv = [(u.MJy,
-                                u.MJy / u.sr,
-                                lambda x: (x / spec.meta['_pixel_scale_factor']),
-                                lambda x: x)]
-                    else:
-                        eqv = u.spectral_density(spectral_values)
-                else:
-                    eqv = u.spectral_density(spectral_values)
+                        eqv += [(u.MJy,
+                                 u.MJy / u.sr,
+                                 lambda x: (x / spec.meta['_pixel_scale_factor']),
+                                 lambda x: x)]
 
         else:  # spectral axis
             eqv = u.spectral() + u.pixel_scale(1*u.pix)
diff --git a/jdaviz/tests/test_app.py b/jdaviz/tests/test_app.py
index 33d6ad4c0f..bca0a62e9b 100644
--- a/jdaviz/tests/test_app.py
+++ b/jdaviz/tests/test_app.py
@@ -272,3 +272,16 @@ def test_to_unit(cubeviz_helper):
                        ([1, 2] * original_units)
                        .to_value(target_units,
                                  equivalencies=u.spectral_density(cube.spectral_axis[0])))
+
+    # Change from Fnu to Flam/sr (with values shape matching spectral axis)
+
+    values = np.ones(3001)
+    original_units = u.MJy
+    target_units = u.erg / u.cm**2 / u.s / u.AA / u.sr
+
+    new_values = uc.to_unit(cubeviz_helper, data, cid, values, original_units, target_units)
+
+    assert np.allclose(new_values,
+                       (values * original_units)
+                       .to_value(target_units,
+                                 equivalencies=u.spectral_density(cube.spectral_axis)))