From 0580dde9dbd0535531c650cb2091e214daa17a22 Mon Sep 17 00:00:00 2001 From: Jack Date: Tue, 26 Nov 2024 07:04:51 -0800 Subject: [PATCH 01/17] added detection catalog template class --- txpipe/ingest/ssi.py | 173 ++++++++++++++++++++++++++++++++----------- 1 file changed, 128 insertions(+), 45 deletions(-) diff --git a/txpipe/ingest/ssi.py b/txpipe/ingest/ssi.py index 546496b3..fa77b230 100644 --- a/txpipe/ingest/ssi.py +++ b/txpipe/ingest/ssi.py @@ -16,7 +16,88 @@ "coord_dec": "dec", } -class TXIngestSSIGCR(PipelineStage): +class TXIngestSSI(PipelineStage): + """ + Base-Class for ingesting SSI catalogs + + These can be GCR catalogs, or precursor surveys + """ + + name = "TXIngestSSI" + + def create_photometry(self, input_name, output_name, column_names, dummy_columns): + """ + Create and save photometry dataset from an input catalog. + + This method sets up the output HDF5 file structure and processes photometry + data from the input catalog. It maps the input columns to the desired output + column names, creates datasets for photometric magnitudes and their errors, + and iterates over chunks of data to save them in the output file. + + Parameters: + ---------- + input_name : str + name of (FITS) input + + output_name : str + name of (HDF5) output + + column_names : dict + A dictionary mapping the input column names to the desired output column names. + Keys are input column names, and values are the corresponding output names. + + dummy_columns : dict + A dictionary of columns to be added to the output with fixed values. + Keys are the names of the dummy columns, and values are the constant values + to fill those columns. + + """ + #get some basic onfo about the input file + f = self.open_input(input_name) + n = f[1].get_nrows() + dtypes = f[1].get_rec_dtype()[0] + f.close() + + chunk_rows = self.config["chunk_rows"] + + cols = list(column_names.keys()) + + #set up the output file columns + output = self.open_output(output_name) + g = output.create_group("photometry") + for col in cols: + dtype = dtypes[col] + + if "_mag" in col: + #per band + dtype = dtype.subdtype[0] + for b in "griz": + g.create_dataset(column_names[col]+f"_{b}", (n,), dtype=dtype) + + #also create an empty array for the mag errors + #TO DO: compute mag errors from flux errors + g.create_dataset(column_names[col]+f"_err_{b}", (n,), dtype=dtype) + else: + g.create_dataset(column_names[col], (n,), dtype=dtype) + + #iterate over the input file and save to the output columns + for (s, e, data) in self.iterate_fits(input_name, 1, cols, chunk_rows): + print(s,e,n) + for col in cols: + if "_mag" in col: + for iband,b in enumerate("griz"): + g[column_names[col]+f"_{b}"][s:e] = data[col][:,iband] + else: + g[column_names[col]][s:e] = data[col] + + # Set up any dummy columns with sentinal values + # that were not in the original files + for col_name in dummy_columns.keys(): + g.create_dataset(col_name, data=np.full(n,dummy_columns[col_name])) + + output.close() + +class TXIngestSSIGCR(TXIngestSSI): """ Class for ingesting SSI catalogs using GCR @@ -278,7 +359,7 @@ def finalize_output(self, outfile, group, ntot): -class TXIngestSSIMatched(PipelineStage): +class TXIngestSSIMatched(TXIngestSSI): """ Base-stage for ingesting a matched SSI catalog @@ -316,16 +397,6 @@ def run(self): """ print('Ingesting DES Balrog matched catalog') - #get some basic onfo about the input file - f = self.open_input("balrog_matched_catalog") - n = f[1].get_nrows() - dtypes = f[1].get_rec_dtype()[0] - f.close() - - print(f'{n} objects in matched catalog') - - chunk_rows = self.config["chunk_rows"] - #we will only load a subset of columns to save space column_names = { "bal_id": "bal_id", # Unique identifier for object (created during balrog process) @@ -339,46 +410,58 @@ def run(self): "meas_EXTENDED_CLASS_SOF": "EXTENDED_CLASS_SOF", # Star galaxy classifier (0,1=star, 2,3=Galaxy) "meas_FLAGS_GOLD_SOF_ONLY": "FLAGS_GOLD", # Measured flags (short version) } - cols = list(column_names.keys()) + dummy_columns = { + "redshift_true":10.0, + } + + self.create_photometry("balrog_matched_catalog", "matched_ssi_photometry_catalog", + column_names, dummy_columns) - #set up the output file columns - output = self.open_output("matched_ssi_photometry_catalog") - g = output.create_group("photometry") - for col in cols: - dtype = dtypes[col] +class TXIngestSSIDetection(PipelineStage): + """ + Base-stage for ingesting an SSI "detection" catalog + list of injected objects and their detected properties (if any) - if "_mag" in col: - #per band - dtype = dtype.subdtype[0] - for b in "griz": - g.create_dataset(column_names[col]+f"_{b}", (n,), dtype=dtype) + This stage will just read in a file in a given format and output to a + HDF file that TXPIPE can use - #also create an empty array for the mag errors - #TO DO: compute mag errors from flux errors - g.create_dataset(column_names[col]+f"_err_{b}", (n,), dtype=dtype) - else: - g.create_dataset(column_names[col], (n,), dtype=dtype) + """ - #iterate over the input file and save to the output columns - for (s, e, data) in self.iterate_fits("balrog_matched_catalog", 1, cols, chunk_rows): - print(s,e,n) - for col in cols: - if "_mag" in col: - for iband,b in enumerate("griz"): - g[column_names[col]+f"_{b}"][s:e] = data[col][:,iband] - else: - g[column_names[col]][s:e] = data[col] + name = "TXIngestSSIDetection" - # Set up any dummy columns with sentinal values - # that were not in the original files - dummy_columns = { - "redshift_true":10.0, - } - for col_name in dummy_columns.keys(): - g.create_dataset(col_name, data=np.full(n,dummy_columns[col_name])) + outputs = [ + ("detection_ssi_photometry_catalog", HDFFile), + ] - output.close() + config_options = { + "chunk_rows": 100_000, + } + +class TXIngestSSIDetectionDESBalrog(TXIngestSSIDetection): + """ + Class for ingesting an "SSI" "detection" catalog from DES (AKA Balrog) + """ + name = "TXIngestSSIDetectionDESBalrog" + inputs = [ + ("balrog_detection_catalog", FitsFile), + ] + + def run(self): + """ + Run the analysis for this stage. + """ + #we will only load a subset of columns to save space + column_names = { + #add column names here + } + dummy_columns = { + #"redshift_true":10.0, + #add any neede dummy columns here + } + + self.create_photometry("balrog_detection_catalog", "balrog_detection_catalog", + column_names, dummy_columns) \ No newline at end of file From 4b292dfe83db2008dc15cbdc89971f43cfb331c1 Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Wed, 27 Nov 2024 13:47:00 -0500 Subject: [PATCH 02/17] updated docs --- txpipe/ingest/ssi.py | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/txpipe/ingest/ssi.py b/txpipe/ingest/ssi.py index fa77b230..490c379e 100644 --- a/txpipe/ingest/ssi.py +++ b/txpipe/ingest/ssi.py @@ -197,17 +197,13 @@ def run(self): output_file.close() - class TXMatchSSI(PipelineStage): """ - Class for ingesting SSI injection and photometry catalogs + Class for matching an SSI injection catalog and a photometry catalog Default inputs are in TXPipe photometry catalog format - Will perform its own matching between the catalogs to output a - matched SSI catalog for further use - - TO DO: make a separate stage to ingest the matched catalog directly + Outpus a matched SSI catalog for further use """ name = "TXMatchSSI" @@ -357,8 +353,6 @@ def finalize_output(self, outfile, group, ntot): outfile.close() return - - class TXIngestSSIMatched(TXIngestSSI): """ Base-stage for ingesting a matched SSI catalog @@ -379,7 +373,6 @@ class TXIngestSSIMatched(TXIngestSSI): "magnification":0, # magnification label } - class TXIngestSSIMatchedDESBalrog(TXIngestSSIMatched): """ Class for ingesting a matched "SSI" catalog from DES (AKA Balrog) From 0471151d95a35c0b43279123638b89b325315cd6 Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Wed, 27 Nov 2024 14:53:10 -0500 Subject: [PATCH 03/17] split ssi detection info to new catalog --- txpipe/ingest/ssi.py | 48 +++++++++++++++++++++++++++++++++----------- 1 file changed, 36 insertions(+), 12 deletions(-) diff --git a/txpipe/ingest/ssi.py b/txpipe/ingest/ssi.py index 490c379e..b9f13bd5 100644 --- a/txpipe/ingest/ssi.py +++ b/txpipe/ingest/ssi.py @@ -410,7 +410,7 @@ def run(self): self.create_photometry("balrog_matched_catalog", "matched_ssi_photometry_catalog", column_names, dummy_columns) -class TXIngestSSIDetection(PipelineStage): +class TXIngestSSIDetection(TXIngestSSI): """ Base-stage for ingesting an SSI "detection" catalog list of injected objects and their detected properties (if any) @@ -430,7 +430,6 @@ class TXIngestSSIDetection(PipelineStage): "chunk_rows": 100_000, } - class TXIngestSSIDetectionDESBalrog(TXIngestSSIDetection): """ Class for ingesting an "SSI" "detection" catalog from DES (AKA Balrog) @@ -441,20 +440,45 @@ class TXIngestSSIDetectionDESBalrog(TXIngestSSIDetection): inputs = [ ("balrog_detection_catalog", FitsFile), ] + + outputs = [ + ("injection_catalog", HDFFile), + ("ssi_detection_catalog", HDFFile), + ] def run(self): """ - Run the analysis for this stage. + We will split the Balrog "detection" catalog into two catalogs + One containing the injected catalog + The other contains info on whether a given injection has been detected """ + ## Extract the injection catalog + column_names_inj = { + 'bal_id': 'bal_id', # Unique identifier for object (created during balrog process) + 'true_ra': 'ra', # *injected* ra + 'true_dec': 'dec', # *injected* dec + 'true_bdf_mag_deredden': 'inj_mag', # true magnitude (de-reddened) + 'meas_FLAGS_GOLD_SOF_ONLY': 'flags', # measured data flags + 'meas_EXTENDED_CLASS_SOF': 'meas_EXTENDED_CLASS_SOF', # star galaxy separator + } + + self.create_photometry("balrog_detection_catalog", "injection_catalog", + column_names_inj, {}) + + # extract the "detection" file #we will only load a subset of columns to save space - column_names = { - #add column names here - } - dummy_columns = { - #"redshift_true":10.0, - #add any neede dummy columns here - } + column_names_det = { + 'bal_id': 'bal_id', # Unique identifier for object (created during balrog process) + 'detected': 'detected', # 0 or 1. Is there a match between this object in injection_catalog and the ssi_photometry_catalog (search radius of 0.5'') + 'match_flag_0.5_asec': 'match_flag_0.5_asec', # 0,1 or 2. Is there a match between this object in injection_catalog and the ssi_uninjected_photometry_catalog (search radius 0.5''). 1 if detection is lower flux than injection, 2 if brighter + 'match_flag_0.75_asec': 'match_flag_0.75_asec', # as above but with search radius of 0.75'' + 'match_flag_1.0_asec': 'match_flag_1.0_asec', # etc + 'match_flag_1.25_asec': 'match_flag_1.25_asec', + 'match_flag_1.5_asec': 'match_flag_1.5_asec', + 'match_flag_1.75_asec': 'match_flag_1.75_asec', + 'match_flag_2.0_asec': 'match_flag_2.0_asec' + } - self.create_photometry("balrog_detection_catalog", "balrog_detection_catalog", - column_names, dummy_columns) \ No newline at end of file + self.create_photometry("balrog_detection_catalog", "ssi_detection_catalog", + column_names_det, {}) \ No newline at end of file From a30059cff1f6515206ca4685a43aeda9bf911128 Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Wed, 27 Nov 2024 14:56:04 -0500 Subject: [PATCH 04/17] black on ingest ssi --- txpipe/ingest/ssi.py | 183 ++++++++++++++++++++++++------------------- 1 file changed, 102 insertions(+), 81 deletions(-) diff --git a/txpipe/ingest/ssi.py b/txpipe/ingest/ssi.py index b9f13bd5..1798e3f2 100644 --- a/txpipe/ingest/ssi.py +++ b/txpipe/ingest/ssi.py @@ -16,22 +16,23 @@ "coord_dec": "dec", } + class TXIngestSSI(PipelineStage): """ - Base-Class for ingesting SSI catalogs - + Base-Class for ingesting SSI catalogs + These can be GCR catalogs, or precursor surveys """ - + name = "TXIngestSSI" - + def create_photometry(self, input_name, output_name, column_names, dummy_columns): """ Create and save photometry dataset from an input catalog. - This method sets up the output HDF5 file structure and processes photometry - data from the input catalog. It maps the input columns to the desired output - column names, creates datasets for photometric magnitudes and their errors, + This method sets up the output HDF5 file structure and processes photometry + data from the input catalog. It maps the input columns to the desired output + column names, creates datasets for photometric magnitudes and their errors, and iterates over chunks of data to save them in the output file. Parameters: @@ -43,60 +44,61 @@ def create_photometry(self, input_name, output_name, column_names, dummy_columns name of (HDF5) output column_names : dict - A dictionary mapping the input column names to the desired output column names. + A dictionary mapping the input column names to the desired output column names. Keys are input column names, and values are the corresponding output names. dummy_columns : dict - A dictionary of columns to be added to the output with fixed values. - Keys are the names of the dummy columns, and values are the constant values + A dictionary of columns to be added to the output with fixed values. + Keys are the names of the dummy columns, and values are the constant values to fill those columns. - + """ - #get some basic onfo about the input file + # get some basic onfo about the input file f = self.open_input(input_name) n = f[1].get_nrows() dtypes = f[1].get_rec_dtype()[0] f.close() - + chunk_rows = self.config["chunk_rows"] - + cols = list(column_names.keys()) - - #set up the output file columns + + # set up the output file columns output = self.open_output(output_name) g = output.create_group("photometry") for col in cols: dtype = dtypes[col] if "_mag" in col: - #per band + # per band dtype = dtype.subdtype[0] for b in "griz": - g.create_dataset(column_names[col]+f"_{b}", (n,), dtype=dtype) + g.create_dataset(column_names[col] + f"_{b}", (n,), dtype=dtype) - #also create an empty array for the mag errors - #TO DO: compute mag errors from flux errors - g.create_dataset(column_names[col]+f"_err_{b}", (n,), dtype=dtype) + # also create an empty array for the mag errors + # TO DO: compute mag errors from flux errors + g.create_dataset(column_names[col] + f"_err_{b}", (n,), dtype=dtype) else: g.create_dataset(column_names[col], (n,), dtype=dtype) - #iterate over the input file and save to the output columns - for (s, e, data) in self.iterate_fits(input_name, 1, cols, chunk_rows): - print(s,e,n) + # iterate over the input file and save to the output columns + for s, e, data in self.iterate_fits(input_name, 1, cols, chunk_rows): + print(s, e, n) for col in cols: if "_mag" in col: - for iband,b in enumerate("griz"): - g[column_names[col]+f"_{b}"][s:e] = data[col][:,iband] + for iband, b in enumerate("griz"): + g[column_names[col] + f"_{b}"][s:e] = data[col][:, iband] else: g[column_names[col]][s:e] = data[col] - # Set up any dummy columns with sentinal values + # Set up any dummy columns with sentinal values # that were not in the original files for col_name in dummy_columns.keys(): - g.create_dataset(col_name, data=np.full(n,dummy_columns[col_name])) + g.create_dataset(col_name, data=np.full(n, dummy_columns[col_name])) output.close() + class TXIngestSSIGCR(TXIngestSSI): """ Class for ingesting SSI catalogs using GCR @@ -138,7 +140,7 @@ def run(self): sys.path.insert(0, self.config["GCRcatalog_path"]) import GCRCatalogs - #attempt to ingest three types of catalog + # attempt to ingest three types of catalog catalogs_labels = [ "injection_catalog", "ssi_photometry_catalog", @@ -175,11 +177,15 @@ def run(self): group[column_names[q]] = group[q] except KeyError: # skip quantities that are missing - warnings.warn(f"quantity {q} was missing from the GCRCatalog object") + warnings.warn( + f"quantity {q} was missing from the GCRCatalog object" + ) continue except TypeError: - warnings.warn(f"Quantity {q} coud not be saved as it has a data type not recognised by hdf5") + warnings.warn( + f"Quantity {q} coud not be saved as it has a data type not recognised by hdf5" + ) # convert fluxes to mags using txpipe/utils/conversion.py bands = "ugrizy" @@ -197,6 +203,7 @@ def run(self): output_file.close() + class TXMatchSSI(PipelineStage): """ Class for matching an SSI injection catalog and a photometry catalog @@ -269,10 +276,14 @@ def match_cats(self): ) out_start = 0 - for ichunk, (in_start,in_end,data) in enumerate(self.iterate_hdf("ssi_photometry_catalog", "photometry", ["ra","dec"], batch_size)): + for ichunk, (in_start, in_end, data) in enumerate( + self.iterate_hdf( + "ssi_photometry_catalog", "photometry", ["ra", "dec"], batch_size + ) + ): phot_coord = SkyCoord( - ra=data["ra"]*u.degree, - dec=data["dec"]*u.degree, + ra=data["ra"] * u.degree, + dec=data["dec"] * u.degree, ) idx, d2d, d3d = phot_coord.match_to_catalog_sky(inj_coord) @@ -353,11 +364,12 @@ def finalize_output(self, outfile, group, ntot): outfile.close() return + class TXIngestSSIMatched(TXIngestSSI): """ Base-stage for ingesting a matched SSI catalog - This stage will just read in a file in a given format and output to a + This stage will just read in a file in a given format and output to a HDF file that TXPIPE can use """ @@ -370,9 +382,10 @@ class TXIngestSSIMatched(TXIngestSSI): config_options = { "chunk_rows": 100_000, - "magnification":0, # magnification label + "magnification": 0, # magnification label } + class TXIngestSSIMatchedDESBalrog(TXIngestSSIMatched): """ Class for ingesting a matched "SSI" catalog from DES (AKA Balrog) @@ -383,39 +396,44 @@ class TXIngestSSIMatchedDESBalrog(TXIngestSSIMatched): inputs = [ ("balrog_matched_catalog", FitsFile), ] - + def run(self): """ Run the analysis for this stage. """ - print('Ingesting DES Balrog matched catalog') + print("Ingesting DES Balrog matched catalog") - #we will only load a subset of columns to save space + # we will only load a subset of columns to save space column_names = { - "bal_id": "bal_id", # Unique identifier for object (created during balrog process) - "true_bdf_mag_deredden": "inj_mag", # Magnitude of the original deep field object, dereddened - "true_id": "inj_id", # Original coadd_obj_id of deep field object - "meas_id": "id", # Coadd_object_id of injection - "meas_ra": "ra", # measured RA of the injection - "meas_dec": "dec", # measured DEC of the injection - "meas_cm_mag_deredden": "mag", # measured magnitude of the injection - "meas_cm_T": "cm_T", # measured size parameter T (x^2+y^2) - "meas_EXTENDED_CLASS_SOF": "EXTENDED_CLASS_SOF", # Star galaxy classifier (0,1=star, 2,3=Galaxy) - "meas_FLAGS_GOLD_SOF_ONLY": "FLAGS_GOLD", # Measured flags (short version) - } + "bal_id": "bal_id", # Unique identifier for object (created during balrog process) + "true_bdf_mag_deredden": "inj_mag", # Magnitude of the original deep field object, dereddened + "true_id": "inj_id", # Original coadd_obj_id of deep field object + "meas_id": "id", # Coadd_object_id of injection + "meas_ra": "ra", # measured RA of the injection + "meas_dec": "dec", # measured DEC of the injection + "meas_cm_mag_deredden": "mag", # measured magnitude of the injection + "meas_cm_T": "cm_T", # measured size parameter T (x^2+y^2) + "meas_EXTENDED_CLASS_SOF": "EXTENDED_CLASS_SOF", # Star galaxy classifier (0,1=star, 2,3=Galaxy) + "meas_FLAGS_GOLD_SOF_ONLY": "FLAGS_GOLD", # Measured flags (short version) + } dummy_columns = { - "redshift_true":10.0, - } - - self.create_photometry("balrog_matched_catalog", "matched_ssi_photometry_catalog", - column_names, dummy_columns) + "redshift_true": 10.0, + } + + self.create_photometry( + "balrog_matched_catalog", + "matched_ssi_photometry_catalog", + column_names, + dummy_columns, + ) + class TXIngestSSIDetection(TXIngestSSI): """ Base-stage for ingesting an SSI "detection" catalog list of injected objects and their detected properties (if any) - This stage will just read in a file in a given format and output to a + This stage will just read in a file in a given format and output to a HDF file that TXPIPE can use """ @@ -429,7 +447,8 @@ class TXIngestSSIDetection(TXIngestSSI): config_options = { "chunk_rows": 100_000, } - + + class TXIngestSSIDetectionDESBalrog(TXIngestSSIDetection): """ Class for ingesting an "SSI" "detection" catalog from DES (AKA Balrog) @@ -445,7 +464,7 @@ class TXIngestSSIDetectionDESBalrog(TXIngestSSIDetection): ("injection_catalog", HDFFile), ("ssi_detection_catalog", HDFFile), ] - + def run(self): """ We will split the Balrog "detection" catalog into two catalogs @@ -455,30 +474,32 @@ def run(self): ## Extract the injection catalog column_names_inj = { - 'bal_id': 'bal_id', # Unique identifier for object (created during balrog process) - 'true_ra': 'ra', # *injected* ra - 'true_dec': 'dec', # *injected* dec - 'true_bdf_mag_deredden': 'inj_mag', # true magnitude (de-reddened) - 'meas_FLAGS_GOLD_SOF_ONLY': 'flags', # measured data flags - 'meas_EXTENDED_CLASS_SOF': 'meas_EXTENDED_CLASS_SOF', # star galaxy separator + "bal_id": "bal_id", # Unique identifier for object (created during balrog process) + "true_ra": "ra", # *injected* ra + "true_dec": "dec", # *injected* dec + "true_bdf_mag_deredden": "inj_mag", # true magnitude (de-reddened) + "meas_FLAGS_GOLD_SOF_ONLY": "flags", # measured data flags + "meas_EXTENDED_CLASS_SOF": "meas_EXTENDED_CLASS_SOF", # star galaxy separator } - self.create_photometry("balrog_detection_catalog", "injection_catalog", - column_names_inj, {}) - - # extract the "detection" file - #we will only load a subset of columns to save space + self.create_photometry( + "balrog_detection_catalog", "injection_catalog", column_names_inj, {} + ) + + # Extract the "detection" file + # We will only load a subset of columns to save space column_names_det = { - 'bal_id': 'bal_id', # Unique identifier for object (created during balrog process) - 'detected': 'detected', # 0 or 1. Is there a match between this object in injection_catalog and the ssi_photometry_catalog (search radius of 0.5'') - 'match_flag_0.5_asec': 'match_flag_0.5_asec', # 0,1 or 2. Is there a match between this object in injection_catalog and the ssi_uninjected_photometry_catalog (search radius 0.5''). 1 if detection is lower flux than injection, 2 if brighter - 'match_flag_0.75_asec': 'match_flag_0.75_asec', # as above but with search radius of 0.75'' - 'match_flag_1.0_asec': 'match_flag_1.0_asec', # etc - 'match_flag_1.25_asec': 'match_flag_1.25_asec', - 'match_flag_1.5_asec': 'match_flag_1.5_asec', - 'match_flag_1.75_asec': 'match_flag_1.75_asec', - 'match_flag_2.0_asec': 'match_flag_2.0_asec' + "bal_id": "bal_id", # Unique identifier for object (created during balrog process) + "detected": "detected", # 0 or 1. Is there a match between this object in injection_catalog and the ssi_photometry_catalog (search radius of 0.5'') + "match_flag_0.5_asec": "match_flag_0.5_asec", # 0,1 or 2. Is there a match between this object in injection_catalog and the ssi_uninjected_photometry_catalog (search radius 0.5''). 1 if detection is lower flux than injection, 2 if brighter + "match_flag_0.75_asec": "match_flag_0.75_asec", # as above but with search radius of 0.75'' + "match_flag_1.0_asec": "match_flag_1.0_asec", # etc + "match_flag_1.25_asec": "match_flag_1.25_asec", + "match_flag_1.5_asec": "match_flag_1.5_asec", + "match_flag_1.75_asec": "match_flag_1.75_asec", + "match_flag_2.0_asec": "match_flag_2.0_asec", } - - self.create_photometry("balrog_detection_catalog", "ssi_detection_catalog", - column_names_det, {}) \ No newline at end of file + + self.create_photometry( + "balrog_detection_catalog", "ssi_detection_catalog", column_names_det, {} + ) From be3409760eaef23d25d0aa66e6708b066ed5802b Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Fri, 29 Nov 2024 12:30:13 -0500 Subject: [PATCH 05/17] re-structure class inheritance for SSI ingest stages --- txpipe/ingest/__init__.py | 2 +- txpipe/ingest/ssi.py | 181 ++++++++++++++++++++++++-------------- 2 files changed, 114 insertions(+), 69 deletions(-) diff --git a/txpipe/ingest/__init__.py b/txpipe/ingest/__init__.py index 64b85659..016c0c5a 100644 --- a/txpipe/ingest/__init__.py +++ b/txpipe/ingest/__init__.py @@ -2,4 +2,4 @@ from .mocks import TXCosmoDC2Mock, TXGaussianSimsMock from .gcr import TXMetacalGCRInput, TXIngestStars from .redmagic import TXIngestRedmagic -from .ssi import TXIngestSSIGCR, TXMatchSSI, TXIngestSSIMatched, TXIngestSSIMatchedDESBalrog \ No newline at end of file +from .ssi import TXIngestSSIGCR, TXMatchSSI, TXIngestSSIDESBalrog, TXIngestSSIMatchedDESBalrog, TXIngestSSIDetectionDESBalrog \ No newline at end of file diff --git a/txpipe/ingest/ssi.py b/txpipe/ingest/ssi.py index 1798e3f2..056c8b8d 100644 --- a/txpipe/ingest/ssi.py +++ b/txpipe/ingest/ssi.py @@ -16,7 +16,6 @@ "coord_dec": "dec", } - class TXIngestSSI(PipelineStage): """ Base-Class for ingesting SSI catalogs @@ -61,43 +60,80 @@ def create_photometry(self, input_name, output_name, column_names, dummy_columns chunk_rows = self.config["chunk_rows"] - cols = list(column_names.keys()) + # set up the output file columns + output, g = self.setup_output(output_name, column_names, dtypes, n) + + # iterate over the input file and save to the output columns + self.add_columns(g, input_name, column_names, chunk_rows, n) + + # Set up any dummy columns with sentinal values + # that were not in the original files + for col_name in dummy_columns.keys(): + g.create_dataset(col_name, data=np.full(n, dummy_columns[col_name])) + + output.close() + + def setup_output(self, output_name, column_names, dtypes, n): + """ + Set up the output HDF5 file structure. - # set up the output file columns + Parameters + ---------- + output_name : str + The name of the output HDF5 file. + + column_names : dict + dict of column names to include in the output file. + + dtypes : dict + A dictionary mapping column names to their corresponding data types. + + n : int + The total number of rows in the dataset. + + Returns + ------- + tuple + A tuple containing the open HDF5 output file object and the "photometry" group. + """ + cols = list(column_names.keys()) output = self.open_output(output_name) g = output.create_group("photometry") for col in cols: dtype = dtypes[col] + g.create_dataset(column_names[col], (n,), dtype=dtype) + + return output, g - if "_mag" in col: - # per band - dtype = dtype.subdtype[0] - for b in "griz": - g.create_dataset(column_names[col] + f"_{b}", (n,), dtype=dtype) + def add_columns(self, g, input_name, column_names, chunk_rows, n): + """ + Add data to the HDF5 output file in chunks. - # also create an empty array for the mag errors - # TO DO: compute mag errors from flux errors - g.create_dataset(column_names[col] + f"_err_{b}", (n,), dtype=dtype) - else: - g.create_dataset(column_names[col], (n,), dtype=dtype) + This method reads chunks of data from the input file and writes them + to the corresponding datasets in the output file. - # iterate over the input file and save to the output columns - for s, e, data in self.iterate_fits(input_name, 1, cols, chunk_rows): - print(s, e, n) - for col in cols: - if "_mag" in col: - for iband, b in enumerate("griz"): - g[column_names[col] + f"_{b}"][s:e] = data[col][:, iband] - else: - g[column_names[col]][s:e] = data[col] + Parameters + ---------- + g : h5py.Group + The HDF5 group where data will be written. - # Set up any dummy columns with sentinal values - # that were not in the original files - for col_name in dummy_columns.keys(): - g.create_dataset(col_name, data=np.full(n, dummy_columns[col_name])) + input_name : str + The name of the input file (e.g., a FITS file). - output.close() + column_names : dict + Dict of column names to read from the input file. + + chunk_rows : int + Number of rows to process in each chunk. + n : int + Total number of rows in the dataset. + """ + cols = list(column_names.keys()) + for s, e, data in self.iterate_fits(input_name, 1, cols, chunk_rows): + print(s, e, n) + for col in cols: + g[column_names[col]][s:e] = data[col] class TXIngestSSIGCR(TXIngestSSI): """ @@ -203,7 +239,6 @@ def run(self): output_file.close() - class TXMatchSSI(PipelineStage): """ Class for matching an SSI injection catalog and a photometry catalog @@ -364,29 +399,56 @@ def finalize_output(self, outfile, group, ntot): outfile.close() return - -class TXIngestSSIMatched(TXIngestSSI): +class TXIngestSSIDESBalrog(TXIngestSSI): """ - Base-stage for ingesting a matched SSI catalog - - This stage will just read in a file in a given format and output to a - HDF file that TXPIPE can use - + Base-stage for ingesting a DES SSI catalog AKA "Balrog" """ - name = "TXIngestSSIMatched" + name = "TXIngestSSIDESBalrog" - outputs = [ - ("matched_ssi_photometry_catalog", HDFFile), - ] + def setup_output(self, output_name, column_names, dtypes, n): + """ + For balrog, we need to include if statements to catch the 2D data entries + and add an extendedness column + """ + cols = list(column_names.keys()) + output = self.open_output(output_name) + g = output.create_group("photometry") + for col in cols: + dtype = dtypes[col] - config_options = { - "chunk_rows": 100_000, - "magnification": 0, # magnification label - } + if dtype.subdtype is not None: #this is a multi-dimentional column + assert dtype.subdtype[1]==(4,) #We are assuming this entry is a 2D array with 4 columns (corresponding to griz) + dtype = dtype.subdtype[0] + for b in "griz": + g.create_dataset(column_names[col] + f"_{b}", (n,), dtype=dtype) + if "_mag" in col: #if the column is also a mag + # also create an empty array for the mag errors + # TO DO: compute mag errors from flux errors + g.create_dataset(column_names[col] + f"_err_{b}", (n,), dtype=dtype) + else: + g.create_dataset(column_names[col], (n,), dtype=dtype) + + return output, g -class TXIngestSSIMatchedDESBalrog(TXIngestSSIMatched): + def add_columns(self, g, input_name, column_names, chunk_rows, n): + """ + For balrog, we need to include if statements to catch the 2D data entries + and add a extendedness column + """ + cols = list(column_names.keys()) + for s, e, data in self.iterate_fits(input_name, 1, cols, chunk_rows): + print(s, e, n) + for col in cols: + if len(data[col].shape) == 2: + assert data[col].shape[1] == 4 + for iband, b in enumerate("griz"): + g[column_names[col] + f"_{b}"][s:e] = data[col][:, iband] + else: + g[column_names[col]][s:e] = data[col] + +class TXIngestSSIMatchedDESBalrog(TXIngestSSIDESBalrog): """ Class for ingesting a matched "SSI" catalog from DES (AKA Balrog) """ @@ -397,6 +459,10 @@ class TXIngestSSIMatchedDESBalrog(TXIngestSSIMatched): ("balrog_matched_catalog", FitsFile), ] + outputs = [ + ("matched_ssi_photometry_catalog", HDFFile), + ] + def run(self): """ Run the analysis for this stage. @@ -412,6 +478,7 @@ def run(self): "meas_ra": "ra", # measured RA of the injection "meas_dec": "dec", # measured DEC of the injection "meas_cm_mag_deredden": "mag", # measured magnitude of the injection + #"meas_cm_max_flux_s2n": "snr", # measured S2N of the injection #TODO: add this column "meas_cm_T": "cm_T", # measured size parameter T (x^2+y^2) "meas_EXTENDED_CLASS_SOF": "EXTENDED_CLASS_SOF", # Star galaxy classifier (0,1=star, 2,3=Galaxy) "meas_FLAGS_GOLD_SOF_ONLY": "FLAGS_GOLD", # Measured flags (short version) @@ -426,30 +493,8 @@ def run(self): column_names, dummy_columns, ) - - -class TXIngestSSIDetection(TXIngestSSI): - """ - Base-stage for ingesting an SSI "detection" catalog - list of injected objects and their detected properties (if any) - - This stage will just read in a file in a given format and output to a - HDF file that TXPIPE can use - - """ - - name = "TXIngestSSIDetection" - - outputs = [ - ("detection_ssi_photometry_catalog", HDFFile), - ] - - config_options = { - "chunk_rows": 100_000, - } - - -class TXIngestSSIDetectionDESBalrog(TXIngestSSIDetection): + +class TXIngestSSIDetectionDESBalrog(TXIngestSSIDESBalrog): """ Class for ingesting an "SSI" "detection" catalog from DES (AKA Balrog) """ From 129e1d3dbf7e49ee28ed3439d4b6817fa86ddebc Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Fri, 29 Nov 2024 12:33:40 -0500 Subject: [PATCH 06/17] removed empty mag_err columns --- txpipe/ingest/ssi.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/txpipe/ingest/ssi.py b/txpipe/ingest/ssi.py index 056c8b8d..b5c56790 100644 --- a/txpipe/ingest/ssi.py +++ b/txpipe/ingest/ssi.py @@ -422,11 +422,6 @@ def setup_output(self, output_name, column_names, dtypes, n): dtype = dtype.subdtype[0] for b in "griz": g.create_dataset(column_names[col] + f"_{b}", (n,), dtype=dtype) - - if "_mag" in col: #if the column is also a mag - # also create an empty array for the mag errors - # TO DO: compute mag errors from flux errors - g.create_dataset(column_names[col] + f"_err_{b}", (n,), dtype=dtype) else: g.create_dataset(column_names[col], (n,), dtype=dtype) @@ -485,6 +480,10 @@ def run(self): } dummy_columns = { "redshift_true": 10.0, + "meas_err_g":-99., + "meas_err_r":-99., + "meas_err_i":-99., + "meas_err_z":-99., } self.create_photometry( From ba51baff14673641c5e6e3461dd57263fea14de1 Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Fri, 29 Nov 2024 12:44:22 -0500 Subject: [PATCH 07/17] added extendedness column to matched SSI catalog --- txpipe/ingest/ssi.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/txpipe/ingest/ssi.py b/txpipe/ingest/ssi.py index b5c56790..ed1c366b 100644 --- a/txpipe/ingest/ssi.py +++ b/txpipe/ingest/ssi.py @@ -409,7 +409,7 @@ class TXIngestSSIDESBalrog(TXIngestSSI): def setup_output(self, output_name, column_names, dtypes, n): """ For balrog, we need to include if statements to catch the 2D data entries - and add an extendedness column + and possibly add an extendedness column """ cols = list(column_names.keys()) output = self.open_output(output_name) @@ -424,13 +424,17 @@ def setup_output(self, output_name, column_names, dtypes, n): g.create_dataset(column_names[col] + f"_{b}", (n,), dtype=dtype) else: g.create_dataset(column_names[col], (n,), dtype=dtype) + + if col == "meas_EXTENDED_CLASS_SOF": + #also create an "extendedness" column + g.create_dataset("extendedness", (n,), dtype=dtype) return output, g def add_columns(self, g, input_name, column_names, chunk_rows, n): """ For balrog, we need to include if statements to catch the 2D data entries - and add a extendedness column + and possibly add a extendedness column """ cols = list(column_names.keys()) for s, e, data in self.iterate_fits(input_name, 1, cols, chunk_rows): @@ -442,6 +446,13 @@ def add_columns(self, g, input_name, column_names, chunk_rows, n): g[column_names[col] + f"_{b}"][s:e] = data[col][:, iband] else: g[column_names[col]][s:e] = data[col] + + if col == "meas_EXTENDED_CLASS_SOF": + # "meas_EXTENDED_CLASS_SOF" is (0 or 1) for star, (2 or 3) for galaxy + # extendedness is 0 for star, 1 for galaxy + extendedness = np.where((data[col] == 2) | (data[col] == 3), 1, 0) + g[column_names[col]][s:e] = extendedness + class TXIngestSSIMatchedDESBalrog(TXIngestSSIDESBalrog): """ From a523c6f35fbe5066f57b6de9a102710510e691ab Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Fri, 29 Nov 2024 12:47:08 -0500 Subject: [PATCH 08/17] added SNR column to matched SSI catalog --- txpipe/ingest/ssi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/ingest/ssi.py b/txpipe/ingest/ssi.py index ed1c366b..b745c028 100644 --- a/txpipe/ingest/ssi.py +++ b/txpipe/ingest/ssi.py @@ -484,7 +484,7 @@ def run(self): "meas_ra": "ra", # measured RA of the injection "meas_dec": "dec", # measured DEC of the injection "meas_cm_mag_deredden": "mag", # measured magnitude of the injection - #"meas_cm_max_flux_s2n": "snr", # measured S2N of the injection #TODO: add this column + "meas_cm_max_flux_s2n": "snr", # measured S2N of the injection #TODO: add this column "meas_cm_T": "cm_T", # measured size parameter T (x^2+y^2) "meas_EXTENDED_CLASS_SOF": "EXTENDED_CLASS_SOF", # Star galaxy classifier (0,1=star, 2,3=Galaxy) "meas_FLAGS_GOLD_SOF_ONLY": "FLAGS_GOLD", # Measured flags (short version) From 85632f66ca51b32b3c29234cc93f28fa2abecba2 Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Wed, 4 Dec 2024 16:45:24 -0500 Subject: [PATCH 09/17] added option to rotate maps 180deg when plotting --- txpipe/data_types/types.py | 11 ++++++++++- txpipe/map_plots.py | 5 +++-- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/txpipe/data_types/types.py b/txpipe/data_types/types.py index c5d6d6bb..c012c9e5 100755 --- a/txpipe/data_types/types.py +++ b/txpipe/data_types/types.py @@ -274,12 +274,20 @@ def write_map(self, map_name, pixel, value, metadata): subgroup.create_dataset("pixel", data=pixel) subgroup.create_dataset("value", data=value) - def plot_healpix(self, map_name, view="cart", **kwargs): + def plot_healpix(self, map_name, view="cart", rot180=False, **kwargs): import healpy import numpy as np m, pix, nside = self.read_healpix(map_name, return_all=True) lon, lat = healpy.pix2ang(nside, pix, lonlat=True) + if rot180: #(optional) rotate 180 degrees in the lon direction + lon += 180 + lon[lon > 360.] -= 360. + pix_rot = healpy.ang2pix(nside, lon, lat, lonlat=True) + m_rot = np.ones(healpy.nside2npix(nside))*healpy.UNSEEN + m_rot[pix_rot] = m[pix] + m = m_rot + pix = pix_rot npix = healpy.nside2npix(nside) if len(pix) == 0: print(f"Empty map {map_name}") @@ -291,6 +299,7 @@ def plot_healpix(self, map_name, view="cart", **kwargs): lon_range = [lon[w].min() - 0.1, lon[w].max() + 0.1] lat_range = [lat[w].min() - 0.1, lat[w].max() + 0.1] lat_range = np.clip(lat_range, -90, 90) + lon_range = np.clip(lon_range, 0, 360.) m[m == 0] = healpy.UNSEEN title = kwargs.pop("title", map_name) if view == "cart": diff --git a/txpipe/map_plots.py b/txpipe/map_plots.py index ee3ac433..08257a9f 100644 --- a/txpipe/map_plots.py +++ b/txpipe/map_plots.py @@ -42,6 +42,7 @@ class TXMapPlots(PipelineStage): config_options = { # can also set moll "projection": "cart", + "rot180": False, "debug": False, } @@ -122,11 +123,11 @@ def aux_lens_plots(self): # Depth plots with self.open_output("depth_map", wrapper=True, figsize=(5, 5)) as fig: - m.plot("depth/depth", view=self.config["projection"]) + m.plot("depth/depth", view=self.config["projection"], rot180=self.config["rot180"]) # Bright objects with self.open_output("bright_object_map", wrapper=True, figsize=(5, 5)) as fig: - m.plot("bright_objects/count", view=self.config["projection"]) + m.plot("bright_objects/count", view=self.config["projection"], rot180=self.config["rot180"]) def source_plots(self): import matplotlib.pyplot as plt From 2b993eaf7217f5f730a15c468b27b1b4b9389e09 Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Wed, 4 Dec 2024 16:46:49 -0500 Subject: [PATCH 10/17] added 180deg rotation option to all map plots --- txpipe/map_plots.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/txpipe/map_plots.py b/txpipe/map_plots.py index 08257a9f..7b8a552d 100644 --- a/txpipe/map_plots.py +++ b/txpipe/map_plots.py @@ -97,7 +97,7 @@ def aux_source_plots(self): for i in range(flag_max): plt.subplot(1, flag_max, i + 1) f = 2**i - m.plot(f"flags/flag_{f}", view=self.config["projection"]) + m.plot(f"flags/flag_{f}", view=self.config["projection"], rot180=self.config["rot180"]) fig.close() # PSF plots - 2 x n, for g1 and g2 @@ -105,9 +105,9 @@ def aux_source_plots(self): _, axes = plt.subplots(2, nbin_source, squeeze=False, num=fig.file.number) for i in range(nbin_source): plt.sca(axes[0, i]) - m.plot(f"psf/g1_{i}", view=self.config["projection"]) + m.plot(f"psf/g1_{i}", view=self.config["projection"], rot180=self.config["rot180"]) plt.sca(axes[1, i]) - m.plot(f"psf/g2_{i}", view=self.config["projection"]) + m.plot(f"psf/g2_{i}", view=self.config["projection"], rot180=self.config["rot180"]) fig.close() def aux_lens_plots(self): @@ -149,11 +149,11 @@ def source_plots(self): for i in range(nbin_source): # g1 plt.sca(axes[0, i]) - m.plot(f"g1_{i}", view=self.config["projection"], min=-0.1, max=0.1) + m.plot(f"g1_{i}", view=self.config["projection"], rot180=self.config["rot180"], min=-0.1, max=0.1) # g2 plt.sca(axes[1, i]) - m.plot(f"g2_{i}", view=self.config["projection"], min=-0.1, max=0.1) + m.plot(f"g2_{i}", view=self.config["projection"], rot180=self.config["rot180"], min=-0.1, max=0.1) fig.close() def lens_plots(self): @@ -174,9 +174,9 @@ def lens_plots(self): for i in range(nbin_lens): plt.sca(axes[0, i]) - m.plot(f"ngal_{i}", view=self.config["projection"]) + m.plot(f"ngal_{i}", view=self.config["projection"], rot180=self.config["rot180"]) plt.sca(axes[1, i]) - rho.plot(f"delta_{i}", view=self.config["projection"]) + rho.plot(f"delta_{i}", view=self.config["projection"], rot180=self.config["rot180"]) fig.close() def mask_plots(self): @@ -185,5 +185,5 @@ def mask_plots(self): m = self.open_input("mask", wrapper=True) fig = self.open_output("mask_map", wrapper=True, figsize=(5, 5)) - m.plot("mask", view=self.config["projection"]) + m.plot("mask", view=self.config["projection"], rot180=self.config["rot180"]) fig.close() From 2d30c7c5af5d6c6d79552f13fbba0475da669cff Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Wed, 4 Dec 2024 16:49:17 -0500 Subject: [PATCH 11/17] added SNR to SSI matched catalog input --- txpipe/ingest/ssi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/ingest/ssi.py b/txpipe/ingest/ssi.py index b745c028..cfad6a1d 100644 --- a/txpipe/ingest/ssi.py +++ b/txpipe/ingest/ssi.py @@ -484,7 +484,7 @@ def run(self): "meas_ra": "ra", # measured RA of the injection "meas_dec": "dec", # measured DEC of the injection "meas_cm_mag_deredden": "mag", # measured magnitude of the injection - "meas_cm_max_flux_s2n": "snr", # measured S2N of the injection #TODO: add this column + "meas_cm_max_flux_s2n": "snr", # measured S2N of the injection "meas_cm_T": "cm_T", # measured size parameter T (x^2+y^2) "meas_EXTENDED_CLASS_SOF": "EXTENDED_CLASS_SOF", # Star galaxy classifier (0,1=star, 2,3=Galaxy) "meas_FLAGS_GOLD_SOF_ONLY": "FLAGS_GOLD", # Measured flags (short version) From 19c797b6da59caec31d037c31ae49673502a4fc8 Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Mon, 9 Dec 2024 12:33:44 -0500 Subject: [PATCH 12/17] added aux ssi maps stage --- txpipe/auxiliary_maps.py | 97 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 97 insertions(+) diff --git a/txpipe/auxiliary_maps.py b/txpipe/auxiliary_maps.py index 4c6afa8f..17391e04 100644 --- a/txpipe/auxiliary_maps.py +++ b/txpipe/auxiliary_maps.py @@ -256,3 +256,100 @@ def run(self): print(len(pix)) print(len(depth[pix])) f.write_map("depth/depth", pix, depth[pix], metadata) + +class TXAuxiliarySSIMaps(TXBaseMaps): + """ + Generate auxiliary maps from SSI catalogs + + This class generates maps of: + - depth (measured magnitude) + - depth (true magnitude) + """ + name = "TXAuxiliarySSIMaps" + dask_parallel = True + inputs = [ + ("injection_catalog", HDFFile), # for injection locations + ("matched_ssi_photometry_catalog", HDFFile), + ] + outputs = [ + ("aux_ssi_maps", MapsFile), + ] + + ################### + ################## + + config_options = { + "block_size": 0, + "depth_band": "i", # Make depth maps for this band + "snr_threshold": 10.0, # The S/N value to generate maps for (e.g. 5 for 5-sigma depth) + "snr_delta": 1.0, # The range threshold +/- delta is used for finding objects at the boundary + } + + def run(self): + # Import dask and alias it as 'da' + _, da = import_dask() + + + # Retrieve configuration parameters + block_size = self.config["block_size"] + if block_size == 0: + block_size = "auto" + band = self.config["depth_band"] + + # Open the input photometry catalog file. + # We can't use a "with" statement because we need to keep the file open + # while we're using dask. + f = self.open_input("matched_ssi_photometry_catalog", wrapper=True) + + # Load photometry data into dask arrays. + # This is lazy in dask, so we're not actually loading the data here. + ra = da.from_array(f.file["photometry/ra"], block_size) + block_size = ra.chunksize + dec = da.from_array(f.file["photometry/dec"], block_size) + snr = da.from_array(f.file[f"photometry/snr_{band}"], block_size) + mag_meas = da.from_array(f.file[f"photometry/mag_{band}"], block_size) + mag_true = da.from_array(f.file[f"photometry/inj_mag_{band}"], block_size) + + # Choose the pixelization scheme based on the configuration. + # Might need to review this to make sure we use the same scheme everywhere + pixel_scheme = choose_pixelization(**self.config) + + # Initialize a dictionary to store the maps. + # To start with this is all lazy too, until we call compute + maps = {} + + # Create depth maps using dask and measured magnitudes + pix2, count_map, depth_map, depth_var = make_dask_depth_map( + ra, dec, mag_meas, snr, self.config["snr_threshold"], self.config["snr_delta"], pixel_scheme) + maps["depth/depth_meas"] = (pix2, depth_map[pix2]) + maps["depth/depth_meas_count"] = (pix2, count_map[pix2]) + maps["depth/depth_meas_var"] = (pix2, depth_var[pix2]) + + + # Create depth maps using daskand true magnitudes + pix2, count_map, depth_map, depth_var = make_dask_depth_map( + ra, dec, mag_true, snr, self.config["snr_threshold"], self.config["snr_delta"], pixel_scheme) + maps["depth/depth_true"] = (pix2, depth_map[pix2]) + maps["depth/depth_true_count"] = (pix2, count_map[pix2]) + maps["depth/depth_true_var"] = (pix2, depth_var[pix2]) + + maps, = da.compute(maps) + + # Prepare metadata for the maps. Copy the pixelization-related + # configuration options only here + metadata = { + key: self.config[key] + for key in map_config_options + if key in self.config + } + # Then add the other configuration options + metadata["depth_band"] = band + metadata["depth_snr_threshold"] = self.config["snr_threshold"] + metadata["depth_snr_delta"] = self.config["snr_delta"] + metadata.update(pixel_scheme.metadata) + + # Write the output maps to the output file + with self.open_output("aux_ssi_maps", wrapper=True) as out: + for map_name, (pix, m) in maps.items(): + out.write_map(map_name, pix, m, metadata) + out.file['maps'].attrs.update(metadata) From 20b588edcff254a1e8b62a346bb04f92d299affa Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Mon, 9 Dec 2024 12:39:57 -0500 Subject: [PATCH 13/17] added aux_ssi_plots to map plotting stage --- txpipe/map_plots.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/txpipe/map_plots.py b/txpipe/map_plots.py index 7b8a552d..c00b5631 100644 --- a/txpipe/map_plots.py +++ b/txpipe/map_plots.py @@ -28,6 +28,7 @@ class TXMapPlots(PipelineStage): ("mask", MapsFile), ("aux_source_maps", MapsFile), ("aux_lens_maps", MapsFile), + ("aux_ssi_maps", MapsFile), ] outputs = [ @@ -38,6 +39,8 @@ class TXMapPlots(PipelineStage): ("psf_map", PNGFile), ("mask_map", PNGFile), ("bright_object_map", PNGFile), + ("depth_ssi_meas_map", PNGFile), + ("depth_ssi_true_map", PNGFile), ] config_options = { # can also set moll @@ -58,6 +61,7 @@ def run(self): methods = [ self.aux_source_plots, self.aux_lens_plots, + self.aux_ssi_plots, self.source_plots, self.lens_plots, self.mask_plots, @@ -129,6 +133,25 @@ def aux_lens_plots(self): with self.open_output("bright_object_map", wrapper=True, figsize=(5, 5)) as fig: m.plot("bright_objects/count", view=self.config["projection"], rot180=self.config["rot180"]) + def aux_ssi_plots(self): + import matplotlib.pyplot as plt + if self.get_input("aux_ssi_maps") == "none": + with self.open_output("depth_ssi_meas_map", wrapper=True) as f: + pass + with self.open_output("depth_ssi_true_map", wrapper=True) as f: + pass + return + + m = self.open_input("aux_ssi_maps", wrapper=True) + + # Depth plots (measured magnitude) + with self.open_output("depth_ssi_meas_map", wrapper=True, figsize=(5, 5)) as fig: + m.plot("depth/depth_meas", view=self.config["projection"], rot180=self.config["rot180"]) + + # Depth plots (true magnitude) + with self.open_output("depth_ssi_true_map", wrapper=True, figsize=(5, 5)) as fig: + m.plot("depth/depthtrue", view=self.config["projection"], rot180=self.config["rot180"]) + def source_plots(self): import matplotlib.pyplot as plt From 0287ef8e7160ee7ee2ca372347455d65793e8e2d Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Mon, 9 Dec 2024 13:01:51 -0500 Subject: [PATCH 14/17] fixed typo in depth_true name --- txpipe/map_plots.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/map_plots.py b/txpipe/map_plots.py index c00b5631..e8fd12a3 100644 --- a/txpipe/map_plots.py +++ b/txpipe/map_plots.py @@ -150,7 +150,7 @@ def aux_ssi_plots(self): # Depth plots (true magnitude) with self.open_output("depth_ssi_true_map", wrapper=True, figsize=(5, 5)) as fig: - m.plot("depth/depthtrue", view=self.config["projection"], rot180=self.config["rot180"]) + m.plot("depth/depth_true", view=self.config["projection"], rot180=self.config["rot180"]) def source_plots(self): import matplotlib.pyplot as plt From 00e5585ef0eaff5ba7e0f2385b29516eaa0cb757 Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Tue, 10 Dec 2024 16:45:20 -0500 Subject: [PATCH 15/17] added sii depth example --- examples/ssi/config_ssi_depth.yml | 21 ++++++++++++ examples/ssi/pipeline_ssi_depth.yml | 51 +++++++++++++++++++++++++++++ 2 files changed, 72 insertions(+) create mode 100755 examples/ssi/config_ssi_depth.yml create mode 100755 examples/ssi/pipeline_ssi_depth.yml diff --git a/examples/ssi/config_ssi_depth.yml b/examples/ssi/config_ssi_depth.yml new file mode 100755 index 00000000..b09b45fe --- /dev/null +++ b/examples/ssi/config_ssi_depth.yml @@ -0,0 +1,21 @@ +# Values in this section are accessible to all the different stages. +# They can be overridden by individual stages though. +global: + # This is read by many stages that read complete + # catalog data, and tells them how many rows to read + # at once + chunk_rows: 100_000 + # These mapping options are also read by a range of stages + pixelization: healpix + nside: 256 + sparse: true # Generate sparse maps - faster if using small areas + + +TXIngestSSIDetectionDESBalrog: + name: TXIngestSSIDetectionDESBalrog + +TXMapPlots: + name: TXMapPlots + projection: cart + rot180: True + diff --git a/examples/ssi/pipeline_ssi_depth.yml b/examples/ssi/pipeline_ssi_depth.yml new file mode 100755 index 00000000..57be64f2 --- /dev/null +++ b/examples/ssi/pipeline_ssi_depth.yml @@ -0,0 +1,51 @@ +# Stages to run +stages: + - name: TXIngestSSIDetectionDESBalrog + - name: TXIngestSSIMatchedDESBalrog + - name: TXAuxiliarySSIMaps + - name: TXMapPlots + +modules: txpipe + +# Where to put outputs +output_dir: data/example/outputs_ssi_des_depth/ + +# How to run the pipeline: mini, parsl, or cwl +launcher: + name: mini + interval: 1.0 + +# Where to run the pipeline: cori-interactive, cori-batch, or local +site: + name: local + max_threads: 2 + +# configuration settings +config: examples/ssi/config_ssi_depth.yml + +inputs: + # See README for paths to download these files + + #NERSC + #balrog_detection_catalog: /global/cfs/cdirs/lsst/groups/WL/users/jelvinpo/balrog_txpipe_example/balrog_detection_catalog_sof_run2_v1.4.fits + + #local + balrog_detection_catalog: /Users/jackelvinpoole/DES/cats/y3/balrog/balrog_detection_catalog_sof_run2a_v1.4.fits + balrog_matched_catalog: /Users/jackelvinpoole/DES/cats/y3/balrog/balrog_matched_catalog_sof_run2a_v1.4.fits + + source_maps: None + lens_maps: None + density_maps: None + mask: None + aux_source_maps: None + aux_lens_maps: none + + fiducial_cosmology: data/fiducial_cosmology.yml + +# if supported by the launcher, restart the pipeline where it left off +# if interrupted +resume: false +# where to put output logs for individual stages +log_dir: data/example/outputs_ssi_des_mag/logs/ +# where to put an overall parsl pipeline log +pipeline_log: data/example/outputs_ssi_des_mag/log.txt From 07877f3446aea9acfa595ef3a693584dd5c2e5e6 Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Tue, 10 Dec 2024 17:12:32 -0500 Subject: [PATCH 16/17] moved SSI map plotting to separate sub-class --- txpipe/map_plots.py | 84 ++++++++++++++++++++++++++++++++------------- 1 file changed, 61 insertions(+), 23 deletions(-) diff --git a/txpipe/map_plots.py b/txpipe/map_plots.py index e8fd12a3..0c30686d 100644 --- a/txpipe/map_plots.py +++ b/txpipe/map_plots.py @@ -28,7 +28,6 @@ class TXMapPlots(PipelineStage): ("mask", MapsFile), ("aux_source_maps", MapsFile), ("aux_lens_maps", MapsFile), - ("aux_ssi_maps", MapsFile), ] outputs = [ @@ -39,8 +38,6 @@ class TXMapPlots(PipelineStage): ("psf_map", PNGFile), ("mask_map", PNGFile), ("bright_object_map", PNGFile), - ("depth_ssi_meas_map", PNGFile), - ("depth_ssi_true_map", PNGFile), ] config_options = { # can also set moll @@ -61,7 +58,6 @@ def run(self): methods = [ self.aux_source_plots, self.aux_lens_plots, - self.aux_ssi_plots, self.source_plots, self.lens_plots, self.mask_plots, @@ -133,25 +129,6 @@ def aux_lens_plots(self): with self.open_output("bright_object_map", wrapper=True, figsize=(5, 5)) as fig: m.plot("bright_objects/count", view=self.config["projection"], rot180=self.config["rot180"]) - def aux_ssi_plots(self): - import matplotlib.pyplot as plt - if self.get_input("aux_ssi_maps") == "none": - with self.open_output("depth_ssi_meas_map", wrapper=True) as f: - pass - with self.open_output("depth_ssi_true_map", wrapper=True) as f: - pass - return - - m = self.open_input("aux_ssi_maps", wrapper=True) - - # Depth plots (measured magnitude) - with self.open_output("depth_ssi_meas_map", wrapper=True, figsize=(5, 5)) as fig: - m.plot("depth/depth_meas", view=self.config["projection"], rot180=self.config["rot180"]) - - # Depth plots (true magnitude) - with self.open_output("depth_ssi_true_map", wrapper=True, figsize=(5, 5)) as fig: - m.plot("depth/depth_true", view=self.config["projection"], rot180=self.config["rot180"]) - def source_plots(self): import matplotlib.pyplot as plt @@ -210,3 +187,64 @@ def mask_plots(self): fig = self.open_output("mask_map", wrapper=True, figsize=(5, 5)) m.plot("mask", view=self.config["projection"], rot180=self.config["rot180"]) fig.close() + +class TXMapPlotsSSI(TXMapPlots): + """ + Make plots of all the available maps that use SSI inputs + + This makes plots of: + - depth (meas mag) + - depth (true mag) + """ + name = "TXMapPlotsSSI" + + inputs = [ + ("aux_ssi_maps", MapsFile), + ] + + outputs = [ + ("depth_ssi_meas_map", PNGFile), + ("depth_ssi_true_map", PNGFile), + ] + + def run(self): + # PSF tests + import matplotlib + + matplotlib.use("agg") + import matplotlib.pyplot as plt + + # Plot from each file separately, just + # to organize this file a bit + methods = [ + self.aux_ssi_plots, + ] + + # We don't want this to fail if some maps are missing. + for m in methods: + try: + m() + except: + if self.config["debug"]: + raise + sys.stderr.write(f"Failed to make maps with method {m.__name__}") + + + def aux_ssi_plots(self): + import matplotlib.pyplot as plt + if self.get_input("aux_ssi_maps") == "none": + with self.open_output("depth_ssi_meas_map", wrapper=True) as f: + pass + with self.open_output("depth_ssi_true_map", wrapper=True) as f: + pass + return + + m = self.open_input("aux_ssi_maps", wrapper=True) + + # Depth plots (measured magnitude) + with self.open_output("depth_ssi_meas_map", wrapper=True, figsize=(5, 5)) as fig: + m.plot("depth/depth_meas", view=self.config["projection"], rot180=self.config["rot180"]) + + # Depth plots (true magnitude) + with self.open_output("depth_ssi_true_map", wrapper=True, figsize=(5, 5)) as fig: + m.plot("depth/depth_true", view=self.config["projection"], rot180=self.config["rot180"]) From 6824a9cbf6f39ae37b0b6bc287698fbff746a7cb Mon Sep 17 00:00:00 2001 From: elvinpoole Date: Tue, 10 Dec 2024 17:17:55 -0500 Subject: [PATCH 17/17] modified ssi depth example to run on nersc --- examples/ssi/pipeline_ssi_depth.yml | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/examples/ssi/pipeline_ssi_depth.yml b/examples/ssi/pipeline_ssi_depth.yml index 57be64f2..15ab0069 100755 --- a/examples/ssi/pipeline_ssi_depth.yml +++ b/examples/ssi/pipeline_ssi_depth.yml @@ -3,7 +3,7 @@ stages: - name: TXIngestSSIDetectionDESBalrog - name: TXIngestSSIMatchedDESBalrog - name: TXAuxiliarySSIMaps - - name: TXMapPlots + - name: TXMapPlotsSSI modules: txpipe @@ -26,19 +26,9 @@ config: examples/ssi/config_ssi_depth.yml inputs: # See README for paths to download these files - #NERSC - #balrog_detection_catalog: /global/cfs/cdirs/lsst/groups/WL/users/jelvinpo/balrog_txpipe_example/balrog_detection_catalog_sof_run2_v1.4.fits - - #local - balrog_detection_catalog: /Users/jackelvinpoole/DES/cats/y3/balrog/balrog_detection_catalog_sof_run2a_v1.4.fits - balrog_matched_catalog: /Users/jackelvinpoole/DES/cats/y3/balrog/balrog_matched_catalog_sof_run2a_v1.4.fits - - source_maps: None - lens_maps: None - density_maps: None - mask: None - aux_source_maps: None - aux_lens_maps: none + #NERSC locations of public DES balrog data + balrog_detection_catalog: /global/cfs/cdirs/lsst/groups/WL/users/jelvinpo/balrog_txpipe_example/balrog_detection_catalog_sof_run2_v1.4.fits + balrog_matched_catalog: /global/cfs/cdirs/lsst/groups/WL/users/jelvinpo/balrog_txpipe_example/balrog_matched_catalog_sof_run2_v1.4.fits fiducial_cosmology: data/fiducial_cosmology.yml