Merge pull request #378 from LSSTDESC/ssi_depth_maps

SSI depth maps

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
+

examples/ssi/pipeline_ssi_depth.yml

@@ -0,0 +1,41 @@
+# Stages to run
+stages:
+  - name: TXIngestSSIDetectionDESBalrog
+  - name: TXIngestSSIMatchedDESBalrog
+  - name: TXAuxiliarySSIMaps
+  - name: TXMapPlotsSSI
+
+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 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
+
+# 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

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)

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":

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
+from .ssi import TXIngestSSIGCR, TXMatchSSI, TXIngestSSIDESBalrog, TXIngestSSIMatchedDESBalrog, TXIngestSSIDetectionDESBalrog