add tiling, mosaic, metadetection

  # type: "uniform"

  # object_density: 37 # arcmin^-2

input_dir: "/public/home/fangyuedong/project/50sqDeg_L1_outputs/50sqDeg_Photo_W2/"

output_dir: "/public/home/fangyuedong/project/demo_csst_injection/test/"

pointing_label_list: ["MSC_0000000"]

chip_label_list: ["07"]

# input_dir: "/public/home/fangyuedong/project/50sqDeg_L1_outputs/50sqDeg_Photo_W2/"

# output_dir: "/public/home/fangyuedong/project/demo_csst_injection/test/"

# pointing_label_list: ["MSC_0000000"]

# chip_label_list: ["07"]

input_dir: "/public/home/fangyuedong/project/50sqDeg_L1_outputs/50sqDeg_Photo_W3/"

output_dir: "/public/home/fangyuedong/project/injected_50sqDeg_L1_outputs/g1_002_g_2002_50sqDeg_Photo_W3/"

pointing_label_list: ["MSC_0000000", "MSC_0000001",

"MSC_0000002", "MSC_0000003", "MSC_0000004", "MSC_0000005", "MSC_0000006", "MSC_0000007", "MSC_0000008", "MSC_0000009"]

chip_label_list: null

###############################################

# PSF setting

  # rotate galaxy ellipticity

  rotateEll: 0. # [degree]

  # if we want to apply constant shear (optional)

  g1_const: 0.02

  g2_const: 0.02

###############################################

# Instrumental effects setting

# Size of the taget Gaussian PSF

PSF_FWHM: 0.187

# Size of the image used to draw the PSF

PSF_IMG_SIZE: 151

# Size of one block [needs to be a square]

# IMG_SIZE: 9216

BLK_SIZE: 1024

# Size of the padding for each block

PAD_SIZE: 64

# Boundary used to avoid edge effects

# Objects for which the centre is within this distance from the edge will be

# masked out.

BOUND_SIZE: 20

OUT_COLUMNS: [flags, s2n, T, psf_T]

METADETECT_CONFIG:

    # Shape measurement method

    # wmom: weighted moments

    model: wmom

    # model: pgauss

    # Size of the weight function for the moments

    weight:

        fwhm: 6  # arcsec

        # fwhm: 1.

    # Metacal settings

    metacal:

        psf: fitgauss

        # Kind of shear applied to the image

        types: [noshear, 1p, 1m, 2p, 2m]

    # sx:

    #     # in sky sigma

    #     # DETECT_THRESH

    #     detect_thresh: 2.5

    #     # Minimum contrast parameter for deblending

    #     # DEBLEND_MINCONT

    #     deblend_cont: 0.00001

    #     # minimum number of pixels above threshold

    #     # DETECT_MINAREA: 6

    #     minarea: 3

    #     filter_type: conv

    #     # 7x7 convolution mask of a gaussian PSF with FWHM = 3.0 pixels.

    #     filter_kernel: [

    #         [0.004963, 0.021388, 0.051328, 0.068707, 0.051328, 0.021388, 0.004963],  # noqa

    #         [0.021388, 0.092163, 0.221178, 0.296069, 0.221178, 0.092163, 0.021388],  # noqa

    #         [0.051328, 0.221178, 0.530797, 0.710525, 0.530797, 0.221178, 0.051328],  # noqa

    #         [0.068707, 0.296069, 0.710525, 0.951108, 0.710525, 0.296069, 0.068707],  # noqa

    #         [0.051328, 0.221178, 0.530797, 0.710525, 0.530797, 0.221178, 0.051328],  # noqa

    #         [0.021388, 0.092163, 0.221178, 0.296069, 0.221178, 0.092163, 0.021388],  # noqa

    #         [0.004963, 0.021388, 0.051328, 0.068707, 0.051328, 0.021388, 0.004963],  # noqa

    #     ]

    sx: null

    # This is for the cutout at each detection

    meds:

        min_box_size: 32

        max_box_size: 128

        box_type: iso_radius

        rad_min: 4

        rad_fac: 2

        box_padding: 2

    # check for an edge hit

    bmask_flags: 2**30

    nodet_flags: 2**0

 No newline at end of file

import os

import galsim

import h5py as h5

import healpy as hp

import numpy as np

import treecorr

import matplotlib.pyplot as plt

# CONSTANTS

NSIDE = 128

# bundle_file_list = ['galaxies_C6_bundle000199.h5', 'galaxies_C6_bundle000200.h5', 'galaxies_C6_bundle000241.h5', 'galaxies_C6_bundle000242.h5', 'galaxies_C6_bundle000287.h5', 'galaxies_C6_bundle000288.h5', 'galaxies_C6_bundle000714.h5', 'galaxies_C6_bundle000715.h5', 'galaxies_C6_bundle000778.h5', 'galaxies_C6_bundle000779.h5', 'galaxies_C6_bundle000842.h5', 'galaxies_C6_bundle000843.h5', 'galaxies_C6_bundle002046.h5', 'galaxies_C6_bundle002110.h5', 'galaxies_C6_bundle002111.h5',

#                     'galaxies_C6_bundle002173.h5', 'galaxies_C6_bundle002174.h5', 'galaxies_C6_bundle002238.h5', 'galaxies_C6_bundle002596.h5', 'galaxies_C6_bundle002597.h5', 'galaxies_C6_bundle002656.h5', 'galaxies_C6_bundle002657.h5', 'galaxies_C6_bundle002711.h5', 'galaxies_C6_bundle002712.h5', 'galaxies_C6_bundle002844.h5', 'galaxies_C6_bundle002845.h5', 'galaxies_C6_bundle002884.h5', 'galaxies_C6_bundle002885.h5', 'galaxies_C6_bundle002921.h5', 'galaxies_C6_bundle002922.h5']

gals_cat_dir = '/public/share/yangxuliu/CSSOSDataProductsSims/data_test2024/galcat_C9/cat2CSSTSim_bundle-50sqDeg'

bundle_file_list = ['galaxies_C6_bundle002656.h5']

# def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7):

#     assert NSIDE == 2**healpixOrder

#     shift = healpixOrder - bundleOrder

#     shift = 2*shift

#     nside_bundle = 2**bundleOrder

#     nside_healpix = 2**healpixOrder

#     healpixID_nest = hp.ring2nest(nside_healpix, healpixID_ring)

#     bundleID_nest = (healpixID_nest >> shift)

#     bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest)

#     return bundleID_ring

def load_catalog_truth(file_path):

    gals_cat = h5.File(file_path, 'r')['galaxies']

    ra_arr = np.array([])

    dec_arr = np.array([])

    magr_arr = np.array([])

    for pix in gals_cat.keys():

        gals = gals_cat[str(pix)]

        ra_arr = np.concatenate((ra_arr, gals['ra']))

        dec_arr = np.concatenate((dec_arr, gals['dec']))

        magr_arr = np.concatenate((magr_arr, gals['mag_csst_r']))

    print(ra_arr.shape)

    # ind = magr_arr < 25.

    # print(ra_arr[ind].shape)

    # return ra_arr[ind], dec_arr[ind], magr_arr[ind]

    return ra_arr, dec_arr, magr_arr

def generate_rand_radec(N, ra_range=(0, 360), dec_range=(-30, 30)):

    """

    Generate N random points uniformly distributed over RA/Dec ranges.

    """

    ra = np.random.uniform(ra_range[0], ra_range[1], N)

    sin_dec_min = np.sin(np.radians(dec_range[0]))

    sin_dec_max = np.sin(np.radians(dec_range[1]))

    sin_dec = np.random.uniform(sin_dec_min, sin_dec_max, N)

    dec = np.degrees(np.arcsin(sin_dec))

    return ra, dec

def compute_correlation(ra, dec, ra_rand, dec_rand, ra_units='deg', dec_units='deg', min_sep=0.01, max_sep=10., nbins=20, sep_units='deg'):

    cat = treecorr.Catalog(ra=ra, dec=dec, ra_units=ra_units, dec_units=dec_units)

    cat_rand = treecorr.Catalog(ra=ra_rand, dec=dec_rand, ra_units=ra_units, dec_units=dec_units)

    nn = treecorr.NNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins, sep_units='deg')

    rr = treecorr.NNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins, sep_units='deg')

    dr = treecorr.NNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins, sep_units='deg')

    nn.process(cat)

    rr.process(cat_rand)

    dr.process(cat, cat_rand)

    return nn, rr, dr

def get_xi(nn, rr, dr):

    xi, varxi = nn.calculateXi(rr=rr, dr=dr)

    r = np.exp(nn.meanlogr)  # in degrees

    sig = np.sqrt(varxi)

    return r, xi, sig, varxi

def plot_wtheta(r, xi, sig, min_sep=0.01, max_sep=10., color='blue', legend=' ', leg_symble=[], leg_label=[]):

    xi, varxi = nn.calculateXi(rr=rr, dr=dr)

    r = np.exp(nn.meanlogr)  # in degrees

    sig = np.sqrt(varxi)

    plt.plot(r, xi, color=color)

    plt.plot(r, -xi, color=color, ls=':')

    plt.errorbar(r[xi>0], xi[xi>0], yerr=sig[xi>0], color=color, lw=0.5, ls='')

    plt.errorbar(r[xi<0], -xi[xi<0], yerr=sig[xi<0], color=color, lw=0.5, ls='')

    leg = plt.errorbar(-r, xi, yerr=sig, color=color)

    plt.xscale('log')

    plt.yscale('log', nonpositive='clip')

    plt.xlabel(r'$\theta$ (degrees)')

    # plt.legend([leg], [r'$w(\theta)$' + legend], loc='lower left')

    leg_symble.append(leg)

    leg_label.append(r'$w(\theta)$' + legend)

    plt.xlim([min_sep, max_sep])

    # plt.savefig('wtheta.png')

# def plot_wtheta(nn, rr, dr, min_sep=0.01, max_sep=10.):

#     xi, varxi = nn.calculateXi(rr=rr, dr=dr)

#     r = np.exp(nn.meanlogr)  # in degrees

#     sig = np.sqrt(varxi)

#     plt.plot(r, xi, color='blue')

#     plt.plot(r, -xi, color='blue', ls=':')

#     plt.errorbar(r[xi>0], xi[xi>0], yerr=sig[xi>0], color='blue', lw=0.1, ls='')

#     plt.errorbar(r[xi<0], -xi[xi<0], yerr=sig[xi<0], color='blue', lw=0.1, ls='')

#     # plt.errorbar(r[xi>0], xi[xi>0], yerr=sig[xi>0], fmt='o-', color='blue')

#     # plt.errorbar(r[xi<0], -xi[xi<0], yerr=sig[xi<0], fmt='o-', color='blue')

#     leg = plt.errorbar(-r, xi, yerr=sig, color='blue')

#     plt.xscale('log')

#     plt.yscale('log', nonpositive='clip')

#     plt.xlabel(r'$\theta$ (degrees)')

#     plt.legend([leg], [r'$w(\theta)$'], loc='lower left')

#     plt.xlim([min_sep, max_sep])

#     plt.savefig('wtheta.png')

if __name__ == "__main__":

    file_path = os.path.join(gals_cat_dir, bundle_file_list[0])

    ra, dec, magr = load_catalog_truth(file_path)

    # ra_rand, dec_rand = generate_rand_radec(N=len(ra)*2, ra_range=(ra.min(), ra.max()), dec_range=(dec.min(), dec.max()))

    # nn, rr, dr= compute_correlation(ra, dec, ra_rand, dec_rand, max_sep=5.)

    # plot_wtheta(nn, rr, dr, max_sep=5.)

    plt.figure()

    leg_symble = []

    leg_label = []

    ind = magr < 21.

    ra_rand, dec_rand = generate_rand_radec(N=len(ra[ind])*5, ra_range=(ra.min(), ra.max()), dec_range=(dec.min(), dec.max()))

    nn, rr, dr= compute_correlation(ra[ind], dec[ind], ra_rand, dec_rand, max_sep=5.)

    r, xi, sig, _ = get_xi(nn, rr, dr)

    plot_wtheta(r, xi, sig, max_sep=5., color='blue', legend=' mag_r < 21.', leg_symble=leg_symble, leg_label=leg_label)

    ind = magr < 24.

    # ra_rand, dec_rand = generate_rand_radec(N=len(ra[ind])*2, ra_range=(ra.min(), ra.max()), dec_range=(dec.min(), dec.max()))

    nn, rr, dr= compute_correlation(ra[ind], dec[ind], ra_rand, dec_rand, max_sep=5.)

    r, xi, sig, _ = get_xi(nn, rr, dr)

    plot_wtheta(r, xi, sig, max_sep=5., color='red', legend=' mag_r < 24.', leg_symble=leg_symble, leg_label=leg_label)

    # ind = magr < 25.

    # ra_rand, dec_rand = generate_rand_radec(N=len(ra[ind])*2, ra_range=(ra.min(), ra.max()), dec_range=(dec.min(), dec.max()))

    # nn, rr, dr= compute_correlation(ra[ind], dec[ind], ra_rand, dec_rand, max_sep=5.)

    # r, xi, sig, _ = get_xi(nn, rr, dr)

    # plot_wtheta(r, xi, sig, max_sep=5., color='green', legend=' mag_r < 25.', leg_symble=leg_symble, leg_label=leg_label)

    plt.legend(leg_symble, leg_label, loc='lower left')

    plt.savefig('wtheta.png')

 No newline at end of file

import numpy as np

import pandas as pd

from astropy.io import fits

from astropy.io.fits import Header

from astropy.wcs import WCS

from reproject import reproject_interp

def read_fits_image_wcs(path, ext=1):

    print(path)

    with fits.open(path, memmap=True) as hdul:

        data = hdul[ext].data.astype(np.float32, copy=False)

        wcs = WCS(hdul[ext].header)

    ny, nx = data.shape

    return data, wcs, nx, ny

def reproject_to_tile(data, wcs_in, tile_wcs, shape_out):

    # Returns (warped_data, footprint)

    warped, footprint = reproject_interp((data, wcs_in), tile_wcs, shape_out=shape_out)

    return warped.astype(np.float32, copy=False), footprint.astype(

        np.float32, copy=False

    )

def load_tile_from_parquet(parquet_path, tile_id):

    df = pd.read_parquet(

        parquet_path,

        columns=[

            "tile_id",

            "wcs_header",

            "nx_total",

            "ny_total",

            "ra_bbox_min",

            "ra_bbox_max",

            "dec_bbox_min",

            "dec_bbox_max",

            "ra_center",

            "dec_center",

        ],

    )

    row = df.loc[df.tile_id == tile_id].iloc[0]

    hdr = Header.fromstring(row["wcs_header"], sep="\n")

    tile_wcs = WCS(hdr)

    shape_out = (int(row["ny_total"]), int(row["nx_total"]))  # (ny, nx)

    tile_bbox = (

        float(row["ra_bbox_min"]),

        float(row["ra_bbox_max"]),

        float(row["dec_bbox_min"]),

        float(row["dec_bbox_max"]),

    )

    tile_center = (float(row["ra_center"]), float(row["dec_center"]))

    return tile_wcs, shape_out, tile_bbox, tile_center

def unwrap_ra(ra_deg, ra_ref_deg):

    ra_deg = np.asarray(ra_deg, dtype=float)

    return (ra_deg - ra_ref_deg + 180.0) % 360.0 - 180.0 + ra_ref_deg

def exposure_radec_bbox(exp_wcs, nx, ny, ra_ref):

    corners = np.array([[0, 0], [nx, 0], [nx, ny], [0, ny]], dtype=float)

    sky = exp_wcs.pixel_to_world(corners[:, 0], corners[:, 1])

    ra = unwrap_ra(sky.ra.deg, ra_ref)

    dec = sky.dec.deg

    return float(ra.min()), float(ra.max()), float(dec.min()), float(dec.max())

def bboxes_overlap(a, b):

    # a,b: (ra_min, ra_max, dec_min, dec_max) in a consistent unwrapped RA frame

    return not (a[1] < b[0] or a[0] > b[1] or a[3] < b[2] or a[2] > b[3])

def coadd_one_tile(

    tile_wcs,

    shape_out,

    tile_bbox,

    tile_center,

    science_paths,

    weight_paths=None,  # optional list aligned with science_paths

    mask_paths=None,  # optional list aligned with science_paths

    sci_ext=1,

    wgt_ext=1,

    msk_ext=1,

    bad_mask_bits=None,  # e.g. integer bitmask of "bad" pixels

    exp_scalefactors=None,  # optional per-exposure photometric scaling to common ZP

):

    """

    Tile-centric coadd.

    - science_paths: list[str]

    - weight_paths: list[str] or None (same length); if None => uniform weights

    - mask_paths: list[str] or None (same length)

    - exp_scalefactors: list[float] or None: multiply science by scale, weight by 1/scale^2

    """

    ra_ref = tile_center[0]  # unwrap frame near tile

    tile_bbox_u = (

        tile_bbox[0],

        tile_bbox[1],

        tile_bbox[2],

        tile_bbox[3],

    )  # already unwrapped near tile center

    print(f"the length of sci_paths: {len(science_paths)}")

    sum_img = np.zeros(shape_out, dtype=np.float32)

    sum_wgt = np.zeros(shape_out, dtype=np.float32)

    nexp = np.zeros(shape_out, dtype=np.uint16)  # coverage count

    out_msk = np.zeros(shape_out, dtype=np.uint32)  # optional combined mask

    if weight_paths is None:

        weight_paths = [None] * len(science_paths)

    if mask_paths is None:

        mask_paths = [None] * len(science_paths)

    if exp_scalefactors is None:

        exp_scalefactors = [1.0] * len(science_paths)

    used = 0

    for sci_path, wgt_path, msk_path, scale in zip(

        science_paths, weight_paths, mask_paths, exp_scalefactors

    ):

        sci, wcs_in, nx, ny = read_fits_image_wcs(sci_path, ext=sci_ext)

        # prefilter: exposure bbox vs tile bbox

        exp_bbox = exposure_radec_bbox(wcs_in, nx, ny, ra_ref)

        if not bboxes_overlap(exp_bbox, tile_bbox_u):

            print(f"{sci_path} not overlap with {tile_bbox}")

            print(f"exp_bbox: {exp_bbox}")

            continue

        else:

            print("Start reprojection...")

        # photometric scaling to common flux system

        # science *= scale ; weight should scale as 1/scale^2

        sci = sci * np.float32(scale)

        warped_sci, footprint = reproject_to_tile(sci, wcs_in, tile_wcs, shape_out)

        # Weight map

        if wgt_path is None:

            warped_wgt = footprint  # simple coverage weight

        else:

            wgt, wcs_w, _, _ = read_fits_image_wcs(wgt_path, ext=wgt_ext)

            # If weight map is in same WCS as science, use wcs_in; else use wcs_w

            warped_wgt, fp_w = reproject_to_tile(wgt, wcs_in, tile_wcs, shape_out)

            warped_wgt = warped_wgt * fp_w

        # apply weight scaling for photometric rescale

        warped_wgt = warped_wgt / (np.float32(scale) ** 2)

        # Mask (optional)

        if msk_path is not None:

            with fits.open(msk_path, memmap=True) as hdul:

                msk = hdul[msk_ext].data

                msk_wcs = WCS(hdul[msk_ext].header)

            # For masks, nearest-neighbor is preferable, but reproject_interp uses interpolation.

            # Practical workaround: treat mask as float and threshold after.

            warped_msk_f, fp_m = reproject_to_tile(

                msk.astype(np.float32), msk_wcs, tile_wcs, shape_out

            )

            warped_msk = (warped_msk_f + 0.5).astype(np.uint32)  # approximate NN

            # update combined mask (bitwise OR where footprint exists)

            good_fp = fp_m > 0

            out_msk[good_fp] |= warped_msk[good_fp]

            if bad_mask_bits is not None:

                bad = (warped_msk & np.uint32(bad_mask_bits)) != 0

                warped_wgt[bad] = 0.0

        # Accumulate (ignore NaNs and zero weight)

        good = np.isfinite(warped_sci) & np.isfinite(warped_wgt) & (warped_wgt > 0)

        sum_img[good] += warped_sci[good] * warped_wgt[good]

        sum_wgt[good] += warped_wgt[good]

        nexp[good] += 1

        used += 1

    coadd = np.full(shape_out, np.nan, dtype=np.float32)

    ok = sum_wgt > 0

    coadd[ok] = sum_img[ok] / sum_wgt[ok]

    return coadd, sum_wgt, nexp, out_msk, used

def write_tile_products(out_prefix, tile_wcs, coadd, wsum, nexp, mask=None):

    hdr = tile_wcs.to_header(relax=True)

    fits.writeto(f"{out_prefix}_coadd.fits", coadd, header=hdr, overwrite=True)

    fits.writeto(f"{out_prefix}_wsum.fits", wsum, header=hdr, overwrite=True)

    fits.writeto(f"{out_prefix}_nexp.fits", nexp, header=hdr, overwrite=True)

    if mask is not None:

        fits.writeto(f"{out_prefix}_mask.fits", mask, header=hdr, overwrite=True)

# -----------------------------

# Example usage

# -----------------------------

if __name__ == "__main__":

    # 1) tile definition

    tile_id = "CSSTT0000123"  # change

    tile_wcs, shape_out, tile_bbox, tile_center = load_tile_from_parquet(

        "csst_tileset.parquet", tile_id

    )

    # 2) your exposure lists

    science_paths = ["exp1_sci.fits", "exp2_sci.fits", "exp3_sci.fits"]

    # optionally:

    weight_paths = [

        "exp1_wgt.fits",

        "exp2_wgt.fits",

        "exp3_wgt.fits",

    ]  # or None

    mask_paths = ["exp1_msk.fits", "exp2_msk.fits", "exp3_msk.fits"]  # or None

    # optional photometric scaling factors to common zeropoint (else all 1.0)

    exp_scalefactors = [1.0, 1.0, 1.0]

    # choose which mask bits are "bad" (example: bits 0, 2, 5)

    bad_mask_bits = (1 << 0) | (1 << 2) | (1 << 5)

    coadd, wsum, nexp, out_msk, used = coadd_one_tile(

        tile_wcs,

        shape_out,

        tile_bbox,

        tile_center,

        science_paths,

        weight_paths=weight_paths,

        mask_paths=mask_paths,

        bad_mask_bits=bad_mask_bits,

        exp_scalefactors=exp_scalefactors,

    )

    print("Used exposures:", used)

    write_tile_products(f"tile_{tile_id}", tile_wcs, coadd, wsum, nexp, mask=out_msk)