-
Notifications
You must be signed in to change notification settings - Fork 85
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Update redshift notebook to include unit conversion in jdaviz (#150)
* Include unit conversion from jdaviz * Fix style errors --------- Co-authored-by: Camilla Pacifici <[email protected]> Co-authored-by: Hatice Karatay <[email protected]>
- Loading branch information
3 people
authored
Oct 10, 2023
1 parent
1cc50ae
commit 9d91f65
Showing
1 changed file
with
74 additions
and
29 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -26,7 +26,8 @@ | |
| " - [Clean up the spectrum](#cleanup)\n", | ||
| " - [Run the cross correlation function](#run_crosscorr)\n", | ||
| "\n", | ||
| "**Author**: Camilla Pacifici ([email protected])" | ||
| "**Author**: Camilla Pacifici ([email protected])<br>\n", | ||
| "**Updated**: September 14, 2023" | ||
| ] | ||
| }, | ||
| { | ||
|
|
@@ -90,7 +91,7 @@ | |
| "# astropy\n", | ||
| "import astropy # again for the version number\n", | ||
| "import astropy.units as u\n", | ||
| "from astropy.io import ascii\n", | ||
| "from astropy.io import ascii, fits\n", | ||
| "from astropy.utils.data import download_file\n", | ||
| "from astropy.modeling.models import Linear1D\n", | ||
| "from astropy.nddata import StdDevUncertainty\n", | ||
|
|
@@ -109,7 +110,7 @@ | |
| "from matplotlib import pyplot as plt\n", | ||
| "\n", | ||
| "# display\n", | ||
| "from IPython.core.display import display, HTML" | ||
| "from IPython.display import display, HTML" | ||
| ] | ||
| }, | ||
| { | ||
|
|
@@ -179,6 +180,32 @@ | |
| "fn_template = download_file('https://stsci.box.com/shared/static/3rkurzwl0l79j70ddemxafhpln7ljle7.dat', cache=True)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "1415a59f", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "Jdaviz will read by default the surface brightness extension (in MJy/sr), but I prefer the flux extension (in Jy). I create a Spectrum1D object reading the file myself." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "4abd95d6", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "hdu = fits.open(f'{data_dir}/{fn}')\n", | ||
| "wave = hdu[1].data['WAVELENGTH'] * u.Unit(hdu[1].header['TUNIT1'])\n", | ||
| "flux = hdu[1].data['FLUX'] * u.Unit(hdu[1].header['TUNIT2'])\n", | ||
| "# std = hdu[1].data['FLUX_ERROR'] * u.Unit(hdu[1].header['TUNIT3']) # These are all nan. Define an artificial uncertainty\n", | ||
| "\n", | ||
| "fluxobs = Spectrum1D(spectral_axis=wave,\n", | ||
| " flux=flux,\n", | ||
| " uncertainty=StdDevUncertainty(0.1*flux))\n", | ||
| "fluxobs" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "68fa01fc", | ||
|
|
@@ -210,7 +237,8 @@ | |
| "outputs": [], | ||
| "source": [ | ||
| "# Load spectrum\n", | ||
| "viz.load_spectrum(f'{data_dir}/{fn}', data_label=\"NIRSpec\")" | ||
| "# viz.load_data(f'{data_dir}/{fn}', data_label=\"NIRSpec\") # To load the surface brightness directly from file\n", | ||
| "viz.load_data(fluxobs, data_label='NIRSpec')" | ||
| ] | ||
| }, | ||
| { | ||
|
|
@@ -258,26 +286,24 @@ | |
| "source": [ | ||
| "<a id='get_template'></a>\n", | ||
| "### Get a template and prepare it for use\n", | ||
| "The template is used for cross correlation, so it can be renormalized for convenience. The units have to match the units of the observed spectrum.\n", | ||
| "\n", | ||
| "**Developer note:** the template spectrum is actually in erg/s/cm2/A. We are taking a shortcut here since the continuum shape and normalization do not matter for the purposes of running the cross-corralation algorithm. However, it would be nice to be able to use the Unit Conversion capabilities of Jdaviz to do the appropriate unit conversion. This capability is currently being refactor. We will update this notebook when the work is completed." | ||
| "The template is used for cross correlation, so it can be renormalized for convenience. The units have to match the units of the observed spectrum. We can do the continuum subtraction in erg/(s cm^2 A) since the continuum is close to linear and then run it by Jdaviz to get the appropriate conversion." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "52a6979a", | ||
| "id": "5d6e4beb", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Define unit\n", | ||
| "spec_unit = u.MJy / u.steradian\n", | ||
| "spec_unit = u.erg / u.s / u.cm**2 / u.AA\n", | ||
| "\n", | ||
| "# Read spectrum with the ascii function\n", | ||
| "template = ascii.read(fn_template)\n", | ||
| "# Create Spectrum1D object\n", | ||
| "template = Spectrum1D(spectral_axis=template['col1']/1E4*u.um, \n", | ||
| " flux=(template['col2']/1E6)*spec_unit)" | ||
| " flux=(template['col2']/1E24)*spec_unit) # Normalize to something close to the observed spectrum" | ||
| ] | ||
| }, | ||
| { | ||
|
|
@@ -287,7 +313,7 @@ | |
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Cut to useful range - template and obs MUST overlap, so we go to 2.4um\n", | ||
| "# Cut to useful range - template and obs MUST overlap, so we go to 2.4um.\n", | ||
| "use_tmp = (template.spectral_axis.value > 0.45) & (template.spectral_axis.value < 2.4)\n", | ||
| "template_cut = Spectrum1D(spectral_axis=template.spectral_axis[use_tmp], flux=template.flux[use_tmp])" | ||
| ] | ||
|
|
@@ -331,7 +357,24 @@ | |
| "# Use fit_generic_continuum\n", | ||
| "fit_temp = fit_generic_continuum(template_forcont, model=Linear1D())\n", | ||
| "cont_temp = fit_temp(template_cut.spectral_axis)\n", | ||
| "template_sub = template_cut - cont_temp" | ||
| "template_sub = template_cut - cont_temp\n", | ||
| "\n", | ||
| "\n", | ||
| "plt.figure(figsize=[10, 6])\n", | ||
| "plt.plot(template_cut.spectral_axis, template_cut.flux)\n", | ||
| "plt.plot(template_cut.spectral_axis, cont_temp)\n", | ||
| "plt.xlabel(\"wavelength ({:latex})\".format(template_sub.spectral_axis.unit))\n", | ||
| "plt.ylabel(\"flux ({:latex})\".format(template_sub.flux.unit))\n", | ||
| "plt.title(\"Plot template and continuum\")\n", | ||
| "plt.show()" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "c0def8ed", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "The diagram shows the template spectrum and the best-fitting continuum." | ||
| ] | ||
| }, | ||
| { | ||
|
|
@@ -387,26 +430,17 @@ | |
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Read the spectrum directly as Spectrum1D\n", | ||
| "spec1d = Spectrum1D.read(f'{data_dir}/{fn}')\n", | ||
| "# The uncertainty seems to be all NaN. Maybe there was a problem with the reduction.\n", | ||
| "# Let's just assume that the uncertainty is 5% of the flux\n", | ||
| "spec1d = Spectrum1D(spectral_axis=spec1d.spectral_axis, flux=spec1d.flux, \n", | ||
| " uncertainty=StdDevUncertainty(0.05*(spec1d.flux)))\n", | ||
| "\n", | ||
| "print(spec1d)\n", | ||
| "\n", | ||
| "# Define Spectral Region\n", | ||
| "region = SpectralRegion(2.0*u.um, 3.0*u.um)\n", | ||
| "# Extract region\n", | ||
| "spec1d_cont = extract_region(spec1d, region)\n", | ||
| "spec1d_cont = extract_region(fluxobs, region)\n", | ||
| "# Run fitting function\n", | ||
| "fit_obs = fit_generic_continuum(spec1d_cont, model=Linear1D(5))\n", | ||
| "# Apply to spectral axis\n", | ||
| "cont_obs = fit_obs(spec1d.spectral_axis)\n", | ||
| "cont_obs = fit_obs(fluxobs.spectral_axis)\n", | ||
| "\n", | ||
| "# Subtract continuum\n", | ||
| "spec1d_sub = spec1d - cont_obs" | ||
| "spec1d_sub = fluxobs - cont_obs" | ||
| ] | ||
| }, | ||
| { | ||
|
|
@@ -496,7 +530,7 @@ | |
| "id": "783799c1", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We visualize the observed and template continuum-subtracted spectra in a new instance of Specviz. Hit the Home button to see the entire wavelength range." | ||
| "We visualize the observed and template continuum-subtracted spectra in a new instance of Specviz. Hit the Home button to see the entire wavelength range. The template spectrum will change unit to match the observed spectrum." | ||
| ] | ||
| }, | ||
| { | ||
|
|
@@ -512,6 +546,17 @@ | |
| "viz3.show()" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "73ef283d", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "template_newunit = viz3.get_data('template', use_display_units=True)\n", | ||
| "template_newunit" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "8195cee8", | ||
|
|
@@ -529,7 +574,7 @@ | |
| "outputs": [], | ||
| "source": [ | ||
| "# Call the function\n", | ||
| "corr, lag = correlation.template_correlate(spec1d_masked, template_sub, lag_units=u.one)\n", | ||
| "corr, lag = correlation.template_correlate(spec1d_masked, template_newunit, lag_units=u.one)\n", | ||
| "\n", | ||
| "# Plot the correlation\n", | ||
| "plt.plot(lag, corr)\n", | ||
|
|
@@ -581,8 +626,8 @@ | |
| "source": [ | ||
| "# Visualize the redshifted template and the observed spectrum\n", | ||
| "viz4 = Specviz()\n", | ||
| "viz4.load_spectrum(spec1d_masked, data_label='Observed spectrum')\n", | ||
| "viz4.load_spectrum(template_sub_z, data_label='Redshifted best template')\n", | ||
| "viz4.load_data(spec1d_masked, data_label='Observed spectrum')\n", | ||
| "viz4.load_data(template_sub_z, data_label='Redshifted best template')\n", | ||
| "viz4.show()" | ||
| ] | ||
| }, | ||
|
|
@@ -611,7 +656,7 @@ | |
| "name": "python", | ||
| "nbconvert_exporter": "python", | ||
| "pygments_lexer": "ipython3", | ||
| "version": "3.11.3" | ||
| "version": "3.9.13" | ||
| } | ||
| }, | ||
| "nbformat": 4, | ||
|
|
||