diff --git a/notebooks/NIRSpec/IFU_cube_continuum_fit/NGC4151_FeII_ContinuumFit.ipynb b/notebooks/NIRSpec/IFU_cube_continuum_fit/NGC4151_FeII_ContinuumFit.ipynb index 0bda20cc5..a03136daf 100644 --- a/notebooks/NIRSpec/IFU_cube_continuum_fit/NGC4151_FeII_ContinuumFit.ipynb +++ b/notebooks/NIRSpec/IFU_cube_continuum_fit/NGC4151_FeII_ContinuumFit.ipynb @@ -34,7 +34,7 @@ " - astropy.utils.data for accessing the data\n", " - specutils.fitting for spectral data fitting\n", " - specutils Spectrum1D for modeling emission lines\n", - " - jdaviz.app to use cubeviz in the notebook\n", + " - jdaviz to use cubeviz in the notebook\n", "\n" ] }, @@ -63,6 +63,16 @@ "from glue.core.roi import XRangeROI" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import jdaviz\n", + "print(jdaviz.__version__)" + ] + }, { "cell_type": "code", "execution_count": null, @@ -246,6 +256,15 @@ " cubeviz.load_data(fn)" ] }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "# Change spectral unit to Angstrom if needed\n", + "unit_conv = cubeviz.plugins['Unit Conversion']\n", + "unit_conv.spectral_unit = 'Angstrom'" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -320,6 +339,36 @@ "spatial_regions" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "By default, the spectra have been extracted with the `sum` function and default options in the spectral extraction plugin. For the purpose of this notebook, we would like to have spectra extracted with the `mean` function instead. We can do that from the interface like in the screenshot here below or running the following cells.
\n", + "\"Open
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "spec_ext = cubeviz.plugins['Spectral Extraction']\n", + "spec_ext.function = 'Mean'\n", + "\n", + "spec_ext.aperture = 'Subset 1'\n", + "spec_ext.add_results = 'Spectrum (Subset 1, mean)'\n", + "spec_ext.extract()\n", + "\n", + "spec_ext.aperture = 'Subset 2'\n", + "spec_ext.add_results = 'Spectrum (Subset 2, mean)'\n", + "spec_ext.extract()\n", + "\n", + "spec_ext.aperture = 'Subset 3'\n", + "spec_ext.add_results = 'Spectrum (Subset 3, mean)'\n", + "spec_ext.extract()" + ] + }, { "cell_type": "code", "execution_count": null, @@ -327,9 +376,9 @@ "outputs": [], "source": [ "# Extract spectra corresponding to the colored regions in cubeviz\n", - "spectrum1 = cubeviz.get_data(\"contents[SCI]\", spatial_subset='Subset 1', function=\"mean\") # AGN Center\n", - "spectrum2 = cubeviz.get_data(\"contents[SCI]\", spatial_subset='Subset 2', function=\"mean\") # Red shifted component\n", - "spectrum3 = cubeviz.get_data(\"contents[SCI]\", spatial_subset='Subset 3', function=\"mean\") # Blue shifted component\n", + "spectrum1 = cubeviz.get_data(\"Spectrum (Subset 1, mean)\") # AGN Center\n", + "spectrum2 = cubeviz.get_data(\"Spectrum (Subset 2, mean)\") # Red shifted component\n", + "spectrum3 = cubeviz.get_data(\"Spectrum (Subset 3, mean)\") # Blue shifted component\n", "spectrum1" ] }, @@ -415,6 +464,7 @@ "# Visualize new subsets\n", "plt.figure()\n", "plt.plot(spec_agn.spectral_axis, spec_agn.flux, color='black')\n", + "plt.title('Spectrum subset 1')\n", "plt.xlabel('wavelength')\n", "plt.ylabel('flux')" ] @@ -429,6 +479,7 @@ "plt.figure()\n", "plt.plot(spec_feii_blue.spectral_axis, spec_feii_blue.flux, color='b')\n", "plt.plot(spec_feii_red.spectral_axis, spec_feii_red.flux, color='r')\n", + "plt.title('Spectra subset 2 and 3')\n", "plt.xlabel('wavelength')\n", "plt.ylabel('flux')" ] @@ -441,7 +492,9 @@ "\n", "Open up Model Fitting Plugin. There are a number of fields to fill in and drop down menus to select from. It is important to keep in mind that the Data menu will provide only spectra to model, while the Spectral Region menu will provide only spectral region subsets to choose. In other words, you can fit the spectra in specific spectral regions. If no spectral region is selected, the entire wavelength array will be fit by the mode.
\n", "\n", - "Select Data Cube: Entire cube
\n", + "We will first fit an individual spectrum to test if our fitting parameters are good, then we will switch to fitting the full cube.
\n", + "**Single spectrum**\n", + "Select Data: Spectrum Subset 1 (mean)
\n", "Select Spectral Region: Subset 5
\n", "Model: Linear1D
\n", "ModelID: L
\n", @@ -451,7 +504,10 @@ "Model Label: LinFitCont
\n", "Click \"Fit model\", which fits the collapsed spectrum.
\n", "View the fit in the spectral viewer and confirm you are happy with it.
\n", - "Toggle \"Cube fit\".\n", + "\n", + "**Full cube**\n", + "Toggle \"Cube fit\" at the top.
\n", + "\n", "Change the name to LinFitCont_cube and click again \"Fit model\".
\n", "\n", "\"Model\n", @@ -481,7 +537,7 @@ " plugin_mf = cubeviz.plugins['Model Fitting']\n", " plugin_mf.open_in_tray()\n", " # Input the appropriate datasets\n", - " plugin_mf.spatial_subset = 'Entire Cube'\n", + " plugin_mf.dataset = 'Spectrum (sum)'\n", " plugin_mf.spectral_subset = 'Subset 5'\n", " # Input model component\n", " plugin_mf.create_model_component(model_component='Linear1D',\n", @@ -499,8 +555,10 @@ " # Open model fitting plugin\n", " plugin_mf = cubeviz.plugins['Model Fitting']\n", " plugin_mf.open_in_tray()\n", + " # Set the fitting to cube\n", + " plugin_mf.cube_fit = True\n", " # Input the appropriate datasets\n", - " plugin_mf.spatial_subset = 'Entire Cube'\n", + " plugin_mf.dataset = 'contents[SCI]'\n", " plugin_mf.spectral_subset = 'Subset 5'\n", " # Input model component\n", " plugin_mf.create_model_component(model_component='Linear1D',\n", @@ -508,7 +566,6 @@ " # Model equation gets populated automatically\n", " plugin_mf.equation = 'L2' \n", " # After we run this, we go to the GUI and check that the fit makes sense\n", - " plugin_mf.cube_fit = True\n", " plugin_mf.add_results.label = 'LinFitCont_cube'\n", " plugin_mf.calculate_fit()" ] @@ -520,6 +577,7 @@ "outputs": [], "source": [ "models = cubeviz.get_models()\n", + "#models\n", "models['LinFitCont_cube (30, 30)']" ] }, @@ -541,7 +599,11 @@ "outputs": [], "source": [ "# List available data\n", - "print(cubeviz.data_labels)" + "print(cubeviz.data_labels)\n", + "\n", + "# Get full original cube out\n", + "scidata = cubeviz.get_data(\"contents[SCI]\")\n", + "scidata" ] }, { @@ -552,20 +614,12 @@ "source": [ "# Extract SCI cube and continuum model from Cubeviz above and make a continuum subtracted cube\n", "if 'LinFitCont_cube' in cubeviz.app.data_collection:\n", - " tsci = cubeviz.get_data(\"contents[SCI]\")\n", - " tcont_psf_cube = cubeviz.get_data(\"LinFitCont_cube\")\n", + " cont_psf_cube = cubeviz.get_data(\"LinFitCont_cube\")\n", " print('Check shape of the objects')\n", - " print(tsci.shape)\n", - " print(tcont_psf_cube.shape)\n", - " print('Transpose because it inverted RA and Dec')\n", - " sciflux = tsci.flux.value.transpose(1, 0, 2)\n", - " cont_psf_cubeflux = tcont_psf_cube.flux.value.transpose(1, 0, 2)\n", - " sci = Spectrum1D(spectral_axis=tsci.spectral_axis,\n", - " flux=sciflux*tsci.flux.unit)\n", - " cont_psf_cube = Spectrum1D(spectral_axis=tcont_psf_cube.spectral_axis,\n", - " flux=cont_psf_cubeflux*tcont_psf_cube.flux.unit)\n", + " print(scidata.shape)\n", + " print(cont_psf_cube.shape)\n", " # Obtain continuum subtracted cube\n", - " sci_contsub = sci-cont_psf_cube\n", + " sci_contsub = scidata-cont_psf_cube\n", " # Save to file\n", " # sci_contsub.write('NGC4151_Hband_ContinuumSubtract.fits', format='wcs1d-fits', overwrite=True)\n", " # cont_psf_cube.write('NGC4151_Hband_ContinuumPSF.fits', format='wcs1d-fits', overwrite=True)" @@ -576,8 +630,7 @@ "metadata": {}, "source": [ "**Developer Note**:
\n", - "- the cube extracted form cubeviz has RA and Dec flipped when reinjested\n", - "- I hit a traceback if I try to save the cubes to file" + "- I hit a traceback if I try to save the cubes to file, because they do not have a proper header" ] }, { @@ -615,8 +668,8 @@ "\n", "start_time = time.time()\n", "\n", - "cont_sub_cube = np.zeros([nz, ny, nx])\n", - "cont_psf_cube = np.zeros([nz, ny, nx])\n", + "sci_contsub_np = np.zeros([nx, ny, nz])\n", + "cont_psf_cube_np = np.zeros([nx, ny, nz])\n", "\n", "for i in range(1, nx-2):\n", " for j in range(1, ny-2):\n", @@ -624,12 +677,14 @@ " cont_fit = np.polyfit(wave[continuummin:continuummax], flux1[continuummin:continuummax], 1)\n", " fitval = np.poly1d(cont_fit)\n", " continuum = fitval(wave) \n", - " cont_sub_cube[:, j, i] = flux1 - continuum\n", - " cont_psf_cube[:, j, i] = continuum \n", + " sci_contsub_np[i, j, :] = flux1 - continuum\n", + " cont_psf_cube_np[i, j, :] = continuum \n", + "\n", + "print(sci_contsub_np.shape)\n", "\n", "del newheader_cube['MODE']\n", - "fits.writeto('NGC4151_Hband_ContinuumSubtract_numpy.fits', cont_sub_cube, newheader_cube, overwrite=True)\n", - "fits.writeto('NGC4151_Hband_ContinuumPSF_numpy.fits', cont_psf_cube, newheader_cube, overwrite=True)\n", + "fits.writeto('NGC4151_Hband_ContinuumSubtract_numpy.fits', sci_contsub_np, newheader_cube, overwrite=True)\n", + "fits.writeto('NGC4151_Hband_ContinuumPSF_numpy.fits', cont_psf_cube_np, newheader_cube, overwrite=True)\n", "print('Continuum subtracted cube saved. PSF continuum cube saved.')" ] }, @@ -655,7 +710,7 @@ "For this example, we recommend setting up a 3 component gaussian model with the following inputs
\n", "Open up Model Fitting Plugin. There are a number of fields to fill in and drop down menus to select from. It is important to keep in mind that the Data menu will provide only spectra to model, while the Spectral Region menu will provide only spectral region subsets to choose. In other words, you can fit the spectra in specific spectral regions. If no spectral region is selected, the entire wavelength array will be fit by the mode.
\n", "\n", - "Data: Entire Cube
\n", + "Data: Spectrum (sum)
\n", "Spectral region: Subset 1
\n", "Model: Three different Gaussians with ModelID's set to G1, G2, and G3
\n", "Model Parameters:
\n", @@ -753,6 +808,13 @@ "print(spec)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Developer note**: currently there is no way to use the spectral extraction plugin on a created model cube. The 1D model should be calculated on the 'mean' extracted spectrum if the 'mean' model is desired." + ] + }, { "cell_type": "code", "execution_count": null, @@ -760,10 +822,17 @@ "outputs": [], "source": [ "# Get gauss model spectrum and model cube\n", - "all_spec = cubeviz2.get_data('Continuum Subtracted[FLUX]', function='mean') # AGN Center Data Cube\n", + "spec_ext2 = cubeviz2.plugins['Spectral Extraction']\n", + "spec_ext2.function = 'Mean'\n", + "spec_ext2.aperture = 'Entire Cube'\n", + "spec_ext2.add_results = 'Spectrum entire (mean)'\n", + "spec_ext2.extract()\n", + "\n", + "# This is used only as spectral axis for the models afterwards\n", + "all_spec = cubeviz2.get_data('Spectrum entire (mean)')\n", "\n", "if 'GaussAll' in alldata:\n", - " gauss_spec = cubeviz2.get_data('GaussAll', function='mean') # AGN Center Model Spec\n", + " gauss_spec = cubeviz2.get_data('GaussAll') # AGN Center Model Spec. This is 'sum', not 'mean'\n", " print('Model spectrum 1D available')\n", " print(gauss_spec)\n", "else:\n", @@ -787,7 +856,8 @@ "if gauss_cube is False:\n", " # Get both the model cube and the continuum cube not created with Cubeviz\n", " fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_cube_continuum_fit/gauss_model_cube.fits', cache=False)\n", - " gauss_cube = fits.getdata(fn)\n", + " tgauss_cube = fits.getdata(fn)\n", + " gauss_cube = tgauss_cube.transpose(2, 1, 0)\n", " print('Shape of downloaded model cube: ', gauss_cube.shape)\n", " fn_continuum = 'NGC4151_Hband_ContinuumPSF_numpy.fits'\n", " continuum_cube = fits.open(fn_continuum, memmap=False)\n", @@ -833,22 +903,41 @@ "outputs": [], "source": [ "# Overwrite gauss model with only 2 of the components of interest\n", - "gauss_cube_2component = gauss_cube*0\n", + "if 'GaussAll_cube' in alldata:\n", + " gauss_cube_2component = gauss_cube.flux * 0.\n", + " model_label = \"GaussAll_cube\"\n", + " specunit = 1.\n", + " ampunit = 1.\n", + "else:\n", + " gauss_cube_2component = gauss_cube * 0.\n", + " model_label = \"GaussAll_3d\"\n", + " specunit = u.Angstrom\n", + " ampunit = u.Unit('count')\n", "\n", - "model_label = \"GaussAll_3d\" # Change this when cube model fitting converges\n", + "print(gauss_cube_2component.shape)\n", "\n", - "nz, ny, nx = gauss_cube_2component.shape\n", + "nx, ny, nz = gauss_cube_2component.shape\n", "for i in range(0, nx-1):\n", " for j in range(0, ny-1):\n", " amp1 = params[model_label]['amplitude_0'][i][j]\n", " amp2 = params[model_label]['amplitude_2'][i][j]\n", - " m1 = params[model_label]['mean_0'][i][j]\n", - " m2 = params[model_label]['mean_2'][i][j]\n", - " stdev1 = params[model_label]['stddev_0'][i][j]\n", - " stdev2 = params[model_label]['stddev_2'][i][j]\n", - " g1 = Gaussian1D(amplitude=amp1*u.Unit('count'), mean=m1*u.m, stddev=stdev1*u.m)\n", - " g2 = Gaussian1D(amplitude=amp2*u.Unit('count'), mean=m2*u.m, stddev=stdev2*u.m)\n", - " gauss_cube_2component[:, i, j] = g1(all_spec.spectral_axis)+g2(all_spec.spectral_axis)" + " m1 = params[model_label]['mean_0'][i][j]*1E10 # Original model was in meters\n", + " m2 = params[model_label]['mean_2'][i][j]*1E10\n", + " stdev1 = params[model_label]['stddev_0'][i][j]*1E10\n", + " stdev2 = params[model_label]['stddev_2'][i][j]*1E10\n", + " g1 = Gaussian1D(amplitude=amp1*ampunit, mean=m1*specunit, stddev=stdev1*specunit)\n", + " g2 = Gaussian1D(amplitude=amp2*ampunit, mean=m2*specunit, stddev=stdev2*specunit)\n", + " gauss_cube_2component[i, j, :] = g1(all_spec.spectral_axis)+g2(all_spec.spectral_axis)\n", + "\n", + "gauss_cube_2component_spec = Spectrum1D(spectral_axis=all_spec.spectral_axis,\n", + " flux=gauss_cube_2component*ampunit)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Note:** Why is the fit not done with all components at once instead of adding the continuum component here?" ] }, { @@ -858,8 +947,8 @@ "outputs": [], "source": [ "# Add the continuum cube to the new model cube\n", - "full_model = gauss_cube_2component + continuum_data\n", - "print(continuum_data.shape)" + "full_model = gauss_cube_2component_spec + continuum_data\n", + "print(full_model.shape)" ] }, { @@ -873,10 +962,14 @@ "cube_file = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_cube_continuum_fit/NGC4151_Hband.fits'\n", "newfinalsub_header = fits.getheader(cube_file)\n", "\n", - "with fits.open(cube_file, memmap=False) as original_cube:\n", - " original_data = original_cube['SCI'].data\n", - " final_sub_cube = original_data - full_model\n", - " print(original_data.shape)" + "# Cube from cubeviz is scidata\n", + "final_sub_cube = scidata.flux.value - full_model.flux.value\n", + "\n", + "final_sub_cube_units = Spectrum1D(spectral_axis=scidata.spectral_axis,\n", + " flux=final_sub_cube*ampunit)\n", + "\n", + "print(final_sub_cube_units.shape)\n", + "print(scidata.shape)" ] }, { @@ -898,15 +991,20 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Developer note:** Fits writeto gives an error. Making cell raw for now." + ] + }, + { + "cell_type": "raw", "metadata": {}, - "outputs": [], "source": [ "del newfinalsub_header['MODE']\n", "\n", "fits.writeto('NGC4151_Hband_ContinuumandBrackettModel.fits', full_model, newfull_header, overwrite=True)\n", - "fits.writeto('NGC4151_Hband_FinalSubtract.fits', final_sub_cube, newfinalsub_header, overwrite=True)\n", + "fits.writeto('NGC4151_Hband_FinalSubtract.fits', final_sub_cube_units, newfinalsub_header, overwrite=True)\n", "print('Continuum subtracted cube saved. PSF continuum cube saved.')" ] }, @@ -920,10 +1018,10 @@ "plt.figure()\n", "plt.xlim([16200, 16650])\n", "plt.ylim([600, 900])\n", - "plt.plot(all_spec.spectral_axis, continuum_data[:, 30, 30], label='Continuum')\n", - "plt.plot(all_spec.spectral_axis, original_data[:, 30, 30], label='Original Data')\n", - "plt.plot(all_spec.spectral_axis, full_model[:, 30, 30], label='2 Component Model')\n", - "plt.plot(all_spec.spectral_axis, final_sub_cube[:, 30, 30]+700, label='Model Subtraction+Offset')\n", + "plt.plot(all_spec.spectral_axis, continuum_data[30, 30, :], label='Continuum')\n", + "plt.plot(all_spec.spectral_axis, scidata.flux[30, 30, :], label='Original Data')\n", + "plt.plot(all_spec.spectral_axis, full_model.flux[30, 30, :], label='2 Component Model')\n", + "plt.plot(all_spec.spectral_axis, final_sub_cube_units.flux[30, 30, :]+700*ampunit, label='Model Subtraction+Offset')\n", "plt.legend()\n", "plt.xlabel('wavelength')\n", "plt.ylabel('flux')\n", @@ -947,7 +1045,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.7" + "version": "3.11.11" }, "toc": { "base_numbering": 1, diff --git a/notebooks/NIRSpec/IFU_cube_continuum_fit/requirements.txt b/notebooks/NIRSpec/IFU_cube_continuum_fit/requirements.txt index 78cb78d6b..1de3f946a 100644 --- a/notebooks/NIRSpec/IFU_cube_continuum_fit/requirements.txt +++ b/notebooks/NIRSpec/IFU_cube_continuum_fit/requirements.txt @@ -1 +1 @@ -jdaviz>=3.8.1 +git+https://github.com/spacetelescope/jdaviz.git