diff --git a/astra.yaml b/astra.yaml index e566a1253..8204587ac 100644 --- a/astra.yaml +++ b/astra.yaml @@ -124,8 +124,7 @@ decisions: label: Postage-stamp size shared by vignets, ngmix stamps and PSF models rationale: >- Every stamp width is 51 px (about 9.5 arcsec): the tile VIGNET (cut - from the tile image by read_ext_sexcat in committed data, by SExtractor - under tile_detection = sextractor), the multi-epoch vignetmaker stamps + by the tile SExtractor run), the multi-epoch vignetmaker stamps ngmix fits, the exposure SExtractor VIGNET PSFEx trains on, and the PSFEx model stamp. Only the first two are coupled: ngmix lays the tile VIGNET's -1e30 neighbour markers over each epoch stamp pixel for pixel, @@ -134,7 +133,6 @@ decisions: Values: default_noimaflags.param#VIGNET = 51; default.param#VIGNET = 51; - READ_EXT_SEXCAT_RUNNER.VIGNET_SIZE = 51; VIGNETMAKER_RUNNER_RUN_1.STAMP_SIZE = 51; VIGNETMAKER_RUNNER_RUN_2.STAMP_SIZE = 51; PSF_SIZE = 51. @@ -157,10 +155,9 @@ decisions: FSCALE before the joint fit, which puts them on zero-point 30 only if FSCALE = 10^(-0.4 (PHOTZP - 30)). Nothing in the repo checks this, but it held to 0.02% on a sampled exposure (2114045p). ngmix magnitudes use - MAG_ZP = 30. In committed data the catalogue's tile MAG_AUTO is the - UNIONS catalogue's own; the fixed tile zero-point of 30 - (ZP_FROM_HEADER=False) applies only under tile_detection = sextractor, - and assumes the MegaPipe stacks are calibrated to 30. + MAG_ZP = 30. The tile SExtractor run uses a fixed zero-point of 30 + (ZP_FROM_HEADER=False), which assumes the MegaPipe stacks are + calibrated to 30. Values: MAG_ZEROPOINT = 30.0; config_tile_Sx.ini#SEXTRACTOR_RUNNER.ZP_FROM_HEADER = False; @@ -174,9 +171,9 @@ decisions: header_everywhere: label: Per-image header zero-points on tiles too description: >- - ZP_FROM_HEADER=True for tile SExtractor (tile_detection = - sextractor only); changes tile magnitudes only where a tile's header - zero-point differs from 30. ngmix MAG_ZP stays 30. + ZP_FROM_HEADER=True for tile SExtractor; changes tile magnitudes + only where a tile's header zero-point differs from 30. ngmix MAG_ZP + stays 30. prior_insights: des_psf_blacklist: @@ -247,9 +244,8 @@ analyses: shape_measurement.epoch_masked_fraction_cut). Tiles have no flag image, and ShapePipe's tile SExtractor run requests no IMAFLAGS_ISO (detection.detection_source_mode). Neighbour pixels are masked - separately, from the tile segmentation map - (detection.catalogue_neighbour_marking, - shape_measurement.blend_handling). Sky-fixed masks never touch + separately, from the -1e30 markers SExtractor writes into the tile + VIGNET (shape_measurement.blend_handling). Sky-fixed masks never touch pixels: a star halo changes no stamp. Values: config_exp_psfex.ini#SEXTRACTOR_RUNNER.FLAG_IMAGE = True; @@ -342,12 +338,11 @@ analyses: detection: description: >- Object detection on r-band tiles (the galaxy sample) and on - single-exposure CCDs (PSF-star candidates, always SExtractor). The tile - sample is Gwyn's UNIONS tile catalogue (committed for data) or - ShapePipe's own SExtractor run (image simulations), per tile_detection. - Every default_tile.sex / config_tile_Sx.ini setting below applies only - to the SExtractor run, which follows the MegaPipe parameters of that - catalogue; exposures keep ShapePipe's stock values. + single-exposure CCDs (PSF-star candidates), both by ShapePipe's own + SExtractor runs. The tile run follows the MegaPipe parameters of Gwyn's + UNIONS tile catalogue, and in committed data its detections are joined + to that catalogue (tile_detection); exposures keep ShapePipe's stock + values. inputs: - id: tile_stack type: data @@ -362,41 +357,67 @@ analyses: [tile_detection, detection_threshold_policy, deblending_policy, background_model, weight_map_usage, zero_weight_interpolation, detection_source_mode, - epoch_membership_ccd_bounds, catalogue_neighbour_marking, - photometry_parameters, + epoch_membership_ccd_bounds, photometry_parameters, spurious_detection_cleaning, blend_photometry_mask_type, saturation_level] decisions: tile_detection: - label: Source of the tile galaxy sample + label: Object list and NUMBER of the tile galaxy sample rationale: >- - Real data takes its tile sample from Gwyn's UNIONS per-tile - catalogue, converted to the sexcat the chain reads with its own - NUMBER kept, so shape, photometry and photo-z catalogues share one - object list and one ID (TILE_UNIQUE_ID). PSF stars still come from - ShapePipe's exposure-level SExtractor run, because external star - catalogues are too shallow. detection_source_mode and the tile - SExtractor settings of the other decisions here govern only the - sextractor option; epoch_membership_ccd_bounds applies to both. - Image simulations, which have no UNIONS catalogue, use sextractor. + Both options run ShapePipe's tile SExtractor, so data and image + simulations share one detection and measurement path: windowed + positions, VIGNET neighbour markers and every SExtractor column. + In real data each detection is then paired with its mutual nearest + neighbour within 1 px in Gwyn's UNIONS per-tile catalogue and takes + its NUMBER, so shape, photometry and photo-z catalogues share one + object list and one ID (TILE_UNIQUE_ID); unpaired detections leave + the catalogue. The catalogue is MegaPipe's SExtractor run on the + same DR6 image with the configuration of default_tile.sex + (detection_threshold_policy), so pairs agree to 4e-5 px at the + median. On tile 186.307 every one of the 41,201 detections pairs; + over eight DR6 tiles spanning a rich cluster, a bright star's halo, + low latitude and the survey edge, at least 99.7% of detections and + 99.2% of catalogue objects pair. The catalogue objects with no + detection (0-0.7%, mostly fainter than MAG_AUTO 24.5) are + deblended children that ShapePipe's SExtractor leaves merged, + concentrated around very large objects; they have photometry but + no shape. The join needs the tile image to be the catalogue's + release: on the DR5 image of a tile only ~95% of detections pair. + The run stops when the unpaired detections, or the unpaired + catalogue objects, exceed both 2% of their side and 20 objects: a + margin over deblending near large objects, a stop for other pixels. PSF stars come from the exposure-level + SExtractor run in both options. Image simulations, which have no + UNIONS catalogue, keep SExtractor's own NUMBER. Values: - input_types.data.tile_detection = unions_catalogue. + input_types.data.tile_detection = unions_catalogue; + config_tile_Sx.ini#SEXTRACTOR_RUNNER.MATCH_RADIUS = 1.0; + config_tile_Sx.ini#SEXTRACTOR_RUNNER.MATCH_MIN_FRACTION = 0.98; + config_tile_Sx.ini#SEXTRACTOR_RUNNER.MATCH_TOLERATED_UNPAIRED = 20. default: unions_catalogue options: unions_catalogue: - label: Gwyn's UNIONS per-tile catalogue, converted in place + label: Tile SExtractor joined to Gwyn's UNIONS per-tile catalogue sextractor: - label: SExtractor on the tile image + label: Tile SExtractor with its own NUMBER detection_threshold_policy: label: Detection significance, minimum area, matched filter rationale: >- - Tile SExtractor (tile_detection = sextractor, i.e. image - simulations) uses the MegaPipe tile-catalogue threshold, minimum + Tile SExtractor uses the MegaPipe tile-catalogue threshold, minimum area and filter, so that its galaxy sample reproduces the UNIONS - catalogue committed data adopts. Run on the DR6 image of tile - 186.307, it reproduces the DR6 catalogue to 99.89% completeness + catalogue whose NUMBER committed data adopts (tile_detection). Run + on the DR6 image of tile 186.307, it reproduces the DR6 catalogue to 99.89% completeness and 100% purity with identical segmentation pixels; the 45 misses - are faint deblended children of bright parents (#929). This + are faint deblended children of bright parents (#929). The pixel + stack (MEMORY_PIXSTACK) holds the largest object on every tested + tile, a bright star's halo included, so no object is truncated. + Unexplained: on tiles with such a large object (a cD galaxy, a + bright star's halo) MAG_AUTO and FLUX_RADIUS differ from the + catalogue by more than 1% for 2-11% of objects, most of them past + the large object in SExtractor's scan order, and on one + low-latitude tile for 2% of objects spread over the tile. The + catalogue's run must differ in a setting not recorded here + (memory, version); shapes use ShapePipe's measurements and + photometry the catalogue's, so for those objects the two differ. This validates the configuration on data pixels, not detection on simulated images. The 7x7 Gaussian of FWHM 3 px is near the CFIS average seeing of 0.65 arcsec (about 3.5 px at @@ -410,6 +431,7 @@ analyses: default_tile.sex#ANALYSIS_THRESH = 1.0; default_tile.sex#DETECT_MINAREA = 3; default_tile.sex#FILTER = Y; + default_tile.sex#MEMORY_PIXSTACK = 3000000; config_tile_Sx.ini#SEXTRACTOR_RUNNER.DOT_CONV_FILE = $SP_CONFIG/gauss_3.0_7x7.conv; config_exp_psfex.ini#SEXTRACTOR_RUNNER.DOT_SEX_FILE = @@ -573,11 +595,9 @@ analyses: absent SExtractor falls back to its built-in level (50000 ADU per its documentation). Delivered exposures carry the card (every CCD of 2114045p: 65535) and split_exp keeps each CCD's header; a sampled - tile (CFIS.186.307) has 9558.7, which applies under tile_detection = - sextractor. The PSF-star masks require FLAGS == 0, so the level can - remove bright stars from the PSF sample. In committed data the final - catalogue's FLAGS column is the UNIONS catalogue's own, not this - pass's. Changing the card or pinning a fixed level requires checking + tile (CFIS.186.307) has 9558.7. The PSF-star masks require FLAGS == + 0, so the level can remove bright stars from the PSF sample. + Changing the card or pinning a fixed level requires checking the resulting bright-star selection. Values: default_exp.sex#SATUR_KEY = SATURATE; @@ -600,9 +620,7 @@ analyses: recorded. On the exposure pass MAG_AUTO is the axis of the star-selection magnitude windows and FLUX_AUTO is PSFEx's photometric normalisation, so a different Kron factor moves stars - across those cuts. The tile values apply under tile_detection = - sextractor; in committed data the catalogue's MAG_AUTO is the UNIONS - catalogue's own. The reported apertures and flux-fraction + across those cuts. The reported apertures and flux-fraction measurements must match the simulated catalogue they are compared with. Values: @@ -624,7 +642,12 @@ analyses: column list: tiles have no instrument flag image and no detection coadd exists, so detection sees every tile pixel and the tile catalogue carries no IMAFLAGS_ISO; final_cat.param, the merge's - exact allow-list, does not request it. + exact allow-list, does not request it. Without a detection image + sextractor_runner calls SExtractor in single-image mode, which on + tile 186.307 reproduces the UNIONS catalogue's MAG_AUTO to 0.01 mag + for all but 0.03% of objects. Dual-image mode with the same image + twice is not equivalent: it changes the background, MAG_AUTO and + FLUX_RADIUS of a fifth of the objects (#936). Values: SEXTRACTOR_RUNNER.DETECTION_IMAGE = False; SEXTRACTOR_RUNNER.FLAG_IMAGE = False. @@ -650,13 +673,10 @@ analyses: admits its upper endpoint, column 2080. [HARDCODED] The comparisons are strict in make_post_process; CCD_SIZE sets only the bounds. A WCS inversion failure skips that CCD, lowering N_EPOCH, which caps - how many exposures enter each galaxy's multi-epoch fit. Both - tile_detection options run this post-processing. + how many exposures enter each galaxy's multi-epoch fit. Values: SEXTRACTOR_RUNNER.CCD_SIZE = 33,2080,1,4612; - SEXTRACTOR_RUNNER.MAKE_POST_PROCESS = True; - READ_EXT_SEXCAT_RUNNER.CCD_SIZE = 33,2080,1,4612; - READ_EXT_SEXCAT_RUNNER.MAKE_POST_PROCESS = True. + SEXTRACTOR_RUNNER.MAKE_POST_PROCESS = True. default: trimmed_bounds_33_2080 options: trimmed_bounds_33_2080: @@ -664,36 +684,6 @@ analyses: inclusive_bounds: label: Same bounds, inclusive description: Not implemented; make_post_process uses strict comparisons. - catalogue_neighbour_marking: - label: Neighbour pixels in the catalogue path's VIGNET - rationale: >- - ngmix masks a neighbour only where the tile VIGNET is -1e30 (flag - 2**10: zero weight, noise-filled under noisefill), and SExtractor - writes -1e30 on neighbours' footprints and off the image. The - unions_catalogue converter reproduces that from the catalogue's own - r-band segmentation map (CFIS..r.seg.fits.fz, on the tile's - pixel grid), so both tile_detection options mask the same way. - Segmentation labels are not the catalogue NUMBER; each object claims - the footprint under its centre pixel (the VIGNET centre pixel), which - on 202.301 matches 36,064 of 36,065 objects with none shared. The map - is relabelled to NUMBER (unclaimed footprints -1; an object on sky or - on a claimed footprint gets a 3-pixel-radius disc of its own), and - VIGNET pixels whose relabelled value is neither 0 nor the object's - NUMBER become -1e30. Without it, catalogue-mode ngmix fits neighbour - light: in pilot smk-g10, NGMIX_MAG_NOSHEAR was over 1 mag brighter - than MAG_AUTO for 6.1% of objects (2.7% with SExtractor, g9). - Values: - READ_EXT_SEXCAT_RUNNER.SEGMENTATION = True. - default: segmentation_map - options: - segmentation_map: - label: -1e30 on other objects' segmentation footprints, as SExtractor - none: - label: Image pixels only; ngmix masks no neighbours - excluded: true - excluded_reason: >- - Leaves neighbour light in the fit, unlike the sextractor option - it replaces; the cause of the smk-g10 magnitude outliers. prior_insights: guinot22_sextractor_params: claim: >- @@ -815,8 +805,9 @@ analyses: object_position_columns: label: Windowed centroids (XWIN/YWIN) define every position rationale: >- - PSF interpolation sites, tile and multi-epoch stamp centres, and the - catalogue position all use SExtractor's windowed centroid: in pixels + PSF interpolation sites, tile and multi-epoch stamp centres, epoch + membership and the catalogue position all use the tile SExtractor + run's windowed centroid, in data and image simulations alike: in pixels for tile stamps and exposure-side PSF validation, in world coordinates for tile-side PSF interpolation and multi-epoch stamps. Windowed, isophotal and model centroids differ systematically for @@ -835,6 +826,15 @@ analyses: options: xwin_windowed: label: Windowed centroids everywhere + isophotal_barycentre: + label: Isophotal barycentre (X_IMAGE/Y_IMAGE) + excluded: true + excluded_reason: >- + Noisier than the windowed centroid and pulled by asymmetric or + deblend-clipped footprints; on the DR6 image of tile 186.307 + it sits 0.16 px (median) and 1.6 px (p99) from it for + 20 < MAG_AUTO < 24.5, and the image simulations measure the + windowed one. stamp_positioning_and_padding: label: Nearest-pixel stamp extraction with zero padding rationale: >- @@ -2047,7 +2047,7 @@ analyses: excluded_reason: >- Not needed in the pipeline: the cut needs only each object's TILE_ID and position, and sp_validation applies it downstream. - Leaving it out keeps every UNIONS catalogue row in the shape + Leaving it out keeps every tile-catalogue row in the shape catalogue. failure_sentinels: label: Objects without shape measurements kept, with sentinel values diff --git a/src/shapepipe/modules/read_ext_sexcat_package/__init__.py b/src/shapepipe/modules/read_ext_sexcat_package/__init__.py deleted file mode 100644 index cc65bc9fd..000000000 --- a/src/shapepipe/modules/read_ext_sexcat_package/__init__.py +++ /dev/null @@ -1 +0,0 @@ -"""READ EXTERNAL SEXCAT PACKAGE.""" diff --git a/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py b/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py deleted file mode 100644 index de9af6511..000000000 --- a/src/shapepipe/modules/read_ext_sexcat_package/read_ext_sexcat.py +++ /dev/null @@ -1,358 +0,0 @@ -"""READ EXTERNAL SEXCAT. - -Convert an ASCII SExtractor-format catalogue to FITS-LDAC for use in -ShapePipe tile processing. - -:Author: Martin Kilbinger - -""" - -import numpy as np -from astropy.io import ascii as asc -from astropy.io import fits - - -def _build_ldac_imhead(img_header): - """Build LDAC_IMHEAD extension from a FITS image header. - - The LDAC_IMHEAD binary table has one row and one column containing - the image header cards as an array of 80-character strings. The - HISTORY cards are the ones used by ``make_post_process`` to identify - the single exposures that contributed to the tile. - - Parameters - ---------- - img_header : astropy.io.fits.Header - Header of the tile image FITS file - - Returns - ------- - astropy.io.fits.BinTableHDU - LDAC_IMHEAD extension - - """ - cards = [str(card).ljust(80)[:80] for card in img_header.cards] - n_cards = len(cards) - header_str = "".join(cards) - - arr = np.array([header_str.encode()], dtype=f"S{len(header_str)}") - col = fits.Column( - name="Field Header Card", - format=f"{len(header_str)}A", - array=arr, - ) - hdu = fits.BinTableHDU.from_columns([col]) - hdu.header["TDIM1"] = f"(80,{n_cards})" - hdu.name = "LDAC_IMHEAD" - return hdu - - -# SExtractor's VIGNET value for pixels that are not the object's: off the -# image and on a neighbour's segmentation footprint. ngmix flags these pixels -# (tile VIGNET == -1e30) and noise-fills them at zero weight. -BIG = -1e30 - -# The label of a footprint no catalogue object claims in the relabelled -# segmentation map. Negative, so it never collides with a NUMBER; VIGNET -# marking only asks "self or not self", so neighbours need no identity. -NEIGHBOUR_LABEL = -1 - -# Radius, in pixels, of the disc of its own NUMBER painted for an object with -# no footprint of its own (on sky, or on a footprint claimed by another). -FALLBACK_RADIUS = 3 - - -def _centre_pixels(x_pos, y_pos): - """0-based (column, row) of the pixel holding each 1-based position. - - @sc [decision:preparation.stamp_positioning_and_padding] - """ - col = np.rint(np.asarray(x_pos, dtype=float)).astype(np.int64) - 1 - row = np.rint(np.asarray(y_pos, dtype=float)).astype(np.int64) - 1 - return col, row - - -def relabel_seg(seg, number, x_image, y_image, fallback_radius=FALLBACK_RADIUS): - """Relabel a segmentation map into a catalogue's NUMBER. - - Each object claims, in NUMBER order and first come first served, the - footprint its centre pixel falls in; that footprint becomes its NUMBER. - Footprints nobody claims become ``NEIGHBOUR_LABEL``; sky stays 0. An - object on sky or on an already-claimed footprint gets a disc of radius - ``fallback_radius`` of its NUMBER, and every object on the image ends - holding its own centre pixel (the lower NUMBER wins a shared pixel). - - Parameters - ---------- - seg : numpy.ndarray - Segmentation map, 0 for sky and one positive label per footprint - number : array_like - Catalogue ``NUMBER`` - x_image, y_image : array_like - 1-based pixel positions on the grid of ``seg`` - fallback_radius : int, optional - Radius of the disc painted for an object without a footprint - - Returns - ------- - numpy.ndarray - Relabelled map, ``int32`` - dict - Counts of the claim outcomes ``matched``, ``unclaimed`` (on sky), - ``shared`` and ``off_image``, which partition the catalogue, plus - ``shared_pixel`` (objects rounding to a lower NUMBER's pixel) - - @sc [decision:detection.catalogue_neighbour_marking] - """ - number = np.asarray(number) - col, row = _centre_pixels(x_image, y_image) - n_row, n_col = seg.shape - inside = (col >= 0) & (col < n_col) & (row >= 0) & (row < n_row) - order = np.argsort(number, kind="stable") - - counts = dict(matched=0, unclaimed=0, shared=0, off_image=0, - shared_pixel=0) - table = np.full(max(int(seg.max()), 0) + 1, NEIGHBOUR_LABEL, np.int32) - table[0] = 0 - claimed, fallback = set(), [] - for i in order: - if not inside[i]: - counts["off_image"] += 1 - continue - label = int(seg[row[i], col[i]]) - if label == 0 or label in claimed: - counts["unclaimed" if label == 0 else "shared"] += 1 - fallback.append(i) - else: - claimed.add(label) - table[label] = number[i] - counts["matched"] += 1 - out = table[np.clip(seg, 0, None)] - - r = np.arange(-fallback_radius, fallback_radius + 1) - disc = np.argwhere(np.add.outer(r**2, r**2) <= fallback_radius**2) - disc -= fallback_radius - for i in fallback: - out[np.clip(row[i] + disc[:, 0], 0, n_row - 1), - np.clip(col[i] + disc[:, 1], 0, n_col - 1)] = number[i] - - # Centres last, so no disc can take an object's centre away. - taken = set() - for i in order[inside[order]]: - pixel = (row[i], col[i]) - if pixel in taken: - counts["shared_pixel"] += 1 - continue - taken.add(pixel) - out[pixel] = number[i] - return out, counts - - -def _extract_vignets(image_data, x_pos, y_pos, stamp_size, seg=None, - number=None): - """Extract postage stamps from a tile image array. - - For each object position, a ``stamp_size x stamp_size`` cutout is - extracted, centred on the pixel holding the position. Pixels off the - image are set to ``BIG``, as SExtractor does. With a segmentation map in - the catalogue's numbering (:func:`relabel_seg`), pixels on any footprint - other than the object's own are set to ``BIG`` too, which is how - SExtractor's VIGNET marks neighbours. - - Parameters - ---------- - image_data : numpy.ndarray - 2-D tile image array, shape ``(ny, nx)`` - x_pos : array_like - X pixel positions, 1-based (SExtractor convention) - y_pos : array_like - Y pixel positions, 1-based (SExtractor convention) - stamp_size : int - Side length of the square postage stamp (should be odd) - seg : numpy.ndarray, optional - Segmentation map on the image grid, labelled with ``number`` - number : array_like, optional - Catalogue ``NUMBER``, required with ``seg`` - - Returns - ------- - numpy.ndarray - Array of shape ``(n_obj, stamp_size, stamp_size)``, dtype float32 - - @sc [decision:detection.catalogue_neighbour_marking] - """ - ny, nx = image_data.shape - half = stamp_size // 2 - col, row = _centre_pixels(x_pos, y_pos) - vignets = np.full((len(col), stamp_size, stamp_size), BIG, np.float32) - seg_stamp = np.zeros((stamp_size, stamp_size), np.int32) - - for i, (xi, yi) in enumerate(zip(col, row)): - x0, y0 = xi - half, yi - half - xc0, xc1 = max(0, x0), min(nx, x0 + stamp_size) - yc0, yc1 = max(0, y0), min(ny, y0 + stamp_size) - if xc0 >= xc1 or yc0 >= yc1: - continue - inner = np.s_[yc0 - y0:yc1 - y0, xc0 - x0:xc1 - x0] - vignets[i][inner] = image_data[yc0:yc1, xc0:xc1] - if seg is not None: - seg_stamp[:] = 0 - seg_stamp[inner] = seg[yc0:yc1, xc0:xc1] - vignets[i][(seg_stamp != 0) & (seg_stamp != number[i])] = BIG - - return vignets - - -def _build_ldac_objects(cat_data, vignets): - """Build LDAC_OBJECTS extension from an astropy table. - - Parameters - ---------- - cat_data : astropy.table.Table - Catalogue data read from the ASCII SExtractor file - vignets : numpy.ndarray - Array of shape ``(n_obj, stamp_size, stamp_size)`` - - Returns - ------- - astropy.io.fits.BinTableHDU - LDAC_OBJECTS extension - - """ - stamp_size = vignets.shape[1] - n_pix = stamp_size * stamp_size - - fits_cols = [] - for colname in cat_data.colnames: - arr = np.array(cat_data[colname]) - kind = arr.dtype.kind - if kind in ("i", "u"): - fmt = "J" - elif kind == "f": - fmt = "E" if arr.dtype.itemsize <= 4 else "D" - else: - fmt = f"{arr.dtype.itemsize}A" - fits_cols.append(fits.Column(name=colname, format=fmt, array=arr)) - - _aliases = { - "XWIN_IMAGE": "X_IMAGE", - "YWIN_IMAGE": "Y_IMAGE", - "XWIN_WORLD": "ALPHA_J2000", - "YWIN_WORLD": "DELTA_J2000", - } - col_map = {c.name: c for c in fits_cols} - for alias, source in _aliases.items(): - if source in col_map and alias not in col_map: - src_col = col_map[source] - fits_cols.append( - fits.Column( - name=alias, - format=src_col.format, - array=np.array(cat_data[source]), - ) - ) - - vignet_col = fits.Column( - name="VIGNET", - format=f"{n_pix}E", - array=vignets.reshape(len(vignets), n_pix), - dim=f"({stamp_size},{stamp_size})", - ) - fits_cols.append(vignet_col) - - hdu = fits.BinTableHDU.from_columns(fits_cols) - hdu.name = "LDAC_OBJECTS" - return hdu - - -def make_ldac_from_ascii( - input_cat_path, - image_path, - output_cat_path, - stamp_size=51, - seg_path=None, - w_log=None, -): - """Convert an external ASCII catalogue to FITS-LDAC format. - - Reads an ASCII catalogue in SExtractor format (column definitions as - ``# N NAME ...`` comment lines followed by data rows), attaches the - tile image header as an ``LDAC_IMHEAD`` extension so that - ``make_post_process`` can recover the contributing single exposures, and - writes a standard FITS-LDAC file compatible with all downstream ShapePipe - modules. - - The input columns, ``NUMBER`` included, are copied unchanged; a - ``VIGNET`` column (postage stamps extracted from the tile image) is - added to ``LDAC_OBJECTS``. Given the catalogue's segmentation map, which - shares the tile's pixel grid, the map is relabelled to the catalogue's - ``NUMBER`` (:func:`relabel_seg`) and used to set neighbours' pixels in each ``VIGNET`` to ``BIG``, as - SExtractor does. - - Parameters - ---------- - input_cat_path : str - Path to input ASCII catalogue in SExtractor format - image_path : str - Path to tile image FITS file - output_cat_path : str - Path to the output FITS-LDAC catalogue - stamp_size : int, optional - Side length of the square postage stamp in pixels, default 51 - seg_path : str, optional - Path to the catalogue's segmentation map (FITS, compressed or not) - w_log : logging.Logger, optional - Pipeline logger - - """ - cat_data = asc.read(input_cat_path, format="sextractor") - n_obj = len(cat_data) - if w_log: - w_log.info(f"Read {n_obj} objects from {input_cat_path}") - - with fits.open(image_path) as hdul: - img_header = hdul[0].header - image_data = hdul[0].data.astype(np.float32) - - seg = None - if seg_path is not None: - with fits.open(seg_path) as hdul: - hdu = next(h for h in hdul if h.data is not None) - seg_raw = hdu.data - if seg_raw.shape != image_data.shape: - raise ValueError( - f"Segmentation map {seg_path} has shape {seg_raw.shape}, the" - + f" image {image_data.shape}; they must share one grid." - ) - seg, counts = relabel_seg( - seg_raw, cat_data["NUMBER"], cat_data["X_IMAGE"], - cat_data["Y_IMAGE"], - ) - del seg_raw - if w_log: - w_log.info( - f"Relabelled {seg_path} to NUMBER: " - + ", ".join(f"{k}={v}" for k, v in counts.items()) - ) - - if w_log: - w_log.info( - f"Extracting {stamp_size}x{stamp_size} vignets from {image_path}" - ) - vignets = _extract_vignets( - image_data, - cat_data["X_IMAGE"], - cat_data["Y_IMAGE"], - stamp_size, - seg=seg, - number=np.asarray(cat_data["NUMBER"]), - ) - - ldac_imhead = _build_ldac_imhead(img_header) - ldac_objects = _build_ldac_objects(cat_data, vignets) - - hdul_out = fits.HDUList([fits.PrimaryHDU(), ldac_imhead, ldac_objects]) - hdul_out.writeto(output_cat_path, overwrite=True) - - if w_log: - w_log.info(f"Written FITS-LDAC catalogue to {output_cat_path}") diff --git a/src/shapepipe/modules/read_ext_sexcat_runner.py b/src/shapepipe/modules/read_ext_sexcat_runner.py deleted file mode 100644 index 17e23163d..000000000 --- a/src/shapepipe/modules/read_ext_sexcat_runner.py +++ /dev/null @@ -1,85 +0,0 @@ -"""READ EXTERNAL SEXCAT RUNNER. - -Module runner for ``read_ext_sexcat``. - -:Author: Martin Kilbinger - -""" - -from shapepipe.modules.module_decorator import module_runner -from shapepipe.modules.read_ext_sexcat_package import read_ext_sexcat as rs -from shapepipe.modules.sextractor_package import sextractor_script as ss - - -@module_runner( - version="1.0", - input_module=["get_images_runner"], - file_pattern=["CFIS_cat", "CFIS_image"], - file_ext=[".cat", ".fits"], - depends=["numpy", "astropy"], -) -def read_ext_sexcat_runner( - input_file_list, - run_dirs, - file_number_string, - config, - module_config_sec, - w_log, -): - """Define the Read External SExtractor Catalogue Runner. - - Reads an external ASCII catalogue (SExtractor format), converts it to - a FITS-LDAC catalogue compatible with downstream ShapePipe modules. - The inputs are the catalogue and the tile image, then the catalogue's - segmentation map if SEGMENTATION = True, then the WCS log if - MAKE_POST_PROCESS = True. With the segmentation map, neighbours' pixels - in VIGNET are set to -1e30; the map itself is not written out; the - only output is the FITS-LDAC ``sexcat.fits``. MAKE_POST_PROCESS runs the - multi-epoch post-processing that adds per-exposure HDUs. - """ - cat_path, image_path, *extra_inputs = input_file_list - use_seg = config.has_option( - module_config_sec, "SEGMENTATION" - ) and config.getboolean(module_config_sec, "SEGMENTATION") - seg_path = extra_inputs.pop(0) if use_seg else None - - if config.has_option(module_config_sec, "SUFFIX"): - suffix = config.get(module_config_sec, "SUFFIX") - else: - suffix = "sexcat" - - output_path = f"{run_dirs['output']}/{suffix}{file_number_string}.fits" - - if config.has_option(module_config_sec, "VIGNET_SIZE"): - stamp_size = config.getint(module_config_sec, "VIGNET_SIZE") - else: - stamp_size = 51 - - w_log.info(f"Reading external catalogue: {cat_path}") - w_log.info(f"Reading image header from: {image_path}") - - rs.make_ldac_from_ascii( - cat_path, - image_path, - output_path, - stamp_size=stamp_size, - seg_path=seg_path, - w_log=w_log, - ) - - if config.getboolean(module_config_sec, "MAKE_POST_PROCESS"): - # The WCS log is supplied as a positional input (via FILE_PATTERN) - # when post-processing is enabled, not by the decorator default. - if not extra_inputs: - raise ValueError( - "MAKE_POST_PROCESS requires the WCS log file as the last" - + " input; add 'log_exp_headers' to FILE_PATTERN and" - + f" FILE_EXT in the [{module_config_sec}] config section." - ) - f_wcs_path = extra_inputs[0] - pos_params = config.getlist(module_config_sec, "WORLD_POSITION") - ccd_size = config.getlist(module_config_sec, "CCD_SIZE") - w_log.info("Running post-processing") - ss.make_post_process(output_path, f_wcs_path, pos_params, ccd_size) - - return None, None diff --git a/src/shapepipe/modules/sextractor_package/match_catalogue.py b/src/shapepipe/modules/sextractor_package/match_catalogue.py new file mode 100644 index 000000000..68bca78f7 --- /dev/null +++ b/src/shapepipe/modules/sextractor_package/match_catalogue.py @@ -0,0 +1,225 @@ +"""MATCH CATALOGUE. + +Join a SExtractor tile catalogue to an external catalogue of the same image, +so the tile's rows carry the external catalogue's object numbers. + +The UNIONS tile catalogue (Stephen Gwyn's MegaPipe SExtractor run on the DR6 +tile) defines the object list and the ``NUMBER`` that the photometry and +photo-z catalogues share. ShapePipe runs its own SExtractor on the same image +with the same detection configuration, which reproduces that catalogue for +all but a few objects in a thousand, and takes membership and ``NUMBER`` from +it here: every +measurement (windowed positions, ``VIGNET`` and its neighbour marking, the +SExtractor columns) stays SExtractor's own. + +:Author: Cail Daley + +""" + +import numpy as np +from astropy.io import ascii as asc +from astropy.io import fits +from scipy.spatial import cKDTree + +# The SEG_VIGNET label of a footprint whose SExtractor row has no partner in +# the external catalogue, and so leaves the catalogue. Negative, so it never +# collides with a NUMBER; UberSeg only asks "self or not self". +UNMATCHED_LABEL = -1 + + +def mutual_nearest(x_a, y_a, x_b, y_b, radius): + """One-to-one pairs of mutual nearest neighbours closer than ``radius``. + + Parameters + ---------- + x_a, y_a : array_like + Positions of the first set + x_b, y_b : array_like + Positions of the second set, in the same units + radius : float + Largest separation of a pair + + Returns + ------- + numpy.ndarray + Indices into the first set of the paired rows, ascending + numpy.ndarray + Index into the second set of each one's partner + numpy.ndarray + Separation of each pair + + """ + a = np.column_stack([x_a, y_a]).astype(np.float64) + b = np.column_stack([x_b, y_b]).astype(np.float64) + if len(a) == 0 or len(b) == 0: + empty = np.zeros(0, np.int64) + return empty, empty, np.zeros(0) + dist, nearest_b = cKDTree(b).query(a) + _, nearest_a = cKDTree(a).query(b) + idx_a = np.arange(len(a)) + paired = (dist < radius) & (nearest_a[nearest_b] == idx_a) + return idx_a[paired], nearest_b[paired], dist[paired] + + +def relabel(seg_vignets, old_number, new_number): + """Map segmentation labels from one numbering to another. + + Parameters + ---------- + seg_vignets : numpy.ndarray + Segmentation stamps labelled with the old numbering, 0 on sky + old_number : array_like of int + Every old label in use (the catalogue's ``NUMBER`` before the join) + new_number : array_like of int + New label of each ``old_number``; ``UNMATCHED_LABEL`` for one that + leaves the catalogue + + Returns + ------- + numpy.ndarray + The stamps in the new numbering, int32; labels outside + ``old_number`` become ``UNMATCHED_LABEL`` + + """ + old_number = np.asarray(old_number, np.int64) + table = np.full(int(max(old_number.max(initial=0), + seg_vignets.max(initial=0))) + 1, + UNMATCHED_LABEL, np.int32) + table[0] = 0 + table[old_number] = new_number + if seg_vignets.min(initial=0) < 0: + raise ValueError("Segmentation stamps hold negative labels.") + return table[seg_vignets] + + +def match_catalogue(cat_path, ext_cat_path, radius=1.0, min_fraction=0.98, + tolerated_unpaired=20, w_log=None): + """Take membership and ``NUMBER`` from an external catalogue. + + Each SExtractor row is paired with its mutual nearest neighbour in the + external catalogue (``X_IMAGE``, ``Y_IMAGE``, same image grid) within + ``radius`` pixels. Paired rows take the external ``NUMBER``; unpaired + rows leave the catalogue. A ``SEG_VIGNET`` column, when present, is + relabelled to the new numbering, with the footprints of rows that left + marked ``UNMATCHED_LABEL``. The catalogue is rewritten in place; every + other HDU and column is kept. + + Both catalogues come from the same pixels, so pairs agree to ~1e-4 + pixel, and on eight DR6 tiles at least 99.2% of each side pairs; the + shortfall is deblending that differs near very large objects (a cD + galaxy, a bright star's halo). On different pixels (a DR5 image against + the DR6 catalogue) only ~95% pair. The guard is there to catch that: the + run stops when the unpaired rows of either side exceed both + ``(1 - min_fraction)`` of that side and ``tolerated_unpaired``, so a + near-empty edge tile with a couple of unpaired children passes. On the + SExtractor side it catches a catalogue from other pixels; on the + external side, a detection configuration that finds fewer objects than + the catalogue's, which would leave its objects without shapes. + + Parameters + ---------- + cat_path : str + Path to the SExtractor FITS-LDAC catalogue, rewritten in place + ext_cat_path : str + Path to the external ASCII catalogue in SExtractor format, with + ``NUMBER``, ``X_IMAGE`` and ``Y_IMAGE`` + radius : float, optional + Largest separation of a pair, in pixels; default 1 + min_fraction : float, optional + Smallest acceptable fraction of SExtractor rows, and of external + objects, with a partner; default 0.98 + tolerated_unpaired : int, optional + Number of unpaired rows on either side accepted whatever the + fraction; default 20 + w_log : logging.Logger, optional + Pipeline logger + + Returns + ------- + dict + ``n_sextractor``, ``n_external``, ``n_paired``, ``n_dropped`` + (SExtractor rows without a partner) and ``n_external_only`` + (external objects without a SExtractor row) + + Raises + ------ + ValueError + If more than ``tolerated_unpaired`` rows of either side, and more + than ``1 - min_fraction`` of it, have no partner + + @sc [decision:detection.tile_detection] + """ + ext = asc.read(ext_cat_path, format="sextractor", + include_names=["NUMBER", "X_IMAGE", "Y_IMAGE"]) + with fits.open(cat_path) as hdul: + hdus = [hdu.copy() for hdu in hdul] + objects = next(h for h in hdus if h.name == "LDAC_OBJECTS") + data = objects.data + + i_sex, i_ext, dist = mutual_nearest( + data["X_IMAGE"], data["Y_IMAGE"], ext["X_IMAGE"], ext["Y_IMAGE"], + radius, + ) + counts = dict( + n_sextractor=len(data), n_external=len(ext), n_paired=len(i_sex), + n_dropped=len(data) - len(i_sex), + n_external_only=len(ext) - len(i_sex), + ) + summary = ", ".join(f"{k}={v}" for k, v in counts.items()) + def too_many_unpaired(n_total): + unpaired = n_total - len(i_sex) + return (unpaired > tolerated_unpaired + and unpaired > (1 - min_fraction) * n_total) + + fraction = len(i_sex) / max(len(data), 1) + if too_many_unpaired(len(data)): + raise ValueError( + f"Only {fraction:.4f} of the SExtractor rows of {cat_path} pair" + + f" with {ext_cat_path} within {radius} px (minimum" + + f" {min_fraction} beyond {tolerated_unpaired} unpaired;" + + f" {summary}). The catalogue must come from the" + + " same pixels as the image: check that the tile image is the" + + " release of the catalogue (DR6)." + ) + fraction_ext = len(i_sex) / max(len(ext), 1) + if too_many_unpaired(len(ext)): + raise ValueError( + f"Only {fraction_ext:.4f} of the objects of {ext_cat_path} pair" + + f" with a row of {cat_path} within {radius} px (minimum" + + f" {min_fraction} beyond {tolerated_unpaired} unpaired;" + + f" {summary}). SExtractor finds fewer objects" + + " than the catalogue: check that the detection configuration" + + " is the catalogue's." + ) + + old_number = np.asarray(data["NUMBER"]) + new_number = np.full(len(data), UNMATCHED_LABEL, np.int64) + new_number[i_sex] = np.asarray(ext["NUMBER"])[i_ext] + + columns = [] + for col in objects.columns: + array = data[col.name][i_sex] + if col.name == "NUMBER": + array = new_number[i_sex].astype(array.dtype) + elif col.name == "SEG_VIGNET": + array = relabel(np.asarray(data[col.name]), old_number, + new_number)[i_sex] + columns.append(fits.Column( + name=col.name, format=col.format, unit=col.unit, dim=col.dim, + array=array, + )) + new = fits.BinTableHDU.from_columns( + columns, header=objects.header, name="LDAC_OBJECTS", + ) + fits.HDUList( + [new if h is objects else h for h in hdus] + ).writeto(cat_path, overwrite=True) + + if w_log: + sep = np.quantile(dist, [0.5, 0.99]) if len(dist) else [np.nan] * 2 + w_log.info( + f"Matched {cat_path} to {ext_cat_path} within {radius} px:" + + f" {summary}; pair separation p50 {sep[0]:.2e} px," + + f" p99 {sep[1]:.2e} px" + ) + return counts diff --git a/src/shapepipe/modules/sextractor_package/sextractor_script.py b/src/shapepipe/modules/sextractor_package/sextractor_script.py index 5b97fb07d..c3716db74 100644 --- a/src/shapepipe/modules/sextractor_package/sextractor_script.py +++ b/src/shapepipe/modules/sextractor_package/sextractor_script.py @@ -384,6 +384,13 @@ def set_input_files( Set up all of the input image files. + Without a detection image, SExtractor runs in single-image mode on + the measurement image: ``img -WEIGHT_IMAGE w``. With one, it runs in + dual-image mode, ``det,img -WEIGHT_IMAGE det_w,w``. Dual-image mode + with the same image twice is not equivalent to single-image mode: it + changes the background and the photometry of a fifth of the objects of + a UNIONS tile (#936). + @sc [decision:detection.weight_map_usage,label:convention] weight-map-on-both-images With a weight file, the same map weights detection and measurement (a separate detection weight only when DETECTION_WEIGHT is set); without @@ -443,26 +450,28 @@ def set_input_files( "DETECTION_WEIGHT cannot be True " + "if WEIGHT_FILE is False" ) - # Check for separate image file for detection and measurement + # A separate detection image switches SExtractor to dual-image mode; + # without one, the measurement image alone (single-image mode). if use_detect_img: - self._detect_img_path = self._all_input_path[extra] + self._image_arg = ( + f"{self._all_input_path[extra]},{self._meas_img_path}" + ) extra += 1 else: - self._detect_img_path = self._meas_img_path + self._image_arg = self._meas_img_path - # Check for separate weight file corresponding to the detection image. - # If False, use measurement weight image. - # Note: This could be changed, and no weight image could be used, but - # this might lead to more user errors. + # In dual-image mode the detection image takes the detection weight + # if given, else the measurement weight. if use_weight: - if use_detect_weight: - detect_weight_path = self._all_input_path[extra] - extra += 1 - else: - detect_weight_path = weight_image - self._cmd_line_extra += ( - f" -WEIGHT_IMAGE {detect_weight_path}" + f",{weight_image}" - ) + weight_arg = weight_image + if use_detect_img: + if use_detect_weight: + detect_weight_path = self._all_input_path[extra] + extra += 1 + else: + detect_weight_path = weight_image + weight_arg = f"{detect_weight_path},{weight_image}" + self._cmd_line_extra += f" -WEIGHT_IMAGE {weight_arg}" else: self._cmd_line_extra += " -WEIGHT_TYPE None" @@ -563,7 +572,7 @@ def make_command_line(self, exec_path): """ # Base arguments for SExtractor command_line_base = ( - f"{exec_path} {self._detect_img_path},{self._meas_img_path} " + f"{exec_path} {self._image_arg} " + f"-c {self._path_dot_sex} " + f"-PARAMETERS_NAME {self._path_dot_param} " + f"-FILTER_NAME {self._path_dot_conv} " diff --git a/src/shapepipe/modules/sextractor_runner.py b/src/shapepipe/modules/sextractor_runner.py index e45885b06..65181712f 100644 --- a/src/shapepipe/modules/sextractor_runner.py +++ b/src/shapepipe/modules/sextractor_runner.py @@ -9,6 +9,7 @@ import re from shapepipe.modules.module_decorator import module_runner +from shapepipe.modules.sextractor_package import match_catalogue as mc from shapepipe.modules.sextractor_package import sextractor_script as ss from shapepipe.pipeline.execute import execute @@ -109,6 +110,28 @@ def sextractor_runner( # Parse SExtractor errors stdout, stderr = ss_inst.parse_errors(stderr, stdout) + # MATCH_CATALOGUE (optional, environment-expanded path; empty for none): + # take membership and NUMBER from that external catalogue of the same + # image, before the post-processing keys the epoch HDUs on NUMBER. + match_path = ( + config.getexpanded(module_config_sec, "MATCH_CATALOGUE") + if config.has_option(module_config_sec, "MATCH_CATALOGUE") + else "" + ) + if match_path: + mc.match_catalogue( + ss_inst.path_output_file, + match_path, + radius=config.getfloat(module_config_sec, "MATCH_RADIUS"), + min_fraction=config.getfloat( + module_config_sec, "MATCH_MIN_FRACTION" + ), + tolerated_unpaired=config.getint( + module_config_sec, "MATCH_TOLERATED_UNPAIRED" + ), + w_log=w_log, + ) + # Run sextractor post processing if config.getboolean(module_config_sec, "MAKE_POST_PROCESS"): pos_params = config.getlist(module_config_sec, "WORLD_POSITION") diff --git a/src/shapepipe/pipeline/config.py b/src/shapepipe/pipeline/config.py index 4e809316b..2c294ae65 100644 --- a/src/shapepipe/pipeline/config.py +++ b/src/shapepipe/pipeline/config.py @@ -11,18 +11,24 @@ from configparser import ConfigParser +# ``$VAR``, ``${VAR}`` or ``${VAR:-default}``. +_VAR_RE = re.compile(r"\$(?:\{(\w+)(?::-([^}]*))?\}|(\w+))") + + def _expandvars_strict(value): """Expand Environment Variables Strictly. Expand environment variables in ``value``, raising an error if any referenced variable is unset instead of silently leaving the literal - ``$VAR`` string in place. + ``$VAR`` string in place. ``${VAR:-default}`` expands to ``default`` + when ``VAR`` is unset or empty, as in the shell. Each value is inserted + as is, not expanded again. Parameters ---------- value : str - Configuration value possibly containing ``$VAR`` or ``${VAR}`` - references + Configuration value possibly containing ``$VAR``, ``${VAR}`` or + ``${VAR:-default}`` references Returns ------- @@ -32,11 +38,30 @@ def _expandvars_strict(value): Raises ------ ValueError - If a referenced environment variable is not set + If a referenced environment variable without a default is not set, + or a ``${`` opens none of the forms above (``${VAR-x}``, + ``${VAR:=x}``, ``${ VAR }``) """ - expanded = os.path.expandvars(value) - unset = re.findall(r"\$\{?(\w+)\}?", expanded) + if "${" in _VAR_RE.sub("", value): + raise ValueError( + f"Config value '{value}' has a '${{' that is not $VAR, ${{VAR}}" + + " or ${VAR:-default}." + ) + unset = [] + + def substitute(match): + braced, default, bare = match.groups() + name = braced or bare + env = os.environ.get(name) + if default is not None: + return env or default + if env is None: + unset.append(name) + return match.group(0) + return env + + expanded = _VAR_RE.sub(substitute, value) if unset: raise ValueError( f"Environment variable(s) {', '.join(sorted(set(unset)))} " diff --git a/tests/module/test_final_cat_columns.py b/tests/module/test_final_cat_columns.py deleted file mode 100644 index 51ea46772..000000000 --- a/tests/module/test_final_cat_columns.py +++ /dev/null @@ -1,141 +0,0 @@ -"""The final-catalogue columns each tile_detection mode requests. - -``final_cat_merge`` asks every tile catalogue for the columns of -``workflow/config/cfis/final_cat.param`` and stops on any one missing. That -list is written against SExtractor-mode tiles; under ``tile_detection: -unions_catalogue`` the detection columns are whatever the UNIONS per-tile -catalogue carries, copied through ``read_ext_sexcat``. So the catalogue-mode -request is the param list less ``merge_final_cat.SEXTRACTOR_ONLY_COLUMNS``. - -Here the catalogue-mode detection columns are not written down but derived: -the converter runs on a catalogue with the UNIONS DR6 header (copied verbatim -from ``vos:cfis/tiles_DR6/CFIS.202.301.r.cat``) and its output is what a tile -can carry. The requested SExtractor columns are those of the tile SExtractor -parameter file. Every other requested column comes from stages that run the -same way in both modes (post-processing, ngmix, make_cat), so the detection -columns are the whole of the difference between them. -""" - -import importlib.util -import re -import sys -from pathlib import Path - -import numpy as np -import pytest -from astropy.io import fits - -from shapepipe.modules.read_ext_sexcat_package import read_ext_sexcat as rs - -REPO_ROOT = Path(__file__).resolve().parents[2] -SCRIPTS = REPO_ROOT / "workflow" / "scripts" -CFIS = REPO_ROOT / "workflow" / "config" / "cfis" - -# The header of a UNIONS DR6 per-tile catalogue and its first object, moved to -# pixel (4, 4) so that its stamp falls on the small test image. -DR6_CATALOGUE = """\ -# 1 NUMBER Running object number -# 2 X_IMAGE Object position along x [pixel] -# 3 Y_IMAGE Object position along y [pixel] -# 4 ALPHA_J2000 Right ascension of barycenter (J2000) [deg] -# 5 DELTA_J2000 Declination of barycenter (J2000) [deg] -# 6 MAG_AUTO Kron-like elliptical aperture magnitude [mag] -# 7 MAGERR_AUTO RMS error for AUTO magnitude [mag] -# 8 MAG_BEST Best of MAG_AUTO and MAG_ISOCOR [mag] -# 9 MAGERR_BEST RMS error for MAG_BEST [mag] -# 10 MAG_APER Fixed aperture magnitude vector [mag] -# 11 MAGERR_APER RMS error vector for fixed aperture mag. [mag] -# 12 A_WORLD Profile RMS along major axis (world units) [deg] -# 13 ERRA_WORLD World RMS position error along major axis [deg] -# 14 B_WORLD Profile RMS along minor axis (world units) [deg] -# 15 ERRB_WORLD World RMS position error along minor axis [deg] -# 16 THETA_J2000 Position angle (east of north) (J2000) [deg] -# 17 ERRTHETA_J2000 J2000 error ellipse pos. angle (east of north) [deg] -# 18 ISOAREA_IMAGE Isophotal area above Analysis threshold [pixel**2] -# 19 MU_MAX Peak surface brightness above background [mag * arcsec**(-2)] -# 20 FLUX_RADIUS Fraction-of-light radii [pixel] -# 21 FLAGS Extraction flags - 1 4.0000 4.0000 205.3679556 +60.2576198 14.8792 0.0002 14.8792 0.0002 14.9100 0.0002 0.000389709 2.80088e-07 0.0003095117 1.974931e-07 -1.48 -0.72 5707 16.3802 4.474 0 -""" - - -def _load_script(name, path): - sys.path.insert(0, str(SCRIPTS)) - try: - spec = importlib.util.spec_from_file_location(f"_{name}", path) - module = importlib.util.module_from_spec(spec) - spec.loader.exec_module(module) - finally: - sys.path.remove(str(SCRIPTS)) - return module - - -merge_final_cat = _load_script("merge_final_cat", - SCRIPTS / "merge_final_cat.py") -create_final_cat = _load_script( - "create_final_cat", - REPO_ROOT / "scripts" / "python" / "create_final_cat.py") - - -@pytest.fixture(scope="module") -def param_list(): - return create_final_cat.read_param_file(str(CFIS / "final_cat.param")) - - -@pytest.fixture(scope="module") -def sextractor_columns(): - """Column names the tile SExtractor run writes (vector sizes stripped).""" - names = set() - for line in (CFIS / "default_noimaflags.param").read_text().splitlines(): - entry = line.split("#")[0].strip() - if entry: - names.add(re.sub(r"\(.*\)$", "", entry)) - return names - - -@pytest.fixture(scope="module") -def catalogue_columns(tmp_path_factory): - """Detection columns of a tile under tile_detection: unions_catalogue.""" - tmp = tmp_path_factory.mktemp("dr6") - cat, img, out = (tmp / "CFIS.202.301.r.cat", tmp / "CFIS.202.301.r.fits", - tmp / "sexcat-202-301.fits") - cat.write_text(DR6_CATALOGUE) - fits.PrimaryHDU(np.zeros((8, 8), dtype=np.float32)).writeto(img) - rs.make_ldac_from_ascii(str(cat), str(img), str(out), stamp_size=3) - with fits.open(out) as hdul: - return set(hdul["LDAC_OBJECTS"].columns.names) - {"VIGNET"} - - -def test_sextractor_only_is_exactly_what_the_catalogue_lacks( - param_list, sextractor_columns, catalogue_columns): - """Every requested SExtractor column the catalogue lacks, and no other.""" - lacking = {c for c in param_list - if c in sextractor_columns and c not in catalogue_columns} - assert set(merge_final_cat.SEXTRACTOR_ONLY_COLUMNS) == lacking - - -def test_catalogue_mode_request_is_carried( - param_list, sextractor_columns, catalogue_columns): - """Each requested detection column is one the converter writes.""" - requested = merge_final_cat.requested_columns(param_list, - "unions_catalogue") - detection = [c for c in requested if c in sextractor_columns] - assert detection and set(detection) <= catalogue_columns - # Only detection columns are dropped, and the order is the param file's. - assert requested == [c for c in param_list - if c not in merge_final_cat.SEXTRACTOR_ONLY_COLUMNS] - assert {"MAG_AUTO", "MAGERR_AUTO", "FLUX_RADIUS"} <= set(requested) - - -def test_sextractor_mode_requests_the_param_file(): - for input_type in ("cfis", "cfis_image_sims"): - params = create_final_cat.read_param_file( - str(REPO_ROOT / "workflow" / "config" / input_type - / "final_cat.param")) - assert merge_final_cat.requested_columns( - params, "sextractor") == params - - -def test_unknown_mode_is_refused(param_list): - with pytest.raises(ValueError): - merge_final_cat.requested_columns(param_list, "sextractr") diff --git a/tests/module/test_match_catalogue.py b/tests/module/test_match_catalogue.py new file mode 100644 index 000000000..6b7d47c46 --- /dev/null +++ b/tests/module/test_match_catalogue.py @@ -0,0 +1,222 @@ +"""Joining the tile SExtractor catalogue to the UNIONS catalogue.""" + +import numpy as np +import pytest +from astropy.io import fits + +from shapepipe.modules.sextractor_package import match_catalogue as mc +from shapepipe.pipeline.config import CustomParser + +HEADER = """\ +# 1 NUMBER Running object number +# 2 X_IMAGE Object position along x [pixel] +# 3 Y_IMAGE Object position along y [pixel] +# 4 MAG_AUTO Kron-like elliptical aperture magnitude [mag] +""" + + +def _external(path, number, x, y): + rows = "".join(f"{n:10d} {xi:11.4f} {yi:11.4f} 24.0000\n" + for n, xi, yi in zip(number, x, y)) + path.write_text(HEADER + rows) + return str(path) + + +def _sexcat(path, x, y, seg_vignets=None): + """A FITS-LDAC sexcat: NUMBER 1..n at (x, y), 3x3 stamps.""" + n = len(x) + number = np.arange(1, n + 1, dtype=np.int32) + cols = [ + fits.Column(name="NUMBER", format="J", array=number), + fits.Column(name="X_IMAGE", format="E", array=np.asarray(x)), + fits.Column(name="Y_IMAGE", format="E", array=np.asarray(y)), + fits.Column(name="XWIN_IMAGE", format="D", + array=np.asarray(x) + 0.01), + fits.Column(name="VIGNET", format="9E", dim="(3,3)", + array=np.arange(n * 9, dtype=np.float32).reshape(n, 9)), + ] + if seg_vignets is not None: + cols.append(fits.Column(name="SEG_VIGNET", format="9J", dim="(3,3)", + array=seg_vignets.reshape(n, 9))) + imhead = fits.BinTableHDU.from_columns( + [fits.Column(name="Field Header Card", format="80A", + array=np.array(["SIMPLE = T"]))], + name="LDAC_IMHEAD") + objects = fits.BinTableHDU.from_columns(cols, name="LDAC_OBJECTS") + fits.HDUList([fits.PrimaryHDU(), imhead, objects]).writeto(path) + return str(path) + + +def _objects(path): + with fits.open(path) as hdul: + assert [h.name for h in hdul] == ["PRIMARY", "LDAC_IMHEAD", + "LDAC_OBJECTS"] + return hdul["LDAC_OBJECTS"].data.copy() + + +X = np.array([10.0, 50.0, 90.0, 130.0, 170.0]) +Y = np.array([20.0, 60.0, 100.0, 140.0, 180.0]) + + +def test_rows_take_the_external_number_one_to_one(tmp_path): + """Every row pairs with its object, whatever the external order.""" + cat = _sexcat(tmp_path / "sexcat.fits", X, Y) + order = [3, 0, 4, 1, 2] + numbers = [1001, 1002, 1003, 1004, 1005] + # Shuffled, offset by 1e-4 px, plus one object SExtractor did not detect. + ext = _external(tmp_path / "ext.cat", numbers + [1006], + list(X[order] + 1e-4) + [300.0], + list(Y[order] - 1e-4) + [300.0]) + counts = mc.match_catalogue(cat, ext, min_fraction=0.8) + data = _objects(cat) + expected = np.empty(5, int) + expected[order] = numbers + assert list(data["NUMBER"]) == list(expected) + assert counts == dict(n_sextractor=5, n_external=6, n_paired=5, + n_dropped=0, n_external_only=1) + # Every measurement is SExtractor's, row for row. + np.testing.assert_array_equal(data["XWIN_IMAGE"], X + 0.01) + np.testing.assert_array_equal(data["VIGNET"][2].ravel(), + np.arange(18, 27)) + + +def test_unpaired_rows_leave_the_catalogue(tmp_path): + """A detection with no object within the radius is dropped.""" + cat = _sexcat(tmp_path / "sexcat.fits", X, Y) + ext = _external(tmp_path / "ext.cat", [7, 8, 9, 10], + X[[0, 1, 2, 4]], Y[[0, 1, 2, 4]] + 0.9) + counts = mc.match_catalogue(cat, ext, radius=1.0, min_fraction=0.5) + data = _objects(cat) + assert list(data["NUMBER"]) == [7, 8, 9, 10] + np.testing.assert_allclose(data["X_IMAGE"], X[[0, 1, 2, 4]]) + assert counts["n_dropped"] == 1 and counts["n_external_only"] == 0 + + +def test_pairs_are_mutual_nearest_neighbours(tmp_path): + """Two detections near one object: only the nearer one takes it.""" + x = np.array([10.0, 10.6, 90.0]) + y = np.array([10.0, 10.0, 90.0]) + cat = _sexcat(tmp_path / "sexcat.fits", x, y) + ext = _external(tmp_path / "ext.cat", [5, 6], [10.1, 90.0], [10.0, 90.0]) + mc.match_catalogue(cat, ext, min_fraction=0.5) + data = _objects(cat) + assert list(data["NUMBER"]) == [5, 6] + np.testing.assert_allclose(data["X_IMAGE"], [10.0, 90.0]) + + +def test_pairs_lie_within_the_radius(tmp_path): + """Mutual nearest neighbours 2 px apart, beyond the radius, do not pair.""" + x = np.append(X, 300.0) + y = np.append(Y, 300.0) + cat = _sexcat(tmp_path / "sexcat.fits", x, y) + ext = _external(tmp_path / "ext.cat", [1, 2, 3, 4, 5, 6], + np.append(X, 302.0), np.append(Y, 300.0)) + counts = mc.match_catalogue(cat, ext, radius=1.0, min_fraction=0.8) + assert list(_objects(cat)["NUMBER"]) == [1, 2, 3, 4, 5] + assert counts["n_dropped"] == 1 and counts["n_external_only"] == 1 + + +def test_too_few_catalogue_objects_paired_stop_the_run(tmp_path): + """Fewer detections than catalogue objects (a drifted detection + configuration) fails, though every detection pairs.""" + cat = _sexcat(tmp_path / "sexcat.fits", X[:4], Y[:4]) + ext = _external(tmp_path / "ext.cat", [1, 2, 3, 4, 5], X, Y) + with pytest.raises(ValueError, match="fewer objects"): + mc.match_catalogue(cat, ext, min_fraction=0.99, tolerated_unpaired=0) + assert list(_objects(cat)["NUMBER"]) == [1, 2, 3, 4] + counts = mc.match_catalogue(cat, ext, min_fraction=0.8) + assert counts["n_paired"] == 4 and counts["n_external_only"] == 1 + + +def test_too_few_pairs_stop_the_run(tmp_path): + """Below MATCH_MIN_FRACTION (other pixels: the DR5 image) the run fails.""" + cat = _sexcat(tmp_path / "sexcat.fits", X, Y) + ext = _external(tmp_path / "ext.cat", [1, 2, 3, 4], X[:4], Y[:4]) + with pytest.raises(ValueError, match="DR6"): + mc.match_catalogue(cat, ext, min_fraction=0.99, tolerated_unpaired=0) + # Nothing was written. + assert list(_objects(cat)["NUMBER"]) == [1, 2, 3, 4, 5] + mc.match_catalogue(cat, ext, min_fraction=0.8) + assert list(_objects(cat)["NUMBER"]) == [1, 2, 3, 4] + + +def _grid(n): + """n well-separated positions, 10 px apart.""" + i = np.arange(n) + return 5.0 + 10.0 * (i % 50), 5.0 + 10.0 * (i // 50) + + +def test_small_tile_with_few_unpaired_children_passes(tmp_path): + """An edge tile: 5 of 100 catalogue objects unpaired (5%) is within the + 20 tolerated, so the run goes on with the default guard.""" + x, y = _grid(105) + cat = _sexcat(tmp_path / "sexcat.fits", x[:100], y[:100]) + ext = _external(tmp_path / "ext.cat", np.arange(1, 106), x, y) + counts = mc.match_catalogue(cat, ext) + assert counts["n_paired"] == 100 and counts["n_external_only"] == 5 + + +def test_large_tile_with_unpaired_children_passes(tmp_path): + """A cluster tile: 30 of 2000 catalogue objects unpaired (1.5%) is more + than 20 but under 2%, so the run goes on with the default guard.""" + x, y = _grid(2000) + cat = _sexcat(tmp_path / "sexcat.fits", x[:1970], y[:1970]) + ext = _external(tmp_path / "ext.cat", np.arange(1, 2001), x, y) + counts = mc.match_catalogue(cat, ext) + assert counts["n_external_only"] == 30 + + +@pytest.mark.parametrize("side", ["sextractor", "external"]) +def test_other_pixels_stop_the_run(tmp_path, side): + """5% unpaired on either side (a DR5 image: ~95% pair) stops the run with + the default guard, far beyond the 20 tolerated.""" + x, y = _grid(2000) + n_sex, n_ext = (2000, 1900) if side == "sextractor" else (1900, 2000) + cat = _sexcat(tmp_path / "sexcat.fits", x[:n_sex], y[:n_sex]) + ext = _external(tmp_path / "ext.cat", np.arange(1, n_ext + 1), + x[:n_ext], y[:n_ext]) + with pytest.raises(ValueError, match="DR6" if side == "sextractor" + else "fewer objects"): + mc.match_catalogue(cat, ext) + + +def test_seg_vignet_is_relabelled(tmp_path): + """SEG_VIGNET labels follow NUMBER; a dropped row's footprint is -1.""" + x, y = X[:3], Y[:3] + seg = np.array([ + [[0, 1, 1], [0, 1, 2], [0, 0, 3]], + [[2, 2, 0], [1, 2, 0], [0, 0, 0]], + [[3, 3, 3], [0, 3, 3], [2, 0, 0]], + ], dtype=np.int32) + cat = _sexcat(tmp_path / "sexcat.fits", x, y, seg) + # Row 2 (NUMBER 2) has no partner; rows 1 and 3 become 40 and 70. + ext = _external(tmp_path / "ext.cat", [70, 40], [x[2], x[0]], + [y[2], y[0]]) + mc.match_catalogue(cat, ext, min_fraction=0.5) + data = _objects(cat) + assert list(data["NUMBER"]) == [40, 70] + lut = {0: 0, 1: 40, 2: mc.UNMATCHED_LABEL, 3: 70} + expected = np.vectorize(lut.get)(seg[[0, 2]]) + np.testing.assert_array_equal( + np.asarray(data["SEG_VIGNET"]).reshape(2, 3, 3), expected) + # Each object's own centre label is its new NUMBER, as UberSeg needs. + centre = np.asarray(data["SEG_VIGNET"]).reshape(2, 3, 3)[:, 1, 1] + assert list(centre) == [40, 70] + + +def test_relabel_rejects_unknown_negative_labels(): + with pytest.raises(ValueError): + mc.relabel(np.array([[-1, 1]]), [1], [5]) + + +def test_match_catalogue_key_expands_to_empty_when_unset(monkeypatch): + """config_tile_Sx.ini's ${SP_MATCH_CATALOGUE:-} means "no join" unset.""" + parser = CustomParser() + parser.read_string("[S]\nMATCH = ${SP_MATCH_CATALOGUE:-}\n" + "BARE = $SP_MATCH_CATALOGUE\n") + monkeypatch.delenv("SP_MATCH_CATALOGUE", raising=False) + assert parser.getexpanded("S", "MATCH") == "" + with pytest.raises(ValueError): + parser.getexpanded("S", "BARE") + monkeypatch.setenv("SP_MATCH_CATALOGUE", "/a/CFIS_cat-186-307.cat") + assert parser.getexpanded("S", "MATCH") == "/a/CFIS_cat-186-307.cat" diff --git a/tests/module/test_read_ext_sexcat.py b/tests/module/test_read_ext_sexcat.py deleted file mode 100644 index fa406b9a6..000000000 --- a/tests/module/test_read_ext_sexcat.py +++ /dev/null @@ -1,424 +0,0 @@ -"""UNIT TESTS FOR MODULE PACKAGE: READ_EXT_SEXCAT. - -Drives ``make_ldac_from_ascii`` on a synthetic ASCII SExtractor-format -catalogue and a synthetic tile image, and checks the FITS-LDAC it writes is -what the tile chain downstream of ``tile_detect`` reads: the LDAC_IMHEAD -extension carrying the tile header, the SExtractor column aliases, one -``VIGNET`` stamp per object cut from the image, and the input ``NUMBER`` -kept as is. It follows the catalogue through -``make_cat.save_sextractor_data``, which builds ``TILE_UNIQUE_ID``. The rest -covers the segmentation map: relabelling it to the catalogue's ``NUMBER`` and -setting neighbours' ``VIGNET`` pixels to -1e30, as SExtractor does, which is -all ngmix reads to mask neighbours. -""" - -from pathlib import Path - -import numpy as np -import numpy.testing as npt -import pytest -from astropy.io import fits - -from shapepipe.modules.make_cat_package import make_cat -from shapepipe.modules.read_ext_sexcat_package import read_ext_sexcat as rs - -NX, NY = 40, 30 -STAMP = 5 -# (NUMBER, X_IMAGE, Y_IMAGE): an interior object, one on the left edge, one -# in the top-right corner. -OBJECTS = [(1, 10.0, 12.0), (2, 1.0, 20.0), (7, 40.0, 30.0)] - - -def _write_ascii_cat(path): - lines = [ - "# 1 NUMBER Running object number", - "# 2 X_IMAGE Object position along x [pixel]", - "# 3 Y_IMAGE Object position along y [pixel]", - "# 4 ALPHA_J2000 Right ascension of barycenter [deg]", - "# 5 DELTA_J2000 Declination of barycenter [deg]", - "# 6 MAG_AUTO Kron-like elliptical aperture magnitude [mag]", - ] - for num, x, y in OBJECTS: - lines.append(f"{num} {x} {y} {150.0 + num} {30.0 + num} {20.0 + num}") - path.write_text("\n".join(lines) + "\n") - - -def _write_image(path): - # pixel value = 1000*row + column (0-based), so every stamp pixel names - # where it came from. - data = (np.arange(NY)[:, None] * 1000 + np.arange(NX)[None, :]).astype( - np.float32 - ) - hdu = fits.PrimaryHDU(data) - hdu.header["HISTORY"] = "input image 2605805p.fits" - hdu.header["TILEKEY"] = "kept" - hdu.writeto(path, overwrite=True) - - -@pytest.fixture -def ldac(tmp_path): - cat = tmp_path / "CFIS_cat-301-279.cat" - img = tmp_path / "CFIS_image-301-279.fits" - out = tmp_path / "sexcat-301-279.fits" - _write_ascii_cat(cat) - _write_image(img) - rs.make_ldac_from_ascii( - str(cat), str(img), str(out), stamp_size=STAMP - ) - return out - - -def test_ldac_layout_and_header(ldac): - with fits.open(ldac) as hdul: - assert [h.name for h in hdul] == ["PRIMARY", "LDAC_IMHEAD", "LDAC_OBJECTS"] - cards = hdul["LDAC_IMHEAD"].data[0][0] - assert isinstance(cards, str) or cards.ndim == 1 - text = "".join(cards) if not isinstance(cards, str) else cards - assert "TILEKEY" in text and "2605805p" in text - - -def test_number_is_kept_and_aliases_added(ldac): - with fits.open(ldac) as hdul: - data = hdul["LDAC_OBJECTS"].data - # The input NUMBER (gapped here: 1, 2, 7) is the object's identity and is - # copied unchanged; the converter builds no ID of its own. - npt.assert_array_equal(data["NUMBER"], [o[0] for o in OBJECTS]) - assert "TILE_UNIQUE_ID" not in data.names - npt.assert_array_equal(data["XWIN_IMAGE"], data["X_IMAGE"]) - npt.assert_array_equal(data["YWIN_IMAGE"], data["Y_IMAGE"]) - npt.assert_array_equal(data["XWIN_WORLD"], data["ALPHA_J2000"]) - npt.assert_array_equal(data["YWIN_WORLD"], data["DELTA_J2000"]) - - -def test_vignets_are_cut_from_the_image_and_padded_as_sextractor(ldac): - with fits.open(ldac) as hdul: - vignets = hdul["LDAC_OBJECTS"].data["VIGNET"] - assert vignets.shape == (len(OBJECTS), STAMP, STAMP) - - # Interior object at (10, 12), 1-based: centre pixel is (row 11, col 9). - assert vignets[0, STAMP // 2, STAMP // 2] == 11 * 1000 + 9 - assert vignets[0, 0, 0] == 9 * 1000 + 7 - - # Left edge, x = 1: the two columns left of the image are -1e30, the value - # SExtractor writes off the image. - assert (vignets[1, :, :2] == rs.BIG).all() - assert vignets[1, STAMP // 2, STAMP // 2] == 19 * 1000 + 0 - - # Top-right corner: only the lower-left quadrant of the stamp is in the - # image. - assert (vignets[2, STAMP // 2 + 1:, :] == rs.BIG).all() - assert (vignets[2, :, STAMP // 2 + 1:] == rs.BIG).all() - assert vignets[2, STAMP // 2, STAMP // 2] == 29 * 1000 + 39 - - -def test_tile_unique_id_reaches_the_final_catalogue(ldac, tmp_path): - """make_cat builds the ID from the tile and the catalogue's own NUMBER.""" - final = make_cat.prepare_final_cat_file(str(tmp_path), "-301-279") - n_obj = make_cat.save_sextractor_data(final, str(ldac)) - assert n_obj == len(OBJECTS) - with fits.open(tmp_path / "final_cat-301-279.fits") as hdul: - data = hdul["RESULTS"].data - assert "VIGNET" not in data.names - npt.assert_array_equal( - data["TILE_UNIQUE_ID"], 301279 * 10**6 + np.array([1, 2, 7]) - ) - npt.assert_allclose(data["TILE_ID"], 301.279) - - -# --- the segmentation map ------------------------------------------------- - - -def _seg_map(): - """Two footprints, labelled 7 and 9, on a 20x20 sky.""" - seg = np.zeros((20, 20), dtype=np.int32) - seg[2:6, 2:6] = 7 - seg[12:18, 12:18] = 9 - return seg - - -@pytest.mark.decision("detection.catalogue_neighbour_marking") -class TestRelabel: - """The relabelled map carries each object's NUMBER on its own footprint. - - Segmentation labels are not the catalogue's NUMBER; the claim is by the - pixel under the object's position, which is all both maps share. - """ - - def test_each_object_owns_the_footprint_it_sits_in(self): - seg = _seg_map() - # Positions in an order that is NOT the label order, so a - # relabelling that merely renumbered would fail. - out, counts = rs.relabel_seg( - seg, number=np.array([1, 2]), - x_image=np.array([15.0, 4.0]), y_image=np.array([15.0, 4.0])) - assert counts["matched"] == 2 - assert set(np.unique(out[seg == 9])) == {1} - assert set(np.unique(out[seg == 7])) == {2} - assert np.all(out[seg == 0] == 0) - - def test_an_unclaimed_footprint_becomes_a_neighbour(self): - seg = _seg_map() - out, _ = rs.relabel_seg( - seg, number=np.array([1]), - x_image=np.array([4.0]), y_image=np.array([4.0])) - assert set(np.unique(out[seg == 9])) == {rs.NEIGHBOUR_LABEL} - assert set(np.unique(out[seg == 7])) == {1} - - def test_an_object_on_sky_gets_a_disc(self): - seg = _seg_map() - out, counts = rs.relabel_seg( - seg, number=np.array([1, 5]), - x_image=np.array([4.0, 10.0]), y_image=np.array([4.0, 10.0]), - fallback_radius=2) - assert counts["unclaimed"] == 1 - assert out[9, 9] == 5 - assert np.count_nonzero(out == 5) == np.count_nonzero( - np.add.outer(np.arange(-2, 3) ** 2, np.arange(-2, 3) ** 2) <= 4) - - def test_two_objects_in_one_footprint_both_keep_a_centre(self): - seg = _seg_map() - out, counts = rs.relabel_seg( - seg, number=np.array([1, 2]), - x_image=np.array([14.0, 16.0]), y_image=np.array([14.0, 16.0]), - fallback_radius=1) - assert counts == dict(matched=1, unclaimed=0, shared=1, off_image=0, - shared_pixel=0) - assert out[13, 13] == 1 and out[15, 15] == 2 - assert np.count_nonzero(out == 1) > np.count_nonzero(out == 2) - - def test_every_object_is_self_somewhere(self): - rng = np.random.default_rng(0) - seg = np.zeros((60, 60), dtype=np.int32) - for label in range(1, 12): - row, col = rng.integers(0, 55, size=2) - seg[row:row + 5, col:col + 5] = label - number = np.arange(1, 31) - x_image = rng.uniform(1, 60, size=30) - y_image = rng.uniform(1, 60, size=30) - out, counts = rs.relabel_seg(seg, number, x_image, y_image) - assert sum(counts[k] for k in - ("matched", "unclaimed", "shared", "off_image")) == 30 - for num in number: - assert np.any(out == num), f"object {num} has no self pixels" - assert set(np.unique(out)) <= set(number) | {0, rs.NEIGHBOUR_LABEL} - - def test_two_objects_on_one_pixel(self): - seg = _seg_map() - out, counts = rs.relabel_seg( - seg, number=np.array([4, 6]), x_image=np.array([10.0, 10.2]), - y_image=np.array([10.0, 10.1]), fallback_radius=1) - assert counts["shared_pixel"] == 1 - assert out[9, 9] == 4 - assert np.any(out == 6) - - @pytest.mark.parametrize("x, y", [(0.4, 5.0), (5.0, 21.0)]) - def test_a_position_off_the_image_is_counted(self, x, y): - out, counts = rs.relabel_seg( - _seg_map(), number=np.array([1]), x_image=np.array([x]), - y_image=np.array([y])) - assert counts["off_image"] == 1 - assert not np.any(out == 1) - - -# Objects on _seg_map(): 1 in footprint 7, 2 on sky between the footprints; -# footprint 9 is claimed by nobody. -SEG_OBJECTS = [(1, 4.0, 4.0), (2, 10.0, 9.0)] -SEG_STAMP = 21 - - -def _marked(number, x, y, seg=None): - image = np.ones(_seg_map().shape, np.float32) - seg = _seg_map() if seg is None else seg - relabelled, _ = rs.relabel_seg(seg, number, x, y) - return rs._extract_vignets(image, x, y, SEG_STAMP, seg=relabelled, - number=number) - - -def _stamp_of(array, x, y, fill): - """The SEG_STAMP stamp of ``array`` centred on 1-based (x, y).""" - half = SEG_STAMP // 2 - padded = np.pad(array, half, constant_values=fill) - row, col = int(round(y)) - 1, int(round(x)) - 1 - return padded[row:row + SEG_STAMP, col:col + SEG_STAMP] - - -@pytest.mark.decision("detection.catalogue_neighbour_marking") -def test_neighbour_footprints_become_big_and_nothing_else(): - number = np.array([o[0] for o in SEG_OBJECTS]) - x = np.array([o[1] for o in SEG_OBJECTS]) - y = np.array([o[2] for o in SEG_OBJECTS]) - vignets = _marked(number, x, y) - seg = _seg_map() - - # Object 1 in footprint 7: footprint 9 (unclaimed) is a neighbour; its - # own footprint and the sky keep their image values. - s = _stamp_of(seg, x[0], y[0], fill=-99) - v = vignets[0] - assert (v[s == 9] == rs.BIG).all() - assert (v[s == 7] == 1).all() - assert (v[s == -99] == rs.BIG).all() # off the image - # Object 2's disc is a catalogue object's pixels, so a neighbour too; the - # rest of the sky is untouched. - relabelled, _ = rs.relabel_seg(seg, number, x, y) - disc = _stamp_of(relabelled, x[0], y[0], fill=-99) == 2 - assert disc.any() and (v[disc] == rs.BIG).all() - assert (v[(s == 0) & ~disc] == 1).all() - - # Object 2 on sky: every footprint in its stamp is a neighbour; the sky, - # its own centre included, keeps its image values. - s = _stamp_of(seg, x[1], y[1], fill=-99) - v = vignets[1] - assert (v[(s == 7) | (s == 9)] == rs.BIG).all() - assert (v[s == 0] == 1).all() - assert v[SEG_STAMP // 2, SEG_STAMP // 2] == 1 - - -@pytest.mark.decision("detection.catalogue_neighbour_marking") -def test_a_claimed_neighbour_is_masked_by_its_number(): - """With both footprints claimed, each object masks the other's.""" - number, x, y = np.array([3, 8]), np.array([4.0, 15.0]), np.array([4.0, 15.0]) - vignets = _marked(number, x, y) - seg = _seg_map() - for i, (own, other) in enumerate([(7, 9), (9, 7)]): - s = _stamp_of(seg, x[i], y[i], fill=-99) - assert (vignets[i][s == other] == rs.BIG).all() - assert (vignets[i][s == own] == 1).all() - - -@pytest.mark.decision("detection.catalogue_neighbour_marking") -def test_converter_relabels_and_marks_from_compressed_map(tmp_path): - """End to end from a compressed map, as fetched from vos.""" - cat = tmp_path / "CFIS_cat-301-279.cat" - img = tmp_path / "CFIS_image-301-279.fits" - seg_in = tmp_path / "CFIS_seg-301-279.fitsfz" - out = tmp_path / "sexcat-301-279.fits" - lines = ["# 1 NUMBER", "# 2 X_IMAGE", "# 3 Y_IMAGE", - "# 4 ALPHA_J2000", "# 5 DELTA_J2000"] - lines += [f"{n} {x} {y} 150.0 30.0" for n, x, y in SEG_OBJECTS] - cat.write_text("\n".join(lines) + "\n") - fits.PrimaryHDU(np.ones((20, 20), np.float32)).writeto(img) - fits.HDUList([fits.PrimaryHDU(), - fits.CompImageHDU(_seg_map())]).writeto(seg_in) - - rs.make_ldac_from_ascii(str(cat), str(img), str(out), stamp_size=SEG_STAMP, - seg_path=str(seg_in)) - - seg = _seg_map() - number, x, y = (np.array(c) for c in zip(*SEG_OBJECTS)) - relabelled, _ = rs.relabel_seg(seg, number, x, y) - assert set(np.unique(relabelled[seg == 7])) == {1} - assert set(np.unique(relabelled[seg == 9])) == {rs.NEIGHBOUR_LABEL} - with fits.open(out) as hdul: - v = hdul["LDAC_OBJECTS"].data["VIGNET"][0] - s = _stamp_of(seg, 4.0, 4.0, fill=-99) - assert (v[s == 9] == rs.BIG).all() and (v[s == 7] == 1).all() - - fits.HDUList([fits.PrimaryHDU(), - fits.CompImageHDU(_seg_map()[:10])]).writeto( - seg_in, overwrite=True) - with pytest.raises(ValueError, match="one grid"): - rs.make_ldac_from_ascii(str(cat), str(img), str(out), - stamp_size=SEG_STAMP, seg_path=str(seg_in)) - - -DR6_PATCH = Path(__file__).parent / "data" / "dr6_202.301_seg_patch.fits" - - -@pytest.mark.decision("detection.catalogue_neighbour_marking") -def test_dr6_marks_every_neighbour_pixel_and_no_own_pixel(): - """On a crowded 200x200 patch of the real 202.301 map and catalogue. - - The raw labels are not NUMBER, so this checks the whole chain on real - data: 1-based positions, the centre-pixel claim, and the stamp geometry. - Every stamp fully on the patch is checked against the RAW map: pixels of - footprints other than the one under the object are all -1e30, its own - footprint and the sky are untouched. (SExtractor's own VIGNET, on a - SExtractor seg map of an image-sim tile, marks 93% of neighbour-label - pixels and 0.14% of own pixels: check_sex_vignet.py.) - """ - with fits.open(DR6_PATCH) as hdul: - seg = hdul["SEG"].data - objects = hdul["OBJECTS"].data - number = np.array(objects["NUMBER"]) - x, y = np.array(objects["X_IMAGE"]), np.array(objects["Y_IMAGE"]) - relabelled, counts = rs.relabel_seg(seg, number, x, y) - assert counts["matched"] == len(number) - - stamp, half = 51, 25 - image = np.ones(seg.shape, np.float32) - vignets = rs._extract_vignets(image, x, y, stamp, seg=relabelled, - number=number) - col, row = np.rint(x).astype(int) - 1, np.rint(y).astype(int) - 1 - full = ((col >= half) & (col < seg.shape[1] - half) - & (row >= half) & (row < seg.shape[0] - half)) - assert full.sum() >= 10 - - marked = {"neighbour": [0, 0], "own": [0, 0], "sky": [0, 0]} - for i in np.flatnonzero(full): - s = seg[row[i] - half:row[i] + half + 1, col[i] - half:col[i] + half + 1] - big = vignets[i] == rs.BIG - own = s[half, half] - for kind, where in (("neighbour", (s != 0) & (s != own)), - ("own", s == own), ("sky", s == 0)): - marked[kind][0] += (big & where).sum() - marked[kind][1] += where.sum() - assert marked["neighbour"][1] > 1000 - assert marked["neighbour"][0] == marked["neighbour"][1] - assert marked["own"][0] == 0 - assert marked["sky"][0] == 0 - - -# --- the runner's output against the completeness table ------------------- - - -def test_runner_output_matches_tile_detect_completeness(tmp_path, monkeypatch): - """The runner writes exactly the files ``tile_detect`` expects. - - Runs the real runner, segmentation map on, into a run dir and checks it - with ``completeness.check_counts`` under ``SP_TILE_DETECTION=unions_catalogue``, - so the table and the converter's outputs cannot drift apart. - """ - import configparser - import importlib.util - import logging - - from shapepipe.modules.read_ext_sexcat_runner import read_ext_sexcat_runner - - scripts = Path(__file__).resolve().parents[2] / "workflow" / "scripts" - spec = importlib.util.spec_from_file_location( - "_completeness", scripts / "completeness.py") - completeness = importlib.util.module_from_spec(spec) - spec.loader.exec_module(completeness) - - cat = tmp_path / "CFIS_cat-301-279.cat" - img = tmp_path / "CFIS_image-301-279.fits" - seg_in = tmp_path / "CFIS_seg-301-279.fitsfz" - lines = ["# 1 NUMBER", "# 2 X_IMAGE", "# 3 Y_IMAGE", - "# 4 ALPHA_J2000", "# 5 DELTA_J2000"] - lines += [f"{n} {x} {y} 150.0 30.0" for n, x, y in SEG_OBJECTS] - cat.write_text("\n".join(lines) + "\n") - fits.PrimaryHDU(np.ones((20, 20), np.float32)).writeto(img) - fits.HDUList([fits.PrimaryHDU(), - fits.CompImageHDU(_seg_map())]).writeto(seg_in) - - run_dir = tmp_path / "run_sp_tile_Rx" - out_dir = run_dir / "read_ext_sexcat_runner" / "output" - out_dir.mkdir(parents=True) - config = configparser.ConfigParser() - config["READ_EXT_SEXCAT_RUNNER"] = { - "SEGMENTATION": "True", "MAKE_POST_PROCESS": "False", - "VIGNET_SIZE": str(SEG_STAMP), - } - read_ext_sexcat_runner( - [str(cat), str(img), str(seg_in)], {"output": str(out_dir)}, - "-301-279", config, "READ_EXT_SEXCAT_RUNNER", - logging.getLogger("test"), - ) - - monkeypatch.setenv("SP_TILE_DETECTION", "unions_catalogue") - table = completeness.COMPLETENESS["tile_detect"]["unions_catalogue"] - expect = table["read_ext_sexcat_runner"]["expect"] - written = sorted(p.name for p in out_dir.iterdir()) - assert len(written) == expect, written - ok, details = completeness.check_counts("tile_detect", run_dir) - assert ok, details diff --git a/tests/module/test_sextractor_command_line.py b/tests/module/test_sextractor_command_line.py new file mode 100644 index 000000000..47488df0d --- /dev/null +++ b/tests/module/test_sextractor_command_line.py @@ -0,0 +1,81 @@ +"""SExtractor command line built by ``SExtractorCaller``. + +Without a detection image the call is single-image (``img -WEIGHT_IMAGE w``); +dual-image mode (``det,img -WEIGHT_IMAGE det_w,w``) only with a detection +image, since dual-image mode on the same image twice changes the measurements. +""" + +import shlex + +import pytest + +from shapepipe.modules.sextractor_package.sextractor_script import ( + SExtractorCaller, +) + + +def _command(inputs, **flags): + options = dict( + use_weight=False, use_flag=False, use_psf=False, + use_detection_image=False, use_detection_weight=False, + use_zero_point=False, use_background=False, + ) + options.update(flags) + caller = SExtractorCaller( + inputs, "/out", "-000-000", "d.sex", "d.param", "d.conv", + check_image=[""], **options, + ) + argv = shlex.split(caller.make_command_line("source-extractor")) + options = {argv[i]: argv[i + 1] for i in range(2, len(argv) - 1) + if argv[i].startswith("-")} + return argv, options + + +def test_tile_single_image_with_weight(): + argv, opt = _command(["img.fits", "w.fits"], use_weight=True) + assert argv[:2] == ["source-extractor", "img.fits"] + assert opt["-WEIGHT_IMAGE"] == "w.fits" + assert "-FLAG_IMAGE" not in opt + assert opt["-CATALOG_NAME"] == "/out/sexcat-000-000.fits" + assert opt["-CHECKIMAGE_TYPE"] == "NONE" + + +def test_exposure_single_image_with_weight_and_flag(): + argv, opt = _command(["img.fits", "w.fits", "f.fits"], + use_weight=True, use_flag=True) + assert argv[1] == "img.fits" + assert opt["-WEIGHT_IMAGE"] == "w.fits" + assert opt["-FLAG_IMAGE"] == "f.fits" + + +def test_single_image_without_weight(): + argv, opt = _command(["img.fits"]) + assert argv[1] == "img.fits" + assert opt["-WEIGHT_TYPE"] == "None" + assert "-WEIGHT_IMAGE" not in opt + + +@pytest.mark.parametrize("detection_weight", [False, True]) +def test_dual_image_with_detection_image(detection_weight): + inputs = ["img.fits", "w.fits", "det.fits"] + if detection_weight: + inputs.append("det_w.fits") + argv, opt = _command(inputs, use_weight=True, use_detection_image=True, + use_detection_weight=detection_weight) + assert argv[1] == "det.fits,img.fits" + det_w = "det_w.fits" if detection_weight else "w.fits" + assert opt["-WEIGHT_IMAGE"] == f"{det_w},w.fits" + + +def test_check_images_named_per_unit(): + caller = SExtractorCaller( + ["img.fits", "w.fits"], "/out", "-001-002", "d.sex", "d.param", + "d.conv", True, False, False, False, False, False, False, + check_image=["BACKGROUND", "SEGMENTATION"], + ) + argv = shlex.split(caller.make_command_line("source-extractor")) + assert argv[1] == "img.fits" + assert argv[argv.index("-CHECKIMAGE_TYPE") + 1] == "BACKGROUND,SEGMENTATION" + assert argv[argv.index("-CHECKIMAGE_NAME") + 1] == ( + "/out/background-001-002.fits,/out/segmentation-001-002.fits" + ) diff --git a/tests/unit/test_config_parse.py b/tests/unit/test_config_parse.py index 4aff1efbf..458448ac3 100644 --- a/tests/unit/test_config_parse.py +++ b/tests/unit/test_config_parse.py @@ -251,3 +251,38 @@ def test_module_sections_check_catches_a_missing_repeated_section(): problems = _missing_required_keys(parser) assert any("VIGNETMAKER_RUNNER_RUN_1" in problem for problem in problems) + + +def _expanded(value): + from shapepipe.pipeline.config import CustomParser + + parser = CustomParser() + parser.read_dict({"S": {"K": value}}) + return parser.getexpanded("S", "K") + + +def test_env_default_applies_when_unset_or_empty(monkeypatch): + """${VAR:-x} is x when VAR is unset and when it is empty, as in the shell.""" + monkeypatch.delenv("SP_TEST_VAR", raising=False) + assert _expanded("${SP_TEST_VAR:-x}") == "x" + monkeypatch.setenv("SP_TEST_VAR", "") + assert _expanded("${SP_TEST_VAR:-x}") == "x" + assert _expanded("${SP_TEST_VAR:-}") == "" + assert _expanded("${SP_TEST_VAR}") == "" + + +def test_env_value_is_inserted_unexpanded(monkeypatch): + """A value holding '$' or '${' is inserted as is, not expanded again.""" + monkeypatch.setenv("SP_TEST_VAR", "/a/$HOME/${X:=y}") + assert _expanded("$SP_TEST_VAR/b") == "/a/$HOME/${X:=y}/b" + assert _expanded("${SP_TEST_VAR:-z}") == "/a/$HOME/${X:=y}" + + +@pytest.mark.parametrize( + "value", ["${SP_TEST_VAR-x}", "${SP_TEST_VAR:=x}", "${ SP_TEST_VAR }", + "/p/${SP_TEST_VAR"]) +def test_malformed_brace_is_rejected(monkeypatch, value): + """A '${' that is not $VAR, ${VAR} or ${VAR:-x} fails, set or not.""" + monkeypatch.setenv("SP_TEST_VAR", "v") + with pytest.raises(ValueError, match="not \\$VAR"): + _expanded(value) diff --git a/tests/unit/test_final_cat_merge_invariants.py b/tests/unit/test_final_cat_merge_invariants.py index 21eb6adad..7dcfc6183 100644 --- a/tests/unit/test_final_cat_merge_invariants.py +++ b/tests/unit/test_final_cat_merge_invariants.py @@ -130,7 +130,7 @@ def _campaign(root: Path, drop=None): argv = [sys.executable, str(SCRIPT), "--products-dir", str(products), "--tile-list", str(tile_list), "--index-db", str(index), "--output", str(output), "--campaign", CAMPAIGN, - "--param-file", str(param), "--tile-detection", "sextractor"] + "--param-file", str(param)] return argv, output, sources diff --git a/tests/unit/test_postage_stamp_size.py b/tests/unit/test_postage_stamp_size.py index 0bd31b568..ca7d4af72 100644 --- a/tests/unit/test_postage_stamp_size.py +++ b/tests/unit/test_postage_stamp_size.py @@ -11,12 +11,6 @@ CFIS_CONFIG = CONFIG / "cfis" -def _ini_int(path, section, option): - parser = configparser.ConfigParser(interpolation=None) - assert parser.read(path) == [str(path)] - return parser.getint(section, option) - - def _multi_epoch_stamp_size(path): """STAMP_SIZE of the one vignetmaker section run in MULTI-EPOCH mode.""" parser = configparser.ConfigParser(interpolation=None) @@ -50,11 +44,6 @@ def test_tile_and_epoch_vignet_sizes_match(): """The tile markers align pixel-for-pixel with each multi-epoch stamp.""" param_file = CFIS_CONFIG / "default_noimaflags.param" sizes = { - "config_tile_Uc.ini#READ_EXT_SEXCAT_RUNNER.VIGNET_SIZE": _ini_int( - CFIS_CONFIG / "config_tile_Uc.ini", - "READ_EXT_SEXCAT_RUNNER", - "VIGNET_SIZE", - ), f"{param_file.name}#VIGNET": _param_vignet_size(param_file), } diff --git a/tests/unit/test_workflow_tile_detection.py b/tests/unit/test_workflow_tile_detection.py index 24775304e..498a015de 100644 --- a/tests/unit/test_workflow_tile_detection.py +++ b/tests/unit/test_workflow_tile_detection.py @@ -1,13 +1,13 @@ -"""The two tile_detect modes agree on where the sexcat is. - -``tile_detection: unions_catalogue`` swaps the SExtractor rule for a fetch plus -a conversion, and everything downstream of ``tile_detect`` is unchanged only -because both modes write ``run_sp_tile_Sx`` and both are checked as the -``tile_detect`` stage. The agreements that make that true are between files -that never see each other at run time: the two inis' RUN_NAMEs, -``completeness.STAGE_DIR``, the flavoured ``COMPLETENESS['tile_detect']`` -table, and ``run_report``'s stage list. They are asserted here, statically. -Container-free: the scripts are stdlib-only. +"""The two tile_detection modes run one tile_detect, and agree on its inputs. + +Both modes run SExtractor (config_tile_Sx.ini) in the ``run_sp_tile_Sx`` stage +dir; ``tile_detection: unions_catalogue`` adds ``tile_get_catalogue`` and +points the ini's MATCH_CATALOGUE at its output through the workflow's +SP_MATCH_CATALOGUE. The agreements that make that work are between files that +never see each other at run time: the inis' RUN_NAMEs and fetch patterns, +``completeness``'s tables, ``run_report``'s stage list and tile.smk's path. +They are asserted here, statically. Container-free: the scripts are +stdlib-only. """ import configparser @@ -44,43 +44,35 @@ def _ini(name): def test_modes_are_the_two_the_rules_know(): assert completeness.TILE_DETECTIONS == ("sextractor", "unions_catalogue") - assert set(completeness.COMPLETENESS["tile_detect"]) == set( - completeness.TILE_DETECTIONS) -def test_both_detection_inis_write_the_stage_dir(): - """config_tile_Sx and config_tile_Uc share the run dir STAGE_DIR names.""" +def test_detection_ini_writes_the_stage_dir(): level, subdir = completeness.STAGE_DIR["tile_detect"] assert level == "tile" - for name in ("config_tile_Sx.ini", "config_tile_Uc.ini"): - assert _ini(name)["DEFAULT"]["RUN_NAME"].strip() == subdir, ( - f"{name} writes a run dir other than {subdir}: unit_pre would clear " - "the wrong directory and the chain downstream would read nothing.") - assert _ini("config_tile_Uc.ini")["DEFAULT"]["RUN_DATETIME"].strip() == "False" + assert _ini("config_tile_Sx.ini")["DEFAULT"]["RUN_NAME"].strip() == subdir -def test_fetch_ini_matches_its_stage_and_feeds_the_converter(): +def test_detection_joins_only_when_the_workflow_says_so(): + """MATCH_CATALOGUE is empty unless SP_MATCH_CATALOGUE is exported.""" + sx = _ini("config_tile_Sx.ini")["SEXTRACTOR_RUNNER"] + assert sx["MATCH_CATALOGUE"].strip() == "${SP_MATCH_CATALOGUE:-}" + assert float(sx["MATCH_RADIUS"]) == 1.0 + assert 0.9 < float(sx["MATCH_MIN_FRACTION"]) <= 1 + # The join renumbers before the post-processing keys epochs on NUMBER. + assert sx["MAKE_POST_PROCESS"].strip() == "True" + + +def test_fetch_ini_writes_the_catalogue_tile_smk_joins(): + """Gic's output is the path tile.smk exports as SP_MATCH_CATALOGUE.""" level, gic = completeness.STAGE_DIR["tile_get_catalogue"] assert level == "tile" assert _ini("config_tile_Gic.ini")["DEFAULT"]["RUN_NAME"].strip() == gic - uc = _ini("config_tile_Uc.ini")["READ_EXT_SEXCAT_RUNNER"] - assert f"$SP_RUN/output/{gic}/get_images_runner/output" in uc["INPUT_DIR"] - assert uc["FILE_PATTERN"].split(",")[0].strip() == "CFIS_cat" - # The segmentation map rides the same fetch and is the converter's third - # input, the position SEGMENTATION = True reads it from. gi = _ini("config_tile_Gic.ini")["GET_IMAGES_RUNNER"] - assert [p.strip() for p in gi["OUTPUT_FILE_PATTERN"].split(",")] == [ - "CFIS_cat-", "CFIS_seg-"] - assert [e.strip() for e in gi["INPUT_FILE_EXT"].split(",")] == [ - ".cat", ".fits.fz"] - third = [[v.strip() for v in uc[k].split(",")][2] - for k in ("INPUT_DIR", "FILE_PATTERN", "FILE_EXT")] - assert third == [f"$SP_RUN/output/{gic}/get_images_runner/output", - "CFIS_seg", ".fitsfz"] - assert uc["SEGMENTATION"].strip() == "True" - # The multi-epoch post-processing is what gives the sexcat its EPOCH_k - # extensions; ngmix_range.py refuses a sexcat without them. - assert uc["MAKE_POST_PROCESS"].strip() == "True" + assert gi["OUTPUT_FILE_PATTERN"].strip() == "CFIS_cat-" + assert gi["INPUT_FILE_EXT"].strip() == ".cat" + smk = (REPO_ROOT / "workflow" / "rules" / "tile.smk").read_text() + assert f'/output/{gic}/get_images_runner/output"' in smk + assert 'f"{gic}/CFIS_cat{unit_num(tile)}.cat"' in smk def _stage_dir(tmp_path, runner, n): @@ -91,31 +83,16 @@ def _stage_dir(tmp_path, runner, n): return tmp_path -@pytest.mark.parametrize("mode, runner, expect", [ - ("sextractor", "sextractor_runner", 2), - ("unions_catalogue", "read_ext_sexcat_runner", 1), +@pytest.mark.parametrize("stage, runner, expect", [ + ("tile_detect", "sextractor_runner", 2), + ("tile_get_catalogue", "get_images_runner", 1), ]) -def test_tile_detect_is_checked_per_mode(tmp_path, monkeypatch, mode, runner, expect): - monkeypatch.setenv("SP_TILE_DETECTION", mode) +def test_stage_counts(tmp_path, stage, runner, expect): ok, details = completeness.check_counts( - "tile_detect", _stage_dir(tmp_path, runner, expect)) + stage, _stage_dir(tmp_path, runner, expect)) assert ok and details == [(runner, expect, expect, False)] -def test_unset_mode_is_sextractor(tmp_path, monkeypatch): - """A prologue without the export checks as before: data runs are unchanged.""" - monkeypatch.delenv("SP_TILE_DETECTION", raising=False) - ok, details = completeness.check_counts( - "tile_detect", _stage_dir(tmp_path, "sextractor_runner", 2)) - assert ok and details[0][0] == "sextractor_runner" - - -def test_invalid_mode_is_fatal(tmp_path, monkeypatch): - monkeypatch.setenv("SP_TILE_DETECTION", "steven") - with pytest.raises(ValueError, match="SP_TILE_DETECTION"): - completeness.check_counts("tile_detect", tmp_path) - - @pytest.mark.parametrize("mode, present", [ ("sextractor", False), ("unions_catalogue", True), (None, False), ]) diff --git a/tests/workflow/params_pin.json b/tests/workflow/params_pin.json index eeb61d950..6de487277 100644 --- a/tests/workflow/params_pin.json +++ b/tests/workflow/params_pin.json @@ -8,10 +8,10 @@ "exp_persist": "302e2837542bc1102430c27c81c600b7cda32e8bddcb5fd60d33950987609fff", "exp_psf": "2c4f6d00f1939ccbf05b4982a0202a0ff92a4373a727f4e00f3aaab4eba03352", "exp_split": "6e954f8f3d06f44d3f164675912ce27d6216648855d04968f9168bd7d0f2c4fa", - "final_cat_merge": "f6222c98be9777cbab4131d8857bd3817719df9302d0a321c6cd99f67d5bdc76", + "final_cat_merge": "e7f46859c4503a2220713d7bb2507555515d0a9632d780b20f14c59e32210023", "prepare_all_tiles": "b8f872a22adf014e25a7fa5198f49b71a6fe9e56042ed82b682bc8763970a844", "star_cat_merge": "6277450958474af5270982fa35360f2f237a29f7533c526ee9265dfd5acc07a0", - "tile_detect": "1b4b0871bf28296f8e5d6855ae47ff2cb84542c93593a5ae5c0ffc3c84c45479", + "tile_detect": "8336b148769e9d43b64f0d945c05e7d163c64dba4e6287e8c1964d1502e3ffe1", "tile_exp_forest": "7447ab4a1049de5f0b5c81e5f9ed2a8644c7bdab85cb0a060bfde81261a89b28", "tile_find_exposures": "8704317871744996c44351c2836fcb222d7986a602e9e046a90d684ba7b3c184", "tile_get_catalogue": "7e3f889a955a14b0b917a015c2a8bc90e433a89b1cc5e00b8ca94695c4c93d3c", @@ -24,7 +24,7 @@ "tile_vignets": "9d4ae0d99c18217f2185f245281312454c8a219ec1628176e08c71a5efc4dc91" }, "schema": 1, - "sha256": "f9772987e517e9827c8d501d310ec323be251a62b7b89aad0725a03f8f4f120e", + "sha256": "d67ee5da728e8fee513e1f01a522b386eb6a86d21c9a36250afd6013a06301d5", "unit_pre": { "exp_get_images": "8dec850af212879f225fcf27a5f1281e1a075264c7b97d38c2214395d360168c", "exp_psf": "f2358ddf7385918dc5033d10b37f6dc97a15d02b071a3ea0a4619a5f7e6f5bec", diff --git a/tests/workflow/test_dag.py b/tests/workflow/test_dag.py index 46ae93fc0..dd5fc7f70 100644 --- a/tests/workflow/test_dag.py +++ b/tests/workflow/test_dag.py @@ -31,6 +31,24 @@ def test_rule_set_matches_input_mode(campaign, dag): assert "merge_final_cats" not in dag.declared_rule_names +def test_tile_detect_joins_the_catalogue_iff_unions(campaign, dag): + """Data's tile_detect waits on the fetch and exports its catalogue as + SP_MATCH_CATALOGUE; the image-simulation prologue exports it empty.""" + for job in dag.jobs_for("tile_detect"): + tile = job.wildcards.tile + inputs = {str(f) for f in job.input} + fetch = str(campaign.tile_manifest(tile, "tile_get_catalogue")) + if campaign.tile_detection == "unions_catalogue": + gic = (campaign.run_dir / "tiles" / tile[:2] / tile / "output" + / "run_sp_tile_Gic" / "get_images_runner" / "output") + cat = gic / f"CFIS_cat-{tile.replace('.', '-')}.cat" + assert fetch in inputs + assert f"export SP_MATCH_CATALOGUE='{cat}'" in job.params.pre + else: + assert fetch not in inputs + assert "export SP_MATCH_CATALOGUE=''" in job.params.pre.split("\n") + + def test_clean_exposure_waits_on_persist_iff_psf(campaign, dag): """Reclamation waits for persistence and exactly its in-scope readers.""" jobs = dag.jobs_for("clean_exposure") diff --git a/universes/committed.yaml b/universes/committed.yaml index bc34d7601..05d88e097 100644 --- a/universes/committed.yaml +++ b/universes/committed.yaml @@ -24,7 +24,6 @@ analyses: photometry_parameters: kron_25_35 detection_source_mode: sx_nomask_single_image epoch_membership_ccd_bounds: trimmed_bounds_33_2080 - catalogue_neighbour_marking: segmentation_map preparation: decisions: astrometric_solution_source: delivered_headers diff --git a/workflow/README.md b/workflow/README.md index a1654fcb4..4f0b135c3 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -28,13 +28,13 @@ uv pip install 'snakemake>=9,<10' 'snakemake-executor-plugin-slurm>=2.7,<3' # `psf_model` is `psfex` for data (`fake` for image sims). `mccd` is refused # until `persist_exp.py` and `merge_star_cat.py` read MCCD products. PSFEx is # exercised by smk-g4 through smk-g6. -# `tile_detection` is `unions_catalogue` (the input_types default for data: -# the UNIONS per-tile catalogue at `inputs.catalogues` is fetched and -# converted in place, keeping its NUMBER, and its segmentation map sets -# neighbours' VIGNET pixels to -1e30 as SExtractor does) or `sextractor` (the -# tile is detected with SExtractor; the default for image sims). Either way -# make_cat writes TILE_UNIQUE_ID = tile_id * 10**6 + NUMBER and ngmix masks -# the same neighbours. +# Both `tile_detection` values run SExtractor on the tile image. Under +# `unions_catalogue` (the input_types default for data) the UNIONS per-tile +# catalogue at `inputs.catalogues` is fetched and the detections are joined +# to it, taking its NUMBER; the tile image must be the catalogue's release +# (DR6), or the join fails. `sextractor` (the default for image sims) keeps +# SExtractor's own NUMBER. Either way make_cat writes +# TILE_UNIQUE_ID = tile_id * 10**6 + NUMBER. # The committed launcher loads apptainer/1.4.5 + the /project venv, so a # fresh shell always has the right state. @@ -240,7 +240,7 @@ workflow/ rules/ prepare.smk tile get_images/uncompress/find_exposures exposure.smk per-exposure: get_images, split, psf, persist (no temp()); campaign star_cat_merge - tile.smk per-tile: exp forest, merge_headers, detect (SExtractor, or fetch + convert the UNIONS catalogue), vignets, ngmix, merge, make_cat; campaign final_cat_merge + tile.smk per-tile: exp forest, merge_headers, detect (SExtractor, joined to the UNIONS catalogue on data), vignets, ngmix, merge, make_cat; campaign final_cat_merge scripts/ build_index.py prepare-phase run_index.sqlite builder (plain script) build_forest.py per-tile exposure symlink forest (group-compatible shell) diff --git a/workflow/Snakefile b/workflow/Snakefile index 6aed81cff..ac9d0760c 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -175,9 +175,9 @@ sys.path.insert(0, str(SCRIPTS)) import build_index # noqa: E402 from completeness import STAGE_DIR, TILE_DETECTIONS # noqa: E402 -# Where the tile's galaxy sample comes from: SExtractor on the tile image, or -# the UNIONS per-tile catalogue fetched and converted in place (tile.smk). -# Both write run_sp_tile_Sx, so nothing downstream of tile_detect branches. +# How the tile's galaxy sample is defined: SExtractor on the tile image, alone +# or joined to the UNIONS per-tile catalogue (tile.smk). Both write +# run_sp_tile_Sx, so nothing downstream of tile_detect branches. TILE_DETECTION = config.get("tile_detection", "sextractor") if TILE_DETECTION not in TILE_DETECTIONS: raise WorkflowError( diff --git a/workflow/config.yaml b/workflow/config.yaml index b6f22a01e..8b36c70fa 100644 --- a/workflow/config.yaml +++ b/workflow/config.yaml @@ -21,9 +21,10 @@ input_type: data # run: smk-g6 # Default settings per input_type, can be overridden by a user-defined run config file. -# - tile_detection: allowed are `unions_catalogue` (fetch official UNIONS catalogue -# from vos; recommended for data); `sextractor` (detect objects via a SExtractor -# module run; required for image_sims) +# - tile_detection: both run SExtractor on the tile. `unions_catalogue` (data) +# also fetches the UNIONS per-tile catalogue and joins the detections to it, +# taking its NUMBER (needs the DR6 tile image); `sextractor` keeps +# SExtractor's own (required for image_sims) # - psf_model: allowed are `psfex` (default), `mccd`; `fake` (image_sims only, # PSF read from psf_dict) input_types: @@ -62,8 +63,7 @@ machines: inputs: tiles: $base_dir/unions-wl/tiles exposures: $base_dir/unions-wl/exposures - # UNIONS per-tile catalogues and segmentation maps (CFIS..r.cat, - # CFIS..r.seg.fits.fz), fetched by vcp in + # UNIONS per-tile catalogues (CFIS..r.cat), fetched by vcp in # the job: needs network from compute nodes and ~/.ssl/cadcproxy.pem. catalogues: vos:cfis/tiles_DR6 outputs: diff --git a/workflow/config/cfis/config_tile_Gic.ini b/workflow/config/cfis/config_tile_Gic.ini index 1e6bdb7f5..a9cb21654 100644 --- a/workflow/config/cfis/config_tile_Gic.ini +++ b/workflow/config/cfis/config_tile_Gic.ini @@ -1,5 +1,4 @@ # ShapePipe configuration file for: get the UNIONS per-tile object catalogue -# and its r-band segmentation map # (tile_detection: unions_catalogue in the workflow run config) @@ -54,7 +53,7 @@ TIMEOUT = 96:00:00 ## Module options -# Get the external tile catalogue and its segmentation map +# Get the UNIONS tile catalogue (CFIS..r.cat) [GET_IMAGES_RUNNER] FILE_PATTERN = tile_numbers @@ -68,19 +67,19 @@ NUMBERING_SCHEME = # Where the catalogues are: the run config's inputs.catalogues, a local # directory or a vos: URL (vos:cfis/tiles_DR6) -INPUT_PATH = $SP_INPUT_CATALOGUES, $SP_INPUT_CATALOGUES +INPUT_PATH = $SP_INPUT_CATALOGUES # Input file pattern including tile number as dummy template -INPUT_FILE_PATTERN = CFIS.000.000.r, CFIS.000.000.r.seg +INPUT_FILE_PATTERN = CFIS.000.000.r # Input file extensions -INPUT_FILE_EXT = .cat, .fits.fz +INPUT_FILE_EXT = .cat # Input numbering scheme, python regexp INPUT_NUMBERING = \d{3}\.\d{3} # Output file pattern without number -OUTPUT_FILE_PATTERN = CFIS_cat-, CFIS_seg- +OUTPUT_FILE_PATTERN = CFIS_cat- # Copy/download method, one in 'vos', 'symlink'; the workflow derives it from # inputs.catalogues (vos: URL -> vos, local directory -> symlink) diff --git a/workflow/config/cfis/config_tile_Sx.ini b/workflow/config/cfis/config_tile_Sx.ini index aa5db396f..43a59c493 100644 --- a/workflow/config/cfis/config_tile_Sx.ini +++ b/workflow/config/cfis/config_tile_Sx.ini @@ -115,6 +115,34 @@ CHECKIMAGE = BACKGROUND, SEGMENTATION # File name suffix for the output sextractor files (optional) SUFFIX = sexcat +## Join to an external catalogue (optional) + +# Path to an ASCII SExtractor catalogue of the same image (NUMBER, X_IMAGE, +# Y_IMAGE): each SExtractor row pairs with its mutual nearest neighbour in it +# and takes its NUMBER; rows without a partner leave the catalogue. Empty for +# no join. The workflow sets SP_MATCH_CATALOGUE to the tile's UNIONS +# catalogue under tile_detection: unions_catalogue and leaves it unset +# otherwise (image simulations, tile_detection: sextractor). +# @sc [decision:detection.tile_detection] +MATCH_CATALOGUE = ${SP_MATCH_CATALOGUE:-} + +# Largest separation of a pair, in pixels. The same detection run on the same +# pixels pairs at ~1e-4 px, so the radius only has to exclude neighbours. +# @sc [decision:detection.tile_detection] +MATCH_RADIUS = 1.0 + +# The run stops when the unpaired rows of either side (SExtractor rows, +# catalogue objects) exceed both (1 - MATCH_MIN_FRACTION) of that side and +# MATCH_TOLERATED_UNPAIRED. The guard catches a catalogue measured on other +# pixels (the DR5 image of a tile against its DR6 catalogue pairs ~95%), not +# deblending near very large objects: DR6 tiles pair at least 99.2% of each +# side, and a near-empty edge tile may leave a few children unpaired. +# @sc [decision:detection.tile_detection] +MATCH_MIN_FRACTION = 0.98 + +# @sc [decision:detection.tile_detection] +MATCH_TOLERATED_UNPAIRED = 20 + ## Post-processing # Necessary for tiles, to enable multi-exposure processing diff --git a/workflow/config/cfis/config_tile_Uc.ini b/workflow/config/cfis/config_tile_Uc.ini deleted file mode 100644 index 906ddef02..000000000 --- a/workflow/config/cfis/config_tile_Uc.ini +++ /dev/null @@ -1,97 +0,0 @@ -# ShapePipe configuration file for tile object selection from the UNIONS -# per-tile catalogue (tile_detection: unions_catalogue in the workflow run -# config). Writes the same run dir as config_tile_Sx.ini, so that the tile -# chain downstream reads one path whichever detection mode produced it. - - -## Default ShapePipe options -[DEFAULT] - -# verbose mode (optional), default: True, print messages on terminal -VERBOSE = True - -# Name of run (optional) default: shapepipe_run -RUN_NAME = run_sp_tile_Sx - -# Add date and time to RUN_NAME, optional, default: True -RUN_DATETIME = False - - -## ShapePipe execution options -[EXECUTION] - -# Module name, single string or comma-separated list of valid module runner names -MODULE = read_ext_sexcat_runner - -# Run mode, SMP or MPI -MODE = SMP - - -## ShapePipe file handling options -[FILE] - -# Log file master name, optional, default: shapepipe -LOG_NAME = log_sp - -# Runner log file name, optional, default: shapepipe_runs -RUN_LOG_NAME = log_run_sp - -# NUMBER_LIST selects this unit; the workflow sets SP_UNIT_NUM to the -# dashed tile ID (e.g. -210-282). -NUMBER_LIST = $SP_UNIT_NUM - -# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN -INPUT_DIR = $SP_RUN/output - -# Output directory -OUTPUT_DIR = $SP_RUN/output - - -## ShapePipe job handling options -[JOB] - -# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial -SMP_BATCH_SIZE = 16 - -# Timeout value (optional), default is None, i.e. no timeout limit applied -TIMEOUT = 96:00:00 - - -## Module options - -[READ_EXT_SEXCAT_RUNNER] - -INPUT_DIR = $SP_RUN/output/run_sp_tile_Gic/get_images_runner/output, $SP_RUN/output/run_sp_tile_Git/get_images_runner/output, $SP_RUN/output/run_sp_tile_Gic/get_images_runner/output, $SP_RUN/output/run_sp_tile_Mh_exp/merge_headers_runner/output - -FILE_PATTERN = CFIS_cat, CFIS_image, CFIS_seg, log_exp_headers - -FILE_EXT = .cat, .fits, .fitsfz, .sqlite - -# NUMBERING_SCHEME (optional) string with numbering pattern for input files -NUMBERING_SCHEME = -000-000 - -# File name suffix for the output sextractor files (optional) -SUFFIX = sexcat - -# Side length of the square postage stamp (vignet) extracted from the tile -# image, in pixels (must be odd). Default: 51 -# @sc [decision:postage_stamp_size] -VIGNET_SIZE = 51 - -# The third input is the catalogue's r-band segmentation map: set neighbours' -# VIGNET pixels to -1e30, as SExtractor does, so ngmix masks them -# @sc [decision:detection.catalogue_neighbour_marking] -SEGMENTATION = True - -## Post-processing - -# Necessary for tiles, to enable multi-exposure processing -# @sc [decision:detection.epoch_membership_ccd_bounds] -MAKE_POST_PROCESS = True - -# World coordinate keywords, SExtractor output. Format: KEY_X,KEY_Y -WORLD_POSITION = ALPHA_J2000,DELTA_J2000 - -# Accepted CCD pixel bounds xmin,xmax,ymin,ymax (0-based, strict) -# @sc [decision:detection.epoch_membership_ccd_bounds] -CCD_SIZE = 33,2080,1,4612 diff --git a/workflow/config/cfis/default_tile.sex b/workflow/config/cfis/default_tile.sex index ca7367c3e..cb6370e4c 100644 --- a/workflow/config/cfis/default_tile.sex +++ b/workflow/config/cfis/default_tile.sex @@ -121,7 +121,10 @@ BACK_FILTTHRESH 0.0 # Threshold above which the background- #--------------------- Memory (change with caution!) ------------------------- MEMORY_OBJSTACK 3000 # number of objects in stack -MEMORY_PIXSTACK 300000 # number of pixels in stack +# Large enough that no object overflows the stack, a bright star's halo +# included (300000 overflows on such tiles and truncates the object). +# @sc [decision:detection.detection_threshold_policy] +MEMORY_PIXSTACK 3000000 # number of pixels in stack MEMORY_BUFSIZE 1024 # number of lines in buffer #------------------------------- ASSOCiation --------------------------------- diff --git a/workflow/config/cfis/final_cat.param b/workflow/config/cfis/final_cat.param index 88a831b95..fb8b3f3f7 100644 --- a/workflow/config/cfis/final_cat.param +++ b/workflow/config/cfis/final_cat.param @@ -8,8 +8,9 @@ YWIN_WORLD # @sc [decision:catalogue_assembly.tile_overlap_handling] TILE_ID -# Object number within the tile: sp_validation's extraction joins on it, and -# the per-tile catalogue always carries it. +# Object number within the tile, the UNIONS catalogue's under tile_detection: +# unions_catalogue: sp_validation's extraction joins on it, and the per-tile +# catalogue always carries it. NUMBER # Survey-wide object ID, tile_id * 10**6 + NUMBER (make_cat, both @@ -226,10 +227,7 @@ NGMIX_FLUX_ERR_2P NGMIX_FLUX_ERR_NOSHEAR # magnitudes -# SExtractor detection columns. Under tile_detection: unions_catalogue the -# UNIONS catalogue supplies them, and the merge drops the ones it lacks -# (merge_final_cat.SEXTRACTOR_ONLY_COLUMNS: the WIN, AUTO-flux, APER and FWHM -# columns and SNR_WIN); MAG_AUTO, MAGERR_AUTO and FLUX_RADIUS are in both. +# SExtractor detection columns, from the tile SExtractor run (tile_detect). MAG_AUTO MAGERR_AUTO MAG_WIN diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk index f105d23b4..107bce4f9 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -37,22 +37,22 @@ Distinct tiles share no edge, so each is one group job per tile: job reaches the widest partition set (plus cpubackfill) — which is why each member's runtime is measured p99 plus margin and not the old ceiling. -The heavy middle (tile_detect) stays out: it is a 16 GB / 8 thread SExtractor -run that the shape chain does not need co-scheduled, and folding it in would add -its runtime to a sum that has no room. +tile_detect stays out: the shape chain does not need it co-scheduled, and +folding it in would add its runtime to a sum that has no room. -``tile_detection: unions_catalogue`` (config.yaml) replaces that SExtractor run -with two rules: tile_get_catalogue fetches the UNIONS per-tile catalogue and -its r-band segmentation map (get_images_runner, config_tile_Gic.ini) and -tile_detect converts it to the FITS-LDAC sexcat SExtractor would have written -(read_ext_sexcat_runner, config_tile_Uc.ini), with the tile image's header, -VIGNET stamps cut from the tile image with neighbours' footprints set to -1e30 -as SExtractor sets them (so ngmix masks neighbours in both modes), and the -multi-epoch post-processing. It keeps the catalogue's own -NUMBER, from which make_cat builds ``TILE_UNIQUE_ID`` exactly as in SExtractor -mode. The converter writes run_sp_tile_Sx/read_ext_sexcat_runner, and the rule links it as -sextractor_runner, the one path every downstream config reads; the manifest is -tile_detect.json in both modes, so tile_vignets onwards is the same DAG. +tile_detect runs SExtractor on the tile image (sextractor_runner, +config_tile_Sx.ini) for data and image simulations alike, so both get the same +detection, windowed positions, VIGNET neighbour marking and columns. +``tile_detection: unions_catalogue`` (config.yaml, the data default) adds one +rule and one step: tile_get_catalogue fetches the UNIONS per-tile catalogue +(get_images_runner, config_tile_Gic.ini), and tile_detect joins its SExtractor +rows to it before the multi-epoch post-processing (MATCH_CATALOGUE, set +through SP_MATCH_CATALOGUE). Each row pairs with its mutual nearest catalogue +object within 1 px and takes its NUMBER, from which make_cat builds +``TILE_UNIQUE_ID``, the key shared with the photometry and photo-z +catalogues; unpaired rows leave the catalogue. The catalogue is SExtractor run +by MegaPipe on the same DR6 image with the same configuration, so the join is +exact; on other pixels (a DR5 image) it fails loudly. There is no `tile_mask` rule, and there will not be one (PR #847). ShapePipe generates no masks: tiles have no instrument flag image of their own, so @@ -472,36 +472,13 @@ rule tile_merge_headers: shell: sp_shell("tile_merge_headers", "config_tile_Mh_exp.ini") -# The tile's galaxy sample. One of two definitions of tile_detect, chosen at -# parse time by the run config; both produce manifests/tile_detect.json and -# run_sp_tile_Sx/sextractor_runner/output/sexcat.fits. -if TILE_DETECTION == "sextractor": +# The UNIONS per-tile catalogue (CFIS..r.cat), from a local mirror or +# from vos, the way tile_get_images fetches the image; tile_detect joins its +# SExtractor detections to it. Reads only tile_numbers.txt, which unit_pre +# writes; the edge on the image manifest is what puts it after the prepare +# phase. +if TILE_DETECTION == "unions_catalogue": - # SExtractor object detection on the tile. - rule tile_detect: - input: - uz = f"{TILE_DIR}/manifests/tile_uncompress.json", - mh = rules.tile_merge_headers.output.manifest, - output: - manifest = f"{TILE_DIR}/manifests/tile_detect.json" - log: - f"{TILE_DIR}/logs/tile_detect.json" - params: - pre = lambda wc: unit_pre("tile_detect", wc.tile), - script_hash = SCRIPT_HASH - threads: 8 - resources: - mem_mb = lambda wc, attempt: 16000 * attempt, - runtime = 180 - shell: - sp_shell("tile_detect", "config_tile_Sx.ini") - -else: - - # Fetch the UNIONS per-tile catalogue (CFIS..r.cat) and segmentation - # map (CFIS..r.seg.fits.fz), from a local mirror or from vos, the way tile_get_images fetches the image. Reads - # only tile_numbers.txt, which unit_pre writes; the edge on the image - # manifest is what puts it after the prepare phase. rule tile_get_catalogue: input: git = f"{TILE_DIR}/manifests/tile_get_images.json", @@ -523,36 +500,63 @@ else: shell: sp_shell("tile_get_catalogue", "config_tile_Gic.ini") - # Convert the catalogue to the FITS-LDAC sexcat the chain expects. The - # completeness check counts read_ext_sexcat_runner's own output; the link - # after it is what makes that output readable at the sextractor_runner path - # of every downstream config, and it is made only on success so a failed - # run leaves nothing at that path. - rule tile_detect: - input: - cat = rules.tile_get_catalogue.output.manifest, - git = f"{TILE_DIR}/manifests/tile_get_images.json", - mh = rules.tile_merge_headers.output.manifest, - output: - manifest = f"{TILE_DIR}/manifests/tile_detect.json" - log: - f"{TILE_DIR}/logs/tile_detect.json" - params: - pre = lambda wc: unit_pre( - "tile_detect", wc.tile, - env={"SP_TILE_DETECTION": TILE_DETECTION}), - script_hash = SCRIPT_HASH - threads: 1 - resources: - # The whole tile image plus one 51x51 float32 stamp per object. - mem_mb = lambda wc, attempt: 8000 * attempt, - runtime = 60 - shell: - sp_shell("tile_detect", "config_tile_Uc.ini", - post="if [ $rc -eq 0 ]; then\n" - ' ln -s read_ext_sexcat_runner ' - '"$SP_RUN/output/run_sp_tile_Sx/sextractor_runner"\n' - "fi\n") + +# tile_detect's inputs: the catalogue manifest too when it joins one. +DETECT_INPUTS = {"uz": f"{TILE_DIR}/manifests/tile_uncompress.json", + "mh": rules.tile_merge_headers.output.manifest} +if TILE_DETECTION == "unions_catalogue": + DETECT_INPUTS["cat"] = rules.tile_get_catalogue.output.manifest + + +def detect_env(tile): + """tile_detect's prologue exports: the catalogue to join, if any. + + Empty under tile_detection: sextractor (image simulations), so that no + value left exported in the submitting shell reaches the job. + """ + if TILE_DETECTION != "unions_catalogue": + return {"SP_MATCH_CATALOGUE": ""} + gic = f"{tile_dir(tile)}/output/run_sp_tile_Gic/get_images_runner/output" + return {"SP_MATCH_CATALOGUE": f"{gic}/CFIS_cat{unit_num(tile)}.cat"} + + +# SExtractor object detection on the tile; under unions_catalogue joined to +# the UNIONS catalogue (config_tile_Sx.ini, MATCH_CATALOGUE). +rule tile_detect: + input: + **DETECT_INPUTS + output: + manifest = f"{TILE_DIR}/manifests/tile_detect.json" + log: + f"{TILE_DIR}/logs/tile_detect.json" + params: + pre = lambda wc: unit_pre("tile_detect", wc.tile, + env=detect_env(wc.tile)), + script_hash = SCRIPT_HASH + # One core: the container's SExtractor is built without threads + # (NTHREADS 4 warns and runs in the same 36 s), default_tile.sex sets + # NTHREADS 1, the join and the post-processing are serial Python, and + # `-b {threads}` only sets the SMP batch over input sets, of which a tile + # is one. + # + # MEASURED on eight DR6 tiles chosen to span conditions (high latitude, + # b = 18 deg, the A2199 cluster, Alioth's halo, two survey-edge tiles at + # 96% zero weight; 1-10 exposures), the step as shapepipe_run runs it on + # candide: wall time 23-160 s (SExtractor + join <= 93 s, post-processing + # <= 97 s), peak RSS 1.84 GiB, set by the join holding the ~430 MB + # catalogue twice. Image simulations run the same step on tiles of the same + # size (10000 x 10000 px); SExtractor alone on one takes 34 s and 0.74 GiB. + # 4000 MB is 2.1x the worst tile, and an OOM retries at 8000. The runtime + # is ~7x the slowest tile, for /scratch I/O on nibi (a 400 MB image and + # weight in, a ~430 MB catalogue out). nibi bills max(cores, mem_GB/4), so + # the job bills 1 core-equivalent. + threads: 1 + retries: 1 + resources: + mem_mb = lambda wc, attempt: 4000 * attempt, + runtime = lambda wc, attempt: 20 * attempt + shell: + sp_shell("tile_detect", "config_tile_Sx.ini") # Configured PSF interpolation to galaxies + vignet postage stamps: the last # stage that reads exposure products, and the bulk intra-tile intermediate. The store it @@ -1040,7 +1044,6 @@ rule final_cat_merge: tile_list = str(config["tile_list"]), index_db = str(INDEX_DB), param_file = str(CONFIG_DIR / "final_cat.param"), - tile_detection = TILE_DETECTION, campaign = CAMPAIGN, snapshot = str(SNAPSHOT_JSON), inputs = unit_fingerprint(TILES_READY), @@ -1065,5 +1068,4 @@ rule final_cat_merge: " --output {output.merged}" " --campaign '{params.campaign}'" " --param-file '{params.param_file}'" - " --tile-detection {params.tile_detection}" " --snapshot-json '{params.snapshot}'" diff --git a/workflow/scripts/completeness.py b/workflow/scripts/completeness.py index 52720b724..b39a951ed 100644 --- a/workflow/scripts/completeness.py +++ b/workflow/scripts/completeness.py @@ -72,15 +72,13 @@ import sys from pathlib import Path -# How the tile's galaxy sample is produced: SExtractor on the tile image, or -# the UNIONS per-tile catalogue converted in place. The run config's -# `tile_detection:` picks one; the rules export it as $SP_TILE_DETECTION only -# when it is not the default, so a SExtractor run's prologue is unchanged. +# How the tile's galaxy sample is defined: SExtractor on the tile image alone, +# or those detections joined to the UNIONS per-tile catalogue, whose NUMBER +# they take. The run config's `tile_detection:` picks one. TILE_DETECTIONS = ("sextractor", "unions_catalogue") # stage -> {runner_subdir: {expect, [warn], [subpath]}} -# exp_psf and tile_vignets are selected by $SP_PSF at check time, tile_detect -# by $SP_TILE_DETECTION. +# exp_psf and tile_vignets are selected by $SP_PSF at check time. # @sc [decision:per_unit_completeness] COMPLETENESS = { # --- tile prepare (phase A) --- @@ -142,13 +140,9 @@ # --- tile post --- "tile_merge_headers": {"merge_headers_runner": dict(expect=1)}, - # The fetched UNIONS catalogue and its r-band segmentation map. - "tile_get_catalogue": {"get_images_runner": dict(expect=2)}, - "tile_detect": { - "sextractor": {"sextractor_runner": dict(expect=2)}, - # The FITS-LDAC sexcat converted from the fetched catalogue. - "unions_catalogue": {"read_ext_sexcat_runner": dict(expect=1)}, - }, + # The fetched UNIONS catalogue. + "tile_get_catalogue": {"get_images_runner": dict(expect=1)}, + "tile_detect": {"sextractor_runner": dict(expect=2)}, "tile_vignets": { "psfex": { "psfex_interp_runner": dict(expect=1), @@ -223,13 +217,6 @@ def check_counts(stage, run_dir): raise ValueError( f"Invalid SP_PSF={psf_model!r}; expected one of {sorted(table)}." ) from exc - elif stage == "tile_detect": - detection = os.environ.get("SP_TILE_DETECTION", TILE_DETECTIONS[0]) - if detection not in TILE_DETECTIONS: - raise ValueError( - f"Invalid SP_TILE_DETECTION={detection!r}; expected one of " - f"{', '.join(TILE_DETECTIONS)}.") - table = table[detection] details, ok = [], True for runner, spec in table.items(): n = count_products(run_dir, runner, spec) @@ -262,8 +249,6 @@ def check_counts(stage, run_dir): "exp_psf": ("exp", "run_sp_exp_SxSePsf"), "tile_merge_headers": ("tile", "run_sp_tile_Mh_exp"), "tile_get_catalogue": ("tile", "run_sp_tile_Gic"), - # Both detection modes write here (config_tile_Sx.ini / config_tile_Uc.ini - # share the RUN_NAME): the chain downstream reads one path. "tile_detect": ("tile", "run_sp_tile_Sx"), "tile_vignets": ("tile", "run_sp_tile_PiViVi"), "tile_ngmix": ("tile", "run_sp_tile_ngmix_Ng${SP_NGMIX_CHUNK}u"), diff --git a/workflow/scripts/merge_final_cat.py b/workflow/scripts/merge_final_cat.py index afa05373a..0d4594ec4 100644 --- a/workflow/scripts/merge_final_cat.py +++ b/workflow/scripts/merge_final_cat.py @@ -52,18 +52,6 @@ argument validator and then falls through to the ordinary walk, so it is not a way to add one tile by hand.) -WHICH COLUMNS DEPEND ON WHERE THE TILE'S DETECTIONS CAME FROM. The param file -lists the columns of a SExtractor-mode tile catalogue. Under ``tile_detection: -unions_catalogue`` the detection columns are instead the UNIONS catalogue's own -(read_ext_sexcat copies them as they are), and that catalogue lacks some of the -SExtractor quantities the param file names; ``SEXTRACTOR_ONLY_COLUMNS`` lists -them, and ``requested_columns`` subtracts exactly that list in that mode. The -merge stays strict on everything else: a column still requested and absent -from a tile stops the merge. tests/module/test_final_cat_columns.py derives the -catalogue-mode columns by running the converter on the UNIONS catalogue's -header and asserts the list is exactly the requested SExtractor columns the -catalogue does not carry, so it cannot silently go stale in either direction. - WHICH TILES — AND WHY THE JOB DERIVES THE SET RATHER THAN BEING TOLD IT. The set is the CAMPAIGN's: every tile both declared in ``tile_list`` and present in the index, which is exactly the Snakefile's TILES_READY, rebuilt here from the same @@ -99,41 +87,12 @@ # Same directory; the rule invokes this file by path, so it is sys.path[0]. import build_index import hdf5_reconcile -from completeness import TILE_DETECTIONS # /scripts/python/create_final_cat.py, from /workflow/scripts/this. CFC_PATH = (Path(__file__).resolve().parents[2] / "scripts" / "python" / "create_final_cat.py") -# Requested SExtractor columns the UNIONS per-tile catalogue (DR6) does not -# carry, so a unions_catalogue-mode tile catalogue cannot have them. None has a -# DR6 column measuring the same quantity: DR6's MAG_APER is a magnitude through -# its own aperture, not FLUX_APER; MAG_AUTO, MAGERR_AUTO and FLUX_RADIUS are -# the same SExtractor measurements in both modes and stay requested. -SEXTRACTOR_ONLY_COLUMNS = ( - "MAG_WIN", "MAGERR_WIN", # Gaussian-windowed magnitude - "FLUX_AUTO", "FLUXERR_AUTO", # DR6 carries the AUTO magnitude only - "FLUX_APER", "FLUXERR_APER", # ShapePipe's aperture, in flux - "SNR_WIN", # Gaussian-windowed SNR - "FWHM_IMAGE", "FWHM_WORLD", # Gaussian-core FWHM -) - - -def requested_columns(param_list: list, tile_detection: str) -> list: - """The columns a tile catalogue of this detection mode must carry. - - The param file's list, in its order, less ``SEXTRACTOR_ONLY_COLUMNS`` when - the detections are the UNIONS catalogue's (see the module docstring). - """ - if tile_detection not in TILE_DETECTIONS: - raise ValueError(f"tile_detection {tile_detection!r} is not one of " - f"{TILE_DETECTIONS}") - if tile_detection == "sextractor": - return list(param_list) - return [c for c in param_list if c not in SEXTRACTOR_ONLY_COLUMNS] - - def spval_group(campaign: str) -> str: """The hdf5 group the campaign's per-tile datasets live under. @@ -192,9 +151,6 @@ def main() -> None: help="names the campaign's group in the output file") p.add_argument("--param-file", required=True, type=Path, help="the input type's final_cat.param — the column list") - p.add_argument("--tile-detection", required=True, choices=TILE_DETECTIONS, - help="where the tile detections came from (config " - "tile_detection); selects the columns requested") p.add_argument("--hdu", type=int, default=1) p.add_argument("--snapshot-json", type=Path, default=None, help="sp run's code snapshot (bin/sp's " @@ -202,9 +158,7 @@ def main() -> None: args = p.parse_args() cfc = load_create_final_cat() - param_list = requested_columns( - cfc.read_param_file(str(args.param_file), verbose=False), - args.tile_detection) + param_list = cfc.read_param_file(str(args.param_file), verbose=False) if not param_list: sys.exit(f"merge_final_cat: no columns read from {args.param_file}") # read_data/copy_data read their knobs out of this dict, exactly as diff --git a/workflow/scripts/run_report.py b/workflow/scripts/run_report.py index a8efe82ae..bd846d34d 100644 --- a/workflow/scripts/run_report.py +++ b/workflow/scripts/run_report.py @@ -60,7 +60,7 @@ TILE_STAGES = ["tile_get_images", "tile_uncompress", "tile_find_exposures", "tile_merge_headers", "tile_detect", "tile_vignets", "tile_ngmix", "tile_merge_cats", "tile_make_cat"] -# The catalogue fetch exists only when the run converts the UNIONS catalogue +# The catalogue fetch exists only when the run joins the UNIONS catalogue # (config.yaml's tile_detection); the callers pass the mode in the environment # so a SExtractor run does not report the stage as not run. if os.environ.get("SP_TILE_DETECTION") == "unions_catalogue":