Manipulating an SED object

SED objects represent SED templates belonging to one of the three possible classes “STAR”, “QSO” (for AGN type of objects), and “GAL” for galaxies. SED templates available with LePhare can be found under the sed directory.

The notebook can be downloaded here.

import os
import lephare as lp
import numpy as np
from matplotlib import pylab as plt
import struct

%matplotlib inline
%load_ext wurlitzer

User defined LEPHAREDIR is set. Code runs depend on all required
auxiliary data being present at /home/docs/.cache/lephare/data.
User defined LEPHAREWORK is set. All intermediate files will be written to:
 /home/docs/.cache/lephare/work
OMP_NUM_THREADS set to 1 (max available CPUs: 2)

keymap = lp.all_types_to_keymap(lp.default_cosmos_config)
# Get the auxiliary files required.
lp.data_retrieval.get_auxiliary_data(keymap=keymap, additional_files=["sed/STAR/PICKLES/f2ii.sed"])

Registry file downloaded and saved as data_registry.txt.

/home/docs/checkouts/readthedocs.org/user_builds/lephare/envs/latest/lib/python3.10/site-packages/lephare/data_retrieval.py:431: UserWarning: /home/docs/.cache/lephare/data/sed/STAR/STAR_MOD_ALL.list is present locally and will not be overwritten.
  warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/lephare/envs/latest/lib/python3.10/site-packages/lephare/data_retrieval.py:431: UserWarning: /home/docs/.cache/lephare/data/sed/GAL/COSMOS_SED/COSMOS_MOD.list is present locally and will not be overwritten.
  warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/lephare/envs/latest/lib/python3.10/site-packages/lephare/data_retrieval.py:431: UserWarning: /home/docs/.cache/lephare/data/sed/QSO/SALVATO09/AGN_MOD.list is present locally and will not be overwritten.
  warnings.warn(
Downloading file 'sed/STAR/PICKLES/f2ii.sed' from 'https://raw.githubusercontent.com/lephare-photoz/lephare-data/main/sed/STAR/PICKLES/f2ii.sed' to '/home/docs/.cache/lephare/data'.

Checking/downloading 442 files...
442 completed.
All files downloaded successfully and are non-empty.
Checking/downloading 1 files...
1 completed.
All files downloaded successfully and are non-empty.

filter_name = "cosmos/u_cfht.lowres"
f1 = lp.flt(0, filter_name, 1, 0)
f1.read(os.path.join(lp.LEPHAREDIR, "filt", filter_name))
f1.plot_filter_curve()

# we can write the config to a file to keep a record
config_file = "./config_file.para"
lp.write_para_config(keymap, config_file)
# before running PhotoZ we must run filt, sedtolib and maggal
filter_driver = lp.Filter(config_file=config_file)
filter_driver.run()

[<lephare._lephare.flt at 0x7432eb0ebd30>,
 <lephare._lephare.flt at 0x7432eb0ebd80>,
 <lephare._lephare.flt at 0x7432eb0ebc40>,
 <lephare._lephare.flt at 0x7432eb0eb6f0>,
 <lephare._lephare.flt at 0x7432eb0e9260>,
 <lephare._lephare.flt at 0x7432eb2b0180>,
 <lephare._lephare.flt at 0x7432eb2b0400>,
 <lephare._lephare.flt at 0x7432eb2b0450>,
 <lephare._lephare.flt at 0x7432eb2b04a0>,
 <lephare._lephare.flt at 0x7432eb2b04f0>,
 <lephare._lephare.flt at 0x7432eb2b0540>,
 <lephare._lephare.flt at 0x7432eb2b0590>,
 <lephare._lephare.flt at 0x7432eb2b05e0>,
 <lephare._lephare.flt at 0x7432eb2b0630>,
 <lephare._lephare.flt at 0x7432eb2b0680>,
 <lephare._lephare.flt at 0x7432eb2b06d0>,
 <lephare._lephare.flt at 0x7432eb2b0720>,
 <lephare._lephare.flt at 0x7432eb2b0770>,
 <lephare._lephare.flt at 0x7432eb2b07c0>,
 <lephare._lephare.flt at 0x7432eb2b0810>,
 <lephare._lephare.flt at 0x7432eb2b0860>,
 <lephare._lephare.flt at 0x7432eb2b08b0>,
 <lephare._lephare.flt at 0x7432eb2b0900>,
 <lephare._lephare.flt at 0x7432eb2b0950>,
 <lephare._lephare.flt at 0x7432eb2b09a0>,
 <lephare._lephare.flt at 0x7432eb2b09f0>,
 <lephare._lephare.flt at 0x7432eb2b0a40>,
 <lephare._lephare.flt at 0x7432eb2b0a90>,
 <lephare._lephare.flt at 0x7432eb2b0ae0>,
 <lephare._lephare.flt at 0x7432eb2b0b30>]

filter_output = os.path.join(
    os.environ["LEPHAREWORK"], "filt", filter_driver.keymap["FILTER_FILE"].value + ".dat"
)
filters = np.loadtxt(filter_output, dtype={"names": ("lamb", "val", "bid"), "formats": (float, float, int)})
plt.loglog(filters["lamb"], filters["val"])
plt.xlabel("wavelength");

!ls $LEPHAREDIR/sed

GAL  QSO  STAR

Let’s start with the template of a star. A SED object is created most easily with 3 arguments: its name, an identifying integer, and the type it belongs to. Then calling the read function reads the ASCII file passed as argument.

star_sed = "f2ii.sed"
sed_filename = os.path.join(lp.LEPHAREDIR, "sed/STAR/PICKLES/", star_sed)
sed = lp.StarSED(star_sed, 1)
sed.read(sed_filename)

The python code exposes the templates vectors through the data method

x = sed.data()[0]
y = sed.data()[1]
plt.plot(x, y)

[<matplotlib.lines.Line2D at 0x7432eaf14a90>]

In the context of LePhare, such a SED object is going to be written as a byte compressed file, and read back downstream to compute expected magnitudes and to perform the fit. Writing an reading this binary stored file goes as follows:

rootname = star_sed.split(".")[0]
sed.writeSED(rootname + ".bin", rootname + ".phys", rootname + ".doc")

We can read it back into a new SED object, and check that the values have been correctly read back

sed2 = lp.StarSED(star_sed, 2)
sed2.readSEDBin(rootname + ".bin")
x2 = sed.data()[0]
y2 = sed.data()[1]
assert np.all(x == x2)
assert np.all(y == y2)
print(sed.name, sed2.name)

f2ii.sed MOD_1

sed = lp.Sedtolib(config_file)
sed.run(typ="STAR", star_sed=f"{lp.LEPHAREDIR}/sed/STAR/STAR_MOD_ALL.list", star_fscale="1")

sed_output = os.path.join(os.environ["LEPHAREWORK"], "lib_bin", keymap["STAR_LIB"].value + ".bin")

sed_output

'/home/docs/.cache/lephare/work/lib_bin/LIB_STAR.bin'

# Code to unpack binary not running
# buf = open(sed_output, "rb").read()
# counter = 0
# off = 0
# while counter <= 0:
#     nrec, jtype, nbw = struct.unpack("iil", buf[off : off + 16])
#     wave = struct.unpack("d" * nbw, buf[off + 16 : off + 16 + 8 * nbw])
#     print(nrec, jtype, nbw)
#     off += 16 + 8 * nbw
#     nrec, jtype, nbw = struct.unpack("iil", buf[off : off + 16])
#     spec = struct.unpack("d" * nbw, buf[off + 16 : off + 16 + 8 * nbw])
#     print(nrec, jtype, nbw)
#     off += 16 + 8 * nbw
#     counter += 1
# plt.plot(wave, spec)
# # plt.xlim(1000,1.e4)
# # plt.ylim(1.e-10,1.e-6)

# n = 6973
# struct.unpack("d", buf[off + 16 + 8 * (n - 1) : off + 16 + 8 * n])

# nbw

# buf[38896 + 16 + 8 : 38896 + 16 + 16]

# struct.unpack("d", buf[38896 + 16 + 8 : 38896 + 16 + 16])

# print(4.2399775704000006e-08 / 1.23542470e01)

# struct.unpack("d" * 1, buf[off + 16 : off + 16 + 8])

# struct.calcsize("l")

# wave[1], spec[1]

d = np.loadtxt(f"{lp.LEPHAREDIR}/sed/STAR/PICKLES/o5v.sed.ext")
plt.plot(d[:, 0], d[:, 1])

[<matplotlib.lines.Line2D at 0x7432eaf5b8e0>]

d[:, 1]

array([0.00000000e+00, 1.23542470e+01, 9.92826800e+00, ...,
       2.15824600e-65, 4.93845048e-66, 1.13000525e-66], shape=(4860,))

stream_buf = open(sed_output, "rb").read()

struct.unpack("iil", stream_buf[0:16])

(1, 4860, 0)

k = lp.keyword("TEST", "0.1,0.1,6")
print(k.split_double("0.03", 3))
print(k.split_double("0.03", 3))

[0.1, 0.1, 6.0]
[0.1, 0.1, 6.0]

mag = lp.MagGal(config_file)
mag.run(typ="STAR", lib_ascii="YES", star_lib_out="ALLSTAR_COSMOS")

sed.run(typ="QSO", qso_sed=f"{lp.LEPHAREDIR}/sed/QSO/SALVATO09/AGN_MOD.list")

sed.run(typ="GAL", gal_sed=f"{lp.LEPHAREDIR}/sed/GAL/COSMOS_SED/COSMOS_MOD.list", gal_lib="LIB_VISTA")

mag = lp.MagGal(config_file)
mag.run(
    typ="QSO",
    lib_ascii="NO",
    mod_extinc="0,1000",
    eb_v="0.,0.1,0.2,0.3",
    extinc_law="SB_calzetti.dat",
    z_step="0.04,0,6",
    verbose=True,
)