Adds reference catalog creation

[egg]

# egg configuration. Each key value pair corresponds to input arguments to egg.

#

# Run

#   $ egg-gencat help

# in order to find out all the available input arguments.

# print verbose output

verbose: 1

# EGG filename

out: data/catalog/test/egg.fits

# area in deg2

area: 0.10

# right ascension of center

ra0: +150.11916667

# declination of center

dec0: +2.20583333

# apparent flux bands

bands: [lsst-u,johnson-B,lsst-g,lsst-r,lsst-i,lsst-z,lsst-y,euclid-VIS,euclid-nisp-Y,euclid-nisp-J,euclid-nisp-H,spitzer-irac1,spitzer-irac2,spitzer-irac3,spitzer-irac4,wise-w1,wise-w2,wise-w3,wise-w4,mock-1000-4000,mock-1450,mock-4400,mock-5007,mock-15um,vista-Y,vista-J,vista-H,vista-Ks]

# absolute flux bands

rfbands: [lsst-u,johnson-B,lsst-g,lsst-r,lsst-i,lsst-z,lsst-y,euclid-VIS,euclid-nisp-Y,euclid-nisp-J,euclid-nisp-H,spitzer-irac1,spitzer-irac2,spitzer-irac3,spitzer-irac4,wise-w1,wise-w2,wise-w3,wise-w4,mock-1000-4000,mock-1450,mock-4400,mock-5007,mock-15um,vista-Y,vista-J,vista-H,vista-Ks]

# input stellar mass function

mass_func: data/egg/share/mass_func_cosmos2020_agn_zou24_loglambda32.fits

# minimum redshift

zmin: 0.21

# maximum redshift

zmax: 5.49

# redshift slice resolution

dz: 0.05

# minimum stellar mass in logMsun (Salpeter IMF)

mmin: 8.74

# NOTE: 0.30 deg should cover most cases i.e. r~0-200 Mpc/h at z~0.20-5.50

clust_r0: 0.30

clust_fclust_him: 0.40

clust_fclust_lom: 0.20

# save SED?

save_sed: 1

# random number seed

seed: 20250911

# filter database

filter_db: data/egg/share/filter-db/db.dat

        # Handles something else

        return None

    def write_reference_catalog(self, filename: str, maglim: float, selection_band: str):

        """

        Write a mock LSST reference catalog in csv format to the given

        'filename'. The reference catalog is cut to brighter than the magnitude

        limit 'maglim' in the band 'selection_band'.

        """

        if os.path.exists(filename):

            logger.info("Reference catalog exists. Will not overwrite.")

            return

        # Cut at some magnitude limit e.g. r < 24 and non-variable sources

        select = util.flux_to_mag(self.get_flux_total(selection_band)) < float(maglim)

        select *= ~self.catalog["is_agn"]

        c = self.catalog[select]

        # Initialize the data to be written

        data = {

            # ID

            "source_id": c["ID"],

            # coordinates

            "ra": c["RA"],

            "dec": c["DEC"],

            # Fluxes

            "u": util.flux_to_mag(c["lsst-u_total"]),

            "g": util.flux_to_mag(c["lsst-g_total"]),

            "r": util.flux_to_mag(c["lsst-r_total"]),

            "i": util.flux_to_mag(c["lsst-i_total"]),

            "z": util.flux_to_mag(c["lsst-z_total"]),

            "y": util.flux_to_mag(c["lsst-y_total"]),

            # NOTE: for pmra see here:

            #   https://gea.esac.esa.int/archive/documentation/GDR2/Gaia_archive/chap_datamodel/sec_dm_main_tables/ssec_dm_gaia_source.html#:~:text=pmra%20%3A%20Proper%20motion%20in%20right,direction%20of%20increasing%20right%20ascension.

            #   https://pipelines.lsst.io/modules/lsst.meas.algorithms/creating-a-reference-catalog.html

            "pmra":           np.where(np.isfinite(c["pmracosd"]), c["pmracosd"], 0),

            "pmdec":          np.where(np.isfinite(c["pmdec"]),    c["pmdec"],    0),

            "epoch":          c["RA"].size * ["2000-01-01 12:00:00"],

            "parallax":       np.where(np.isfinite(c["Z"]), 0.0, util.distance_modulus_to_parallax(c["D"])),

            # Errors

            "ra_error":       np.zeros(c.size),

            "dec_error":      np.zeros(c.size),

            "pmra_error":     np.zeros(c.size),

            "pmdec_error":    np.zeros(c.size),

            "parallax_error": np.zeros(c.size),

            # Correlations/covariances

            "ra_dec_corr":         np.zeros(c.size),

            "ra_parallax_corr":    np.zeros(c.size),

            "ra_pmra_corr":        np.zeros(c.size),

            "ra_pmdec_corr":       np.zeros(c.size),

            "dec_parallax_corr":   np.zeros(c.size),

            "dec_pmra_corr":       np.zeros(c.size),

            "dec_pmdec_corr":      np.zeros(c.size),

            "parallax_pmra_corr":  np.zeros(c.size),

            "parallax_pmdec_corr": np.zeros(c.size),

            "pmra_pmdec_corr":     np.zeros(c.size),

        }

        logger.info("Writing the reference catalog...")

        with open(filename, 'w') as f:

            to_write = ','.join(data.keys())

            f.write(f"{to_write}\n")

            for i in range(c["RA"].size):

                to_write = ','.join([f"{v[i]}" for v in data.values()])

                f.write(f"{to_write}\n")

        # Run the LSST tools ontop

        dirname = f"{os.path.dirname(filename)}/reference_catalog"

        util.create_directory(dirname)

        os.system(f"convertReferenceCatalog {dirname} {ROOT}/src/python/config_reference_catalog.py {filename}")

    def __getitem__(self, key):

        return self.catalog_combined[key]

from lsst.meas.algorithms import convertRefcatManager

#config.manager.retarget(convertRefcatManager.ConvertRefcatManager)

config.manager.retarget(convertRefcatManager.ConvertAgileManager)

config.dataset_config.ref_dataset_name = "my-refcat"

#!/usr/bin/env python3

# Create the truth summary entry in the butler repository

import os

import numpy as np

import pandas

from catalog_galaxy_agn import get_catalog_from_config

from lsst.daf import butler

from lsst.geom import Angle, SpherePoint

import util

dirname = "data/catalog/dr1_new_new"

###

# initializations

#c = get_catalog_from_config("etc/config_dr1.ini")

c = get_catalog_from_config("etc/config_dr1_new_new.ini")

my_butler = butler.Butler(f"{dirname}/repo", collections="skymaps")

my_skymap = my_butler.get("skyMap")

# Set nan values to some sufficiently small magnitude

def get_mag(b):

    ret = util.flux_to_mag(c[f"lsst-{b}_total"])

    return np.where(np.isfinite(ret), ret, 99.0)

umag, gmag, rmag, imag, zmag, ymag = [get_mag(b) for b in "ugrizy"]

def get_flux(b):

    return c[f"lsst-{b}_total"] * 1000.0

uflux, gflux, rflux, iflux, zflux, yflux = [get_flux(b) for b in "ugrizy"]

# get truth_type

#   0: agn type 1

#   1: agn type 2

#   2: galaxy

#   3: star

is_star = ~(c["Z"] > 0)

is_agn = c["is_agn"]

is_agn_type1 = is_agn * ~c["is_optical_type2"]

is_agn_type2 = is_agn *  c["is_optical_type2"]

is_galaxy = ~is_star * ~is_agn

truth_type = (

    0 * is_agn_type1 +

    1 * is_agn_type2 +

    2 * is_galaxy +

    3 * is_star

)

files = {}

for i in range(c["RA"].size):

    print(i, end='\r')

    longitude = Angle(c["RA"][i] * np.pi / 180.0)

    latitude = Angle(c["DEC"][i] * np.pi / 180.0)

    coords = SpherePoint(longitude, latitude)

    for ret in my_skymap.findTractPatchList([coords]):

        # Gets the tract/patch

        tract = ret[0]._id

        patch = ret[1][0]._sequentialIndex

        # Opens the files for writing

        #filename = f"data/dr1/extra/truth/{tract}/{patch}/truth_catalog.csv"

        #filename = f"data/dr1/extra/truth/{tract}/truth_catalog.csv"

        filename = f"{dirname}/extra/truth/{tract}/truth_catalog.csv"

        if not os.path.exists(os.path.dirname(filename)):

            os.makedirs(os.path.dirname(filename))

        if not filename in files:

            files[filename] = open(filename, 'w')

            print("id,ra,dec,flux_u,flux_g,flux_r,flux_i,flux_z,flux_y,truth_type,is_unique_truth_entry,is_variable,is_pointsource,tract,patch,skymap", file=files[filename])

        print(

            ','.join([

                str(c["ID"][i]),                            # id

                str(c["RA"][i]),                            # ra

                str(c["DEC"][i]),                           # dec

                str(uflux[i]),                              # u

                str(gflux[i]),                              # g

                str(rflux[i]),                              # r

                str(iflux[i]),                              # i

                str(zflux[i]),                              # z

                str(yflux[i]),                              # y

                str(truth_type[i]),                         # truth_type

                "True",                                     # is_unique_truth_entry

                "True" if truth_type[i] <= 1 else "False",  # is_variable

                "True" if truth_type[i] == 3 else "False",  # is_pointsource

                str(tract),                                 # tract

                str(patch),                                 # patch

                "DC2",                                      # skymap

            ]),

            file=files[filename]

        )

###

# Creates the truth summary

cwd = os.getcwd() + '/'

#filename = "data/dr1/extra/truth/truth_summary.csv"

filename = f"{dirname}/extra/truth/truth_summary.csv"

print("Writing", filename)

with open(filename, 'w') as f:

    print("filename,tract,skymap", file=f)

    for filename_csv, v in files.items():

        ###

        # Converts .csv to .parq

        v.close()

        filename_parq = filename_csv.replace(".csv", ".parq")

        pandas.read_csv(filename_csv).to_parquet(filename_parq)

        print(filename_parq)

        tract = filename_parq.split('/')[-2]

        patch = filename_parq.split('/')[-1]

        #print(','.join([cwd + filename_csv, tract, patch, "DC2"]), file=f)

        print(','.join([cwd + filename_parq, tract, "DC2"]), file=f)

print(f"===============================================================================")

print(f"Run the following commands next")

print(f"===============================================================================")

print(f"butler remove-runs {dirname}/repo truth_summary")

print(f"butler ingest-files --transfer direct {dirname}/repo truth_summary truth_summary {dirname}/extra/truth/truth_summary.csv")

print(f"butler collection-chain --mode=extend {dirname}/repo u/viitanen/output truth_summary")

print(f"===============================================================================")

# https://rtn-029.lsst.io/#creating-and-populating-the-repository

#   butler register-dataset-type $REPO 'truth_summary' DataFrame skymap tract

os.system(f"butler remove-runs {dirname}/repo truth_summary")

os.system(f"butler ingest-files --transfer direct {dirname}/repo truth_summary truth_summary {dirname}/extra/truth/truth_summary.csv")

os.system(f"butler collection-chain --mode=extend {dirname}/repo u/viitanen/output truth_summary")

#!/bin/bash -x

if [ $# -lt 2 ] ; then

	echo "Usage: $0 [repo] (raws)"

	exit 1

fi

repo=$1

shift

raws="$@"

butler ingest-raws --transfer direct ${repo} ${raws}

butler define-visits --collections LSSTCam/raw/all ${repo} LSSTCam

# The following needs only to be done once, but the LSSTCam/raw/all

# collection must already exist and that only happens after the

# initial ingest of raw data. Additional raw data can be ingested

# with this script if this line is commented out.

butler collection-chain ${repo} LSSTCam/defaults 'LSSTCam/raw/all,LSSTCam/calib,u/jchiang/calib_imSim2.0_w_2023_04,refcats,skymaps'