some small fixes and edits to work with modap

frb/data/analysis/CIGALE/VLT_HAWKI_H.dat

@@ -1,4 +1,5 @@
 #VLT_HAWKI_H
+#photon
 #VLT HAWKI H band http://www.eso.org/sci/facilities/paranal/instruments/hawki/inst.html
 4992.5 0.000316136
 5012.47 0.00026036

frb/data/analysis/CIGALE/VLT_HAWKI_J.dat

@@ -1,4 +1,5 @@
 #VLT_HAWKI_J
+#photon
 #VLT HAWKI J band http://www.eso.org/sci/facilities/paranal/instruments/hawki/inst.html
 5014.98 -0.000258185
 5034.96 0.000143311

frb/data/analysis/CIGALE/VLT_HAWKI_Ks.dat

@@ -1,4 +1,5 @@
 #VLT_HAWKI_Ks
+#photon
 #VLT HAWKI Ks band http://www.eso.org/sci/facilities/paranal/instruments/hawki/inst.html
 4970.0 0.000266825
 4989.88 -0.00022829

frb/galaxies/cigale.py

@@ -181,6 +181,16 @@ def gen_cigale_in(photometry_table, zcol, idcol=None, infile="cigale_in.fits",
         'Pan-STARRS_y': 'panstarrs.ps1.y',
         'LRISr_I': 'LRIS_I',
         'LRISb_V': 'LRIS_V',
+        'LRISb_G': 'LRIS_G',
+        'LRISr_R': 'LRIS_R',
+        'GMOS_N_g': "GMOS_N_g",
+        'GMOS_N_r': "GMOS_N_r",
+        'GMOS_N_i': "GMOS_N_i",
+        'GMOS_N_z': "GMOS_N_z",
+        'GMOS_S_g': "GMOS_S_g",
+        'GMOS_S_r': "GMOS_S_r",
+        'GMOS_S_i': "GMOS_S_i",
+        'GMOS_S_z': "GMOS_S_z",
         'WFC3_F160W': 'hst.wfc3.ir.F160W',
         'WFC3_F300X': 'WFC3_F300X', 
         'Spitzer_3.6': 'spitzer.irac.l1',
@@ -392,10 +402,17 @@ def run(photometry_table, zcol,
             series = ['stellar_attenuated', 'stellar_unattenuated', 'dust', 'agn', 'model']
             SED(config = cigconf, sed_type = "mJy", nologo = True,
                 xrange = (False, False), yrange =  (False, False),
-                series = series, format =  "pdf", outdir =  Path("out"))
+                series = series, format =  "png", outdir =  Path("out"))
             # Set back to a GUI
+            # import matplotlib
+            # matplotlib.use('TkAgg')
             import matplotlib
-            matplotlib.use('TkAgg')
+            if matplotlib.get_backend().lower() != 'tkagg':
+                try:
+                    matplotlib.use('TkAgg')
+                except ImportError:
+                    pass
+
 
     # Rename the default output directory?
     if outdir != 'out':
@@ -452,3 +469,54 @@ def host_run(host, cut_photom=None, cigale_file=None):
         sfh_file = cigale_file.replace('CIGALE', 'CIGALE_SFH')
         os.system('mv {:s}/{:s}_SFH.fits {:s}'.format(host.name, host.name, sfh_file))
     return
+
+def add_text_SED(host, cigale_results, out_name=None):
+    """
+    Add a text file with the SED results
+    to the host object.
+    Args
+    ----
+    host (FRBGalaxy): A host galaxy.
+    """
+    try:
+        from pcigale_plots.plot_types.sed import SED
+    except ImportError:
+        console.print(f"{ERROR} This wrapper is compatible with CIGALE v. 2025. and later. Please update your version.")
+        pass
+
+    # Get unique surveys to put on plot
+    surveys = set()
+    for key in host.photom.keys():
+        survey = key.split("_")[0]  # take everything before the first underscore
+        if survey == "Pan-STARRS":
+            survey = "PS1"
+        surveys.add(survey)
+    surveys = sorted(list(surveys))
+
+    from PIL import Image, ImageDraw, ImageFont
+
+    # Open the image
+    img = Image.open(f"{host.name}/{host.name}_best_model.png")
+
+    # Create a drawing object
+    draw = ImageDraw.Draw(img)
+
+    # Define font (use a system font or specify a .ttf file)
+    font = ImageFont.load_default(size=22)
+
+    plot_text = f"Mstar = {cigale_results['Mstar']:.2e} Msun \n"
+    plot_text += f"SFR = {cigale_results['SFR_photom']:.2f} Msun/yr\n"
+    plot_text += "Photometry: \n  + " + "\n  + ".join(surveys)
+
+    # Add text (x, y coordinates, text, color, font)
+    draw.text((20, 15), plot_text, fill="black", font=font)
+
+    # Save result
+    if out_name is not None:
+        img.save(f"{host.name}/{out_name}")
+        print(f"Saved {host.name}/{out_name}")
+    else:
+        img.save(f"{host.name}/{host.name}_best_model_text.png")
+        print(f"Saved {host.name}/{host.name}_best_model_text.png")
+
+    return

frb/galaxies/eazy.py

@@ -598,4 +598,4 @@ def getEazyPz(idx, MAIN_OUTPUT_FILE='photz', OUTPUT_DIRECTORY='./OUTPUT', CACHE_
     if get_prior:
         return tempfilt['zgrid'], pzi, prior
     else:
-        return tempfilt['zgrid'], pzi
+        return tempfilt['zgrid'], pzi

frb/galaxies/photom.py

@@ -187,9 +187,13 @@ def extinction_correction(filt, EBV, RV=3.1, max_wave=None, required=True):
 
     source_flux = 1.
     #calculate linear correction
-    delta = np.trapezoid(throughput * source_flux * 
-                     10 ** (-0.4 * Alambda), x=wave) / np.trapezoid(
-                         throughput * source_flux, x=wave)
+    if getattr(np, "trapz", None) is not None:
+        delta = np.trapz(throughput * source_flux * 
+                        10 ** (-0.4 * Alambda), wave) / np.trapz(
+                            throughput * source_flux, wave)
+    else:
+        delta = np.trapezoid(throughput * source_flux * 10 ** (-0.4 * Alambda), wave) / np.trapezoid(
+                            throughput * source_flux, wave)
 
     correction = 1./delta
 
@@ -212,7 +216,7 @@ def correct_photom_table(photom, EBV, name, max_wave=None, required=True):
             Required keys: 'Name', filters
         EBV (float):
             E(B-V) (can get from frb.galaxies.nebular.get_ebv which uses IRSA Dust extinction query
-        name (str):\
+        name (str):
             Name of the object to correct
         required (bool, optional):
             Crash out if the transmission curve is not present
@@ -250,12 +254,13 @@ def correct_photom_table(photom, EBV, name, max_wave=None, required=True):
                 continue
         except:
             embed(header='187 in photom')
-        # SDSS
-        if 'SDSS' in filt:
-            if 'extinction_{}'.format(filt[-1]) in photom.keys():
-                print("Appying SDSS-provided extinction correction")
-                cut_photom[key] -= cut_photom['extinction_{}'.format(filt[-1])]
-                continue
+        print('Correcting filter {} for Galactic extinction'.format(filt))
+        # SDSS ### removing this because SDSS correction is sometimes different for very dusty sightlines, but TBD
+        # if 'SDSS' in filt:
+        #     if 'extinction_{}'.format(filt[-1]) in photom.keys():
+        #         print("Appying SDSS-provided extinction correction")
+        #         cut_photom[key] -= cut_photom['extinction_{}'.format(filt[-1])]
+        #         continue
         # Hack for LRIS
         if 'LRIS' in filt:
             _filter = 'LRIS_{}'.format(filt[-1])
@@ -268,6 +273,7 @@ def correct_photom_table(photom, EBV, name, max_wave=None, required=True):
                                              required=required)
         mag_dust = 2.5 * np.log10(1. / dust_correct)
         cut_photom[key] += mag_dust
+
     # Add it back in
     photom[idx] = cut_photom
 

frb/halos/hmf.py

@@ -49,8 +49,8 @@ def init_hmf():
 #  Needs to come after def init_hmf()
 try:
     import hmf_emulator
-except:
-    print("If you wish to use the halos modules, you need the Aemulus HMF emulator. Please install it: github.com/AemulusProject/hmf_emulator")
+except ImportError:
+    warnings.warn("If you wish to use the halos modules, you need the Aemulus HMF emulator. Please install it: github.com/AemulusProject/hmf_emulator", ImportWarning, stacklevel=2)
 else:
     hmfe = init_hmf()
 

frb/halos/models.py

@@ -55,8 +55,8 @@ def init_hmf():
 # Storing for use
 try:
     import hmf_emulator
-except:
-    print("If you wish to use the halos modules, you need the Aemulus HMF emulator. Please install it: github.com/AemulusProject/hmf_emulator")
+except ImportError:
+    warnings.warn("If you wish to use the halos modules, you need the Aemulus HMF emulator. Please install it: github.com/AemulusProject/hmf_emulator", ImportWarning, stacklevel=2)
 else:
     hmfe = init_hmf()
 

frb/scripts/r_vs_dm.py

@@ -165,4 +165,4 @@ def r_vs_dm(outfile='fig_r_vs_z.png',
 ''' Run me '''
 # Command line execution
 if __name__ == '__main__':
-    r_vs_dm()
+    r_vs_dm()

frb/surveys/catalog_utils.py

@@ -100,7 +100,11 @@ def match_ids(IDs, match_IDs, require_in_match=True):
     """
     rows = -1 * np.ones_like(IDs).astype(int)
     # Find which IDs are in match_IDs
-    in_match = np.in1d(IDs, match_IDs)
+    try:
+        in_match = np.isin(IDs, match_IDs)
+    except AttributeError:
+        # Older NumPy versions
+        in_match = np.in1d(IDs, match_IDs)
     if require_in_match:
         if np.sum(~in_match) > 0:
             raise IOError("qcat.match_ids: One or more input IDs not in match_IDs")
@@ -213,8 +217,8 @@ def xmatch_catalogs(cat1:Table, cat2:Table, dist:units.Quantity = 5*units.arcsec
         dist1 = None
         dist2 = None
     # Get corodinates
-    cat1_coord = SkyCoord(cat1[RACol1]*units.deg, cat1[DecCol1]*units.deg, distance=dist1)
-    cat2_coord = SkyCoord(cat2[RACol2]*units.deg, cat2[DecCol2]*units.deg, distance=dist2)
+    cat1_coord = SkyCoord(cat1[RACol1], cat1[DecCol1], unit='deg', distance=dist1)
+    cat2_coord = SkyCoord(cat2[RACol2], cat2[DecCol2], unit='deg', distance=dist2)
 
     # Match 2D
     idx, d2d, d3d = cat1_coord.match_to_catalog_sky(cat2_coord)

frb/surveys/delve.py

@@ -20,14 +20,19 @@
 photom = {}
 photom['DELVE'] = {}
 photom['DELVE']['DELVE_ID'] = 'quick_object_id'
+# photom['DELVE']['DELVE_ID'] = 'coadd_object_id'
 photom['DELVE']['ra'] = 'ra'
 photom['DELVE']['dec'] = 'dec'
 photom['DELVE']['ebv'] = 'ebv' # Schegel, Finkbeiner, Davis (1998)
+# photom['DELVE']['ebv'] = '_sfd98'
 DELVE_bands = ['g', 'r', 'i', 'z']
 for band in DELVE_bands:
     photom['DELVE'][f'DELVE_{band}'] = f'mag_auto_{band}' #mag
     photom['DELVE'][f'DELVE_{band}_err'] = f'magerr_auto_{band}' #magerr
     photom['DELVE'][f'class_star_{band}'] = f'class_star_{band}' #morphology class
+    ## 'bdf_mag_' 'bdf_mag_err_'  'gap_mag_' 'gap_mag_err_'  'mag_detmodel_' 'magerr_detmodel_'
+    # photom['DELVE'][f'DELVE_{band}'] = 'bdf_mag_'+band
+    # photom['DELVE'][f'DELVE_{band}_err'] = 'bdf_mag_err_'+band
 
 
 class DELVE_Survey(dlsurvey.DL_Survey):
@@ -48,8 +53,10 @@ def __init__(self, coord, radius, **kwargs):
         self.survey = 'DELVE'
         self.bands = DELVE_bands
         self.svc = sia.SIAService("https://datalab.noao.edu/sia/delve_dr2")
+        # self.svc = sia.SIAService("https://datalab.noao.edu/sia/delve_dr3")
         self.qc_profile = "default"
         self.database = "delve_dr2.objects"
+        # self.database = "delve_dr3.coadd_objects"
         self.default_query_fields = list(photom['DELVE'].values())
 
     def get_catalog(self, query=None, query_fields=None, print_query=False,**kwargs):

frb/surveys/galex.py

@@ -63,7 +63,7 @@ def get_catalog(self,query_fields=None, print_query=False):
             query_fields = _DEFAULT_query_fields
         else:
             query_fields = _DEFAULT_query_fields+query_fields
-
+        import astropy.units as u
         data = {}
         data['ra'] = self.coord.ra.value
         data['dec'] = self.coord.dec.value

frb/surveys/nedlvs.py

@@ -35,7 +35,9 @@ def __init__(self, coord, radius=90.*u.deg, cosmo=None, **kwargs):
 
         # Set redshift distances using the cosmology of choice
         redshift_dist_sources = self.datatab['DistMpc_method']=='Redshift'
-        self.datatab['DistMpc'][redshift_dist_sources] = self.cosmo.luminosity_distance(self.datatab['z'][redshift_dist_sources]).to('Mpc').value
+        # self.datatab['DistMpc'][redshift_dist_sources] = self.cosmo.luminosity_distance(self.datatab['z'][redshift_dist_sources]).to('Mpc').value
+        z_vals = np.asarray(self.datatab['z'][redshift_dist_sources])
+        self.datatab['DistMpc'][redshift_dist_sources] = self.cosmo.luminosity_distance(z_vals).to('Mpc').value
         self.datatab['phys_sep'] = self.datatab['DistMpc']*u.Mpc*np.sin(self.datatab['ang_sep'].to('rad').value)
 
     def get_column_names(self):

frb/surveys/panstarrs.py

@@ -13,6 +13,8 @@
 from ..galaxies.defs import PanSTARRS_bands
 import importlib_resources
 
+from .images import grab_from_url
+
 import warnings
 import requests
 

frb/surveys/psrcat.py

@@ -6,8 +6,11 @@
 from astropy.coordinates import SkyCoord
 from astropy import units
 
+import warnings
 try:
-    from pulsars import io as pio
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore", category=UserWarning, message="pkg_resources")
+        from pulsars import io as pio
 except ImportError:
     print("Warning:  You need FRB/pulsars installed to use PSRCat")
 

frb/surveys/survey_utils.py

@@ -270,13 +270,64 @@ def search_all_surveys(coord:SkyCoord, radius:u.Quantity, include_radio:bool=Fal
         # Fill in any empty separations and sort them.
         combined_cat['separation'] = coord.separation(SkyCoord(combined_cat['ra'], combined_cat['dec'], unit='deg')).to(u.arcmin)
         combined_cat.sort('separation')
-
+        combined_cat = pick_best_row_by_phot(combined_cat)
         # Make the ra, dec, separation the first columns
         colnames = combined_cat.colnames
         other_cols = np.setdiff1d(colnames, ['ra', 'dec', 'separation'])
         combined_cat = combined_cat[['ra', 'dec', 'separation']+other_cols.tolist()]
 
     return combined_cat
+
+def pick_best_row_by_phot(cat, mag_cols=None, max_sep=None):
+    """
+    From a catalog sorted by separation, pick the best row
+    based on having any good photometry in the specified
+    magnitude columns.
+    Args:
+        cat (Table): Astropy table sorted by separation.
+        mag_cols (list, optional): List of magnitude column names
+            to consider for "good photometry". If None,
+            will use all valid filters from frb.galaxies.defs.
+        max_sep (Quantity, optional): Maximum separation
+            to consider. If None, no maximum separation
+            is applied.
+
+    Returns:
+        best_row (Table): A 1-row table with the best row.
+    """
+    if len(cat) == 0:
+        return cat
+
+    if mag_cols is None:
+        from frb.galaxies.defs import valid_filters
+        mag_cols = valid_filters
+
+    # require separation column already computed in arcmin/arcsec/whatever
+    sep = cat['separation']
+    if max_sep is not None:
+        inrad = sep <= max_sep
+    else:
+        inrad = np.ones(len(cat), dtype=bool)
+
+    # “has any good mag”
+    has_good = np.zeros(len(cat), dtype=bool)
+    for c in mag_cols:
+        if c in cat.colnames:
+            m = np.array(cat[c])
+            good = np.isfinite(m) & (m < 900) & (m > -10)
+            has_good |= good
+
+    ok = inrad & has_good
+    if not np.any(ok):
+        # fall back: just nearest, but you’ll know it’s junk
+        i = np.argmin(sep)
+    else:
+        # nearest among rows with any good photometry
+        i = np.argmin(sep[ok])
+        i = np.where(ok)[0][i]
+
+    return cat[i:i+1]  # returns a 1-row table
+
 
 def PS1_tile(coord:SkyCoord, side:u.Quantity=1*u.deg, **kwargs)->Table:
     """
@@ -331,4 +382,4 @@ def PS1_tile(coord:SkyCoord, side:u.Quantity=1*u.deg, **kwargs)->Table:
             else:
                 continue
     combined_table = remove_duplicates(combined_table, 'Pan-STARRS_ID')
-    return combined_table
+    return combined_table

frb/surveys/twomass.py

@@ -16,8 +16,8 @@
 photom = {}
 photom['2MASS'] = {}
 for band in MASS_bands:
-    photom["2MASS"]["2MASS"+'_{:s}'.format(band)] = '{:s}_m'.format(band.lower()) # Many options for apertures, tbd
-    photom["2MASS"]["2MASS"+'_{:s}_err'.format(band)] = '{:s}_msig'.format(band.lower())
+    photom["2MASS"]["2MASS"+'_{:s}'.format(band)] = '{:s}_m'.format(band) # default aperture queried like this is k20fe
+    photom["2MASS"]["2MASS"+'_{:s}_err'.format(band)] = '{:s}_msig'.format(band)
     photom["2MASS"]["2MASS_ID"] = 'designation'
 photom["2MASS"]['ra'] = 'ra'
 photom["2MASS"]['dec'] = 'dec'