From 12f839e9e2a0ed6d15b233d5c3edad0ffed62fd9 Mon Sep 17 00:00:00 2001
From: Andreas Faisst <anfaisst@gmail.com>
Date: Fri, 14 Jun 2024 20:56:09 +0000
Subject: [PATCH] fixed bug with specutils reading in wrong units; added proper
 error to HST

---
 spectroscopy/code_src/mast_functions.py | 15 +++++++++++----
 spectroscopy/spectra_generator.md       | 22 +++++++++++++---------
 2 files changed, 24 insertions(+), 13 deletions(-)

diff --git a/spectroscopy/code_src/mast_functions.py b/spectroscopy/code_src/mast_functions.py
index d8a9c204..500f303d 100644
--- a/spectroscopy/code_src/mast_functions.py
+++ b/spectroscopy/code_src/mast_functions.py
@@ -151,11 +151,17 @@ def JWST_get_spec_helper(sample_table, search_radius_arcsec, datadir, verbose):
                     file_idx = np.where( [ tab["productFilename"][jj] in download_results["Local Path"][iii] for iii in range(len(download_results))] )[0]
                     
                     # open spectrum
+                    # Note that specutils returns the wrong units. Use Table.read() instead.
                     filepath = download_results["Local Path"][file_idx[0]]
+                    spec1d = Table.read(filepath , hdu=1)
+                    
                     #print(filepath)
-                    spec1d = Spectrum1D.read(filepath)
+                    #spec1d = Spectrum1D.read(filepath)
                     
-                    dfsingle = pd.DataFrame(dict(wave=[spec1d.spectral_axis] , flux=[spec1d.flux], err=[np.repeat(0,len(spec1d.flux))],
+                    dfsingle = pd.DataFrame(dict(#wave=[spec1d.spectral_axis] , flux=[spec1d.flux], err=[np.repeat(0,len(spec1d.flux))],,
+                                                 wave=[spec1d["WAVELENGTH"].data * spec1d["WAVELENGTH"].unit] ,
+                                                 flux=[spec1d["FLUX"].data * spec1d["FLUX"].unit],
+                                                 err=[spec1d["FLUX_ERROR"].data * spec1d["FLUX_ERROR"].unit],
                                                  label=[stab["label"]],
                                                  objectid=[stab["objectid"]],
                                                  #objID=[tab["objID"][jj]], # REMOVE
@@ -205,7 +211,7 @@ def JWST_group_spectra(df, verbose, quickplot):
     objects_unique = np.unique(tab["label"])
 
     for obj in objects_unique:
-        if verbose: print("++Processing Object {} ++".format(obj))
+        print("Grouping object {}".format(obj))
 
         ## Get filters
         filters_unique = np.unique(tab["filter"])
@@ -349,7 +355,8 @@ def HST_get_spec(sample_table, search_radius_arcsec, datadir, verbose):
                     spec1d = Spectrum1D.read(filepath)
 
                     # Note: this should be in erg/s/cm2/A and any wavelength unit.
-                    dfsingle = pd.DataFrame(dict(wave=[spec1d.spectral_axis] , flux=[spec1d.flux], err=[np.repeat(0,len(spec1d.flux))],
+                    dfsingle = pd.DataFrame(dict(#wave=[spec1d.spectral_axis] , flux=[spec1d.flux], err=[np.repeat(0,len(spec1d.flux))],
+                                                 wave=[spec1d.spectral_axis] , flux=[spec1d.flux], err=[spec1d.uncertainty.array * spec1d.uncertainty.unit],
                                                  label=[stab["label"]],
                                                  objectid=[stab["objectid"]],
                                                  #objID=[tab["objID"][jj]],
diff --git a/spectroscopy/spectra_generator.md b/spectroscopy/spectra_generator.md
index 746b5e72..ceb999a6 100644
--- a/spectroscopy/spectra_generator.md
+++ b/spectroscopy/spectra_generator.md
@@ -134,14 +134,14 @@ labels = []
 coords.append(SkyCoord("{} {}".format("09 54 49.40" , "+09 16 15.9"), unit=(u.hourangle, u.deg) ))
 labels.append("NGC3049")
 
-coords.append(SkyCoord("{} {}".format("12 45 17.44 " , "27 07 31.8"), unit=(u.hourangle, u.deg) ))
-labels.append("NGC4670")
+#coords.append(SkyCoord("{} {}".format("12 45 17.44 " , "27 07 31.8"), unit=(u.hourangle, u.deg) ))
+#labels.append("NGC4670")
 
-coords.append(SkyCoord("{} {}".format("14 01 19.92" , "−33 04 10.7"), unit=(u.hourangle, u.deg) ))
-labels.append("Tol_89")
+#coords.append(SkyCoord("{} {}".format("14 01 19.92" , "−33 04 10.7"), unit=(u.hourangle, u.deg) ))
+#labels.append("Tol_89")
 
-coords.append(SkyCoord(233.73856 , 23.50321, unit=u.deg ))
-labels.append("Arp220")
+#coords.append(SkyCoord(233.73856 , 23.50321, unit=u.deg ))
+#labels.append("Arp220")
 
 coords.append(SkyCoord( 150.091 , 2.2745833, unit=u.deg ))
 labels.append("COSMOS1")
@@ -149,8 +149,8 @@ labels.append("COSMOS1")
 coords.append(SkyCoord( 150.1024475 , 2.2815559, unit=u.deg ))
 labels.append("COSMOS2")
 
-coords.append(SkyCoord("{} {}".format("150.000" , "+2.00"), unit=(u.deg, u.deg) ))
-labels.append("COSMOS3")
+#coords.append(SkyCoord("{} {}".format("150.000" , "+2.00"), unit=(u.deg, u.deg) ))
+#labels.append("COSMOS3")
 
 coords.append(SkyCoord("{} {}".format("+53.15508" , "-27.80178"), unit=(u.deg, u.deg) ))
 labels.append("JADESGS-z7-01-QU")
@@ -232,7 +232,7 @@ This archive includes spectra taken by
 ```python
 %%time
 ## Get Spectra for HST
-df_spec_HST = HST_get_spec(sample_table , search_radius_arcsec = 0.5, datadir = "./data/", verbose = False)
+df_spec_HST = HST_get_spec(sample_table , search_radius_arcsec = 0.5, datadir = "./data/", verbose = True)
 df_spec.append(df_spec_HST)
 ```
 
@@ -243,6 +243,10 @@ df_jwst = JWST_get_spec(sample_table , search_radius_arcsec = 0.5, datadir = "./
 df_spec.append(df_jwst)
 ```
 
+```python
+df_spec.data
+```
+
 ### 2.3 SDSS Archive
 
 This includes SDSS spectra.