Skip to content

Commit

Permalink
Finish fixing
Browse files Browse the repository at this point in the history
  • Loading branch information
Camilla Pacifici committed Dec 13, 2024
1 parent 50b64a9 commit 36281c2
Show file tree
Hide file tree
Showing 2 changed files with 97 additions and 59 deletions.
154 changes: 96 additions & 58 deletions notebooks/NIRSpec/IFU_cube_continuum_fit/NGC4151_FeII_ContinuumFit.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -257,12 +257,10 @@
]
},
{
"cell_type": "code",
"execution_count": null,
"cell_type": "raw",
"metadata": {},
"outputs": [],
"source": [
"# Change spectral unit to Angstrom\n",
"# Change spectral unit to Angstrom if needed\n",
"unit_conv = cubeviz.plugins['Unit Conversion']\n",
"unit_conv.spectral_unit = 'Angstrom'"
]
Expand Down Expand Up @@ -345,7 +343,7 @@
"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",
"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. <br>\n",
"<img src='./cubeviz_spectral_extraction_plugin.png' alt=\"Open the spectral extraction plugin, select a spatial subset, and the extraction function mean.\" width=\"500\"/><br>"
]
},
Expand Down Expand Up @@ -560,7 +558,7 @@
" # Set the fitting to cube\n",
" plugin_mf.cube_fit = True\n",
" # Input the appropriate datasets\n",
" plugin_mf.dataset = 'contents[FLUX]'\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",
Expand All @@ -579,8 +577,8 @@
"outputs": [],
"source": [
"models = cubeviz.get_models()\n",
"models\n",
"#models['LinFitCont_cube (30, 30)']"
"#models\n",
"models['LinFitCont_cube (30, 30)']"
]
},
{
Expand All @@ -601,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"
]
},
{
Expand All @@ -612,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)"
Expand All @@ -636,8 +630,7 @@
"metadata": {},
"source": [
"**Developer Note**:<br>\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"
]
},
{
Expand Down Expand Up @@ -675,21 +668,23 @@
"\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",
" flux1 = cube[:, j, i] \n",
" 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.')"
]
},
Expand All @@ -715,7 +710,7 @@
"For this example, we recommend setting up a 3 component gaussian model with the following inputs<br />\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. <br />\n",
"\n",
"Data: Entire Cube<br />\n",
"Data: Spectrum (sum)<br />\n",
"Spectral region: Subset 1<br />\n",
"Model: Three different Gaussians with ModelID's set to G1, G2, and G3<br />\n",
"Model Parameters: <br />\n",
Expand Down Expand Up @@ -813,17 +808,31 @@
"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,
"metadata": {},
"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",
Expand All @@ -847,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",
Expand Down Expand Up @@ -893,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?"
]
},
{
Expand All @@ -918,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)"
]
},
{
Expand All @@ -933,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)"
]
},
{
Expand All @@ -958,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.')"
]
},
Expand All @@ -980,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",
Expand All @@ -1007,7 +1045,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.10"
"version": "3.11.11"
},
"toc": {
"base_numbering": 1,
Expand Down
2 changes: 1 addition & 1 deletion notebooks/NIRSpec/IFU_cube_continuum_fit/requirements.txt
View file
Original file line number Diff line number Diff line change
@@ -1 +1 @@
jdaviz>=4.0.0
git+https://github.com/spacetelescope/jdaviz.git

0 comments on commit 36281c2

Please sign in to comment.