diff --git a/notebooks/mos_spectroscopy_advanced/MOSspec_advanced.ipynb b/notebooks/mos_spectroscopy_advanced/MOSspec_advanced.ipynb
index 6a340854d..eacfe8a06 100644
--- a/notebooks/mos_spectroscopy_advanced/MOSspec_advanced.ipynb
+++ b/notebooks/mos_spectroscopy_advanced/MOSspec_advanced.ipynb
@@ -85,7 +85,7 @@
     "\n",
     "# jdaviz\n",
     "import jdaviz\n",
-    "from jdaviz import Mosviz, Specviz2d, Specviz\n",
+    "from jdaviz import Mosviz, Specviz2d, Specviz  # noqa\n",
     "\n",
     "# glue\n",
     "from glue.core.roi import XRangeROI\n",
@@ -112,10 +112,10 @@
    "source": [
     "## Check versions. Should be:\n",
     "\n",
-    "Numpy:  1.24.2<br>\n",
-    "Astropy:  5.2.1<br>\n",
-    "Specutils:  1.9.1<br>\n",
-    "Jdaviz: 3.4.0"
+    "Numpy:  1.25.2<br>\n",
+    "Astropy:  5.3.2<br>\n",
+    "Specutils:  1.11.0<br>\n",
+    "Jdaviz: 3.7"
    ]
   },
   {
@@ -181,6 +181,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# pathtodata = Path('/home/shared/preloaded-fits/jdaviz_data/mos_ceers_data')\n",
     "pathtodata = Path('./data')"
    ]
   },
@@ -199,18 +200,21 @@
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
+   "cell_type": "raw",
+   "metadata": {},
    "source": [
     "mosviz = Mosviz()\n",
     "mosviz.load_data(directory=pathtodata, instrument=\"nirspec\")\n",
     "mosviz.show()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<img src='./mosviz_screenshot.png' alt=\"Mosviz with its 2D spectrum, 1D spectrum, and table viewer.\" width=\"900\"/>"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -247,9 +251,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "specviz2d = Specviz2d()\n",
@@ -266,13 +268,6 @@
     "This is being worked on."
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Specviz2d cannot yet work with the WCS of the 2D spectrum, so we use this workaround to get the spectral axis from the 1D spectrum extracted by the pipeline (the x1d file)."
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -280,7 +275,7 @@
    "outputs": [],
    "source": [
     "# Get out extracted spectrum from Specviz2d\n",
-    "spectra = specviz2d.app.get_data_from_viewer('spectrum-viewer', 'Spectrum 1D')\n",
+    "spectra = specviz2d.get_data('Spectrum 1D')\n",
     "spectra"
    ]
   },
@@ -290,23 +285,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Get wavelength array from x1d spectrum\n",
-    "hdu = fits.open(file1d)\n",
-    "hdu[1].data['WAVELENGTH']"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Combine to get the extracted spectrum with WCS\n",
     "# Include some fake uncertainty for now\n",
-    "spec1d_wcs = Spectrum1D(spectral_axis=hdu[1].data['WAVELENGTH']*u.um,\n",
-    "                        flux=spectra.flux,\n",
-    "                        uncertainty=StdDevUncertainty(0.1*spectra.flux))\n",
-    "spec1d_wcs"
+    "spec1d = Spectrum1D(spectral_axis=spectra.spectral_axis,\n",
+    "                    flux=spectra.flux,\n",
+    "                    uncertainty=StdDevUncertainty((np.zeros(len(spectra.flux)) + 1E-13) * spectra.unit))\n",
+    "spec1d"
    ]
   },
   {
@@ -317,7 +300,7 @@
    "source": [
     "# And open in Specviz\n",
     "specviz = Specviz()\n",
-    "specviz.load_data(spec1d_wcs, data_label='spec1d calibrated')\n",
+    "specviz.load_data(spec1d, data_label='spec1d calibrated')\n",
     "specviz.show()"
    ]
   },
@@ -325,7 +308,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "There are still some artifacts, but we can select a subset masking the artifacts and get out a spectrum without unwanted spikes. We can do so using the tool to select a subset with the \"add\" option (in the top bar) to select multiple regions as part of a single subset."
+    "There are still some artifacts in the data, but we can select a subset masking the artifacts and get out a spectrum without unwanted spikes. We can do so using the tool to select a subset with the \"add\" option (in the top bar) to select multiple regions as part of a single subset."
    ]
   },
   {
@@ -336,7 +319,7 @@
    "source": [
     "# Create a subset in the area of interest if it has not been created manually\n",
     "try:\n",
-    "    region1 = specviz.get_spectra(data_label='spec1d calibrated', subset_to_apply='Subset 1')\n",
+    "    region1 = specviz.get_data(data_label='spec1d calibrated', spectral_subset='Subset 1')\n",
     "    print(region1)\n",
     "    region1_exists = True\n",
     "except Exception:\n",
@@ -357,8 +340,8 @@
    "outputs": [],
    "source": [
     "# Get spectrum out with mask\n",
-    "spec1d_wcs_region = specviz.get_spectral_regions()\n",
-    "spec1d_wcs_masked = extract_region(spec1d_wcs, spec1d_wcs_region['Subset 1'], return_single_spectrum=True)"
+    "spec1d_region = specviz.get_spectral_regions()\n",
+    "spec1d_masked = extract_region(spec1d, spec1d_region['Subset 1'], return_single_spectrum=True)"
    ]
   },
   {
@@ -368,7 +351,7 @@
    "outputs": [],
    "source": [
     "# Load in specviz\n",
-    "specviz.load_data(spec1d_wcs_masked, data_label='spec1d masked')"
+    "specviz.load_data(spec1d_masked, data_label='spec1d masked')"
    ]
   },
   {
@@ -378,8 +361,8 @@
    "outputs": [],
    "source": [
     "# Write the extracted spectrum to a fits file\n",
-    "file_extracted = Path('./extracted _spectrum.fits')\n",
-    "spec1d_wcs_masked.write(file_extracted, overwrite=True)"
+    "file_extracted = Path('./extracted_spectrum.fits')\n",
+    "spec1d_masked.write(file_extracted, overwrite=True)"
    ]
   },
   {
@@ -445,7 +428,7 @@
    "source": [
     "# Region with some continuum for gaussian fit\n",
     "sv.toolbar_active_subset.selected = []\n",
-    "sv.apply_roi(XRangeROI(1.11, 1.15))  "
+    "sv.apply_roi(XRangeROI(1.05, 1.25))  "
    ]
   },
   {
@@ -524,7 +507,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Open line list plugin. Select the SDSS IV line list. Load the Oxygen lines and Hb. Go back to line analysis plugin and associate the Oxygen II line with the line we just analyzed."
+    "Open line list plugin. Select the SDSS IV line list. Load the Oxygen II lines and Hb. Go back to line analysis plugin and associate the Oxygen II line with the line we just analyzed."
    ]
   },
   {
@@ -572,6 +555,16 @@
     "plugin_mf.get_model_component('G')"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plugin_mf.set_model_component('G', 'stddev', 0.001)\n",
+    "plugin_mf.set_model_component('G', 'mean', 1.128)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -579,7 +572,7 @@
    "outputs": [],
    "source": [
     "# Model equation gets populated automatically\n",
-    "plugin_mf.equation"
+    "plugin_mf.equation = 'P2+G'"
    ]
   },
   {
@@ -622,8 +615,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cont_spec1d = fit_generic_continuum(spec1d_wcs_masked)\n",
-    "cont_fit = cont_spec1d(spec1d_wcs_masked.spectral_axis)"
+    "cont_spec1d = fit_generic_continuum(spec1d_masked)\n",
+    "cont_fit = cont_spec1d(spec1d_masked.spectral_axis)"
    ]
   },
   {
@@ -633,18 +626,18 @@
    "outputs": [],
    "source": [
     "plt.figure(figsize=[10, 6])\n",
-    "plt.plot(spec1d_wcs_masked.spectral_axis, spec1d_wcs_masked.flux, label=\"data\")\n",
-    "plt.plot(spec1d_wcs_masked.spectral_axis, cont_fit, label=\"modeled continuum\")\n",
-    "plt.xlabel(\"wavelength ({:latex})\".format(spec1d_wcs_masked.spectral_axis.unit))\n",
-    "plt.ylabel(\"flux ({:latex})\".format(spec1d_wcs_masked.flux.unit))\n",
+    "plt.plot(spec1d_masked.spectral_axis, spec1d_masked.flux, label=\"data\")\n",
+    "plt.plot(spec1d_masked.spectral_axis, cont_fit, label=\"modeled continuum\")\n",
+    "plt.xlabel(\"wavelength ({:latex})\".format(spec1d_masked.spectral_axis.unit))\n",
+    "plt.ylabel(\"flux ({:latex})\".format(spec1d_masked.flux.unit))\n",
     "plt.legend()\n",
     "plt.title(\"Observed spectrum and fitted continuum\")\n",
     "plt.show()\n",
     "\n",
     "plt.figure(figsize=[10, 6])\n",
-    "plt.plot(spec1d_wcs_masked.spectral_axis, spec1d_wcs_masked.uncertainty.array, label=\"data\")\n",
-    "plt.xlabel(\"wavelength ({:latex})\".format(spec1d_wcs_masked.spectral_axis.unit))\n",
-    "plt.ylabel(\"uncertainty ({:latex})\".format(spec1d_wcs_masked.uncertainty.unit))\n",
+    "plt.plot(spec1d_masked.spectral_axis, spec1d_masked.uncertainty.array, label=\"data\")\n",
+    "plt.xlabel(\"wavelength ({:latex})\".format(spec1d_masked.spectral_axis.unit))\n",
+    "plt.ylabel(\"uncertainty ({:latex})\".format(spec1d_masked.uncertainty.unit))\n",
     "plt.legend()\n",
     "plt.title(\"Uncertianty of observed spectrum\")\n",
     "plt.show()"
@@ -664,7 +657,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "spec1d_sub = spec1d_wcs_masked - cont_fit\n",
+    "spec1d_sub = spec1d_masked - cont_fit\n",
     "spec1d_sub"
    ]
   },
@@ -874,7 +867,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "spec1d_norm = spec1d_wcs_masked / cont_fit"
+    "spec1d_norm = spec1d_masked / cont_fit"
    ]
   },
   {
@@ -941,7 +934,7 @@
    "outputs": [],
    "source": [
     "# Redshift taken from the Specviz analysis\n",
-    "zz = 2.0269\n",
+    "zz = 2.0256\n",
     "\n",
     "f_lamb_units = u.erg / u.s / (u.cm**2) / u.AA\n",
     "\n",
@@ -961,12 +954,12 @@
    "outputs": [],
    "source": [
     "# Change the units of the observed spectrum to match the template\n",
-    "spec1d_wcs_masked_flamb = spec1d_wcs_masked.new_flux_unit(f_lamb_units)\n",
+    "spec1d_masked_flamb = spec1d_masked.new_flux_unit(f_lamb_units)\n",
     "# The new_flux_unit function does not change the uncertainty and specviz complains that there is a mismatch\n",
-    "# so we re-add the uncertainty as a percentage of the new flux\n",
-    "spec1d_wcs_masked_flamb_unc = Spectrum1D(spectral_axis=spec1d_wcs_masked_flamb.spectral_axis,\n",
-    "                                         flux=spec1d_wcs_masked_flamb.flux,\n",
-    "                                         uncertainty=StdDevUncertainty(0.1*spec1d_wcs_masked_flamb.flux))"
+    "# so we re-add the uncertainty like we did a few cells above\n",
+    "spec1d_masked_flamb_unc = Spectrum1D(spectral_axis=spec1d_masked_flamb.spectral_axis,\n",
+    "                                     flux=spec1d_masked_flamb.flux,\n",
+    "                                     uncertainty=StdDevUncertainty((np.zeros(len(spec1d_masked_flamb.flux)) + 1E-20) * spec1d_masked_flamb.unit))"
    ]
   },
   {
@@ -976,18 +969,18 @@
    "outputs": [],
    "source": [
     "# Take a look at the observed spectrum and one of the templates at the correct redshift\n",
-    "mean_obs = np.mean(spec1d_wcs_masked_flamb_unc.flux)\n",
+    "mean_obs = np.mean(spec1d_masked_flamb_unc.flux)\n",
     "mean_temp = np.mean(templatelist[0].flux)\n",
     "temp_for_plot = Spectrum1D(spectral_axis=templatelist[0].spectral_axis * (1.+zz),\n",
     "                           flux=templatelist[0].flux*mean_obs/mean_temp)\n",
     "\n",
     "plt.figure(figsize=[10, 6])\n",
-    "plt.plot(spec1d_wcs_masked_flamb_unc.spectral_axis, spec1d_wcs_masked_flamb_unc.flux, label='data')\n",
+    "plt.plot(spec1d_masked_flamb_unc.spectral_axis, spec1d_masked_flamb_unc.flux, label='data')\n",
     "plt.plot(temp_for_plot.spectral_axis, temp_for_plot.flux, label='model', alpha=0.6)\n",
-    "plt.xlabel(\"wavelength ({:latex})\".format(spec1d_wcs_masked_flamb_unc.spectral_axis.unit))\n",
+    "plt.xlabel(\"wavelength ({:latex})\".format(spec1d_masked_flamb_unc.spectral_axis.unit))\n",
     "plt.ylabel(\"flux (normalized)\")\n",
     "plt.xlim(1.1, 1.7)\n",
-    "plt.ylim(0, 4e-18)\n",
+    "plt.ylim(0, 2e-18)\n",
     "plt.legend()\n",
     "plt.title(\"Observed spectrum compared to one template at correct redshift\")\n",
     "plt.show()"
@@ -999,21 +992,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "tm_results = template_comparison.template_match(observed_spectrum=spec1d_wcs_masked_flamb_unc, \n",
+    "tm_results = template_comparison.template_match(observed_spectrum=spec1d_masked_flamb_unc, \n",
     "                                                spectral_templates=templatelist, \n",
     "                                                resample_method=\"flux_conserving\", \n",
     "                                                redshift=zz)\n",
     "tm_results[0]"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "**Developer note**<br>\n",
-    "The fit returns nan and I do not understand why. Maybe the observed spectrum is not properly flux calibrated, but there should still be a best chi^2 from the pile of templates."
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -1021,16 +1006,37 @@
    "outputs": [],
    "source": [
     "plt.figure(figsize=[10, 6])\n",
-    "plt.plot(spec1d_wcs_masked_flamb_unc.spectral_axis, spec1d_wcs_masked_flamb_unc.flux, label=\"data\")\n",
+    "plt.plot(spec1d_masked_flamb_unc.spectral_axis, spec1d_masked_flamb_unc.flux, label=\"data\")\n",
     "plt.plot(tm_results[0].spectral_axis, tm_results[0].flux, color='r', alpha=0.5, label='model')\n",
     "plt.xlim(1.0, 1.7)\n",
-    "plt.xlabel(\"wavelength ({:latex})\".format(spec1d_wcs_masked_flamb_unc.spectral_axis.unit))\n",
-    "plt.ylabel(\"flux ({:latex})\".format(spec1d_wcs_masked_flamb_unc.flux.unit))\n",
+    "plt.ylim(0, 5e-19)\n",
+    "plt.xlabel(\"wavelength ({:latex})\".format(spec1d_masked_flamb_unc.spectral_axis.unit))\n",
+    "plt.ylabel(\"flux ({:latex})\".format(spec1d_masked_flamb_unc.flux.unit))\n",
     "plt.legend()\n",
     "plt.title(\"Observed spectrum and best-fitting model template\")\n",
     "plt.show()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## New instance of Specviz with spectrum and template\n",
+    "Passing the spectra with different (but compatible) units. Specviz adopts the first and converts the second spectrum appropriately."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "specviz_2 = Specviz()\n",
+    "specviz_2.load_data(spec1d_masked, data_label='observed') # This is in MJy\n",
+    "specviz_2.load_data(tm_results[0], data_label='model') # This is in erg/(s cm^2 A)\n",
+    "specviz_2.show()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -1042,7 +1048,8 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Notebook created by Camilla Pacifici (cpacifici@stsci.edu)"
+    "Notebook created by Camilla Pacifici (cpacifici@stsci.edu)<br>\n",
+    "Updated on September 14, 2023"
    ]
   }
  ],
@@ -1062,9 +1069,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.9"
+   "version": "3.9.13"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }
diff --git a/notebooks/mos_spectroscopy_advanced/mosviz_screenshot.png b/notebooks/mos_spectroscopy_advanced/mosviz_screenshot.png
new file mode 100644
index 000000000..85a31545e
Binary files /dev/null and b/notebooks/mos_spectroscopy_advanced/mosviz_screenshot.png differ