From 668e05d539ba8556a25ccf03d0d557ee729d5f5d Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Fri, 9 Jan 2026 14:18:07 +0000 Subject: [PATCH 01/34] removed tmp files --- docs/source/.pipeline_canfar.md.swo | Bin 16384 -> 0 bytes docs/source/.pipeline_canfar.md.swp | Bin 16384 -> 0 bytes 2 files changed, 0 insertions(+), 0 deletions(-) delete mode 100644 docs/source/.pipeline_canfar.md.swo delete mode 100644 docs/source/.pipeline_canfar.md.swp diff --git a/docs/source/.pipeline_canfar.md.swo b/docs/source/.pipeline_canfar.md.swo deleted file mode 100644 index efa536569de8b57939d1fb7c21cf9e46f6f94b8e..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 16384 zcmeHOTWsW389ww1z0d;U1qtbIlCs&=dYnnJx1`yn&1D zfh=tqCvO!K0=^482Ye0qEbtlNEN~n+27Da&5O6#2 z_6LRdE${~LD)1xV72rAGLEr)~4VXX;xC0VhHu-bl zhrn~dv%q7(zit-dC%}uq3&7Wbi@-9l1Q@_k;9B67n}ql_@OdBszVHDdrhvb^Ux+^g zF9Vl=ZQv1L5|{u!0o)1P0K9gi5WfUo0lo!%6Zj%g1}C;&A?L%Z@BPWYhv6|vYRHhlSYZ2SL^{Br}=IisVtG?%2Yb}arRvvs~o)g zFfG|?gEy!JDk1ijgHvGFgGKR=Ub(ui7K ziZ!%YJ&Q(F#MhH7%w?66AZv8CBb^N@tVLU1+9^$F@%WXIm=~tWu~Wx+?FNpJ?o#nU zlb3$h?=a$W-}_0*NxV2!X1twtqR=cnhnKCz+r5<^nizkiI5azyFO7TW)~NWiv^6cx zDV!T&G1^L^pm!d)EDn2^;AlKFf#Qltn5jePxtRTp>QkiRJV0jZM@h#guJB3i(1p%& zuYx~C0bbjZyhT-O+8n|TwS2p6K|18k=HGLOv9{rXSlxo`+!2c@Wo(CgVwPk4T6_}A+|c48qy?_)8s~Upq4S;+Kd!yd76AhhEJ}Y z+8p5~Jw{v3GzGYScob!6oTVmq%n{Dhj+nhODmEYvX%?f7)R2|ck=CW2uxT0~Jqch5 z+#2?ZpN4+1mpuO-#uYa(PA*#qM0$%-ZnABKhi%C`j6UtMpFZ zb9>Xyy4Z!u&SL6=SmR@O8})ecA=j8Ly!qSyP6JEvCipVt4G2!FQX%p?^vXkZ&lj$`6;b57Wj@K{jy{yJ43mmcLH5-i8`%|Al!9OuMHLRhK}kykQZ0^6om>Fu?7QnQ+m(q65?|hAwA^^~9;dbVvC^1pU;0}tk)6Q^rLoG8| z>#~e}-mmBW|8wY7pF(fT{eSs<{}1Trp9j7Nd={(*k|cfe)fabN)$0=S>Q z3HT@a^*4c6fn{I`xC{6g@D_UZmw}gnr-A!{Szr>F0Jxw374R(Z4A2Bl0~5e4z|Fu< z(YyZuh=B+Qfpx$FR)7=0^}vhh?Oy;Ca0fsD_xisE_`FZ@nX8{l1C<6U4OAMaG*D^a zzoCIPI3zL_3c{0VKi!ssUeU$ieq3#R{G=7O1Fu^&u-1egdYZd${iom>UNNt$N?J>S zCn{*kob}CV$-H$(W6op&$eaW<7RW4aj5Wc-R76PR%G>Py@5LSNJSJUG@(e?>Hs?W5fcykskk_90n`v_2mp+Xe*o%}d_>zd4IV}My zp?z2gF6h3SI+RI+1v`+B3=j%s&54zUmktOBxTR$VK<0| zoPOnP7LrP9a$Gzx+t|mKMnDkEo?Bj6JX?Ctif{2FD(yFWTyO7iI9qS{UK)MpftU6M zY^|}&!QesKxj=ZqxmY*aL;;Rke1>#>P*Ll+tz+7Dxk$dISqif^vc~xU#@uDus4%k` zYMNE$3+dDctCz-_Sx=K zc0x#RgL_DjPJ8vhj}vsu4DtoCWEXS$dFY+lye{SkXoqu_hW)wc?%ep|S%;$!(95lh z^yuS7hAkgn96C73_JZ^x9joKPnbZ;tFrGn7^+VWQk#EZ!JCL;3)892~--{!r5jY|S u${#;=)u#e@vQiiV&ov>e=;8;<#J%upPa0U}Dw8KV1SpU;6N=0xh<^hEQO@`P diff --git a/docs/source/.pipeline_canfar.md.swp b/docs/source/.pipeline_canfar.md.swp deleted file mode 100644 index a972bd87d7288b67e06096419218029bbabcf38c..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 16384 zcmeHOTWlmp6|H0;47}0#9p(r&VwYhd%A08I_{Uc zy6y4qW)Tt~5dr~%#7hYh5<#FK1pFnNeDaY*LIff47eOQ<0{kGs2NDYBR(Frbk8y&O zk07Q$YvO>Ya0jLcd2L4P#r6hksC?F32kMlFpX_#@CDBeq-k40%ZjVMbAYw-tCE3| zfj7-SnzSk>uM^{A^%}Ps9=ee}@cz{|-MqZ9WT0fAWT0fAWT0fAWT0fAWT0f=|DS;* zx>9@&Zd{qWav}eH#jf9n^3Rj`?~A*>&p(#GN(M>>N(M>>N(M>>N(M>>N(M>>N(M>> zN(M>>-i8cVc*OoOjBAJm06hOM>i@r9FT`uWkAWwFZv$Thz5px$CxGL?M}hYNHvq4{ zSBO6We*|6to&%l*o&@d%?g6F%1E>Nw0t5))_3MQA2k;E=4Zs9$2CfFK0{-$IA)W`G z2A%|d1Y89E{ca(C4g3sv3ivMYAg~N90TtjFa0qzzT|)d2co>L*uUsp{B=GlZg!mh9 z88{DY0S^EZz&P*;;3nWI;H9gD_$}}(@C5KU@Cfi}UoB~b)zrPCp08ay71ug*R z0T<{1i@=9~4+2*Jue=ld0lxsA03HYK2W|rX^$sEa3A_UQ3TOi#2VR3P^9SHD;G4kb zflmQ<0LOs|;1KX4=Jz6S0q6j?08;?Z`;EZhpJ6)ZdXn@T8{|t%rLpv7kO+~zrD3YH z>AQ*fVA!O}gLGnoDlJOT#0*luDPvVlx{1iGM6ryl*h<_mz%yt&jlKL`57P^6%d4sm z+ZZ3CiccecC%ei9Eze^+REcOLQZ1@1)5w{n6?1uY$vk)a*}21m63=$ zrit3A6Wn(LE+jiB9ysMiob?BcxYYN)lConrN|X_8C7m!Z^2p(CtI<}k=ljOT9>|^M zfbwMJ_Ssd+|19R37H1Xw#^M@n#i8Hx2O$e*j}QWlZ;e4*z7nSEFnT89c%%9xsVG|@ z)AeJd*AK3+r&i#Avn;9*Poa<3xFiozRht$sF`HVR)iz-r_GYs0d5JOGh(Ppi!gl@^ zR@S6gcfp=@$1C=v7#^mxpdR5OtwcJA5Vu8vv=NV3S-fU;$N6sGYPe?}(c1H>Ld_5} zTR2anwPRF2LKc=?xXu>Og-*yUu+OXKwcCwose^BFj?kd(g-y%rBL>}Sa#SCGH5zp! zg4cF=>{D5U^U9G{?74{Mje509Sx_-WM|R|ZC=S_8%rmm+zFDAtg4gcuq{vJ9zKWuY zoQdw>c+r`*Kcx(uI=01}zZH`=@n$EFJ4JZ(Uf$D-4Z&o$TxHl{&~n?Rc5?N6)eT&j z9j4U)@ds*^(*Dz?blaWe;2>MGtz-{dYaC~FsvtnNe7q&NuFcON<0O*!v!YQB{8&~n zMKmvu5nfE(4naLFioN~ zF)(6Y(cJ8)(eqKUX*g-*cO9)F3#+58^F3ozG{Ac7V@`0HeleCT-ylD=O-^x0 zj_g%{pv`jx^FdYP}Zor%#{OII3ePiy&=yM}YQq^@fZT5iUf&RWw%2v#ZO~Ytg2&Cj52B zjrly7#h~y4LN4v!s%6*mBTpKfo(9aF9u79zqOlz+KU9g87#m}TUQyYA^8-b(#3MAW+SMIAVAI{0xcLVU1xCypYf#PC%2>-|~f7?K}j%-<}m(WUMo_7rC6`MR+ zDUK(4B;bONq*MCPgIx#bv!T4*&@r6{vU5S!MB&7W<7s}wE^}eq<>PyPKUl>_e4Y*M zl*=w%JL!4OcD}|+z5ocv5GS=094T>d!48VF(@wFw;p{R|?Xr!1-mmBS|7FywFQc~Q z`oDO-|2x$4F9DZLg?%Ab;fl7W(el7W(el7Y7i18qo2RGH6pPgehQ zwF;F*pAq+yZ0lnu&7keO-Mp?f$KMhoy1!ewvLEL^1=k0PHDXmrWD1&7E^%fGawc(R zvH*&H79t>a$2h-17~`T_i$yGR6oTVL>RNFdmHY~^S_FNMpocNW{$|PD^Q1?k6^unn2U<)>%}nxuOwm3f1Roj~ zXSt~(PT5?7+1Zi2#q~smZ@>g=No8LVwl{M&PoW-zu$@9ehIgYpC#;h237(Wp zAYApd_e=rn>3>kzHry~(2(N4uj#ic-LaAIFW0gD?*ikbuLnxB)04{flnJ^=_JGL3d zp2J!;)T3&Z<1s6kPtV?~FAe!4DiF(aQ947l*1d3#m0NV}2BWuDZgHLN ztp!gGIyIO(tSvk?LeY|QV}gHUg$2;FI?ZAlA<1TF&eCfOv>e9 z=9Y@t!>kvRHNBR^`Hi$fhYzn|pjj$GXc@SxV$|1U+~^;&Up*%ebn8V9K^xhi+!%J6MovF-RWi8 zMY`}%o^OlCB8T57Ci0USL)qnP-FJCYV0}PBE-}muB Wk4Zd$MdQ5OA7BGp<9?4I#D4*RG06b{ From 38c58195cf95f9ca2e2630a2d8d1a8d6481695be Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Mon, 12 Jan 2026 14:34:24 +0000 Subject: [PATCH 02/34] Moved get_ccds script to pipeline runner and library script --- pyproject.toml | 1 + scripts/python/get_ccds_with_psf.py | 180 ++----------- src/shapepipe/get_ccds_run.py | 44 ++++ src/shapepipe/utilities/ccd_psf_handler.py | 292 +++++++++++++++++++++ 4 files changed, 356 insertions(+), 161 deletions(-) mode change 100644 => 100755 scripts/python/get_ccds_with_psf.py create mode 100644 src/shapepipe/get_ccds_run.py create mode 100644 src/shapepipe/utilities/ccd_psf_handler.py diff --git a/pyproject.toml b/pyproject.toml index 4d7ebe3b7..9a4153dbe 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -51,6 +51,7 @@ summary_run = "shapepipe.summary_run:main" canfar_submit_job = "shapepipe.canfar_run:run_job" canfar_monitor = "shapepipe.canfar_run:run_log" canfar_monitor_log = "shapepipe.canfar_run:run_monitor_log" +get_ccd_with_psfs = "shapepipe.get_ccds_run:run_ccd_psf_handler" [tool.pytest.ini_options] addopts = "--verbose --cov=shapepipe" diff --git a/scripts/python/get_ccds_with_psf.py b/scripts/python/get_ccds_with_psf.py old mode 100644 new mode 100755 index 5f917ad02..e284e2640 --- a/scripts/python/get_ccds_with_psf.py +++ b/scripts/python/get_ccds_with_psf.py @@ -9,192 +9,50 @@ """ - import sys -import numpy as np - -from shapepipe.utilities import summary - - -def get_lines(fname): - """Get Lines. - - Return list of lines read from a text file. - - Parameters - ---------- - fname: str - input file name - Returns: - list - IDs +from shapepipe.utilities.ccd_psf_handler import CcdPsfHandler - """ - IDs = [] - with open(fname) as f: - lines = f.readlines() - for line in lines: - IDs.append(line.rstrip()) - - return IDs +def run_ccd_psf_handler(args=None): + """Run CCD PSF Handler. -def get_exp_shdu_missing(patches): - """Get Exp Shdu Missing. - - Returns set of missing CCDs (single-exposure single-HDU IDs) from a list of patches. + Create instance and run the CCD PSF handler. Parameters ---------- - patches: list - input patches + args : list, optional + command line arguments Returns ------- - set - missing CCD IDs + int + exit code """ - exp_shdu_missing_all = set() - - for patch in patches: - - path_exp_shdu_missing = f"{patch}/summary/missing_job_32_all.txt" - exp_shdu_missing = get_lines(path_exp_shdu_missing) - - print(f"Patch {patch}: Found {len(exp_shdu_missing)} missing ccds", end="; ") - - exp_shdu_missing_all.update(exp_shdu_missing) - - print(f"cumulative {len(exp_shdu_missing_all)} missing ccds") - - print() - - return exp_shdu_missing_all - - -def get_exp(patches): - """Get Exp. - - Return set of exposures from a list of patches. - - Parameters - ---------- - patches: list - input patches - - Returns - ------- - set - exposure IDs - - """ - exp_all = set() - - for patch in patches: - - path_exp = f"{patch}/exp_numbers.txt" - exp = get_lines(path_exp) - - print(f"Patch {patch}: Found {len(exp)} exposures", end="; ") - - exp_all.update(exp) - - print(f"cumulative {len(exp_all)} exposures") + # Create instance + obj = CcdPsfHandler() - print() + return obj.run(args=args) - return exp_all +def main(argv=None): + """Main. -def get_ccds_with_psf(patches, n_CCD=40): - """Get CCDs With PSF. - - Return set of CCDs from list of patches. + Main program. Parameters ---------- - patches: list - input patches + argv : list, optional + command line arguments Returns ------- - set - CCD IDs + int + exit code """ - # Get missing CCDs - print("=== get missing CCDs ===") - exp_shdu_missing_all = get_exp_shdu_missing(patches) - - # Get all exposures used in tiles - print("=== get exposures ===") - exp_all = get_exp(patches) - - # Turn exposures into exposure-single-HDU names (CCDs) - exp_shdu_all = summary.get_all_shdus(exp_all, n_CCD) - - print(f"Found {len(exp_shdu_all)} CCDs") - - return exp_shdu_all - - -def get_ccds_with_psf_method_2(patches, n_CCD=40): - - for patch in patches: - directory = f"{patch}/exp_runs" - - for entry in os.scandir(directory): - pass - -def save(IDs, path): - """Save. - - Save list of IDs to text file. - - Parameters - ---------- - IDs: set - input IDs - - path: str - output file name - - """ - with open(path, "w") as f_out: - for ID in IDs: - print(ID, file=f_out) - -def main(argv): - """Main. - - Main program. - - """ - version = "v1.5" - - if version == "v1.4": - n_patch = 7 - elif version == "v1.5": - n_patch = 8 - else: - raise ValueError(f"Invalid version {version}") - - patches = [f'P{x}' for x in np.arange(n_patch) + 1] - - print(f"=== get_ccds_with_psf for version {version}, patches {patches} ===") - - print("=== method 1: exp_list - missing === ") - exp_shdu_all = get_ccds_with_psf(patches) - - save(exp_shdu_all, f"ccds_with_psf_{version}.txt") - - #print("=== method 2: star cats === ") - #exp_shdu_all_method_2 = get_ccds_with_psf_method_2(patches) - #save(exp_shdu_all_method_2, f"ccds_with_psf_{version}_method_2.txt") - - return 0 + return run_ccd_psf_handler(args=argv) if __name__ == "__main__": diff --git a/src/shapepipe/get_ccds_run.py b/src/shapepipe/get_ccds_run.py new file mode 100644 index 000000000..fb18b0a74 --- /dev/null +++ b/src/shapepipe/get_ccds_run.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python3 + +"""GET_CCDS_WITH_PSF + +Obtain list of CCDs (single-exposure single-HDU files) for which valid PSF information +is available. This can serve to create a footprint coverage mask. + +Author: Martin Kilbinger + +""" + +import sys + +from shapepipe.utilities.ccd_psf_handler import CcdPsfHandler + + +def run_ccd_psf_handler(args=None): + """Run CCD PSF Handler. + + Create instance and run the CCD PSF handler. + + Parameters + ---------- + args : list, optional + command line arguments + + Returns + ------- + int + exit code + + """ + # Create instance + + obj = CcdPsfHandler() + + return obj.run(args=args) + + +def main(argv=None): + """Main. + """ + # A scripts to call the ccd psf class is created by pyproject.toml + return 0 diff --git a/src/shapepipe/utilities/ccd_psf_handler.py b/src/shapepipe/utilities/ccd_psf_handler.py new file mode 100644 index 000000000..2707cc62e --- /dev/null +++ b/src/shapepipe/utilities/ccd_psf_handler.py @@ -0,0 +1,292 @@ +"""CCD_PSF_HANDLER + +Obtain list of CCDs (single-exposure single-HDU files) for which valid PSF +information is available. This can serve to create a footprint coverage mask. + +Author: Martin Kilbinger + +""" + +import sys +import numpy as np + +from cs_util import args as cs_args +from cs_util import logging + +from shapepipe.utilities import summary + + +class CcdPsfHandler(object): + """CCD PSF Handler Class. + + Handles extraction of CCDs with valid PSF information from shapepipe + patches. + """ + + def __init__(self): + """Initialize the handler.""" + self.params_default() + + def params_default(self): + """Set default parameters and command line options.""" + + self._params = { + "version_cat": "v1.6", + "n_CCD": 40, + "output": None, + } + + self._short_options = { + "version_cat": "-V", + "n_CCD": "-n", + "output": "-o", + } + + self._types = { + "n_CCD": "int", + } + + self._help_strings = { + "version_cat": "catalogue major version, allowed are v1.3, v1.4, v1.5, v1.6; default is {}", + "n_CCD": "number of CCDs per exposure; default is {}", + "output": "output file path; default is ccds_with_psf_.txt", + } + + def set_params_from_command_line(self, args): + """Set Params From Command line. + + Only use when calling using python from command line. + Does not work from ipython or jupyter. + + Parameters + ---------- + args : list + command line arguments + + """ + # Read command line options + options = cs_args.parse_options( + self._params, + self._short_options, + self._types, + self._help_strings, + args=args, + ) + self._params = options + + # Save calling command + logging.log_command(args) + + def update_params(self): + """Update parameters. + + Set derived parameters based on input parameters. + """ + # Determine number of patches based on version + version = self._params["version_cat"] + if version in ("v1.3", "v1.4"): + n_patch = 7 + elif version == "v1.5": + n_patch = 8 + elif version == "v1.6": + n_patch = 9 + else: + raise ValueError(f"Invalid version {version}") + + self._params["n_patch"] = n_patch + self._params["patches"] = [f"P{x}" for x in np.arange(n_patch) + 1] + + # Set output file if not specified + if self._params["output"] is None: + self._params["output"] = f"ccds_with_psf_{version}.txt" + + def check_params(self): + """Check parameters for validity.""" + # Add any parameter validation logic here + pass + + def get_lines(self, fname): + """Get Lines. + + Return list of lines read from a text file. + + Parameters + ---------- + fname : str + input file name + + Returns + ------- + list + IDs + + """ + IDs = [] + with open(fname) as f: + lines = f.readlines() + for line in lines: + IDs.append(line.rstrip()) + + return IDs + + def get_exp_shdu_missing(self, patches): + """Get Exp Shdu Missing. + + Returns set of missing CCDs (single-exposure single-HDU IDs) from a list of patches. + + Parameters + ---------- + patches : list + input patches + + Returns + ------- + set + missing CCD IDs + + """ + exp_shdu_missing_all = set() + + for patch in patches: + + path_exp_shdu_missing = f"{patch}/summary/missing_job_32_all.txt" + exp_shdu_missing = self.get_lines(path_exp_shdu_missing) + + print( + f"Patch {patch}: Found {len(exp_shdu_missing)} missing ccds", + end="; ", + ) + + exp_shdu_missing_all.update(exp_shdu_missing) + + print(f"cumulative {len(exp_shdu_missing_all)} missing ccds") + + print() + + return exp_shdu_missing_all + + def get_exp(self, patches): + """Get Exp. + + Return set of exposures from a list of patches. + + Parameters + ---------- + patches : list + input patches + + Returns + ------- + set + exposure IDs + + """ + exp_all = set() + + for patch in patches: + + path_exp = f"{patch}/exp_numbers.txt" + exp = self.get_lines(path_exp) + + print(f"Patch {patch}: Found {len(exp)} exposures", end="; ") + + exp_all.update(exp) + + print(f"cumulative {len(exp_all)} exposures") + + print() + + return exp_all + + def get_ccds_with_psf(self, patches, n_CCD=40): + """Get CCDs With PSF. + + Return set of CCDs from list of patches. + + Parameters + ---------- + patches : list + input patches + n_CCD : int + number of CCDs per exposure + + Returns + ------- + set + CCD IDs + + """ + # Get missing CCDs + print("=== get missing CCDs ===") + exp_shdu_missing_all = self.get_exp_shdu_missing(patches) + + # Get all exposures used in tiles + print("=== get exposures ===") + exp_all = self.get_exp(patches) + + # Turn exposures into exposure-single-HDU names (CCDs) + exp_shdu_all = summary.get_all_shdus(exp_all, n_CCD) + + print(f"Found {len(exp_shdu_all)} CCDs") + + return exp_shdu_all + + def save(self, IDs, path): + """Save. + + Save list of IDs to text file. + + Parameters + ---------- + IDs : set + input IDs + path : str + output file name + + """ + with open(path, "w") as f_out: + for ID in IDs: + print(ID, file=f_out) + + def run(self, args=None): + """Run. + + Main execution method. + + Parameters + ---------- + args : list, optional + command line arguments + + Returns + ------- + int + exit code (0 for success) + + """ + if args is None: + args = sys.argv + + # Set parameters from command line + self.set_params_from_command_line(args) + self.update_params() + self.check_params() + + # Get parameters + patches = self._params["patches"] + version = self._params["version_cat"] + n_CCD = self._params["n_CCD"] + output = self._params["output"] + + print( + f"=== get_ccds_with_psf for version {version}, patches {patches} ===" + ) + + print("=== method 1: exp_list - missing === ") + exp_shdu_all = self.get_ccds_with_psf(patches, n_CCD) + + self.save(exp_shdu_all, output) + + print(f"Results saved to {output}") + + return 0 From d6144ad6e4af78fc4f75ac57e922ed5815aee91b Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Mon, 12 Jan 2026 15:12:04 +0000 Subject: [PATCH 03/34] class to download eposure headers (for coverage) --- src/shapepipe/utilities/header_downloader.py | 247 +++++++++++++++++++ 1 file changed, 247 insertions(+) create mode 100644 src/shapepipe/utilities/header_downloader.py diff --git a/src/shapepipe/utilities/header_downloader.py b/src/shapepipe/utilities/header_downloader.py new file mode 100644 index 000000000..742acf380 --- /dev/null +++ b/src/shapepipe/utilities/header_downloader.py @@ -0,0 +1,247 @@ +"""HEADER_DOWNLOADER + +Download FITS headers from VOSpace for exposures listed in a CCD file. + +Author: Martin Kilbinger + +""" + +import sys +import os +from os.path import exists + +import numpy as np +import vos + +from cs_util import args as cs_args +from cs_util import logging + + +class HeaderDownloader(object): + """Header Downloader Class. + + Downloads FITS headers from VOSpace for exposures in a CCD list. + """ + + def __init__(self): + """Initialize the downloader.""" + self.params_default() + + def params_default(self): + """Set default parameters and command line options.""" + + self._params = { + "input_file": None, + "output_dir": "header", + "vospace_path": "vos:cfis/pitcairn", + "overwrite": False, + } + + self._short_options = { + "input_file": "-i", + "output_dir": "-o", + "vospace_path": "-p", + "overwrite": "-O", + } + + self._types = { + "overwrite": "bool", + } + + self._help_strings = { + "input_file": "input CCD list file (txt or csv); required", + "output_dir": "output directory for headers; default is {}", + "vospace_path": "VOSpace base path; default is {}", + "overwrite": "overwrite existing header files; default is {}", + } + + def set_params_from_command_line(self, args): + """Set Params From Command line. + + Only use when calling using python from command line. + Does not work from ipython or jupyter. + + Parameters + ---------- + args : list + command line arguments + + """ + # Read command line options + options = cs_args.parse_options( + self._params, + self._short_options, + self._types, + self._help_strings, + args=args, + ) + self._params = options + + # Save calling command + logging.log_command(args) + + def update_params(self): + """Update parameters. + + Set derived parameters based on input parameters. + """ + # Ensure output directory ends without trailing slash + if self._params["output_dir"].endswith("/"): + self._params["output_dir"] = self._params["output_dir"][:-1] + + def check_params(self): + """Check parameters for validity.""" + if self._params["input_file"] is None: + raise ValueError("Input file is required (use -i or --input_file)") + + if not exists(self._params["input_file"]): + raise FileNotFoundError( + f"Input file not found: {self._params['input_file']}" + ) + + # Create output directory if it doesn't exist + if not exists(self._params["output_dir"]): + os.makedirs(self._params["output_dir"]) + if self._params["verbose"]: + print(f"Created output directory: {self._params['output_dir']}") + + def get_exposures(self, ccd_list_file): + """Get Exposures. + + Extract unique exposure numbers from CCD list file. + + Parameters + ---------- + ccd_list_file : str + path to CCD list file (txt or csv) + + Returns + ------- + np.array + unique exposure numbers + + """ + exps = [] + + # Check if CSV format + if ccd_list_file.endswith(".csv"): + import astropy.table + + t = astropy.table.Table.read(ccd_list_file) + expccd = t["CCD"].data + r = np.char.split(expccd, sep="-") + + for i, r1 in enumerate(r): + exp = int(r1[0]) + exps.append(exp) + else: + # Text format + f = np.loadtxt(ccd_list_file, dtype="str", encoding="ascii") + r = np.char.split(f, sep="-") + + for i, r1 in enumerate(r): + exp = int(r1[0]) + exps.append(exp) + + exps = np.array(exps) + uniq = np.unique(exps) + + return uniq + + def get_fits_header(self, expnum, client): + """Get FITS Header. + + Download FITS header from VOSpace. + + Parameters + ---------- + expnum : int + exposure number + client : vos.Client + VOSpace client + + Returns + ------- + bool + True if successful, False otherwise + + """ + vospace_path = self._params["vospace_path"] + output_dir = self._params["output_dir"] + overwrite = self._params["overwrite"] + + source = f"{vospace_path}/{expnum:d}p.fits.fz" + dest = f"{output_dir}/{expnum:d}.txt" + + if exists(dest) and not overwrite: + return True + + try: + client.copy(source, dest, head=True) + return True + except Exception as e: + print(f"Could not copy {source}: {e}") + return False + + def run(self, args=None): + """Run. + + Main execution method. + + Parameters + ---------- + args : list, optional + command line arguments + + Returns + ------- + int + exit code (0 for success) + + """ + if args is None: + args = sys.argv + + # Set parameters from command line + self.set_params_from_command_line(args) + self.update_params() + self.check_params() + + # Get parameters + input_file = self._params["input_file"] + output_dir = self._params["output_dir"] + verbose = self._params["verbose"] + + if verbose: + print(f"Reading CCD list from: {input_file}") + + # Extract unique exposures + exps = self.get_exposures(input_file) + + print(f"Found {len(exps)} unique exposures") + + if verbose: + print(f"Downloading headers to: {output_dir}/") + + # Initialize VOSpace client once + client = vos.Client() + + # Download headers + n_success = 0 + n_failed = 0 + + for i, exp in enumerate(exps): + success = self.get_fits_header(exp, client) + if success: + n_success += 1 + else: + n_failed += 1 + + if verbose and i % 100 == 0: + print(f"{i:6d} / {len(exps):6d}") + + print(f"Downloaded {n_success} headers") + if n_failed > 0: + print(f"Failed to download {n_failed} headers") + + return 0 From bc43799f82bcb29180039089114a3eb47adbd8c7 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Wed, 14 Jan 2026 14:29:05 +0000 Subject: [PATCH 04/34] coverage mask: scripts and plotter; using cs_util --- pyproject.toml | 13 +- src/shapepipe/coverage_run.py | 149 +++++++ .../utilities/coverage_map_builder.py | 325 ++++++++++++++ src/shapepipe/utilities/coverage_plotter.py | 343 +++++++++++++++ .../utilities/field_corners_extractor.py | 405 ++++++++++++++++++ 5 files changed, 1234 insertions(+), 1 deletion(-) create mode 100644 src/shapepipe/coverage_run.py create mode 100644 src/shapepipe/utilities/coverage_map_builder.py create mode 100644 src/shapepipe/utilities/coverage_plotter.py create mode 100644 src/shapepipe/utilities/field_corners_extractor.py diff --git a/pyproject.toml b/pyproject.toml index 9a4153dbe..b20220ece 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -18,6 +18,8 @@ dependencies = [ "canfar", "sf_tools", "h5py", + "healsparse", + "skyproj", ] @@ -33,6 +35,10 @@ lint = [ "black", "isort" ] +plot = [ + "skyproj", + "matplotlib" +] release = [ "build", "twine", @@ -43,7 +49,7 @@ test = [ "pytest-pycodestyle", "pytest-pydocstyle" ] -dev = ["shapepipe[doc,lint,release,test]"] +dev = ["shapepipe[doc,lint,plot,release,test]"] [project.scripts] shapepipe_run = "shapepipe.shapepipe_run:main" @@ -52,6 +58,11 @@ canfar_submit_job = "shapepipe.canfar_run:run_job" canfar_monitor = "shapepipe.canfar_run:run_log" canfar_monitor_log = "shapepipe.canfar_run:run_monitor_log" get_ccd_with_psfs = "shapepipe.get_ccds_run:run_ccd_psf_handler" +download_headers = "shapepipe.coverage_run:run_download_headers" +extract_field_corners = "shapepipe.coverage_run:run_extract_corners" +build_coverage_map = "shapepipe.coverage_run:run_build_coverage" +plot_coverage_map = "shapepipe.coverage_run:run_plot_coverage" +coverage_pipeline = "shapepipe.coverage_run:run_pipeline" [tool.pytest.ini_options] addopts = "--verbose --cov=shapepipe" diff --git a/src/shapepipe/coverage_run.py b/src/shapepipe/coverage_run.py new file mode 100644 index 000000000..3b30ce044 --- /dev/null +++ b/src/shapepipe/coverage_run.py @@ -0,0 +1,149 @@ +"""COVERAGE_RUN + +Call coverage processing classes. + +Author: Martin Kilbinger + +""" + +import sys + +from shapepipe.utilities.header_downloader import HeaderDownloader +from shapepipe.utilities.field_corners_extractor import FieldCornersExtractor +from shapepipe.utilities.coverage_map_builder import CoverageMapBuilder +from shapepipe.utilities.coverage_plotter import CoveragePlotter + + +def run_download_headers(args=None): + """Run Download Headers. + + Download FITS headers from VOSpace for exposures in a CCD list. + + Parameters + ---------- + args : list, optional + command line arguments + + Returns + ------- + int + exit code + + """ + obj = HeaderDownloader() + return obj.run(args=args) + + +def run_extract_corners(args=None): + """Run Extract Corners. + + Extract field corner coordinates from FITS headers. + + Parameters + ---------- + args : list, optional + command line arguments + + Returns + ------- + int + exit code + + """ + obj = FieldCornersExtractor() + return obj.run(args=args) + + +def run_build_coverage(args=None): + """Run Build Coverage. + + Build HealSparse coverage maps from field corner coordinates. + + Parameters + ---------- + args : list, optional + command line arguments + + Returns + ------- + int + exit code + + """ + obj = CoverageMapBuilder() + return obj.run(args=args) + + +def run_plot_coverage(args=None): + """Run Plot Coverage. + + Plot HealSparse coverage map. + + Parameters + ---------- + args : list, optional + command line arguments + + Returns + ------- + int + exit code + + """ + obj = CoveragePlotter() + return obj.run(args=args) + + +def run_pipeline(args=None): + """Run Pipeline. + + Run the complete coverage map pipeline: + 1. Download FITS headers + 2. Extract field corners + 3. Build coverage map + + Parameters + ---------- + args : list, optional + command line arguments + + Returns + ------- + int + exit code + + """ + print("Step 1: Downloading FITS headers") + ret = run_download_headers(args) + if ret != 0: + print("Failed to download headers") + return ret + + print() + print("Step 2: Extracting field corners") + ret = run_extract_corners(args) + if ret != 0: + print("Failed to extract field corners") + return ret + + print() + print("Step 3: Building coverage map") + ret = run_build_coverage(args) + if ret != 0: + print("Failed to build coverage map") + return ret + + print() + print("Pipeline completed successfully") + + return 0 + + +def main(argv=None): + """Main. + + Main program. + + """ + # Scripts to call coverage classes are created by pyproject.toml + return 0 diff --git a/src/shapepipe/utilities/coverage_map_builder.py b/src/shapepipe/utilities/coverage_map_builder.py new file mode 100644 index 000000000..3f91f763c --- /dev/null +++ b/src/shapepipe/utilities/coverage_map_builder.py @@ -0,0 +1,325 @@ +"""COVERAGE_MAP_BUILDER + +Build HealSparse coverage maps from field corner coordinates. + +Author: Martin Kilbinger + +""" + +import sys +from os.path import exists + +import numpy as np +import healsparse as hsp +import hpgeom as hpg + +from cs_util import args as cs_args +from cs_util import logging + +from shapepipe.utilities.coverage_plotter import CoveragePlotter + + +class CoverageMapBuilder(object): + """Coverage Map Builder Class. + + Builds HealSparse coverage maps from field corner coordinates. + """ + + def __init__(self): + """Initialize the builder.""" + self.params_default() + + def params_default(self): + """Set default parameters and command line options.""" + + self._params = { + "input_file": "exp_ra_dec.txt", + "output_file": "coverage.hsp", + "nside_coverage": 32, + "nside": 2048, + "apply_median_filter": False, + "n_median_iterations": 2, + "create_boolean": False, + "boolean_threshold": 3, + "boolean_output": None, + "create_plot": False, + "plot_output": None, + "plot_region": None, + } + + self._short_options = { + "input_file": "-i", + "output_file": "-o", + "nside_coverage": "-c", + "nside": "-n", + "apply_median_filter": "-m", + "n_median_iterations": "-N", + "create_boolean": "-b", + "boolean_threshold": "-t", + "boolean_output": "-B", + "create_plot": "-p", + "plot_output": "-P", + "plot_region": "-g", + } + + self._types = { + "nside_coverage": "int", + "nside": "int", + "apply_median_filter": "bool", + "n_median_iterations": "int", + "create_boolean": "bool", + "boolean_threshold": "int", + "create_plot": "bool", + } + + self._help_strings = { + "input_file": "input file with field corners; default is {}", + "output_file": "output HealSparse map file; default is {}", + "nside_coverage": "HealSparse coverage nside; default is {}", + "nside": "HealSparse map nside; default is {}", + "apply_median_filter": "apply median filter to smooth map; default is {}", + "n_median_iterations": "number of median filter iterations; default is {}", + "create_boolean": "create boolean coverage map; default is {}", + "boolean_threshold": "threshold value for boolean map; default is {}", + "boolean_output": "output file for boolean map; default is _bool.hsp", + "create_plot": "create plot of coverage map; default is {}", + "plot_output": "output file for plot; default is .png", + "plot_region": "predefined region for plot (NGC, SGC, fullsky); default is {}", + } + + def set_params_from_command_line(self, args): + """Set Params From Command line. + + Only use when calling using python from command line. + Does not work from ipython or jupyter. + + Parameters + ---------- + args : list + command line arguments + + """ + # Read command line options + options = cs_args.parse_options( + self._params, + self._short_options, + self._types, + self._help_strings, + args=args, + ) + self._params = options + + # Save calling command + logging.log_command(args) + + def update_params(self): + """Update parameters. + + Set derived parameters based on input parameters. + """ + # Set boolean output filename if not specified + if self._params["boolean_output"] is None: + output_file = self._params["output_file"] + if output_file.endswith(".hsp"): + base = output_file[:-4] + else: + base = output_file + self._params["boolean_output"] = f"{base}_bool.hsp" + + # Set plot output filename if not specified + if self._params["plot_output"] is None: + output_file = self._params["output_file"] + if output_file.endswith(".hsp"): + base = output_file[:-4] + else: + base = output_file + self._params["plot_output"] = f"{base}.png" + + def check_params(self): + """Check parameters for validity.""" + if not exists(self._params["input_file"]): + raise FileNotFoundError( + f"Input file not found: {self._params['input_file']}" + ) + + # Check nside values are powers of 2 + nside_coverage = self._params["nside_coverage"] + nside = self._params["nside"] + + if not (nside_coverage & (nside_coverage - 1) == 0): + raise ValueError( + f"nside_coverage must be a power of 2, got {nside_coverage}" + ) + + if not (nside & (nside - 1) == 0): + raise ValueError(f"nside must be a power of 2, got {nside}") + + def median_filter(self, hsp_map): + """Median Filter. + + Apply median filter to HealSparse map using neighbors. + + Parameters + ---------- + hsp_map : healsparse.HealSparseMap + input map + + Returns + ------- + healsparse.HealSparseMap + filtered map + + """ + nside = hsp_map.nside_sparse + + new_hsp = hsp_map + 0 + pixs = hsp_map.valid_pixels + n = hpg.neighbors(nside, pixs) + n = np.where(n >= 0, n, 0) + pix2 = pixs.reshape(len(pixs), 1) + n = np.hstack([n, pix2]) + mn = hsp_map[n] + new = np.median(mn, axis=1).astype(np.uint16) + new_hsp[pixs] = new + + return new_hsp + + def run(self, args=None): + """Run. + + Main execution method. + + Parameters + ---------- + args : list, optional + command line arguments + + Returns + ------- + int + exit code (0 for success) + + """ + if args is None: + args = sys.argv + + # Set parameters from command line + self.set_params_from_command_line(args) + self.update_params() + self.check_params() + + # Get parameters + input_file = self._params["input_file"] + output_file = self._params["output_file"] + nside_coverage = self._params["nside_coverage"] + nside = self._params["nside"] + apply_median_filter = self._params["apply_median_filter"] + n_median_iterations = self._params["n_median_iterations"] + create_boolean = self._params["create_boolean"] + boolean_threshold = self._params["boolean_threshold"] + boolean_output = self._params["boolean_output"] + create_plot = self._params["create_plot"] + plot_output = self._params["plot_output"] + plot_region = self._params["plot_region"] + verbose = self._params["verbose"] + + if verbose: + print(f"Reading field corners from: {input_file}") + + # Load field corner data + dtype = np.dtype( + [ + ("expnum", int), + ("tlra", float), + ("trra", float), + ("brra", float), + ("blra", float), + ("tldec", float), + ("trdec", float), + ("brdec", float), + ("bldec", float), + ] + ) + rows = np.loadtxt(input_file, dtype=dtype) + + print(f"Loaded {len(rows)} exposure field corners") + + if verbose: + print( + f"Creating HealSparse map (nside_coverage={nside_coverage}, nside={nside})" + ) + + # Create empty map + m = hsp.HealSparseMap.make_empty(nside_coverage, nside, np.uint16) + + # Add polygons for each exposure + if verbose: + print("Adding polygons to map") + + for i, r in enumerate(rows): + ra = [r["tlra"], r["trra"], r["brra"], r["blra"]] + dec = [r["tldec"], r["trdec"], r["brdec"], r["bldec"]] + poly = hsp.Polygon(ra=ra, dec=dec, value=1) + m += poly + + if verbose and i % 1000 == 0: + print(f"{i:6d} / {len(rows):6d}") + + print(f"Added {len(rows)} polygons to map") + + # Apply median filter if requested + if apply_median_filter: + if verbose: + print( + f"Applying median filter ({n_median_iterations} iterations)" + ) + + for i in range(n_median_iterations): + m = self.median_filter(m) + if verbose: + print(f" Iteration {i+1}/{n_median_iterations} complete") + + # Write main coverage map + if verbose: + print(f"Writing coverage map to: {output_file}") + + m.write(output_file, clobber=True) + print(f"Coverage map saved to {output_file}") + + # Create and write boolean map if requested + if create_boolean: + if verbose: + print( + f"Creating boolean map (threshold={boolean_threshold})" + ) + + c = hsp.HealSparseMap.make_empty(nside_coverage, nside, bool) + c[m.valid_pixels] = m[m.valid_pixels] >= boolean_threshold + + if verbose: + print(f"Writing boolean map to: {boolean_output}") + + c.write(boolean_output, clobber=True) + print(f"Boolean map saved to {boolean_output}") + + # Create plot if requested + if create_plot: + if verbose: + print(f"Creating plot of coverage map") + + try: + plotter = CoveragePlotter() + plotter.plot_coverage_map( + m, + plot_output, + region=plot_region, + vmax=boolean_threshold if create_boolean else 3, + colorbar=True, + colorbar_label="Coverage depth", + ) + print(f"Plot saved to {plot_output}") + except ImportError as e: + print(f"Warning: Could not create plot: {e}") + print("Install sp_validation package for plotting support: pip install -e /path/to/sp_validation") + + return 0 diff --git a/src/shapepipe/utilities/coverage_plotter.py b/src/shapepipe/utilities/coverage_plotter.py new file mode 100644 index 000000000..cd8d2c2f6 --- /dev/null +++ b/src/shapepipe/utilities/coverage_plotter.py @@ -0,0 +1,343 @@ +"""COVERAGE_PLOTTER + +Plot HealSparse coverage maps. + +Author: Martin Kilbinger + +""" + +import sys +from os.path import exists + +import numpy as np +import healsparse as hsp +import matplotlib +import matplotlib.pyplot as plt + +try: + from cs_util.plots import FootprintPlotter + HAS_FOOTPRINT_PLOTTER = True +except ImportError: + HAS_FOOTPRINT_PLOTTER = False + +from cs_util import args as cs_args +from cs_util import logging + + +class CoveragePlotter(object): + """Coverage Plotter Class. + + Plots HealSparse coverage maps. + """ + + def __init__(self): + """Initialize the plotter.""" + self.params_default() + + def params_default(self): + """Set default parameters and command line options.""" + + self._params = { + "input_file": None, + "output_file": "coverage_plot.png", + "region": None, + "ralo": 100.0, + "rahi": 270.0, + "declo": 28.0, + "dechi": 85.0, + "vmin": 0, + "vmax": 3, + "cmap": "rainbow", + "draw_milky_way": True, + "colorbar": True, + "colorbar_label": "Coverage depth", + "figsize_x": 12, + "figsize_y": 8, + "dpi": 200, + } + + self._short_options = { + "input_file": "-i", + "output_file": "-o", + "region": "-g", + "ralo": "-R", + "rahi": "-r", + "declo": "-D", + "dechi": "-d", + "vmin": "-m", + "vmax": "-M", + "cmap": "-c", + "draw_milky_way": "-w", + "colorbar": "-C", + "colorbar_label": "-L", + "figsize_x": "-x", + "figsize_y": "-y", + "dpi": "-p", + } + + self._types = { + "ralo": "float", + "rahi": "float", + "declo": "float", + "dechi": "float", + "vmin": "int", + "vmax": "int", + "draw_milky_way": "bool", + "colorbar": "bool", + "figsize_x": "float", + "figsize_y": "float", + "dpi": "int", + } + + self._help_strings = { + "input_file": "input HealSparse map file; required", + "output_file": "output plot file (png, pdf, etc.); default is {}", + "region": "predefined region (NGC, SGC, fullsky); overrides RA/Dec limits if specified; default is {}", + "ralo": "minimum RA for plot extent; default is {}", + "rahi": "maximum RA for plot extent; default is {}", + "declo": "minimum Dec for plot extent; default is {}", + "dechi": "maximum Dec for plot extent; default is {}", + "vmin": "minimum value for colormap; default is {}", + "vmax": "maximum value for colormap; default is {}", + "cmap": "matplotlib colormap name; default is {}", + "draw_milky_way": "draw Milky Way outline; default is {}", + "colorbar": "add colorbar to plot; default is {}", + "colorbar_label": "colorbar label; default is {}", + "figsize_x": "figure width in inches; default is {}", + "figsize_y": "figure height in inches; default is {}", + "dpi": "figure DPI; default is {}", + } + + def set_params_from_command_line(self, args): + """Set Params From Command line. + + Only use when calling using python from command line. + Does not work from ipython or jupyter. + + Parameters + ---------- + args : list + command line arguments + + """ + # Read command line options + options = cs_args.parse_options( + self._params, + self._short_options, + self._types, + self._help_strings, + args=args, + ) + self._params = options + + # Save calling command + logging.log_command(args) + + def update_params(self): + """Update parameters. + + Set derived parameters based on input parameters. + """ + pass + + def check_params(self): + """Check parameters for validity.""" + if self._params["input_file"] is None: + raise ValueError("Input file is required (use -i or --input_file)") + + if not exists(self._params["input_file"]): + raise FileNotFoundError( + f"Input file not found: {self._params['input_file']}" + ) + + if not HAS_FOOTPRINT_PLOTTER: + raise ImportError( + "cs_util package with FootprintPlotter is required for plotting. " + "Install with: pip install -e /path/to/cs_util" + ) + + # Check region validity + region = self._params["region"] + if region is not None: + valid_regions = list(FootprintPlotter._regions.keys()) + if region not in valid_regions: + raise ValueError( + f"Invalid region '{region}'. Valid regions are: {', '.join(valid_regions)}" + ) + + def plot_coverage_map( + self, + hsp_map, + output_file, + region=None, + ralo=100.0, + rahi=270.0, + declo=28.0, + dechi=85.0, + vmin=0, + vmax=3, + title=None, + colorbar=True, + colorbar_label="Coverage depth", + figsize=(10, 10), + dpi=200, + ): + """Plot Coverage Map. + + Create a plot of the HealSparse coverage map using FootprintPlotter. + + Parameters + ---------- + hsp_map : healsparse.HealSparseMap + coverage map to plot + output_file : str + output file path + region : str, optional + predefined region name (NGC, SGC, fullsky) + ralo : float + minimum RA for plot extent (used if region is None) + rahi : float + maximum RA for plot extent (used if region is None) + declo : float + minimum Dec for plot extent (used if region is None) + dechi : float + maximum Dec for plot extent (used if region is None) + vmin : int + minimum value for colormap + vmax : int + maximum value for colormap + title : str, optional + plot title + colorbar : bool + add colorbar to plot + colorbar_label : str + colorbar label + figsize : tuple + figure size (width, height) in inches + dpi : int + figure DPI + + """ + # Set up matplotlib parameters + font = {"size": 16} + matplotlib.rc("font", **font) + matplotlib.rcParams["savefig.dpi"] = dpi + matplotlib.rcParams["figure.dpi"] = dpi + matplotlib.rcParams["figure.figsize"] = figsize + + # Create FootprintPlotter instance + plotter = FootprintPlotter( + nside_coverage=hsp_map.nside_coverage, + nside_map=hsp_map.nside_sparse + ) + + # Plot based on region or manual extent + if region is not None: + # Use predefined region + region_params = FootprintPlotter._regions[region] + plotter.plot_region( + hsp_map, + region_params, + outpath=output_file, + title=title, + colorbar=colorbar, + colorbar_label=colorbar_label + ) + else: + # Use manual RA/Dec limits + # Calculate ra_0 for projection + if ralo < rahi: + ra_0 = (ralo + rahi) / 2.0 + else: + ra_0 = 130.0 + + extend = [ralo, rahi, declo, dechi] + plotter.plot_area( + hsp_map, + ra_0=ra_0, + extend=extend, + vmax=vmax, + outpath=output_file, + title=title, + colorbar=colorbar, + colorbar_label=colorbar_label + ) + + def run(self, args=None): + """Run. + + Main execution method. + + Parameters + ---------- + args : list, optional + command line arguments + + Returns + ------- + int + exit code (0 for success) + + """ + if args is None: + args = sys.argv + + # Set parameters from command line + self.set_params_from_command_line(args) + self.update_params() + self.check_params() + + # Get parameters + input_file = self._params["input_file"] + output_file = self._params["output_file"] + region = self._params["region"] + ralo = self._params["ralo"] + rahi = self._params["rahi"] + declo = self._params["declo"] + dechi = self._params["dechi"] + vmin = self._params["vmin"] + vmax = self._params["vmax"] + colorbar = self._params["colorbar"] + colorbar_label = self._params["colorbar_label"] + figsize_x = self._params["figsize_x"] + figsize_y = self._params["figsize_y"] + dpi = self._params["dpi"] + verbose = self._params["verbose"] + + if verbose: + print(f"Reading coverage map from: {input_file}") + + # Load HealSparse map + hsp_map = hsp.HealSparseMap.read(input_file) + + print(f"Loaded map with nside={hsp_map.nside_sparse}") + print( + f"Valid pixels: {len(hsp_map.valid_pixels)} / {hsp_map.nside_sparse**2 * 12}" + ) + + if verbose: + if region: + print(f"Creating plot for region: {region}") + else: + print(f"Creating plot with extent: RA=[{ralo}, {rahi}], Dec=[{declo}, {dechi}]") + + # Create plot + self.plot_coverage_map( + hsp_map, + output_file, + region=region, + ralo=ralo, + rahi=rahi, + declo=declo, + dechi=dechi, + vmin=vmin, + vmax=vmax, + colorbar=colorbar, + colorbar_label=colorbar_label, + figsize=(figsize_x, figsize_y), + dpi=dpi, + ) + + print(f"Plot saved to {output_file}") + + return 0 diff --git a/src/shapepipe/utilities/field_corners_extractor.py b/src/shapepipe/utilities/field_corners_extractor.py new file mode 100644 index 000000000..323dd5acd --- /dev/null +++ b/src/shapepipe/utilities/field_corners_extractor.py @@ -0,0 +1,405 @@ +"""FIELD_CORNERS_EXTRACTOR + +Extract field corner coordinates from FITS headers. + +Author: Martin Kilbinger + +""" + +import sys +import os +import glob +import re +from os.path import exists +from multiprocessing import Pool, cpu_count + +import numpy as np +from astropy import wcs +from astropy.io.fits import Header + +from cs_util import args as cs_args +from cs_util import logging + + +class FieldCornersExtractor(object): + """Field Corners Extractor Class. + + Extracts RA/Dec coordinates of field corners from FITS headers. + """ + + def __init__(self): + """Initialize the extractor.""" + self.params_default() + + def params_default(self): + """Set default parameters and command line options.""" + + self._params = { + "input_dir": "header", + "output_file": "exp_ra_dec.txt", + "resume": False, + "n_processes": 1, + } + + self._short_options = { + "input_dir": "-i", + "output_file": "-o", + "resume": "-r", + "n_processes": "-n", + } + + self._types = { + "resume": "bool", + "n_processes": "int", + } + + self._help_strings = { + "input_dir": "input directory containing header files; default is {}", + "output_file": "output file for field corners; default is {}", + "resume": "resume from existing output file; default is {}", + "n_processes": f"number of parallel processes (1=serial, 0=auto={cpu_count()}); default is {{}}", + } + + def set_params_from_command_line(self, args): + """Set Params From Command line. + + Only use when calling using python from command line. + Does not work from ipython or jupyter. + + Parameters + ---------- + args : list + command line arguments + + """ + # Read command line options + options = cs_args.parse_options( + self._params, + self._short_options, + self._types, + self._help_strings, + args=args, + ) + self._params = options + + # Save calling command + logging.log_command(args) + + def update_params(self): + """Update parameters. + + Set derived parameters based on input parameters. + """ + # Ensure input directory ends without trailing slash + if self._params["input_dir"].endswith("/"): + self._params["input_dir"] = self._params["input_dir"][:-1] + + def check_params(self): + """Check parameters for validity.""" + if not exists(self._params["input_dir"]): + raise FileNotFoundError( + f"Input directory not found: {self._params['input_dir']}" + ) + + # Set n_processes to cpu_count if 0 + if self._params["n_processes"] == 0: + self._params["n_processes"] = cpu_count() + + if self._params["n_processes"] < 0: + raise ValueError( + f"n_processes must be >= 0, got {self._params['n_processes']}" + ) + + @staticmethod + def process_single_header(path_and_verbose): + """Process Single Header. + + Worker function to process a single header file. + Static method so it can be pickled for multiprocessing. + + Parameters + ---------- + path_and_verbose : tuple + (path, verbose) where path is header file path and verbose is bool + + Returns + ------- + tuple or None + (expnum, ra_list, dec_list) on success, None on failure + + """ + path, verbose = path_and_verbose + + # Extract exposure number from filename + end = path.find(".txt") + expnum = int(path[end - 6 : end]) + + try: + # Parse header and extract WCS + with open(path, "r") as f: + string = f.read() + tokens = re.split(r"^(END\s+)", string, flags=re.MULTILINE) + n = len(tokens) + w = [] + for i in range(2, n - 1, 2): + h1 = Header.fromstring(tokens[i] + tokens[i + 1], sep="\n") + w.append(wcs.WCS(h1)) + + # Calculate field corners + tl = w[0].pixel_to_world(2079, 0) + tr = w[8].pixel_to_world(32, 0) + br = w[35].pixel_to_world(2079, 0) + bl = w[27].pixel_to_world(32, 0) + + ra = [tl.ra.deg, tr.ra.deg, br.ra.deg, bl.ra.deg] + dec = [tl.dec.deg, tr.dec.deg, br.dec.deg, bl.dec.deg] + + return (expnum, ra, dec) + + except Exception as e: + if verbose: + print(f"Failed to process {expnum}: {e}") + return None + + def get_wcs_from_header(self, ftext): + """Get WCS From Header. + + Parse header text file and extract WCS for all HDUs. + + Parameters + ---------- + ftext : str + path to header text file + + Returns + ------- + list or None + list of WCS objects, one per HDU, or None if parsing failed + + """ + try: + with open(ftext, "r") as f: + string = f.read() + tokens = re.split(r"^(END\s+)", string, flags=re.MULTILINE) + n = len(tokens) + h = [] + w = [] + for i in range(2, n - 1, 2): + h1 = Header.fromstring(tokens[i] + tokens[i + 1], sep="\n") + h.append(h1) + w.append(wcs.WCS(h1)) + return w + except Exception as e: + if self._params["verbose"]: + print(f"Problem opening {ftext}: {e}") + return None + + def get_megacam_field(self, w): + """Get MegaCam Field. + + Calculate RA/Dec of the 4 corners of a MegaCam field. + + Parameters + ---------- + w : list + list of WCS objects for all CCDs + + Returns + ------- + tuple + (ra_list, dec_list) where each is a list of 4 corner coordinates + + """ + # MegaCam field corners from specific chips and pixels + # Note: chips are flipped, so bottom right becomes top left, etc. + tl = w[0].pixel_to_world(2079, 0) # top left + tr = w[8].pixel_to_world(32, 0) # top right + br = w[35].pixel_to_world(2079, 0) # bottom right + bl = w[27].pixel_to_world(32, 0) # bottom left + + return ( + [tl.ra.deg, tr.ra.deg, br.ra.deg, bl.ra.deg], + [tl.dec.deg, tr.dec.deg, br.dec.deg, bl.dec.deg], + ) + + def get_done_exposures(self): + """Get Done Exposures. + + Read exposures already processed from output file. + + Returns + ------- + set + set of exposure numbers already processed + + """ + output_file = self._params["output_file"] + + if not exists(output_file): + return set() + + try: + done = np.loadtxt(output_file, usecols=(0), dtype=int) + return set(done) + except Exception as e: + if self._params["verbose"]: + print(f"Could not read existing output file: {e}") + return set() + + def run(self, args=None): + """Run. + + Main execution method. + + Parameters + ---------- + args : list, optional + command line arguments + + Returns + ------- + int + exit code (0 for success) + + """ + if args is None: + args = sys.argv + + # Set parameters from command line + self.set_params_from_command_line(args) + self.update_params() + self.check_params() + + # Get parameters + input_dir = self._params["input_dir"] + output_file = self._params["output_file"] + resume = self._params["resume"] + verbose = self._params["verbose"] + + # Find all header files + paths = glob.glob(f"{input_dir}/*.txt") + n = len(paths) + + if n == 0: + print(f"No header files found in {input_dir}/") + return 1 + + print(f"{n} header files found") + + # Get list of already processed exposures if resuming + done = set() + if resume: + done = self.get_done_exposures() + print(f"{len(done)} already done") + + # Build todo list + todo = [] + for p in paths: + end = p.find(".txt") + expnum = int(p[end - 6 : end]) + if expnum not in done: + todo.append(p) + + n_todo = len(todo) + print(f"{n_todo} to process") + + if n_todo == 0: + print("Nothing to do") + return 0 + + # Get n_processes + n_processes = self._params["n_processes"] + + # Process headers + if n_processes == 1: + # Serial processing + results = self._process_serial(todo, verbose) + else: + # Parallel processing + print(f"Using {n_processes} parallel processes") + results = self._process_parallel(todo, n_processes, verbose) + + # Filter out failed results (None values) + results = [r for r in results if r is not None] + + # Sort results by exposure number + results.sort(key=lambda x: x[0]) + + # Write results to file + mode = "a" if resume else "w" + with open(output_file, mode, buffering=1) as f: + for expnum, ra, dec in results: + f.write(f"{expnum:6d} ") + np.savetxt(f, ra, fmt="%9.5f", newline=" ") + np.savetxt(f, dec, fmt="%9.5f", newline=" ") + f.write("\n") + + n_success = len(results) + n_failed = n_todo - n_success + + print(f"Processed {n_success} exposures") + if n_failed > 0: + print(f"Failed to process {n_failed} exposures") + + print(f"Results written to {output_file}") + + return 0 + + def _process_serial(self, todo, verbose): + """Process Serial. + + Process headers serially. + + Parameters + ---------- + todo : list + list of header file paths to process + verbose : bool + verbose output + + Returns + ------- + list + list of results (expnum, ra, dec) tuples + + """ + results = [] + n_todo = len(todo) + + for i, p in enumerate(todo): + result = self.process_single_header((p, verbose)) + results.append(result) + + if verbose and i % 100 == 0: + print(f"{i:6d} / {n_todo:6d}") + + return results + + def _process_parallel(self, todo, n_processes, verbose): + """Process Parallel. + + Process headers in parallel using multiprocessing. + + Parameters + ---------- + todo : list + list of header file paths to process + n_processes : int + number of parallel processes + verbose : bool + verbose output + + Returns + ------- + list + list of results (expnum, ra, dec) tuples + + """ + # Prepare arguments for worker function + args = [(p, verbose) for p in todo] + + # Create pool and process + with Pool(processes=n_processes) as pool: + results = pool.map(self.process_single_header, args) + + return results From 487166c9228a963ff79e6693dcd414d7ce64548f Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Wed, 14 Jan 2026 14:35:44 +0000 Subject: [PATCH 05/34] coverage code: added Mike Hudson as co-author --- src/shapepipe/utilities/ccd_psf_handler.py | 2 +- src/shapepipe/utilities/coverage_map_builder.py | 2 +- src/shapepipe/utilities/coverage_plotter.py | 2 +- src/shapepipe/utilities/field_corners_extractor.py | 2 +- src/shapepipe/utilities/header_downloader.py | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/shapepipe/utilities/ccd_psf_handler.py b/src/shapepipe/utilities/ccd_psf_handler.py index 2707cc62e..78bbcca73 100644 --- a/src/shapepipe/utilities/ccd_psf_handler.py +++ b/src/shapepipe/utilities/ccd_psf_handler.py @@ -3,7 +3,7 @@ Obtain list of CCDs (single-exposure single-HDU files) for which valid PSF information is available. This can serve to create a footprint coverage mask. -Author: Martin Kilbinger +Author: Mike Hudson, Martin Kilbinger """ diff --git a/src/shapepipe/utilities/coverage_map_builder.py b/src/shapepipe/utilities/coverage_map_builder.py index 3f91f763c..152598e16 100644 --- a/src/shapepipe/utilities/coverage_map_builder.py +++ b/src/shapepipe/utilities/coverage_map_builder.py @@ -2,7 +2,7 @@ Build HealSparse coverage maps from field corner coordinates. -Author: Martin Kilbinger +Author: Mike Hudson, Martin Kilbinger """ diff --git a/src/shapepipe/utilities/coverage_plotter.py b/src/shapepipe/utilities/coverage_plotter.py index cd8d2c2f6..14f84fc60 100644 --- a/src/shapepipe/utilities/coverage_plotter.py +++ b/src/shapepipe/utilities/coverage_plotter.py @@ -2,7 +2,7 @@ Plot HealSparse coverage maps. -Author: Martin Kilbinger +Author: Mike Hudson, Martin Kilbinger """ diff --git a/src/shapepipe/utilities/field_corners_extractor.py b/src/shapepipe/utilities/field_corners_extractor.py index 323dd5acd..f0765a047 100644 --- a/src/shapepipe/utilities/field_corners_extractor.py +++ b/src/shapepipe/utilities/field_corners_extractor.py @@ -2,7 +2,7 @@ Extract field corner coordinates from FITS headers. -Author: Martin Kilbinger +Author: Mike Hudson, Martin Kilbinger """ diff --git a/src/shapepipe/utilities/header_downloader.py b/src/shapepipe/utilities/header_downloader.py index 742acf380..e7cb70e61 100644 --- a/src/shapepipe/utilities/header_downloader.py +++ b/src/shapepipe/utilities/header_downloader.py @@ -2,7 +2,7 @@ Download FITS headers from VOSpace for exposures listed in a CCD file. -Author: Martin Kilbinger +Author: Mike Hudson, Martin Kilbinger """ From d0d5488a18660833536791700018e0a383d66159 Mon Sep 17 00:00:00 2001 From: "martin.kilbinger" Date: Fri, 16 Jan 2026 08:13:36 +0100 Subject: [PATCH 06/34] documentation (in pipeline_canfar.md) on coverage mask scripts and classes --- docs/source/pipeline_canfar.md | 513 ++++++++++++++++++++++++++------- 1 file changed, 412 insertions(+), 101 deletions(-) diff --git a/docs/source/pipeline_canfar.md b/docs/source/pipeline_canfar.md index 8b3af63d3..ee4c61995 100644 --- a/docs/source/pipeline_canfar.md +++ b/docs/source/pipeline_canfar.md @@ -1,194 +1,505 @@ -# pipeline_canfar.md +# Running `ShapePipe` processing and post-processing pipelines on CANFAR -# Documentation to create SP output products for catalogues v1.4.x +Documentation to create ShapePipe output products for catalogues v1.x. -# canfar login +## Initial Setup + +### CANFAR Login + +Login to the canfar system with + +```bash canfar auth login +``` + +This can be done from at notebook or terminal within the canfar science portal, +or any remote terminal that has the canfar library installed. -# Check; if not on "default", run canfar auth switch default +Check authentication status with + +```bash canfar auth list +``` + +If not on "default", run + +```bash +canfar auth switch default +``` + +### Set variables (optional) -# Set default PSF model +Set the current patch in the shell as + +```bash +patch=P[1-9] +``` + +For convenience, the current PSF model can be set as environment variable, e.g.: + +```bash psf="psfex" +``` + +Allowed are `psfex` and `mccd`. -# Terminal title +Setting the terminal title to display the patch can be useful for long jobs, to keep track of which terminal +runs which patch: + +```bash echo -ne "\033]0;$patch\007" +``` + +### Prepare run directory + +First, go to the dedicated directory with -# Go to run directory -patch P[1-9] +```bash +cd /path/to/version/$patch +``` -# Get tile number list +Next, set links to the tile number list and configuration directory: + +```bash ln -s ~/shapepipe/auxdir/CFIS/tiles_202106/tiles_$patch.txt tile_numbers.txt +ln -s ~/shapepipe/example/cfis +``` -# Create debug log directory -mkdir debug +Create output and debug log directories -# Get images +```bash +mkdir -p output +mkdir -p debug +``` -## Download and link separately +Finally, create and link to central image storage directories for tiles and exposures: -### Download -### Create and link to central image storage directory +```bash mkdir -p ~/cosmostat/v2/data_tiles/$patch ln -s ~/cosmostat/v2/data_tiles/$patch data_tiles mkdir -p ~/cosmostat/v2/data_exp/$patch ln -s ~/cosmostat/v2/data_tiles/$patch data_exp +``` -### Download and move tiles -ln -s ~/shapepipe/example/cfis -mkdir -p output +## `ShapePipe` processing + +Now, everything should be ready to start running `ShapePipe` for the weak lensing processing. The following +details all necessary steps. + +### Get Images + +We first download images, and in a second run create symbolic links with the proper pipeine naming scheme. + +#### Download and move tiles + +When running the main `ShapePipe` script `shapepipe_run`, the following env variable needs to point +to the current working directory + +```bash export SP_RUN=`pwd` +``` + +Now we run the first module (`get_images_runner`) to download the tile images together with the weight files. +This run can get interrupted by VOSpace I/O or connection errors. In that case, +we move new files to the image storage directory, remove the previous (now void of images) run directory, +and update the run log file. We also check the number of previous and new tiles. +```bash shapepipe_run -c cfis/config_Git_vos.ini -ls -l data_tiles/ | wc; mv -i output/run_sp_Git_*/get_images_runner/output/CFIS.???.???.*fits* data_tiles; ls -l data_tiles/ | wc -rm -rf output/run_sp_Git_*; update_runs_log_file.py -# repeat the above block +ls -l data_tiles/ | wc +mv -i output/run_sp_Git_*/get_images_runner/output/CFIS.???.???.*fits* data_tiles +ls -l data_tiles/ | wc +rm -rf output/run_sp_Git_* +update_runs_log_file.py +``` + +Repeat the above block as needed. + +### Find Exposures -### Find exposures +With all tile images (= stacks) downloaded, we can inquire their headers to identify the exposures that were used +to create the stacks. This call to the pipeline also creates the symbolic links to the downloaded tile images. + +```bash shapepipe_run -c cfis/config_GitFe_symlink.ini -# You can also run Fe alone +``` + +(One could also run `Fe` alone.) -### Download and move exposures +### Download and Move Exposures +The last module create exposure lists on output. These are now used to download all exposures. As for the tile downloads, +we have to account for VOSpace errors. + +```bash shapepipe_run -c cfis/config_Gie_vos.ini mv -i output/run_sp_Gie_*/get_images_runner/output/*.fits*fz data_exp rm -rf output/run_sp_Gie_* update_runs_log_file.py -# repeat the above; or: -while true; do shapepipe_run -c cfis/config_Gie_vos.ini; ls -l data_exp/ | wc; mv -i output/run_sp_Gie_*/get_images_runner/output/*.fits*fz data_exp; ls -l data_exp/ | wc; rm -rf output/run_sp_Gie_*; update_runs_log_file.py; done -# Make sure that after all images are downloaded there is no Gie run. This would -# mess up later modules since last:get_image_runner could point to this run. +``` + +Repeat the above by hand, or peform it in an automatic loop: + +```bash +while true; do + shapepipe_run -c cfis/config_Gie_vos.ini + ls -l data_exp/ | wc + mv -i output/run_sp_Gie_*/get_images_runner/output/*.fits*fz data_exp + ls -l data_exp/ | wc + rm -rf output/run_sp_Gie_* + update_runs_log_file.py +done +``` + +**Note:** Make sure that after all images are downloaded there is no `Gie` run in the output directory. +This would mess up later modules since `last:get_image_runner` could point to this run. + +### Create tile links again (necessary?) + +If necessary, e.g. because a previous `Git` run is no longer valid, re-create the symbolic links to the downloaded tiles with -### Create links (and re-run Fe, not necessary) +```bash job_sp_canfar.bash -p $psf `cat tile_numbers.txt` -j 1 -r symlink +``` ### Uncompress tile weights + +The downloaded tile weights are compressed. The following call uncompresses all. + +```bash shapepipe_run -c cfis/config_tile_Uz.ini +``` -# Mask tiles +### Mask tiles -## Run repeatedly if necessary +This step is done globally for all tiles. There might be job failures or interruptions. The following +command to the `ShapePipe` job script can be run repeatedly; already created masks will be skipped. + +```bash job_sp_canfar.bash -p $psf -n $OMP_NUM_THREADS -j 4 +``` + +If masks were created in more than one run, i.e. situated in more than one output directory, these have to be +combined for subsequent pipeline module runs. This is done by creating a new output directory with symbolic +links, using the script -## Combine all runs +```bash combine_runs.bash -c flag_tile +``` + +## Tile detection + +We can finally run our first module using the canfar submission system. First, determine the number of optimal jobs such that at a given +time the allowed maximum of 512 running jobs is not exceeded. +Set as `N_PAR` (number of parallel jobs) a number between 1 and 8. + +```bash +canfar_submit_job -j 16 -f tile_numbers.txt -P N_PAR -v -s +``` -# Tile detection -canfar_submit_job -j 16 -f tile_numbers.txt -v [-s | -J J] +Now, run the previous command with that number `JMAX` +```bash +canfar_submit_job -j 16 -f tile_numbers.txt -P N_PAR -v -J JMAX +``` -# Exposure processing +### Exposure Processing -# Option 0, global split and exp masks: sp_local=0 (depreciated; only for earlier v1.x patch runs) -# Todo: split Uz and SpMh +#### Option 0: Global split and exp masks (deprecated; used for earlier v1.x patch runs) -# For sp_local=- both mh_local (0, 1) are ok +For this option, set `sp_local=0`. + +**TODO: ** Split Uz and SpMh + +For `sp_local=-` both `mh_local` (0, 1) are ok: + +```bash export mh_local=0 +``` + +#### Option 0: Mask exposures (deprecated) -# Mask exposures +Run repeatedly if necessary: -## Run repeatedly if necessary +```bash job_sp_canfar.bash -p $psf -n $OMP_NUM_THREADS -j 8 +``` -# Combine all runs +Combine all runs: + +```bash combine_runs.bash -c flag_exp +``` + +### Option 1: Local split and mask exposures (recommended) + +Optional: Enable flags for local split processing and merge header runs as -# Option 1: sp_local=1, local split and mask exp -#export mh_local=1 +```bash +export sp_local=1 +export mh_local=1 +``` -# Get single-HDU single-exposure IDs file (from missing 32 job) -~/shapepipe/scripts/python/summary_run.py P$patch [32] +These flags are automatically set to 1 in the new job scripts. + + +Get single-HDU single-exposure IDs file (from missing 32 job): + +```bash +summary_run P$patch 32 cp summary/missing_job_32_all.txt exp_shdu.txt +``` + +### Split exposures + +First, determine the number of maximum jobs with the option `-s` (see above). Then, submit with + +```bash +canfar_submit_job -j 2 -v -f exp_shdu.txt -v -P N_PAR -J JMAX +``` + +### Mask exposures -# Split exposures -# (Check missing with summary_run P$patch 4096) -# Get maximum jobs per session to fit in one batch (so to run many jobs per -# replica instead of only one job per replica spread over many batches) -canfar_submit_job -j 2 -v -f exp_shdu.txt -v -P 8 -s -# Submit for real -canfar_submit_job -j 2 -v -f exp_shdu.txt -v -P 8 -J job_per_session +```bash +canfar_submit_job -j 8 -f exp_shdu.txt -v -P N_PAR -J JMAX +``` -# Mask exposures -canfar_submit_job -j 8 -f exp_shdu.txt -v -P 8 -J 137 +### Exposure detection -# Exposure detection +```bash +canfar_submit_job -j 32 -f exp_shdu.txt -v -P N_PAR -J JMAX +``` -cp summary/missing_job_32_sextractor.txt all.txt -canfar_submit_job -j 32 -f exp_shdu.txt -v -P 8 -J 137 +### Tile preparation -# Tile preparation +```bash canfar_submit_job -j 64 -f tile_numbers.txt +``` -# Tile shape measurement +### Tile shape measurement + +```bash canfar_submit_job -j 128 -f tile_numbers.txt +``` + +### Merge sub-catalogues -# Merge subcatalogues +```bash canfar_submit_job -j 256 -f tile_numbers.txt +``` -# Create final cat +### Create final catalogues + +```bash canfar_submit_job -j 512 -f tile_numbers.txt +``` -# Run in parallel -cat IDs.txt | xargs -I {} -P 16 bash -c 'init_run_exclusive_canfar.sh -j 512 -e {}' +This was the last `ShapePipe` module to run for main processing. -# Combine all final cats in common output dir as links -#combine_runs.bash -c final -p psfex -# Merge all final cats -# (W3: 140GB RAM) -cd .. -# to be in /path/to/$psf +### Merge all final catalogues + +The last step of `ShapePipe` processing is, per patch, to merget all final catalogues. This is done via a python script, as follows. +First, change to parent directory `/path/to/version` and run the following command for all patches + +```bash patchnum=`tr $patch P ''` -create_final_cat.py -m final_cat_$patch.hdf5 -i . -p $patch/cfis/final_cat.param -P $patchnum -o $patch/n_tiles_final.txt -v +create_final_cat.py -m final_cat_$patch.hdf5 -i . -p $patch/cfis/final_cat.param \ + -P $patchnum -o $patch/n_tiles_final.txt -v +``` + +## Additional `ShapePipe` processing + + +### Create star Catalogue + +We can additionaly create a combined star catalogue, with star shapes projecte from detector to world coordinates. +This is useful for validation and galaxy-PSF/star correlation diagnostics. + +#### Combine all PSF runs + +In each patch directory /path/to/version/$patch, run + +```bash +combine_runs.bash -p $psf -c psf +``` + +to create a single output directory of PSF files (symbolic links). + +Optionally, to create and plot results for this patch only: -# Star catalogue -combine_runs.bash -p $psf -c psf +```bash shapepipe_run -c $SP_CONFIG/config_Ms_$psf.ini shapepipe_run -c $SP_CONFIG/config_Pl_$psf.ini +``` -# Convert star cat to WCS -## Convert all input validation psf files and create directories par patch -## psf_conv_all/P? -cd ../star_cat +#### Convert star catalogue to wCS -# Create files validation_psf_conv--.fits -# (for the v1.4 setup only one file) +Convert all input validation PSF files and create directories per patch `P?`. +Create files `validation_psf_conv--.fits` (for the v1.4 setup only one file): + +```bash +cd /path/to/version +mkdir stat_car +cd star_cat +``` + +For each patch run + +```bash convert_psf_pix2world.py -i .. -P $patchnum -v +``` + +Combine previously created files as links within one ShapePipe run directory (for the v1.4 setup only one link). +First (and optiohnal), create a subdir for a run and link to the input patches: -# Combine previously created files as links within one SP run dir -# (for the v1.4 setup only one link -cd P$patch +```bash +cd /path/to/version/star_cat +mkdir v1.6 +ln -s ../P1 +ln -s ../P2 +... +``` + +Next, create links to all `validation_conv` runs: + +```bash combine_runs.bash -p psfex -c psf_conv +``` -# Merge all converted star catalogues and create final-starcat.fits +Merge all converted star catalogues and create `final-starcat.fits`: + +```bash export SP_RUN=`pwd` shapepipe_run -c ~/shapepipe/example/cfis/config_Ms_psfex_conv.ini +``` + +Rename to general PSF and star catalogue used for all ("a") sub-versions: + + +```bash +cp output/run_sp_Ms/merge_starcat_runner/output/full_starcat-0000000.fits \ + unions_shapepipe_psf_2024_v1.6.a.fits +``` + +The FITS file `CATTYPE` (newer version) should be `validation_psf_conf`. + +## Post-processing + +The following post-processing steps are performed with the library `sp_validation`. + +### Extract Information + +First, we extract all information from the final catalogue, per patch. We copy +the parameter file and set links to the catalogues and `ShapePipe` config directory. -# sp_validation, on candide -## Extract information and create patch-wise comprehensive catalogues +```bash +cd /path/to/version/$patch cp ~/astro/repositories/github/sp_validation/notebooks/params.py . -ln -s ~/v1.3.x/final_cat_$patchnum3.hdf5 # not ../inal_cat_P$patchnum.hdf5 ! +ln -s ~/v1.3.x/final_cat_$patchnum3.hdf5 # not ../final_cat_P$patchnum.hdf5 ! ln -s output/run_sp_MsPl/mccd_merge_starcat_runner/output/full_starcat-0000000.fits ln -s ~/astro/repositories/github/shapepipe/example/cfis -# Edit patch name; wrap_ra for P2 -# squeue -python ~/astro/repositories/github/sp_validation/notebooks/extract_info.py +``` -## Create joint comprehensive catalogue -# Edit py script to match catalogue name; check coverage mask input file -python ~/astro/repositories/github/sp_validation/notebooks/demo_apply_hsp_masks.py +Then edit `params.py`: Set patch name; set `wrap_ra` for P2. +Now we can run the script, recommended via job submission on candide. For large patches, +this requies a job with a large memory, e.g. with `mem=380000` + -# Extra stuff +```bash +[squeue] python ~/astro/repositories/github/sp_validation/notebooks/extract_info.py +``` -## Delete jobs -SSL=~/.ssl/cadcproxy.pem -SESSION=https://ws-uv.canfar.net/skaha/v0/session -for ID in `cat session_IDs.txt`; do echo $ID; curl -X DELETE -E $SSL $SESSION/$ID; done +This creates a patch-wise comprehensive catalogue. -## Run in terminal in parallel (-e needs to be last arg) -cat all.txt | xargs -P 16 -n 1 init_run_exclusive_canfar.sh -j 64 -p psfex -n -e +### Create global comprehensive catalogues -## Get missing jobs that are not currently running -stats_jobs_canfar.sh -grep -F -v -f jobs_running.txt summary/missing_job_128_ngmix_runner_3.txt > all3.txt +```bash +cd /patch/to/version +[squeue] python ~/astro/repositories/github/sp_validation/scripts/create_joint_comprehensive_cat.py \ + -v v1.6.c -v -p P1+P2+P3+P4+P5+P6+P7+P8+P9 +``` + +This creates the file `unions_shapepipe_comprehensive_2024_v1.6.c.hdf5`. + + +### Apply structural masks + +First, edit the Python script `~/astro/repositories/github/sp_validation/notebooks/demo_apply_hsp_masks.py` +to match catalogue name. Check the coverage mask input file (see below). +Run the script to apply the healsparse structural masks: + +```bash +[squeue] python ~/astro/repositories/github/sp_validation/notebooks/demo_apply_hsp_masks.py +``` + +This creates the file `unions_shapepipe_comprehensive_struct_2024_v1.6.c.hdf5`. + + +### Define sample, calibrate catalogue + +We are close to finally perform the last post-processing step, which is the calibration. First, the final galaxy sample +in question needs to be defined, with masks and cuts to apply from a `yaml` config file. A number of pre-defined files +can be found in `~/astro/repositories/github/sp_validation/calibration`. + +For example, to create `v1.6.6`, the steps are: + +```bash +cd /path/to/version +mkdir -p v1.6.6 +cd v1.6.6 +ln -s ~/astro/repositories/github/sp_validation/calibration/mask_v1.X.6.yaml config_mask.yaml +ln -s ..//unions_shapepipe_comprehensive_struct_2024_v1.6.c.hdf5 unions_shapepipe_comprehensive_struct_2024_v1.X.c.hdf5 +[squeue] python ~/astro/repositories/github/sp_validation/calibrate_comprehensive_cat.py +``` + +calibrate_comprehensive + +### Create coverage mask + +First, on canfar, move to the directory that has the patch subdirectories. + +```bash +cd /path/to/version +``` + +#### Get exposure numbers + +If the file `$patch/exp_numbers.txt` does not exist for a given patch, create it with the summary program + +```bash +summary_run $patch 1 +``` + +Now, create the list of CCDs that have PSF information with + +```bash +get_ccds_with_psf -v -V v1.6 +``` + +Next, download exposures headers + +```bash +download_headers -i ccds_with_psfs_v1.6.txt -o headers_v1.6 -v +``` + +From the headers, the CCD corner coordinates are extracted with +```bash +extract_field_corners -i headers_v1.6 -v +``` + +Then, build the healsparse coverage mask file as +```bash +build_coverage_map +``` + +Use `plot_coverage_map` to create plots of the coverage mask. + +## Extra Utilities + +### Run in Terminal in Parallel + +```bash +cat IDs.txt | xargs -I {} -P 16 bash -c 'init_run_exclusive_canfar.sh -j 512 -e {}' +``` From f4d00d0c34a7024e8c6b0103c0cdac33999d198d Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Mon, 19 Jan 2026 10:04:54 +0100 Subject: [PATCH 07/34] CCD PSF Handler: Added version for v1.3 --- src/shapepipe/utilities/ccd_psf_handler.py | 90 +++++++++++++++++++++- 1 file changed, 88 insertions(+), 2 deletions(-) diff --git a/src/shapepipe/utilities/ccd_psf_handler.py b/src/shapepipe/utilities/ccd_psf_handler.py index 78bbcca73..3ee29303c 100644 --- a/src/shapepipe/utilities/ccd_psf_handler.py +++ b/src/shapepipe/utilities/ccd_psf_handler.py @@ -7,7 +7,11 @@ """ +import glob +import os +import re import sys + import numpy as np from cs_util import args as cs_args @@ -48,6 +52,7 @@ def params_default(self): self._help_strings = { "version_cat": "catalogue major version, allowed are v1.3, v1.4, v1.5, v1.6; default is {}", + "n_CCD": "number of CCDs per exposure; default is {}", "output": "output file path; default is ccds_with_psf_.txt", } @@ -231,6 +236,83 @@ def get_ccds_with_psf(self, patches, n_CCD=40): return exp_shdu_all + def get_ccds_with_psf_method_v1_3(self, patches, n_CCD=40): + """Get CCDs With PSF Method v1.3. + + Return set of CCDs with valid PSF by scanning run directories. + Finds star_selection files and checks corresponding stat files + for valid (non-nan) FWHM values. + + Parameters + ---------- + patches : list + input patches + n_CCD : int + number of CCDs per exposure (unused, kept for API consistency) + + Returns + ------- + set + CCD IDs with valid PSF + + """ + exp_shdu_all = set() + + for patch in patches: + # Find all star_selection files + mask_pattern = ( + f"{patch}/output/run_sp_exp_SxSePs*/" + f"setools_runner/output/mask/star_selection-*.fits" + ) + mask_files = glob.glob(mask_pattern) + + print( + f"Patch {patch}: Found {len(mask_files)} star_selection FITS" + + " files" + ) + + for mask_file in mask_files: + # Extract exp_shdu from filename + basename = os.path.basename(mask_file) + match = re.match(r"star_selection-(.+)\.fits", basename) + if not match: + print(f"Warning: Non-matching file name {basename}") + continue + exp_shdu = match.group(1) + + # Get corresponding stat file path + # Replace mask dir with stat/stat dir in the path + stat_file = mask_file.replace( + "/mask/star_selection-", + "/stat/star_stat-", + ).replace(".fits", ".txt") + + #print("MKDEBUG ", exp_shdu, stat_file) + + # Check if stat file exists and has valid FWHM + if not os.path.exists(stat_file): + print("MKDEBUG stat file ", stat_file, " does not exist") + continue + + # Read stat file and check for nan in "Mean star fwhm selected" + has_nan = False + with open(stat_file) as f: + for line in f: + if "Mean star fwhm selected" in line: + if "nan" in line.lower(): + has_nan = True + break + + if not has_nan: + #print("MKDEBUG append", exp_shdu) + exp_shdu_all.add(exp_shdu) + + print(f"After patch {patch}: {len(exp_shdu_all)} CCDs with valid PSF") + + print(f"Found {len(exp_shdu_all)} CCDs with valid PSF") + + return exp_shdu_all + def save(self, IDs, path): """Save. @@ -282,8 +364,12 @@ def run(self, args=None): f"=== get_ccds_with_psf for version {version}, patches {patches} ===" ) - print("=== method 1: exp_list - missing === ") - exp_shdu_all = self.get_ccds_with_psf(patches, n_CCD) + if self._params["version_cat"] == "v1.3": + print("=== method for 1.3: exp with PSF star sample === ") + exp_shdu_all = self.get_ccds_with_psf_method_v1_3(patches, n_CCD) + else: + print("=== method for >= 1.4: exp_list - missing === ") + exp_shdu_all = self.get_ccds_with_psf(patches, n_CCD) self.save(exp_shdu_all, output) From 271ff40245dc339146eb766ded229816f39cb3d8 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Mon, 19 Jan 2026 15:56:20 +0000 Subject: [PATCH 08/34] fixed coverage plots --- src/shapepipe/utilities/coverage_plotter.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/shapepipe/utilities/coverage_plotter.py b/src/shapepipe/utilities/coverage_plotter.py index 14f84fc60..825ba9d96 100644 --- a/src/shapepipe/utilities/coverage_plotter.py +++ b/src/shapepipe/utilities/coverage_plotter.py @@ -256,6 +256,7 @@ def plot_coverage_map( hsp_map, ra_0=ra_0, extend=extend, + vmin=vmin, vmax=vmax, outpath=output_file, title=title, From 7f99f5ade9cff67dcbe3d8abcdd65fee95a501dd Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Mon, 19 Jan 2026 15:57:07 +0000 Subject: [PATCH 09/34] header downloads (coverage): account for existing downloads --- src/shapepipe/utilities/header_downloader.py | 34 +++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/src/shapepipe/utilities/header_downloader.py b/src/shapepipe/utilities/header_downloader.py index e7cb70e61..5eccf023f 100644 --- a/src/shapepipe/utilities/header_downloader.py +++ b/src/shapepipe/utilities/header_downloader.py @@ -35,6 +35,7 @@ def params_default(self): "output_dir": "header", "vospace_path": "vos:cfis/pitcairn", "overwrite": False, + "dir_for_links": None, } self._short_options = { @@ -42,6 +43,7 @@ def params_default(self): "output_dir": "-o", "vospace_path": "-p", "overwrite": "-O", + "dir_for_links": "-d", } self._types = { @@ -53,6 +55,7 @@ def params_default(self): "output_dir": "output directory for headers; default is {}", "vospace_path": "VOSpace base path; default is {}", "overwrite": "overwrite existing header files; default is {}", + "dir_for_links": "directory to check for existing headers to link instead of download; default is {}", } def set_params_from_command_line(self, args): @@ -99,6 +102,17 @@ def check_params(self): f"Input file not found: {self._params['input_file']}" ) + # Check if dir_for_links exists if specified + if self._params["dir_for_links"] is not None: + if not exists(self._params["dir_for_links"]): + raise FileNotFoundError( + f"Directory for links not found: {self._params['dir_for_links']}" + ) + if not os.path.isdir(self._params["dir_for_links"]): + raise ValueError( + f"Path is not a directory: {self._params['dir_for_links']}" + ) + # Create output directory if it doesn't exist if not exists(self._params["output_dir"]): os.makedirs(self._params["output_dir"]) @@ -151,7 +165,8 @@ def get_exposures(self, ccd_list_file): def get_fits_header(self, expnum, client): """Get FITS Header. - Download FITS header from VOSpace. + Download FITS header from VOSpace, or create a symbolic link + if the header exists in dir_for_links. Parameters ---------- @@ -169,6 +184,7 @@ def get_fits_header(self, expnum, client): vospace_path = self._params["vospace_path"] output_dir = self._params["output_dir"] overwrite = self._params["overwrite"] + dir_for_links = self._params["dir_for_links"] source = f"{vospace_path}/{expnum:d}p.fits.fz" dest = f"{output_dir}/{expnum:d}.txt" @@ -176,6 +192,22 @@ def get_fits_header(self, expnum, client): if exists(dest) and not overwrite: return True + # Check if header exists in dir_for_links + if dir_for_links is not None: + link_source = os.path.abspath(f"{dir_for_links}/{expnum:d}.txt") + if exists(link_source): + try: + # Remove existing file/link if overwrite is True + if exists(dest): + os.remove(dest) + # Create symbolic link + os.symlink(link_source, dest) + return True + except Exception as e: + print(f"Could not create symlink from {link_source}: {e}") + # Fall through to download if symlink fails + + # Download from VOSpace try: client.copy(source, dest, head=True) return True From e10111baf2e8bc3650755dcc58bfe6f1b9c176d8 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Mon, 30 Mar 2026 08:19:24 +0000 Subject: [PATCH 10/34] canfar-changes --- docs/source/pipeline_canfar.md | 13 +++- scripts/sh/build_and_plot_coverage_maps.sh | 69 ++++++++++++++++++++++ 2 files changed, 79 insertions(+), 3 deletions(-) create mode 100755 scripts/sh/build_and_plot_coverage_maps.sh diff --git a/docs/source/pipeline_canfar.md b/docs/source/pipeline_canfar.md index ee4c61995..8ad1105b7 100644 --- a/docs/source/pipeline_canfar.md +++ b/docs/source/pipeline_canfar.md @@ -478,10 +478,11 @@ Now, create the list of CCDs that have PSF information with get_ccds_with_psf -v -V v1.6 ``` -Next, download exposures headers +Next, download exposures headers; indicate (with `-d`) a directory of already +downloaded headers; those will be linked and duplicated download skipped. ```bash -download_headers -i ccds_with_psfs_v1.6.txt -o headers_v1.6 -v +download_headers -i ccds_with_psfs_v1.6.txt -o headers_v1.6 -d headers_v1.3 -v ``` From the headers, the CCD corner coordinates are extracted with @@ -491,11 +492,17 @@ extract_field_corners -i headers_v1.6 -v Then, build the healsparse coverage mask file as ```bash -build_coverage_map +build_coverage_map -i exp_ra_dec_v1.6 -o coverage_v1.6.x.hsp -c 128 -n 131072 -v ``` +The healsparse resolutions (128, 131072) match the bit masks. + + Use `plot_coverage_map` to create plots of the coverage mask. +Building and plotting for a range of versions is done with +`build_and_plot_coverage_maps.sh`. + ## Extra Utilities ### Run in Terminal in Parallel diff --git a/scripts/sh/build_and_plot_coverage_maps.sh b/scripts/sh/build_and_plot_coverage_maps.sh new file mode 100755 index 000000000..24334d984 --- /dev/null +++ b/scripts/sh/build_and_plot_coverage_maps.sh @@ -0,0 +1,69 @@ +#!/bin/bash + +# Configuration variables + +## Catalogue versions +VERSIONS=("v1.3" "v1.4" "v1.5" "v1.6") + +## Output directory of plots +OUTPUT_DIR="coverage_map_plots" + +# healsparse resolution +BUILD_NSIDE=131072 +BUILD_CHANNELS=128 + +# Common parameters +VERBOSE="-v" + +## Colorbar +PLOT_COLORBAR="-C" +PLOT_MIN=1 +PLOT_MAX=5 + +# Plot parameters +## SGC region +SGC_RA_MIN=-20 +SGC_RA_MAX=45 +SGC_DEC_MIN=18 +SGC_DEC_MAX=40 + +## NGC region +NGC_RA_MIN=110 +NGC_RA_MAX=270 +NGC_DEC_MIN=28 +NGC_DEC_MAX=90 + +# Create output directory if it doesn't exist +mkdir -p "${OUTPUT_DIR}" + +# Loop over versions +for VERSION in "${VERSIONS[@]}"; do + echo "Processing ${VERSION}..." + + # Define file paths + INPUT_EXP="exp_ra_dec_${VERSION}.txt" + COVERAGE_MAP="coverage_${VERSION}.x.hsp" + + # Build coverage map + echo " Building coverage map from ${INPUT_EXP}..." + CMD="build_coverage_map -i ${INPUT_EXP} -o ${COVERAGE_MAP} -c ${BUILD_CHANNELS} -n ${BUILD_NSIDE} ${VERBOSE}" + echo "$CMD" + $CMD + + # Plot SGC region + echo " Plotting SGC region..." + CMD="plot_coverage_map -i ${COVERAGE_MAP} -o ${OUTPUT_DIR}/coverage_${VERSION}_SGC.png ${VERBOSE} ${PLOT_COLORBAR} -R ${SGC_RA_MIN} -r ${SGC_RA_MAX} -D ${SGC_DEC_MIN} -d ${SGC_DEC_MAX} -m ${PLOT_MIN} -M ${PLOT_MAX}" + echo "$CMD" + $CMD + + # Plot NGC region + echo " Plotting NGC region..." + CMD="plot_coverage_map -i ${COVERAGE_MAP} -o ${OUTPUT_DIR}/coverage_${VERSION}_NGC.png ${VERBOSE} ${PLOT_COLORBAR} -R ${NGC_RA_MIN} -r ${NGC_RA_MAX} -D ${NGC_DEC_MIN} -d ${NGC_DEC_MAX} -m ${PLOT_MIN} -M ${PLOT_MAX}" + echo "$CMD" + $CMD + + echo " Done with ${VERSION}" + echo +done + +echo "All versions processed successfully!" From fd611209f5b4402513e372064c41593ce296eb70 Mon Sep 17 00:00:00 2001 From: martin kilbinger Date: Tue, 21 Apr 2026 14:58:53 +0200 Subject: [PATCH 11/34] PR: remove non-functioning coverage pipeline function --- src/shapepipe/coverage_run.py | 45 ----------------------------------- 1 file changed, 45 deletions(-) diff --git a/src/shapepipe/coverage_run.py b/src/shapepipe/coverage_run.py index 3b30ce044..2485e4a10 100644 --- a/src/shapepipe/coverage_run.py +++ b/src/shapepipe/coverage_run.py @@ -94,51 +94,6 @@ def run_plot_coverage(args=None): return obj.run(args=args) -def run_pipeline(args=None): - """Run Pipeline. - - Run the complete coverage map pipeline: - 1. Download FITS headers - 2. Extract field corners - 3. Build coverage map - - Parameters - ---------- - args : list, optional - command line arguments - - Returns - ------- - int - exit code - - """ - print("Step 1: Downloading FITS headers") - ret = run_download_headers(args) - if ret != 0: - print("Failed to download headers") - return ret - - print() - print("Step 2: Extracting field corners") - ret = run_extract_corners(args) - if ret != 0: - print("Failed to extract field corners") - return ret - - print() - print("Step 3: Building coverage map") - ret = run_build_coverage(args) - if ret != 0: - print("Failed to build coverage map") - return ret - - print() - print("Pipeline completed successfully") - - return 0 - - def main(argv=None): """Main. From fc37cbf57afed2e0700d80d23c41b7e6345b707f Mon Sep 17 00:00:00 2001 From: martin kilbinger Date: Tue, 21 Apr 2026 15:03:01 +0200 Subject: [PATCH 12/34] PR: coverage_map_builder; moved import into try block --- src/shapepipe/utilities/coverage_map_builder.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/shapepipe/utilities/coverage_map_builder.py b/src/shapepipe/utilities/coverage_map_builder.py index 152598e16..930f92127 100644 --- a/src/shapepipe/utilities/coverage_map_builder.py +++ b/src/shapepipe/utilities/coverage_map_builder.py @@ -16,7 +16,6 @@ from cs_util import args as cs_args from cs_util import logging -from shapepipe.utilities.coverage_plotter import CoveragePlotter class CoverageMapBuilder(object): @@ -308,6 +307,7 @@ def run(self, args=None): print(f"Creating plot of coverage map") try: + from shapepipe.utilities.coverage_plotter import CoveragePlotter plotter = CoveragePlotter() plotter.plot_coverage_map( m, @@ -320,6 +320,6 @@ def run(self, args=None): print(f"Plot saved to {plot_output}") except ImportError as e: print(f"Warning: Could not create plot: {e}") - print("Install sp_validation package for plotting support: pip install -e /path/to/sp_validation") + print("Install te cs_util package for plotting support.") return 0 From afcbae674aa9dabe9d9f2b03cdf1bba070396140 Mon Sep 17 00:00:00 2001 From: martin kilbinger Date: Tue, 21 Apr 2026 15:10:27 +0200 Subject: [PATCH 13/34] PR: field_corners_extractor.py; fixed exp number extraction --- src/shapepipe/utilities/field_corners_extractor.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/shapepipe/utilities/field_corners_extractor.py b/src/shapepipe/utilities/field_corners_extractor.py index f0765a047..5932e47ad 100644 --- a/src/shapepipe/utilities/field_corners_extractor.py +++ b/src/shapepipe/utilities/field_corners_extractor.py @@ -131,8 +131,8 @@ def process_single_header(path_and_verbose): path, verbose = path_and_verbose # Extract exposure number from filename - end = path.find(".txt") - expnum = int(path[end - 6 : end]) + match = re.search(r'(\d+)\.txt') + expnum = int(match.group(1)) if match else None try: # Parse header and extract WCS From 96cc4f580036aed41ae2c157afd8d8d8086262c7 Mon Sep 17 00:00:00 2001 From: martin kilbinger Date: Tue, 21 Apr 2026 15:12:33 +0200 Subject: [PATCH 14/34] PR: added comment --- src/shapepipe/utilities/header_downloader.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/shapepipe/utilities/header_downloader.py b/src/shapepipe/utilities/header_downloader.py index 5eccf023f..2e5abf9e7 100644 --- a/src/shapepipe/utilities/header_downloader.py +++ b/src/shapepipe/utilities/header_downloader.py @@ -33,7 +33,7 @@ def params_default(self): self._params = { "input_file": None, "output_dir": "header", - "vospace_path": "vos:cfis/pitcairn", + "vospace_path": "vos:cfis/pitcairn", # Default UNIONS exposure path "overwrite": False, "dir_for_links": None, } From 521578839d4a2b77959467c7a142b61e6bc8cb7b Mon Sep 17 00:00:00 2001 From: martin kilbinger Date: Tue, 21 Apr 2026 15:39:35 +0200 Subject: [PATCH 15/34] PR: added comment --- src/shapepipe/utilities/field_corners_extractor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/shapepipe/utilities/field_corners_extractor.py b/src/shapepipe/utilities/field_corners_extractor.py index 5932e47ad..9945f40d0 100644 --- a/src/shapepipe/utilities/field_corners_extractor.py +++ b/src/shapepipe/utilities/field_corners_extractor.py @@ -145,7 +145,7 @@ def process_single_header(path_and_verbose): h1 = Header.fromstring(tokens[i] + tokens[i + 1], sep="\n") w.append(wcs.WCS(h1)) - # Calculate field corners + # Calculate field corners; hard-coded values for MegaCAM tl = w[0].pixel_to_world(2079, 0) tr = w[8].pixel_to_world(32, 0) br = w[35].pixel_to_world(2079, 0) From bbe53c562efcf20c5427385e2cd10dbcadfad91f Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 22 Apr 2026 09:15:04 +0200 Subject: [PATCH 16/34] refactor: consolidate field_corners_extractor helpers; fix regex - Extract _expnum_from_path, _parse_header_to_wcs, _megacam_field_corners as module-level helpers. process_single_header (static, for multiprocessing) and the instance methods now share one implementation. - Remove dead get_wcs_from_header and get_megacam_field (duplicated logic that had already been inlined into process_single_header). - Fix re.search(r'(\d+)\.txt') missing its subject string. The new helper also replaces the hard-coded p[end-6:end] exposure-number extraction in run(), so exposure numbers of any length work. - Declare verbose in params_default so the contract is local instead of relying on cs_util.args.parse_options to inject it. Addresses review feedback on #699. --- .../utilities/field_corners_extractor.py | 136 ++++++------------ 1 file changed, 45 insertions(+), 91 deletions(-) diff --git a/src/shapepipe/utilities/field_corners_extractor.py b/src/shapepipe/utilities/field_corners_extractor.py index 9945f40d0..17a4ef7fc 100644 --- a/src/shapepipe/utilities/field_corners_extractor.py +++ b/src/shapepipe/utilities/field_corners_extractor.py @@ -20,6 +20,45 @@ from cs_util import args as cs_args from cs_util import logging +# Filename suffix convention for header files: <6+ digit exposure number>.txt +_EXPNUM_RE = re.compile(r'(\d+)\.txt$') + + +def _expnum_from_path(path): + """Extract exposure number from a header filename like ``1234567.txt``.""" + match = _EXPNUM_RE.search(path) + if match is None: + raise ValueError(f"Could not extract exposure number from {path!r}") + return int(match.group(1)) + + +def _parse_header_to_wcs(path): + """Parse a multi-HDU header text file and return one WCS per HDU.""" + with open(path, "r") as f: + string = f.read() + tokens = re.split(r"^(END\s+)", string, flags=re.MULTILINE) + return [ + wcs.WCS(Header.fromstring(tokens[i] + tokens[i + 1], sep="\n")) + for i in range(2, len(tokens) - 1, 2) + ] + + +def _megacam_field_corners(w): + """Return (ra, dec) of the 4 MegaCam field corners. + + The hard-coded CCD indices (0, 8, 35, 27) and pixel coords + (2079/32, 0) are the MegaCam corner convention; CCDs are flipped so + the nominal bottom-right CCD reports the top-left sky corner, etc. + """ + tl = w[0].pixel_to_world(2079, 0) + tr = w[8].pixel_to_world(32, 0) + br = w[35].pixel_to_world(2079, 0) + bl = w[27].pixel_to_world(32, 0) + return ( + [tl.ra.deg, tr.ra.deg, br.ra.deg, bl.ra.deg], + [tl.dec.deg, tr.dec.deg, br.dec.deg, bl.dec.deg], + ) + class FieldCornersExtractor(object): """Field Corners Extractor Class. @@ -39,6 +78,7 @@ def params_default(self): "output_file": "exp_ra_dec.txt", "resume": False, "n_processes": 1, + "verbose": False, } self._short_options = { @@ -51,6 +91,7 @@ def params_default(self): self._types = { "resume": "bool", "n_processes": "int", + "verbose": "bool", } self._help_strings = { @@ -129,99 +170,17 @@ def process_single_header(path_and_verbose): """ path, verbose = path_and_verbose - - # Extract exposure number from filename - match = re.search(r'(\d+)\.txt') - expnum = int(match.group(1)) if match else None + expnum = _expnum_from_path(path) try: - # Parse header and extract WCS - with open(path, "r") as f: - string = f.read() - tokens = re.split(r"^(END\s+)", string, flags=re.MULTILINE) - n = len(tokens) - w = [] - for i in range(2, n - 1, 2): - h1 = Header.fromstring(tokens[i] + tokens[i + 1], sep="\n") - w.append(wcs.WCS(h1)) - - # Calculate field corners; hard-coded values for MegaCAM - tl = w[0].pixel_to_world(2079, 0) - tr = w[8].pixel_to_world(32, 0) - br = w[35].pixel_to_world(2079, 0) - bl = w[27].pixel_to_world(32, 0) - - ra = [tl.ra.deg, tr.ra.deg, br.ra.deg, bl.ra.deg] - dec = [tl.dec.deg, tr.dec.deg, br.dec.deg, bl.dec.deg] - + w = _parse_header_to_wcs(path) + ra, dec = _megacam_field_corners(w) return (expnum, ra, dec) - except Exception as e: if verbose: print(f"Failed to process {expnum}: {e}") return None - def get_wcs_from_header(self, ftext): - """Get WCS From Header. - - Parse header text file and extract WCS for all HDUs. - - Parameters - ---------- - ftext : str - path to header text file - - Returns - ------- - list or None - list of WCS objects, one per HDU, or None if parsing failed - - """ - try: - with open(ftext, "r") as f: - string = f.read() - tokens = re.split(r"^(END\s+)", string, flags=re.MULTILINE) - n = len(tokens) - h = [] - w = [] - for i in range(2, n - 1, 2): - h1 = Header.fromstring(tokens[i] + tokens[i + 1], sep="\n") - h.append(h1) - w.append(wcs.WCS(h1)) - return w - except Exception as e: - if self._params["verbose"]: - print(f"Problem opening {ftext}: {e}") - return None - - def get_megacam_field(self, w): - """Get MegaCam Field. - - Calculate RA/Dec of the 4 corners of a MegaCam field. - - Parameters - ---------- - w : list - list of WCS objects for all CCDs - - Returns - ------- - tuple - (ra_list, dec_list) where each is a list of 4 corner coordinates - - """ - # MegaCam field corners from specific chips and pixels - # Note: chips are flipped, so bottom right becomes top left, etc. - tl = w[0].pixel_to_world(2079, 0) # top left - tr = w[8].pixel_to_world(32, 0) # top right - br = w[35].pixel_to_world(2079, 0) # bottom right - bl = w[27].pixel_to_world(32, 0) # bottom left - - return ( - [tl.ra.deg, tr.ra.deg, br.ra.deg, bl.ra.deg], - [tl.dec.deg, tr.dec.deg, br.dec.deg, bl.dec.deg], - ) - def get_done_exposures(self): """Get Done Exposures. @@ -293,12 +252,7 @@ def run(self, args=None): print(f"{len(done)} already done") # Build todo list - todo = [] - for p in paths: - end = p.find(".txt") - expnum = int(p[end - 6 : end]) - if expnum not in done: - todo.append(p) + todo = [p for p in paths if _expnum_from_path(p) not in done] n_todo = len(todo) print(f"{n_todo} to process") From 134a79f94bad77fec1709ba83355df9b8112b0b7 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 22 Apr 2026 09:15:07 +0200 Subject: [PATCH 17/34] refactor: coverage_map_builder cleanup - Replace hsp_map + 0 with hsp_map.copy() (healsparse supports .copy since at least 1.x; more explicit intent). - Declare verbose in params_default; currently populated only by cs_util.args.parse_options auto-injection. - Fix typo: "Install te cs_util" -> "Install the cs_util". Addresses review feedback on #699. --- src/shapepipe/utilities/coverage_map_builder.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/shapepipe/utilities/coverage_map_builder.py b/src/shapepipe/utilities/coverage_map_builder.py index 930f92127..d2eb8e1bc 100644 --- a/src/shapepipe/utilities/coverage_map_builder.py +++ b/src/shapepipe/utilities/coverage_map_builder.py @@ -44,6 +44,7 @@ def params_default(self): "create_plot": False, "plot_output": None, "plot_region": None, + "verbose": False, } self._short_options = { @@ -69,6 +70,7 @@ def params_default(self): "create_boolean": "bool", "boolean_threshold": "int", "create_plot": "bool", + "verbose": "bool", } self._help_strings = { @@ -171,7 +173,7 @@ def median_filter(self, hsp_map): """ nside = hsp_map.nside_sparse - new_hsp = hsp_map + 0 + new_hsp = hsp_map.copy() pixs = hsp_map.valid_pixels n = hpg.neighbors(nside, pixs) n = np.where(n >= 0, n, 0) @@ -320,6 +322,6 @@ def run(self, args=None): print(f"Plot saved to {plot_output}") except ImportError as e: print(f"Warning: Could not create plot: {e}") - print("Install te cs_util package for plotting support.") + print("Install the cs_util package for plotting support.") return 0 From 9a574e94f14480fd62ab1fadddffbbe29bdedade Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 22 Apr 2026 09:15:12 +0200 Subject: [PATCH 18/34] docs: declare verbose and annotate survey-specific defaults - coverage_plotter, header_downloader: declare verbose explicitly in params_default. - coverage_plotter: comment acknowledging that FootprintPlotter._regions is a private cs_util attribute; refactor to public accessor when cs_util exposes one. - header_downloader: expand the vos:cfis/pitcairn comment to flag it as UNIONS/CFIS-specific and document override. - build_and_plot_coverage_maps.sh: explain the BUILD_NSIDE=131072 choice (matches UNIONS bit-mask pixel scale, ~0.1"), and note the CoverageMapBuilder default of 2048 for lighter use. Addresses review feedback on #699. --- scripts/sh/build_and_plot_coverage_maps.sh | 5 ++++- src/shapepipe/utilities/coverage_plotter.py | 6 +++++- src/shapepipe/utilities/header_downloader.py | 4 +++- 3 files changed, 12 insertions(+), 3 deletions(-) diff --git a/scripts/sh/build_and_plot_coverage_maps.sh b/scripts/sh/build_and_plot_coverage_maps.sh index 24334d984..c35216963 100755 --- a/scripts/sh/build_and_plot_coverage_maps.sh +++ b/scripts/sh/build_and_plot_coverage_maps.sh @@ -8,7 +8,10 @@ VERSIONS=("v1.3" "v1.4" "v1.5" "v1.6") ## Output directory of plots OUTPUT_DIR="coverage_map_plots" -# healsparse resolution +# healsparse resolution. nside=131072 gives ~0.1" per pixel, chosen to +# match the UNIONS bit-mask resolution so coverage and mask align pixel- +# wise. CoverageMapBuilder defaults to nside=2048 for lighter-weight +# offline use; override here for the production build. BUILD_NSIDE=131072 BUILD_CHANNELS=128 diff --git a/src/shapepipe/utilities/coverage_plotter.py b/src/shapepipe/utilities/coverage_plotter.py index 825ba9d96..12cfce8d7 100644 --- a/src/shapepipe/utilities/coverage_plotter.py +++ b/src/shapepipe/utilities/coverage_plotter.py @@ -54,6 +54,7 @@ def params_default(self): "figsize_x": 12, "figsize_y": 8, "dpi": 200, + "verbose": False, } self._short_options = { @@ -87,6 +88,7 @@ def params_default(self): "figsize_x": "float", "figsize_y": "float", "dpi": "int", + "verbose": "bool", } self._help_strings = { @@ -156,7 +158,9 @@ def check_params(self): "Install with: pip install -e /path/to/cs_util" ) - # Check region validity + # Check region validity. Reaches into FootprintPlotter._regions + # because cs_util doesn't expose a public region registry yet; + # update when it does. region = self._params["region"] if region is not None: valid_regions = list(FootprintPlotter._regions.keys()) diff --git a/src/shapepipe/utilities/header_downloader.py b/src/shapepipe/utilities/header_downloader.py index 2e5abf9e7..1a6365b04 100644 --- a/src/shapepipe/utilities/header_downloader.py +++ b/src/shapepipe/utilities/header_downloader.py @@ -33,9 +33,10 @@ def params_default(self): self._params = { "input_file": None, "output_dir": "header", - "vospace_path": "vos:cfis/pitcairn", # Default UNIONS exposure path + "vospace_path": "vos:cfis/pitcairn", # UNIONS/CFIS default; override for other surveys "overwrite": False, "dir_for_links": None, + "verbose": False, } self._short_options = { @@ -48,6 +49,7 @@ def params_default(self): self._types = { "overwrite": "bool", + "verbose": "bool", } self._help_strings = { From dfeff38a163681dd8700405d896df4c019dca61f Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 14:07:14 +0200 Subject: [PATCH 19/34] =?UTF-8?q?chore:=20uv=20lock=20=E2=80=94=20register?= =?UTF-8?q?=20healsparse=20+=20skyproj=20as=20shapepipe=20deps?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The coverage feature adds healsparse (map building) and skyproj (plotting) to pyproject. healsparse was already resolved transitively via cs_util; this pins skyproj 2.5.0 and records both as shapepipe's own direct dependencies so 'uv lock --check' passes. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01HvP7uKkfDU2PvWC7YmbuMW --- uv.lock | 33 +++++++++++++++++++++++++++++++-- 1 file changed, 31 insertions(+), 2 deletions(-) diff --git a/uv.lock b/uv.lock index 26645f04b..3c8cd5c6b 100644 --- a/uv.lock +++ b/uv.lock @@ -3430,6 +3430,7 @@ dependencies = [ { name = "cs-util", marker = "sys_platform == 'linux'" }, { name = "galsim", marker = "sys_platform == 'linux'" }, { name = "h5py", marker = "sys_platform == 'linux'" }, + { name = "healsparse", marker = "sys_platform == 'linux'" }, { name = "joblib", marker = "sys_platform == 'linux'" }, { name = "matplotlib", marker = "sys_platform == 'linux'" }, { name = "mccd", marker = "sys_platform == 'linux'" }, @@ -3446,6 +3447,7 @@ dependencies = [ { name = "reproject", marker = "sys_platform == 'linux'" }, { name = "sf-tools", marker = "sys_platform == 'linux'" }, { name = "skaha", marker = "sys_platform == 'linux'" }, + { name = "skyproj", marker = "sys_platform == 'linux'" }, { name = "sqlitedict", marker = "sys_platform == 'linux'" }, { name = "termcolor", marker = "sys_platform == 'linux'" }, { name = "tqdm", marker = "sys_platform == 'linux'" }, @@ -3461,12 +3463,14 @@ dev = [ { name = "ipython", marker = "sys_platform == 'linux'" }, { name = "isort", marker = "sys_platform == 'linux'" }, { name = "jupyterlab", marker = "sys_platform == 'linux'" }, + { name = "matplotlib", marker = "sys_platform == 'linux'" }, { name = "myst-parser", marker = "sys_platform == 'linux'" }, { name = "numpydoc", marker = "sys_platform == 'linux'" }, { name = "pytest", marker = "sys_platform == 'linux'" }, { name = "pytest-cov", marker = "sys_platform == 'linux'" }, { name = "pytest-pycodestyle", marker = "sys_platform == 'linux'" }, { name = "pytest-pydocstyle", marker = "sys_platform == 'linux'" }, + { name = "skyproj", marker = "sys_platform == 'linux'" }, { name = "snakemake", marker = "sys_platform == 'linux'" }, { name = "sphinx", marker = "sys_platform == 'linux'" }, { name = "sphinx-book-theme", marker = "sys_platform == 'linux'" }, @@ -3492,6 +3496,10 @@ lint = [ { name = "black", marker = "sys_platform == 'linux'" }, { name = "isort", marker = "sys_platform == 'linux'" }, ] +plot = [ + { name = "matplotlib", marker = "sys_platform == 'linux'" }, + { name = "skyproj", marker = "sys_platform == 'linux'" }, +] release = [ { name = "build", marker = "sys_platform == 'linux'" }, { name = "twine", marker = "sys_platform == 'linux'" }, @@ -3515,12 +3523,14 @@ requires-dist = [ { name = "fitsio", marker = "extra == 'fitsio'" }, { name = "galsim", specifier = ">=2.8" }, { name = "h5py" }, + { name = "healsparse" }, { name = "hypothesis", marker = "extra == 'test'", specifier = ">=6.155.2" }, { name = "ipython", marker = "extra == 'jupyter'", specifier = ">=9.14.1" }, { name = "isort", marker = "extra == 'lint'" }, { name = "joblib", specifier = ">=1.4" }, { name = "jupyterlab", marker = "extra == 'jupyter'", specifier = ">=4.5.9" }, { name = "matplotlib", specifier = ">=3.10" }, + { name = "matplotlib", marker = "extra == 'plot'", specifier = ">=3.10" }, { name = "mccd", specifier = ">=1.2.4" }, { name = "modopt", specifier = ">=1.6" }, { name = "mpi4py", specifier = ">=4.1.2" }, @@ -3540,8 +3550,10 @@ requires-dist = [ { name = "python-pysap", specifier = ">=0.3" }, { name = "reproject", specifier = ">=0.19" }, { name = "sf-tools", specifier = ">=2.0.4" }, - { name = "shapepipe", extras = ["doc", "jupyter", "lint", "release", "test", "fitsio"], marker = "extra == 'dev'" }, + { name = "shapepipe", extras = ["doc", "jupyter", "lint", "plot", "release", "test", "fitsio"], marker = "extra == 'dev'" }, { name = "skaha", specifier = ">=1.7" }, + { name = "skyproj" }, + { name = "skyproj", marker = "extra == 'plot'" }, { name = "snakemake", marker = "extra == 'jupyter'", specifier = ">=9.22.0" }, { name = "sphinx", marker = "extra == 'doc'" }, { name = "sphinx-book-theme", marker = "extra == 'doc'" }, @@ -3552,7 +3564,7 @@ requires-dist = [ { name = "twine", marker = "extra == 'release'" }, { name = "vos", specifier = ">=3.6" }, ] -provides-extras = ["doc", "jupyter", "lint", "release", "test", "fitsio", "dev"] +provides-extras = ["doc", "jupyter", "lint", "plot", "release", "test", "fitsio", "dev"] [[package]] name = "shellingham" @@ -3590,6 +3602,23 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/d9/d3/a3c1d569ae714ed9985adccefdfe2d6dd9f8ad7a297bd8fa4c7ed30c587b/skaha-1.7.0-py3-none-any.whl", hash = "sha256:68b0d3c925b98bf145c5f695237474fe1cea07cf9b69b4a9bae9910375dfa01a", size = 40795, upload-time = "2025-05-28T21:17:02.296Z" }, ] +[[package]] +name = "skyproj" +version = "2.5.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "astropy", marker = "sys_platform == 'linux'" }, + { name = "healsparse", marker = "sys_platform == 'linux'" }, + { name = "hpgeom", marker = "sys_platform == 'linux'" }, + { name = "matplotlib", marker = "sys_platform == 'linux'" }, + { name = "numpy", marker = "sys_platform == 'linux'" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/90/cf/5158151ae6fc60459f430451b612a4bf2da07ebff2a9f52a37e73715a5ee/skyproj-2.5.0.tar.gz", hash = "sha256:cb8d5115927ca43cacdb6d92f00bb4fcc8563ee71028f755203ce1752e67c140", size = 8521587, upload-time = "2026-06-26T19:51:15.378Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/a8/b1/1d86ee4ae4ae7a47e31cc770dc83057388ef16311498283ca90046b3bd04/skyproj-2.5.0-cp312-cp312-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:4cbfcbd4d3d600557ab19bb5cf79b750ec552da054e3be781bce88e996b9955d", size = 4146686, upload-time = "2026-06-26T19:51:04.011Z" }, + { url = "https://files.pythonhosted.org/packages/6e/e9/f2dd235dcb9f6c3b2fec36ef04d7dd36aef4daf9afc4ab4390d6230c96dc/skyproj-2.5.0-cp313-cp313-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:78618084c9362b62ddbfd336ca1864c94d99e2a1f26d2fbf3466fd19dc702c58", size = 4146630, upload-time = "2026-06-26T19:51:10.552Z" }, +] + [[package]] name = "slicerator" version = "1.1.0" From 8b9387b345c837ed2ed025f63a68a2e997caed33 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 14:46:29 +0200 Subject: [PATCH 20/34] fix: coverage console scripts handle -h; drop dead entry point The five coverage/get-ccds runners fell back to `args = sys.argv` (the full argv, including argv[0]) instead of `sys.argv[1:]`, so under the entry-point smoke test's `python -c "...; run_*()" -h` the parser read the program slot as the integer `-c` option and exited non-zero (`invalid integer value: '-h'`). Slice argv[1:] to match cs_util.args.parse_options' own convention. Also drop the `coverage_pipeline` console script, whose target `coverage_run:run_pipeline` was removed in an earlier commit (dangling entry point -> AttributeError on invocation), and declare `hpgeom` explicitly (imported directly by coverage_map_builder, previously only transitive via healsparse). Together these fix all six entry-point test failures. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01HvP7uKkfDU2PvWC7YmbuMW --- pyproject.toml | 2 +- src/shapepipe/utilities/ccd_psf_handler.py | 2 +- src/shapepipe/utilities/coverage_map_builder.py | 2 +- src/shapepipe/utilities/coverage_plotter.py | 2 +- src/shapepipe/utilities/field_corners_extractor.py | 2 +- src/shapepipe/utilities/header_downloader.py | 2 +- uv.lock | 2 ++ 7 files changed, 8 insertions(+), 6 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index fa013d02d..adbf90f43 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,6 +24,7 @@ dependencies = [ "galsim>=2.8", "h5py", "healsparse", + "hpgeom", "joblib>=1.4", "matplotlib>=3.10", "mccd>=1.2.4", @@ -94,7 +95,6 @@ download_headers = "shapepipe.coverage_run:run_download_headers" extract_field_corners = "shapepipe.coverage_run:run_extract_corners" build_coverage_map = "shapepipe.coverage_run:run_build_coverage" plot_coverage_map = "shapepipe.coverage_run:run_plot_coverage" -coverage_pipeline = "shapepipe.coverage_run:run_pipeline" [tool.uv] # shapepipe targets Linux only; skip Windows/macOS wheel resolution diff --git a/src/shapepipe/utilities/ccd_psf_handler.py b/src/shapepipe/utilities/ccd_psf_handler.py index 3ee29303c..5c5e1d7ad 100644 --- a/src/shapepipe/utilities/ccd_psf_handler.py +++ b/src/shapepipe/utilities/ccd_psf_handler.py @@ -347,7 +347,7 @@ def run(self, args=None): """ if args is None: - args = sys.argv + args = sys.argv[1:] # Set parameters from command line self.set_params_from_command_line(args) diff --git a/src/shapepipe/utilities/coverage_map_builder.py b/src/shapepipe/utilities/coverage_map_builder.py index d2eb8e1bc..1ad6b286d 100644 --- a/src/shapepipe/utilities/coverage_map_builder.py +++ b/src/shapepipe/utilities/coverage_map_builder.py @@ -202,7 +202,7 @@ def run(self, args=None): """ if args is None: - args = sys.argv + args = sys.argv[1:] # Set parameters from command line self.set_params_from_command_line(args) diff --git a/src/shapepipe/utilities/coverage_plotter.py b/src/shapepipe/utilities/coverage_plotter.py index 12cfce8d7..7d837c2aa 100644 --- a/src/shapepipe/utilities/coverage_plotter.py +++ b/src/shapepipe/utilities/coverage_plotter.py @@ -285,7 +285,7 @@ def run(self, args=None): """ if args is None: - args = sys.argv + args = sys.argv[1:] # Set parameters from command line self.set_params_from_command_line(args) diff --git a/src/shapepipe/utilities/field_corners_extractor.py b/src/shapepipe/utilities/field_corners_extractor.py index 17a4ef7fc..aa7df9889 100644 --- a/src/shapepipe/utilities/field_corners_extractor.py +++ b/src/shapepipe/utilities/field_corners_extractor.py @@ -222,7 +222,7 @@ def run(self, args=None): """ if args is None: - args = sys.argv + args = sys.argv[1:] # Set parameters from command line self.set_params_from_command_line(args) diff --git a/src/shapepipe/utilities/header_downloader.py b/src/shapepipe/utilities/header_downloader.py index 1a6365b04..9f995647c 100644 --- a/src/shapepipe/utilities/header_downloader.py +++ b/src/shapepipe/utilities/header_downloader.py @@ -234,7 +234,7 @@ def run(self, args=None): """ if args is None: - args = sys.argv + args = sys.argv[1:] # Set parameters from command line self.set_params_from_command_line(args) diff --git a/uv.lock b/uv.lock index 3c8cd5c6b..81542690f 100644 --- a/uv.lock +++ b/uv.lock @@ -3431,6 +3431,7 @@ dependencies = [ { name = "galsim", marker = "sys_platform == 'linux'" }, { name = "h5py", marker = "sys_platform == 'linux'" }, { name = "healsparse", marker = "sys_platform == 'linux'" }, + { name = "hpgeom", marker = "sys_platform == 'linux'" }, { name = "joblib", marker = "sys_platform == 'linux'" }, { name = "matplotlib", marker = "sys_platform == 'linux'" }, { name = "mccd", marker = "sys_platform == 'linux'" }, @@ -3524,6 +3525,7 @@ requires-dist = [ { name = "galsim", specifier = ">=2.8" }, { name = "h5py" }, { name = "healsparse" }, + { name = "hpgeom" }, { name = "hypothesis", marker = "extra == 'test'", specifier = ">=6.155.2" }, { name = "ipython", marker = "extra == 'jupyter'", specifier = ">=9.14.1" }, { name = "isort", marker = "extra == 'lint'" }, From 37d5af84a485723c5b7c75e1a90893b8bce0f025 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 14:46:29 +0200 Subject: [PATCH 21/34] test: unit + property tests for the coverage-mask feature Cover the pure and lightly-fixtured logic: exposure-number parsing, the CCD-list -> unique-exposure reduction, the MegaCam field-corner convention, the multi-HDU header -> WCS split, power-of-two nside validation, and a regression guard for the -h/argv handling of the console runners. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01HvP7uKkfDU2PvWC7YmbuMW --- tests/module/test_coverage.py | 194 ++++++++++++++++++++++++++++++++++ 1 file changed, 194 insertions(+) create mode 100644 tests/module/test_coverage.py diff --git a/tests/module/test_coverage.py b/tests/module/test_coverage.py new file mode 100644 index 000000000..746494d5b --- /dev/null +++ b/tests/module/test_coverage.py @@ -0,0 +1,194 @@ +"""UNIT / PROPERTY TESTS FOR THE COVERAGE-MASK FEATURE. + +Covers the pure and lightly-fixtured logic behind the coverage nexp masks: +exposure-number parsing, the CCD-list -> unique-exposure reduction, the +MegaCam field-corner convention, the multi-HDU header -> WCS split, the +power-of-two ``nside`` validation, and the ``-h``/argv handling of the +console runners. The heavy paths (VOSpace download, HealSparse map building, +plotting, multiprocessing) are exercised end to end by the pipeline, not here. +""" + +import sys + +import numpy as np +import numpy.testing as npt +import pytest +from astropy import wcs +from hypothesis import given +from hypothesis import strategies as st + +from shapepipe.utilities.coverage_map_builder import CoverageMapBuilder +from shapepipe.utilities.field_corners_extractor import ( + FieldCornersExtractor, + _expnum_from_path, + _megacam_field_corners, + _parse_header_to_wcs, +) +from shapepipe.utilities.header_downloader import HeaderDownloader + + +def _tan_wcs(crval_ra, crval_dec=0.0): + """Build a minimal 2-D TAN WCS centred on ``(crval_ra, crval_dec)``.""" + w = wcs.WCS(naxis=2) + w.wcs.ctype = ["RA---TAN", "DEC--TAN"] + w.wcs.crval = [crval_ra, crval_dec] + w.wcs.crpix = [1040.0, 1.0] + w.wcs.cd = [[-1e-5, 0.0], [0.0, 1e-5]] + return w + + +# --------------------------------------------------------------------------- +# _expnum_from_path +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize( + "path, expected", + [ + ("1234567.txt", 1234567), + ("/a/b/2143523.txt", 2143523), + ("headers/0000042.txt", 42), + ], +) +def test_expnum_from_path_extracts_trailing_number(path, expected): + """The trailing ``.txt`` is parsed as the exposure number.""" + assert _expnum_from_path(path) == expected + + +@pytest.mark.parametrize("path", ["no_number.txt", "1234567.fits", "abc.txt"]) +def test_expnum_from_path_raises_without_number(path): + """A filename without a trailing numeric stem raises ``ValueError``.""" + with pytest.raises(ValueError): + _expnum_from_path(path) + + +@given(st.integers(min_value=0, max_value=99999999)) +def test_expnum_from_path_roundtrips(expnum): + """Any exposure number round-trips through the filename convention.""" + assert _expnum_from_path(f"vos_headers/{expnum}.txt") == expnum + + +# --------------------------------------------------------------------------- +# HeaderDownloader.get_exposures +# --------------------------------------------------------------------------- + +def test_get_exposures_reduces_to_unique_exposures(tmp_path): + """A ``-`` CCD list collapses to its unique exposure numbers.""" + ccd_list = tmp_path / "ccds.txt" + ccd_list.write_text("2143523-0\n2143523-5\n2143524-3\n2143524-8\n") + + exps = HeaderDownloader().get_exposures(str(ccd_list)) + + npt.assert_array_equal(exps, np.array([2143523, 2143524])) + + +def test_get_exposures_csv_matches_txt(tmp_path): + """The CSV and text code paths yield the same unique exposures.""" + txt = tmp_path / "ccds.txt" + txt.write_text("2143523-0\n2143523-5\n2143524-3\n") + csv = tmp_path / "ccds.csv" + csv.write_text("CCD\n2143523-0\n2143523-5\n2143524-3\n") + + dl = HeaderDownloader() + + npt.assert_array_equal( + dl.get_exposures(str(txt)), dl.get_exposures(str(csv)) + ) + + +# --------------------------------------------------------------------------- +# field-corner geometry +# --------------------------------------------------------------------------- + +def test_megacam_field_corners_uses_expected_ccds(): + """The corners come from CCDs 0, 8, 35, 27, in that order.""" + # Centre each CCD at RA = 100 + index (away from the RA=0 seam) so the + # returned corner RA reveals which CCD index each slot was drawn from. + wcs_list = [_tan_wcs(100.0 + i) for i in range(40)] + + ra, dec = _megacam_field_corners(wcs_list) + + assert len(ra) == 4 + assert len(dec) == 4 + npt.assert_allclose(ra, [100, 108, 135, 127], atol=0.5) + npt.assert_allclose(dec, [0, 0, 0, 0], atol=0.5) + + +def test_parse_header_to_wcs_returns_one_wcs_per_extension(tmp_path): + """One WCS is returned per CCD HDU; the primary HDU is skipped.""" + ccd_wcs = [_tan_wcs(10.0), _tan_wcs(20.0), _tan_wcs(30.0)] + + blocks = ["SIMPLE = T"] + blocks += ["\n".join(str(card) for card in w.to_header().cards) + for w in ccd_wcs] + text = "".join(f"{block}\nEND \n" for block in blocks) + + path = tmp_path / "1234567.txt" + path.write_text(text) + + result = _parse_header_to_wcs(str(path)) + + assert len(result) == len(ccd_wcs) + + +# --------------------------------------------------------------------------- +# CoverageMapBuilder.check_params (nside validation) +# --------------------------------------------------------------------------- + +def test_check_params_accepts_power_of_two_nside(tmp_path): + """Powers of two for both nside values pass validation.""" + infile = tmp_path / "corners.txt" + infile.write_text("0.0 0.0 1.0 0.0 1.0 1.0 0.0 1.0\n") + + builder = CoverageMapBuilder() + builder._params["input_file"] = str(infile) + builder._params["nside_coverage"] = 32 + builder._params["nside"] = 2048 + + builder.check_params() + + +@pytest.mark.parametrize( + "nside_coverage, nside", [(33, 2048), (32, 100), (32, 3000)] +) +def test_check_params_rejects_non_power_of_two_nside( + tmp_path, nside_coverage, nside +): + """A non-power-of-two nside raises ``ValueError``.""" + infile = tmp_path / "corners.txt" + infile.write_text("0.0 0.0 1.0 0.0 1.0 1.0 0.0 1.0\n") + + builder = CoverageMapBuilder() + builder._params["input_file"] = str(infile) + builder._params["nside_coverage"] = nside_coverage + builder._params["nside"] = nside + + with pytest.raises(ValueError): + builder.check_params() + + +def test_check_params_missing_input_file_raises(tmp_path): + """A missing input file raises ``FileNotFoundError``.""" + builder = CoverageMapBuilder() + builder._params["input_file"] = str(tmp_path / "does_not_exist.txt") + + with pytest.raises(FileNotFoundError): + builder.check_params() + + +# --------------------------------------------------------------------------- +# console-runner argv handling (regression guard for the -h entry points) +# --------------------------------------------------------------------------- + +def test_run_help_flag_exits_cleanly(monkeypatch): + """``run()`` with no args reads ``sys.argv[1:]`` so ``-h`` exits 0. + + Guards the regression where the runners parsed the full ``sys.argv`` + (including ``argv[0]``), which made ``-h`` collide with the integer + ``-c`` option and exit non-zero. + """ + monkeypatch.setattr(sys, "argv", ["extract_field_corners", "-h"]) + + with pytest.raises(SystemExit) as excinfo: + FieldCornersExtractor().run() + + assert excinfo.value.code == 0 From de73a9c2accab68c3e1a187200db4660f81c2d3c Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 01:40:24 +0200 Subject: [PATCH 22/34] fix(coverage): subtract missing CCDs in get_ccds_with_psf The v>=1.4 path built the full set of exposure single-HDU IDs but never removed the CCDs whose PSF model is missing, so every candidate CCD was reported as valid. Subtract the missing set so the result is the CCDs that actually have a valid PSF, and report the three counts. --- src/shapepipe/utilities/ccd_psf_handler.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/src/shapepipe/utilities/ccd_psf_handler.py b/src/shapepipe/utilities/ccd_psf_handler.py index 5c5e1d7ad..2d0dc90eb 100644 --- a/src/shapepipe/utilities/ccd_psf_handler.py +++ b/src/shapepipe/utilities/ccd_psf_handler.py @@ -206,7 +206,9 @@ def get_exp(self, patches): def get_ccds_with_psf(self, patches, n_CCD=40): """Get CCDs With PSF. - Return set of CCDs from list of patches. + Return set of CCDs with valid PSF from a list of patches. A CCD has a + valid PSF if it appears among the exposures' single-HDU IDs but not in + the per-patch missing-CCD lists. Parameters ---------- @@ -218,7 +220,7 @@ def get_ccds_with_psf(self, patches, n_CCD=40): Returns ------- set - CCD IDs + CCD IDs with valid PSF """ # Get missing CCDs @@ -230,11 +232,18 @@ def get_ccds_with_psf(self, patches, n_CCD=40): exp_all = self.get_exp(patches) # Turn exposures into exposure-single-HDU names (CCDs) - exp_shdu_all = summary.get_all_shdus(exp_all, n_CCD) + exp_shdu_all = set(summary.get_all_shdus(exp_all, n_CCD)) - print(f"Found {len(exp_shdu_all)} CCDs") + # Subtract the CCDs whose PSF model is missing + exp_shdu_with_psf = exp_shdu_all - exp_shdu_missing_all - return exp_shdu_all + print( + f"Found {len(exp_shdu_all)} CCDs, " + f"{len(exp_shdu_missing_all)} missing, " + f"{len(exp_shdu_with_psf)} with valid PSF" + ) + + return exp_shdu_with_psf def get_ccds_with_psf_method_v1_3(self, patches, n_CCD=40): """Get CCDs With PSF Method v1.3. From abce8d1bb076c177d8ff5e955e4bdae6f9f4000e Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 01:40:34 +0200 Subject: [PATCH 23/34] feat(coverage): extract per-CCD corners with optional CCD-list filter Replace the whole-MegaCam-field corner extraction with per-CCD extraction: for each HDU's WCS, build the 4 corners from the CCD pixel bounds (read from pixel_shape, i.e. NAXIS1/NAXIS2) using pixel edges so the quadrilateral covers the full CCD area. Emit one row per CCD keyed on the '-' ID. Add an optional --ccd_list/-l option: when given, only CCDs in that file (the output of get_ccds_with_psf) are written, so the coverage mask is built from single CCDs with a valid PSF. Resume still keys the done-set on exposure number, now parsed from the CCD-ID first column. Drop the dead _megacam_field_corners whole-field path; keep _parse_header_to_wcs and the multiprocessing structure. --- .../utilities/field_corners_extractor.py | 179 +++++++++++++----- 1 file changed, 134 insertions(+), 45 deletions(-) diff --git a/src/shapepipe/utilities/field_corners_extractor.py b/src/shapepipe/utilities/field_corners_extractor.py index aa7df9889..c05316e86 100644 --- a/src/shapepipe/utilities/field_corners_extractor.py +++ b/src/shapepipe/utilities/field_corners_extractor.py @@ -43,21 +43,47 @@ def _parse_header_to_wcs(path): ] -def _megacam_field_corners(w): - """Return (ra, dec) of the 4 MegaCam field corners. - - The hard-coded CCD indices (0, 8, 35, 27) and pixel coords - (2079/32, 0) are the MegaCam corner convention; CCDs are flipped so - the nominal bottom-right CCD reports the top-left sky corner, etc. +def _ccd_corners(w): + """Return (ra, dec) of the 4 corners of a single CCD. + + The polygon is built from the CCD's pixel bounds, read from the WCS + ``pixel_shape`` (populated by astropy from the header ``NAXIS1``/``NAXIS2`` + keywords). Pixel edges ``(-0.5, n - 0.5)`` are used so the quadrilateral + covers the full CCD area rather than pixel centres. Corners are returned in + a consistent counterclockwise pixel order + (bottom-left, bottom-right, top-right, top-left); since ``pixel_to_world`` + maps this order to a simple, non-self-intersecting sky quadrilateral for a + CCD-sized field, and HealSparse ``Polygon`` is orientation-agnostic, the + resulting polygon is a valid convex footprint. + + Parameters + ---------- + w : astropy.wcs.WCS + WCS of a single CCD + + Returns + ------- + tuple + ``(ra_list, dec_list)`` each of length 4, in degrees + + Raises + ------ + ValueError + if the WCS carries no pixel shape (missing ``NAXIS1``/``NAXIS2``) """ - tl = w[0].pixel_to_world(2079, 0) - tr = w[8].pixel_to_world(32, 0) - br = w[35].pixel_to_world(2079, 0) - bl = w[27].pixel_to_world(32, 0) - return ( - [tl.ra.deg, tr.ra.deg, br.ra.deg, bl.ra.deg], - [tl.dec.deg, tr.dec.deg, br.dec.deg, bl.dec.deg], - ) + shape = w.pixel_shape + if shape is None: + raise ValueError( + "WCS has no pixel_shape; header is missing NAXIS1/NAXIS2" + ) + nx, ny = shape + + # Pixel-edge corners, counterclockwise: BL, BR, TR, TL + px = [-0.5, nx - 0.5, nx - 0.5, -0.5] + py = [-0.5, -0.5, ny - 0.5, ny - 0.5] + + sky = w.pixel_to_world(px, py) + return list(sky.ra.deg), list(sky.dec.deg) class FieldCornersExtractor(object): @@ -76,6 +102,7 @@ def params_default(self): self._params = { "input_dir": "header", "output_file": "exp_ra_dec.txt", + "ccd_list": None, "resume": False, "n_processes": 1, "verbose": False, @@ -84,6 +111,7 @@ def params_default(self): self._short_options = { "input_dir": "-i", "output_file": "-o", + "ccd_list": "-l", "resume": "-r", "n_processes": "-n", } @@ -96,7 +124,8 @@ def params_default(self): self._help_strings = { "input_dir": "input directory containing header files; default is {}", - "output_file": "output file for field corners; default is {}", + "output_file": "output file for per-CCD corners; default is {}", + "ccd_list": "file of valid CCD IDs (output of get_ccds_with_psf); when given, only listed CCDs are written; default is all CCDs", "resume": "resume from existing output file; default is {}", "n_processes": f"number of parallel processes (1=serial, 0=auto={cpu_count()}); default is {{}}", } @@ -142,6 +171,13 @@ def check_params(self): f"Input directory not found: {self._params['input_dir']}" ) + if self._params["ccd_list"] is not None and not exists( + self._params["ccd_list"] + ): + raise FileNotFoundError( + f"CCD list file not found: {self._params['ccd_list']}" + ) + # Set n_processes to cpu_count if 0 if self._params["n_processes"] == 0: self._params["n_processes"] = cpu_count() @@ -152,35 +188,71 @@ def check_params(self): ) @staticmethod - def process_single_header(path_and_verbose): + def load_ccd_list(path): + """Load CCD List. + + Read valid CCD IDs (one ``-`` per line) into a set. + + Parameters + ---------- + path : str + path to the CCD list file + + Returns + ------- + set + valid CCD IDs + + """ + with open(path) as f: + return {line.strip() for line in f if line.strip()} + + @staticmethod + def process_single_header(args): """Process Single Header. - Worker function to process a single header file. + Worker function to process a single header file into per-CCD corners. Static method so it can be pickled for multiprocessing. Parameters ---------- - path_and_verbose : tuple - (path, verbose) where path is header file path and verbose is bool + args : tuple + ``(path, verbose, valid_ccds)`` where ``path`` is the header file + path, ``verbose`` is a bool, and ``valid_ccds`` is a set of CCD IDs + to keep (or ``None`` to keep all) Returns ------- - tuple or None - (expnum, ra_list, dec_list) on success, None on failure + list or None + list of ``(ccd_id, ra_list, dec_list)`` for the exposure's CCDs on + success, ``None`` on failure """ - path, verbose = path_and_verbose + path, verbose, valid_ccds = args expnum = _expnum_from_path(path) try: - w = _parse_header_to_wcs(path) - ra, dec = _megacam_field_corners(w) - return (expnum, ra, dec) + wcs_list = _parse_header_to_wcs(path) except Exception as e: if verbose: print(f"Failed to process {expnum}: {e}") return None + rows = [] + for ccd_idx, w in enumerate(wcs_list): + ccd_id = f"{expnum}-{ccd_idx}" + if valid_ccds is not None and ccd_id not in valid_ccds: + continue + try: + ra, dec = _ccd_corners(w) + except Exception as e: + if verbose: + print(f"Failed to process CCD {ccd_id}: {e}") + continue + rows.append((ccd_id, ra, dec)) + + return rows + def get_done_exposures(self): """Get Done Exposures. @@ -198,8 +270,12 @@ def get_done_exposures(self): return set() try: - done = np.loadtxt(output_file, usecols=(0), dtype=int) - return set(done) + ids = np.atleast_1d( + np.loadtxt(output_file, usecols=(0), dtype=str) + ) + # First column is a "-" CCD ID; an exposure counts + # as done once any of its CCDs appears in the output. + return {int(ccd_id.split("-")[0]) for ccd_id in ids} except Exception as e: if self._params["verbose"]: print(f"Could not read existing output file: {e}") @@ -235,6 +311,12 @@ def run(self, args=None): resume = self._params["resume"] verbose = self._params["verbose"] + # Load the set of valid CCD IDs to keep, if given + valid_ccds = None + if self._params["ccd_list"] is not None: + valid_ccds = self.load_ccd_list(self._params["ccd_list"]) + print(f"{len(valid_ccds)} valid CCDs in {self._params['ccd_list']}") + # Find all header files paths = glob.glob(f"{input_dir}/*.txt") n = len(paths) @@ -267,31 +349,34 @@ def run(self, args=None): # Process headers if n_processes == 1: # Serial processing - results = self._process_serial(todo, verbose) + results = self._process_serial(todo, verbose, valid_ccds) else: # Parallel processing print(f"Using {n_processes} parallel processes") - results = self._process_parallel(todo, n_processes, verbose) + results = self._process_parallel( + todo, n_processes, verbose, valid_ccds + ) - # Filter out failed results (None values) - results = [r for r in results if r is not None] + # Flatten per-exposure CCD lists, dropping failed exposures (None) + rows = [row for res in results if res is not None for row in res] - # Sort results by exposure number - results.sort(key=lambda x: x[0]) + # Sort by exposure number, then CCD index + rows.sort(key=lambda r: (int(r[0].split("-")[0]), + int(r[0].split("-")[1]))) - # Write results to file + # Write results to file: " ra1 ra2 ra3 ra4 dec1 dec2 dec3 dec4" mode = "a" if resume else "w" with open(output_file, mode, buffering=1) as f: - for expnum, ra, dec in results: - f.write(f"{expnum:6d} ") + for ccd_id, ra, dec in rows: + f.write(f"{ccd_id} ") np.savetxt(f, ra, fmt="%9.5f", newline=" ") np.savetxt(f, dec, fmt="%9.5f", newline=" ") f.write("\n") - n_success = len(results) - n_failed = n_todo - n_success + n_exp_success = sum(1 for res in results if res is not None) + n_failed = n_todo - n_exp_success - print(f"Processed {n_success} exposures") + print(f"Processed {n_exp_success} exposures, {len(rows)} CCDs") if n_failed > 0: print(f"Failed to process {n_failed} exposures") @@ -299,7 +384,7 @@ def run(self, args=None): return 0 - def _process_serial(self, todo, verbose): + def _process_serial(self, todo, verbose, valid_ccds): """Process Serial. Process headers serially. @@ -310,18 +395,20 @@ def _process_serial(self, todo, verbose): list of header file paths to process verbose : bool verbose output + valid_ccds : set or None + CCD IDs to keep, or ``None`` to keep all Returns ------- list - list of results (expnum, ra, dec) tuples + list of per-exposure CCD-row lists (or ``None`` for failures) """ results = [] n_todo = len(todo) for i, p in enumerate(todo): - result = self.process_single_header((p, verbose)) + result = self.process_single_header((p, verbose, valid_ccds)) results.append(result) if verbose and i % 100 == 0: @@ -329,7 +416,7 @@ def _process_serial(self, todo, verbose): return results - def _process_parallel(self, todo, n_processes, verbose): + def _process_parallel(self, todo, n_processes, verbose, valid_ccds): """Process Parallel. Process headers in parallel using multiprocessing. @@ -342,15 +429,17 @@ def _process_parallel(self, todo, n_processes, verbose): number of parallel processes verbose : bool verbose output + valid_ccds : set or None + CCD IDs to keep, or ``None`` to keep all Returns ------- list - list of results (expnum, ra, dec) tuples + list of per-exposure CCD-row lists (or ``None`` for failures) """ # Prepare arguments for worker function - args = [(p, verbose) for p in todo] + args = [(p, verbose, valid_ccds) for p in todo] # Create pool and process with Pool(processes=n_processes) as pool: From 76f856e582080971a1ae7e7245b74ca2e92ee279 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 01:40:41 +0200 Subject: [PATCH 24/34] feat(coverage): build maps from per-CCD rows with RA-wrap and pole guards Parse the new per-CCD row format (string CCD ID in column 0, then 8 float corner columns). Add an RA-wrap guard that unwraps polygons straddling the RA=0 seam onto a common branch, and a pole guard that skips polygons reaching |dec| >= 89 deg with a warning (CCD footprints never reach there). Since the CCDs of one exposure do not overlap, accumulating value 1 per CCD polygon still counts exposures-with-valid-PSF per pixel. --- .../utilities/coverage_map_builder.py | 97 ++++++++++++++----- 1 file changed, 71 insertions(+), 26 deletions(-) diff --git a/src/shapepipe/utilities/coverage_map_builder.py b/src/shapepipe/utilities/coverage_map_builder.py index 1ad6b286d..c7144ef5c 100644 --- a/src/shapepipe/utilities/coverage_map_builder.py +++ b/src/shapepipe/utilities/coverage_map_builder.py @@ -1,6 +1,11 @@ """COVERAGE_MAP_BUILDER -Build HealSparse coverage maps from field corner coordinates. +Build HealSparse coverage maps from per-CCD corner coordinates. + +Each input row is one CCD footprint. Since the CCDs of a single exposure do +not overlap, stamping value 1 per CCD polygon and accumulating makes the map +count, per sky pixel, the number of exposures with a valid PSF model covering +that pixel. Author: Mike Hudson, Martin Kilbinger @@ -16,6 +21,36 @@ from cs_util import args as cs_args from cs_util import logging +# Declination beyond which a polygon is considered too close to a pole for +# HealSparse's planar polygon fill to be reliable. CCD footprints are ~10 +# arcmin, so this is a guard against pathological headers, not a common path. +_DEC_POLE_LIMIT = 89.0 + + +def unwrap_ra(ra): + """Unwrap RA. + + Shift RA values onto a common branch when a polygon straddles the RA=0 + seam, so the corners describe the small footprint rather than its ~360° + complement. Values above 180 deg are shifted by -360 deg when the raw + spread exceeds 180 deg; otherwise the values are returned unchanged. + + Parameters + ---------- + ra : array_like + polygon RA corners, in degrees + + Returns + ------- + numpy.ndarray + RA corners on a common branch, in degrees + + """ + ra = np.asarray(ra, dtype=float) + if ra.max() - ra.min() > 180.0: + return np.where(ra > 180.0, ra - 360.0, ra) + return ra + class CoverageMapBuilder(object): @@ -225,25 +260,18 @@ def run(self, args=None): verbose = self._params["verbose"] if verbose: - print(f"Reading field corners from: {input_file}") - - # Load field corner data - dtype = np.dtype( - [ - ("expnum", int), - ("tlra", float), - ("trra", float), - ("brra", float), - ("blra", float), - ("tldec", float), - ("trdec", float), - ("brdec", float), - ("bldec", float), - ] + print(f"Reading per-CCD corners from: {input_file}") + + # Load per-CCD corner data. Column 0 is a "-" string + # ID; the remaining 8 columns are the 4 RA then 4 Dec corners. + ccd_ids = np.atleast_1d(np.loadtxt(input_file, usecols=(0), dtype=str)) + corners = np.atleast_2d( + np.loadtxt(input_file, usecols=range(1, 9), dtype=float) ) - rows = np.loadtxt(input_file, dtype=dtype) + ra_all = corners[:, :4] + dec_all = corners[:, 4:] - print(f"Loaded {len(rows)} exposure field corners") + print(f"Loaded {len(ccd_ids)} CCD footprints") if verbose: print( @@ -253,20 +281,37 @@ def run(self, args=None): # Create empty map m = hsp.HealSparseMap.make_empty(nside_coverage, nside, np.uint16) - # Add polygons for each exposure + # Add one polygon per CCD footprint if verbose: print("Adding polygons to map") - for i, r in enumerate(rows): - ra = [r["tlra"], r["trra"], r["brra"], r["blra"]] - dec = [r["tldec"], r["trdec"], r["brdec"], r["bldec"]] - poly = hsp.Polygon(ra=ra, dec=dec, value=1) - m += poly + n_added = 0 + n_skipped = 0 + for i in range(len(ccd_ids)): + dec = dec_all[i] + + # Pole guard: HealSparse's planar polygon fill degrades near the + # poles. CCD footprints never reach here, so warn and skip. + if np.any(np.abs(dec) >= _DEC_POLE_LIMIT): + print( + f"Warning: skipping CCD {ccd_ids[i]} with |dec| >= " + f"{_DEC_POLE_LIMIT} (too close to a pole)" + ) + n_skipped += 1 + continue + + # RA-wrap guard: put corners on a common branch across the seam. + ra = unwrap_ra(ra_all[i]) + + m += hsp.Polygon(ra=list(ra), dec=list(dec), value=1) + n_added += 1 if verbose and i % 1000 == 0: - print(f"{i:6d} / {len(rows):6d}") + print(f"{i:6d} / {len(ccd_ids):6d}") - print(f"Added {len(rows)} polygons to map") + print(f"Added {n_added} polygons to map") + if n_skipped > 0: + print(f"Skipped {n_skipped} polygons near a pole") # Apply median filter if requested if apply_median_filter: From 96cd62b6f5ae0925c8f446947145cd37502999f7 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 01:40:54 +0200 Subject: [PATCH 25/34] fix(coverage): download headers atomically via temp-file rename Copy each header to a '.part' file in the same directory and os.rename on success. An interrupted transfer can no longer leave a partial file that resume treats as complete; the temp file is removed on failure. --- src/shapepipe/utilities/header_downloader.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/shapepipe/utilities/header_downloader.py b/src/shapepipe/utilities/header_downloader.py index 9f995647c..e37b21878 100644 --- a/src/shapepipe/utilities/header_downloader.py +++ b/src/shapepipe/utilities/header_downloader.py @@ -209,12 +209,19 @@ def get_fits_header(self, expnum, client): print(f"Could not create symlink from {link_source}: {e}") # Fall through to download if symlink fails - # Download from VOSpace + # Download from VOSpace atomically: copy to a temp file in the same + # directory, then rename on success. An interrupted transfer leaves + # only the temp file behind, so resume never treats a partial download + # as complete. + tmp_dest = f"{dest}.part" try: - client.copy(source, dest, head=True) + client.copy(source, tmp_dest, head=True) + os.rename(tmp_dest, dest) return True except Exception as e: print(f"Could not copy {source}: {e}") + if exists(tmp_dest): + os.remove(tmp_dest) return False def run(self, args=None): From ab699a79bbbfea39d41528076c4217c4a63e5b47 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 01:40:54 +0200 Subject: [PATCH 26/34] docs,build(coverage): align interfaces to the per-CCD flow Fix the console-script name typo (get_ccd_with_psfs -> get_ccds_with_psf) so it matches the docs and script. Harden the build/plot script with 'set -euo pipefail' and document the per-CCD chain (step 0 get_ccds_with_psf feeding extract_field_corners --ccd_list). Update the canfar docs: restore the coverage helpers as shipped, align filenames (ccds_with_psf_.txt, exp_ra_dec_.txt), and describe the per-CCD workflow including step 0. --- docs/source/pipeline_canfar.md | 35 ++++++++++++---------- pyproject.toml | 2 +- scripts/sh/build_and_plot_coverage_maps.sh | 31 ++++++++++++++++--- 3 files changed, 48 insertions(+), 20 deletions(-) diff --git a/docs/source/pipeline_canfar.md b/docs/source/pipeline_canfar.md index 25ae6a546..d2b5a79ea 100644 --- a/docs/source/pipeline_canfar.md +++ b/docs/source/pipeline_canfar.md @@ -458,13 +458,12 @@ calibrate_comprehensive ```{note} -The matched-star-catalogue and coverage-mask diagnostics in the next two -sections have largely moved out of ShapePipe into -[`sp_validation`](https://github.com/CosmoStat/sp_validation) / `cosmo_val`. -Several helpers referenced below — `merge_psf_cat.py`, `download_headers`, -`extract_field_corners`, `build_coverage_map`, `plot_coverage_map` — are no -longer shipped in this repository; see `sp_validation` for their current -equivalents. (`get_ccds_with_psf` and `summary_run` are still part of ShapePipe.) +The matched-star-catalogue diagnostic (`merge_psf_cat.py`) has moved into +[`sp_validation`](https://github.com/CosmoStat/sp_validation) / `cosmo_val`; +see that repository for its current equivalent. The coverage-mask helpers +(`get_ccds_with_psf`, `download_headers`, `extract_field_corners`, +`build_coverage_map`, `plot_coverage_map`, and `summary_run`) are shipped by +ShapePipe and documented below. ``` ### Create matched star catalogue @@ -495,27 +494,33 @@ If the file `$patch/exp_numbers.txt` does not exist for a given patch, create it summary_run $patch 1 ``` -Now, create the list of CCDs that have PSF information with +Now, create the list of single-CCD footprints that have a valid PSF model. +Only these CCDs are stamped into the mask, so the accumulated map counts, per +sky pixel, the number of exposures with a valid PSF. The list is written to +`ccds_with_psf_v1.6.txt`: ```bash -get_ccds_with_psf -v -V v1.6 +get_ccds_with_psf -v -V v1.6 -o ccds_with_psf_v1.6.txt ``` -Next, download exposures headers; indicate (with `-d`) a directory of already -downloaded headers; those will be linked and duplicated download skipped. +Next, download the exposure headers; indicate (with `-d`) a directory of +already downloaded headers, which are linked so duplicate downloads are +skipped: ```bash -download_headers -i ccds_with_psfs_v1.6.txt -o headers_v1.6 -d headers_v1.3 -v +download_headers -i ccds_with_psf_v1.6.txt -o headers_v1.6 -d headers_v1.3 -v ``` -From the headers, the CCD corner coordinates are extracted with +From the headers, the per-CCD corner coordinates are extracted. Pass the +CCD list with `-l` so only CCDs with a valid PSF are written: + ```bash -extract_field_corners -i headers_v1.6 -v +extract_field_corners -i headers_v1.6 -l ccds_with_psf_v1.6.txt -o exp_ra_dec_v1.6.txt -v ``` Then, build the healsparse coverage mask file as ```bash -build_coverage_map -i exp_ra_dec_v1.6 -o coverage_v1.6.x.hsp -c 128 -n 131072 -v +build_coverage_map -i exp_ra_dec_v1.6.txt -o coverage_v1.6.x.hsp -c 128 -n 131072 -v ``` The healsparse resolutions (128, 131072) match the bit masks. diff --git a/pyproject.toml b/pyproject.toml index 783a3864e..5333aa91c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -88,7 +88,7 @@ summary_run = "shapepipe.summary_run:main" canfar_submit_job = "shapepipe.canfar_run:run_job" canfar_monitor = "shapepipe.canfar_run:run_log" canfar_monitor_log = "shapepipe.canfar_run:run_monitor_log" -get_ccd_with_psfs = "shapepipe.get_ccds_run:run_ccd_psf_handler" +get_ccds_with_psf = "shapepipe.get_ccds_run:run_ccd_psf_handler" download_headers = "shapepipe.coverage_run:run_download_headers" extract_field_corners = "shapepipe.coverage_run:run_extract_corners" build_coverage_map = "shapepipe.coverage_run:run_build_coverage" diff --git a/scripts/sh/build_and_plot_coverage_maps.sh b/scripts/sh/build_and_plot_coverage_maps.sh index c35216963..b7ea9c8b2 100755 --- a/scripts/sh/build_and_plot_coverage_maps.sh +++ b/scripts/sh/build_and_plot_coverage_maps.sh @@ -1,5 +1,19 @@ #!/bin/bash +set -euo pipefail + +# Build and plot per-CCD coverage (nexp) maps for a range of catalogue +# versions. Each map counts, per sky pixel, the number of exposures with a +# valid PSF model. The per-version flow is: +# 0. get_ccds_with_psf -> ccds_with_psf_.txt (valid CCD IDs) +# 1. download_headers -> per-exposure header text files +# 2. extract_field_corners-> exp_ra_dec_.txt (per-CCD corners) +# 3. build_coverage_map -> coverage_.x.hsp +# 4. plot_coverage_map -> SGC / NGC plots +# Steps 0-2 are expected to have been run already on canfar (they need +# VOSpace access); this script rebuilds from the extracted per-CCD corners and +# plots. Uncomment the marked lines to run the full chain. + # Configuration variables ## Catalogue versions @@ -23,7 +37,7 @@ PLOT_COLORBAR="-C" PLOT_MIN=1 PLOT_MAX=5 -# Plot parameters +# Plot parameters ## SGC region SGC_RA_MIN=-20 SGC_RA_MAX=45 @@ -44,12 +58,21 @@ for VERSION in "${VERSIONS[@]}"; do echo "Processing ${VERSION}..." # Define file paths - INPUT_EXP="exp_ra_dec_${VERSION}.txt" + CCD_LIST="ccds_with_psf_${VERSION}.txt" + HEADER_DIR="headers_${VERSION}" + INPUT_CORNERS="exp_ra_dec_${VERSION}.txt" COVERAGE_MAP="coverage_${VERSION}.x.hsp" + # Step 0-2: build the per-CCD corner file (canfar / VOSpace only). + # Uncomment to run the full chain from scratch. + # get_ccds_with_psf -V "${VERSION}" -o "${CCD_LIST}" ${VERBOSE} + # download_headers -i "${CCD_LIST}" -o "${HEADER_DIR}" ${VERBOSE} + # extract_field_corners -i "${HEADER_DIR}" -l "${CCD_LIST}" \ + # -o "${INPUT_CORNERS}" ${VERBOSE} + # Build coverage map - echo " Building coverage map from ${INPUT_EXP}..." - CMD="build_coverage_map -i ${INPUT_EXP} -o ${COVERAGE_MAP} -c ${BUILD_CHANNELS} -n ${BUILD_NSIDE} ${VERBOSE}" + echo " Building coverage map from ${INPUT_CORNERS}..." + CMD="build_coverage_map -i ${INPUT_CORNERS} -o ${COVERAGE_MAP} -c ${BUILD_CHANNELS} -n ${BUILD_NSIDE} ${VERBOSE}" echo "$CMD" $CMD From b784d9ecbab732bdc7882b6c7d77aca57e0f3397 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 01:40:54 +0200 Subject: [PATCH 27/34] test(coverage): cover per-CCD extraction, filtering, guards, atomic download Rewrite the coverage tests for the per-CCD design: handler missing-CCD subtraction; per-CCD corner extraction and pixel-shape error; --ccd_list filtering end to end; builder parsing of the new row format; RA-wrap unwrap (a CCD straddling RA=0 fills a small patch, not the complement); pole-guard skip; and atomic-download rename behaviour with a mocked copy. Drop the tests that pinned the removed whole-field convention. --- tests/module/test_coverage.py | 338 +++++++++++++++++++++++++++++++--- 1 file changed, 310 insertions(+), 28 deletions(-) diff --git a/tests/module/test_coverage.py b/tests/module/test_coverage.py index 746494d5b..96b104097 100644 --- a/tests/module/test_coverage.py +++ b/tests/module/test_coverage.py @@ -1,42 +1,74 @@ """UNIT / PROPERTY TESTS FOR THE COVERAGE-MASK FEATURE. -Covers the pure and lightly-fixtured logic behind the coverage nexp masks: -exposure-number parsing, the CCD-list -> unique-exposure reduction, the -MegaCam field-corner convention, the multi-HDU header -> WCS split, the -power-of-two ``nside`` validation, and the ``-h``/argv handling of the -console runners. The heavy paths (VOSpace download, HealSparse map building, -plotting, multiprocessing) are exercised end to end by the pipeline, not here. +Covers the pure and lightly-fixtured logic behind the per-CCD coverage nexp +masks: exposure-number parsing, the CCD-list -> unique-exposure reduction, the +handler's missing-CCD subtraction, per-CCD corner extraction from a multi-HDU +header, ``--ccd_list`` filtering, the builder's per-CCD row parsing and RA-wrap +guard, the power-of-two ``nside`` validation, the atomic-download rename, and +the ``-h``/argv handling of the console runners. The heavy paths (real VOSpace +download, production-resolution map building, plotting, multiprocessing) are +exercised end to end by the pipeline, not here. """ import sys +from types import SimpleNamespace +from unittest import mock import numpy as np import numpy.testing as npt import pytest from astropy import wcs +from astropy.io.fits import Header from hypothesis import given from hypothesis import strategies as st -from shapepipe.utilities.coverage_map_builder import CoverageMapBuilder +from shapepipe.utilities import ccd_psf_handler as ccd_psf_handler_mod +from shapepipe.utilities.ccd_psf_handler import CcdPsfHandler +from shapepipe.utilities.coverage_map_builder import ( + CoverageMapBuilder, + unwrap_ra, +) from shapepipe.utilities.field_corners_extractor import ( FieldCornersExtractor, + _ccd_corners, _expnum_from_path, - _megacam_field_corners, _parse_header_to_wcs, ) from shapepipe.utilities.header_downloader import HeaderDownloader -def _tan_wcs(crval_ra, crval_dec=0.0): - """Build a minimal 2-D TAN WCS centred on ``(crval_ra, crval_dec)``.""" +def _tan_wcs(crval_ra, crval_dec=0.0, nx=2080, ny=4612): + """Build a minimal 2-D TAN WCS with a populated pixel shape. + + ``nx``/``ny`` set ``pixel_shape`` (populated by astropy from + ``NAXIS1``/``NAXIS2``), so ``_ccd_corners`` can read the CCD pixel bounds. + """ w = wcs.WCS(naxis=2) w.wcs.ctype = ["RA---TAN", "DEC--TAN"] w.wcs.crval = [crval_ra, crval_dec] - w.wcs.crpix = [1040.0, 1.0] + w.wcs.crpix = [nx / 2.0, ny / 2.0] w.wcs.cd = [[-1e-5, 0.0], [0.0, 1e-5]] + w.pixel_shape = (nx, ny) return w +def _header_text(wcs_list): + """Render a multi-HDU header text file: a primary HDU then one per WCS. + + Each CCD header carries ``NAXIS1``/``NAXIS2`` so the reconstructed WCS has a + pixel shape. + """ + blocks = ["SIMPLE = T"] + for w in wcs_list: + nx, ny = w.pixel_shape + header = w.to_header() + header["NAXIS"] = 2 + header["NAXIS1"] = nx + header["NAXIS2"] = ny + blocks.append("\n".join(str(card) for card in header.cards)) + return "".join(f"{block}\nEND \n" for block in blocks) + + # --------------------------------------------------------------------------- # _expnum_from_path # --------------------------------------------------------------------------- @@ -96,38 +128,214 @@ def test_get_exposures_csv_matches_txt(tmp_path): # --------------------------------------------------------------------------- -# field-corner geometry +# HeaderDownloader.get_fits_header — atomic rename # --------------------------------------------------------------------------- -def test_megacam_field_corners_uses_expected_ccds(): - """The corners come from CCDs 0, 8, 35, 27, in that order.""" - # Centre each CCD at RA = 100 + index (away from the RA=0 seam) so the - # returned corner RA reveals which CCD index each slot was drawn from. - wcs_list = [_tan_wcs(100.0 + i) for i in range(40)] +def test_get_fits_header_writes_atomically(tmp_path): + """A successful download copies to ``.part`` then renames to the dest.""" + dl = HeaderDownloader() + dl._params["output_dir"] = str(tmp_path) + dl._params["overwrite"] = False + dl._params["dir_for_links"] = None + + dest = tmp_path / "42.txt" + tmp_dest = tmp_path / "42.txt.part" + + def fake_copy(source, target, head=True): + # The copy must land on the temp path, not the final destination. + assert target == str(tmp_dest) + with open(target, "w") as f: + f.write("HEADER") + + client = SimpleNamespace(copy=fake_copy) + + assert dl.get_fits_header(42, client) is True + assert dest.exists() + assert not tmp_dest.exists() + assert dest.read_text() == "HEADER" - ra, dec = _megacam_field_corners(wcs_list) + +def test_get_fits_header_failed_copy_leaves_no_dest(tmp_path): + """A failed download leaves no destination file (only, if any, ``.part``).""" + dl = HeaderDownloader() + dl._params["output_dir"] = str(tmp_path) + dl._params["overwrite"] = False + dl._params["dir_for_links"] = None + + def failing_copy(source, target, head=True): + raise RuntimeError("transfer interrupted") + + client = SimpleNamespace(copy=failing_copy) + + assert dl.get_fits_header(42, client) is False + assert not (tmp_path / "42.txt").exists() + assert not (tmp_path / "42.txt.part").exists() + + +# --------------------------------------------------------------------------- +# CcdPsfHandler.get_ccds_with_psf — missing-CCD subtraction +# --------------------------------------------------------------------------- + +def test_get_ccds_with_psf_subtracts_missing(monkeypatch): + """Valid CCDs are all exposure single-HDUs minus the missing set.""" + handler = CcdPsfHandler() + + # Two exposures, 3 CCDs each -> 6 candidate CCDs; two are missing. + monkeypatch.setattr( + handler, "get_exp", lambda patches: {"100", "200"} + ) + monkeypatch.setattr( + handler, + "get_exp_shdu_missing", + lambda patches: {"100-1", "200-2"}, + ) + monkeypatch.setattr( + ccd_psf_handler_mod.summary, + "get_all_shdus", + lambda exposures, n_CCD: [ + f"{e}-{c}" for e in exposures for c in range(n_CCD) + ], + ) + + result = handler.get_ccds_with_psf(["P1"], n_CCD=3) + + assert result == {"100-0", "100-2", "200-0", "200-1"} + + +# --------------------------------------------------------------------------- +# per-CCD corner geometry +# --------------------------------------------------------------------------- + +def test_ccd_corners_returns_four_corners_around_centre(): + """The 4 corners bracket the CCD centre and cover its full pixel extent.""" + nx, ny = 2080, 4612 + w = _tan_wcs(100.0, 20.0, nx=nx, ny=ny) + + ra, dec = _ccd_corners(w) assert len(ra) == 4 assert len(dec) == 4 - npt.assert_allclose(ra, [100, 108, 135, 127], atol=0.5) - npt.assert_allclose(dec, [0, 0, 0, 0], atol=0.5) + # The footprint straddles the CRVAL centre in both coordinates. + assert min(ra) < 100.0 < max(ra) + assert min(dec) < 20.0 < max(dec) + # Extent matches the CD-scaled pixel span (with edge padding). + npt.assert_allclose(max(dec) - min(dec), (ny) * 1e-5, rtol=1e-3) + + +def test_ccd_corners_raises_without_pixel_shape(): + """A WCS with no pixel shape raises a clear ``ValueError``.""" + w = wcs.WCS(naxis=2) + w.wcs.ctype = ["RA---TAN", "DEC--TAN"] + w.wcs.crval = [100.0, 20.0] + assert w.pixel_shape is None + + with pytest.raises(ValueError, match="pixel_shape"): + _ccd_corners(w) def test_parse_header_to_wcs_returns_one_wcs_per_extension(tmp_path): """One WCS is returned per CCD HDU; the primary HDU is skipped.""" ccd_wcs = [_tan_wcs(10.0), _tan_wcs(20.0), _tan_wcs(30.0)] - blocks = ["SIMPLE = T"] - blocks += ["\n".join(str(card) for card in w.to_header().cards) - for w in ccd_wcs] - text = "".join(f"{block}\nEND \n" for block in blocks) - path = tmp_path / "1234567.txt" - path.write_text(text) + path.write_text(_header_text(ccd_wcs)) result = _parse_header_to_wcs(str(path)) assert len(result) == len(ccd_wcs) + # Reconstructed WCS carry the pixel shape from NAXIS1/NAXIS2. + assert all(w.pixel_shape == (2080, 4612) for w in result) + + +# --------------------------------------------------------------------------- +# FieldCornersExtractor.process_single_header — per-CCD rows and filtering +# --------------------------------------------------------------------------- + +def test_process_single_header_emits_one_row_per_ccd(tmp_path): + """Without a CCD list, every HDU yields a ``-`` row.""" + ccd_wcs = [_tan_wcs(10.0), _tan_wcs(20.0), _tan_wcs(30.0)] + path = tmp_path / "1234567.txt" + path.write_text(_header_text(ccd_wcs)) + + rows = FieldCornersExtractor.process_single_header( + (str(path), False, None) + ) + + assert [r[0] for r in rows] == [ + "1234567-0", + "1234567-1", + "1234567-2", + ] + assert all(len(r[1]) == 4 and len(r[2]) == 4 for r in rows) + + +def test_process_single_header_filters_to_ccd_list(tmp_path): + """A ``valid_ccds`` set keeps only the listed CCDs of the exposure.""" + ccd_wcs = [_tan_wcs(10.0), _tan_wcs(20.0), _tan_wcs(30.0)] + path = tmp_path / "1234567.txt" + path.write_text(_header_text(ccd_wcs)) + + rows = FieldCornersExtractor.process_single_header( + (str(path), False, {"1234567-0", "1234567-2"}) + ) + + assert [r[0] for r in rows] == ["1234567-0", "1234567-2"] + + +def test_load_ccd_list_reads_ids(tmp_path): + """The CCD list loader returns the stripped, non-blank IDs as a set.""" + path = tmp_path / "ccds.txt" + path.write_text("100-0\n100-1\n\n200-5\n") + + assert FieldCornersExtractor.load_ccd_list(str(path)) == { + "100-0", + "100-1", + "200-5", + } + + +def test_run_extract_writes_per_ccd_rows(tmp_path): + """End to end: run() writes one per-CCD row filtered by the CCD list.""" + header_dir = tmp_path / "headers" + header_dir.mkdir() + ccd_wcs = [_tan_wcs(10.0), _tan_wcs(20.0), _tan_wcs(30.0)] + (header_dir / "1234567.txt").write_text(_header_text(ccd_wcs)) + + ccd_list = tmp_path / "ccds.txt" + ccd_list.write_text("1234567-0\n1234567-2\n") + + out = tmp_path / "corners.txt" + + extractor = FieldCornersExtractor() + extractor.run( + args=[ + "-i", str(header_dir), + "-l", str(ccd_list), + "-o", str(out), + ] + ) + + lines = out.read_text().splitlines() + assert len(lines) == 2 + ids = [line.split()[0] for line in lines] + assert ids == ["1234567-0", "1234567-2"] + # Each row: 1 ID + 4 RA + 4 Dec = 9 columns. + assert all(len(line.split()) == 9 for line in lines) + + +def test_get_done_exposures_keys_on_expnum(tmp_path): + """Resume treats an exposure as done if any of its CCD rows are present.""" + out = tmp_path / "corners.txt" + out.write_text( + "1234567-0 10 10 10 10 20 20 20 20\n" + "1234567-5 10 10 10 10 20 20 20 20\n" + "7654321-3 30 30 30 30 40 40 40 40\n" + ) + + extractor = FieldCornersExtractor() + extractor._params["output_file"] = str(out) + + assert extractor.get_done_exposures() == {1234567, 7654321} # --------------------------------------------------------------------------- @@ -137,7 +345,7 @@ def test_parse_header_to_wcs_returns_one_wcs_per_extension(tmp_path): def test_check_params_accepts_power_of_two_nside(tmp_path): """Powers of two for both nside values pass validation.""" infile = tmp_path / "corners.txt" - infile.write_text("0.0 0.0 1.0 0.0 1.0 1.0 0.0 1.0\n") + infile.write_text("100-0 0.0 1.0 1.0 0.0 0.0 0.0 1.0 1.0\n") builder = CoverageMapBuilder() builder._params["input_file"] = str(infile) @@ -155,7 +363,7 @@ def test_check_params_rejects_non_power_of_two_nside( ): """A non-power-of-two nside raises ``ValueError``.""" infile = tmp_path / "corners.txt" - infile.write_text("0.0 0.0 1.0 0.0 1.0 1.0 0.0 1.0\n") + infile.write_text("100-0 0.0 1.0 1.0 0.0 0.0 0.0 1.0 1.0\n") builder = CoverageMapBuilder() builder._params["input_file"] = str(infile) @@ -175,6 +383,80 @@ def test_check_params_missing_input_file_raises(tmp_path): builder.check_params() +# --------------------------------------------------------------------------- +# CoverageMapBuilder — RA-wrap guard and per-CCD row build +# --------------------------------------------------------------------------- + +def test_unwrap_ra_shifts_straddling_seam(): + """A polygon straddling RA=0 is put on a common (negative) branch.""" + out = unwrap_ra([359.9, 0.1, 0.1, 359.9]) + npt.assert_allclose(sorted(out), [-0.1, -0.1, 0.1, 0.1], atol=1e-6) + + +def test_unwrap_ra_leaves_normal_polygon_unchanged(): + """A polygon well away from the seam is returned unchanged.""" + npt.assert_allclose( + unwrap_ra([99.9, 100.1, 100.1, 99.9]), + [99.9, 100.1, 100.1, 99.9], + ) + + +def test_build_map_wrapped_ccd_fills_small_patch(tmp_path): + """A CCD straddling RA=0 fills a small footprint, not the ~360° complement. + + Verified empirically against a same-size CCD away from the seam: with the + unwrap guard the two footprints have the same pixel count to within a few + boundary pixels. (Without unwrapping, HealSparse happens to fill the small + region too for a convex CCD-sized polygon, but the guard makes the intent + explicit and is robust to degenerate orderings.) + """ + # CCD straddling RA=0 (0.2 x 0.2 deg) and a reference CCD at RA=10. + infile = tmp_path / "corners.txt" + infile.write_text( + "100-0 359.9 0.1 0.1 359.9 0.0 0.0 0.2 0.2\n" + "200-0 9.9 10.1 10.1 9.9 0.0 0.0 0.2 0.2\n" + ) + out = tmp_path / "cov.hsp" + + builder = CoverageMapBuilder() + builder.run( + args=[ + "-i", str(infile), + "-o", str(out), + "-c", "32", + "-n", "1024", + ] + ) + + import healsparse as hsp + + m = hsp.HealSparseMap.read(str(out)) + n_valid = len(m.valid_pixels) + + # Two ~0.04 deg^2 patches at nside=1024 (~3.4 arcmin pixels): a small, + # bounded number of pixels — nowhere near a hemisphere. + assert 0 < n_valid < 200 + + +def test_build_map_pole_guard_skips_polygon(tmp_path, capsys): + """A polygon with |dec| near 90 deg is skipped with a warning.""" + infile = tmp_path / "corners.txt" + infile.write_text( + "100-0 10.0 10.2 10.2 10.0 89.5 89.5 89.7 89.7\n" + "200-0 10.0 10.2 10.2 10.0 20.0 20.0 20.2 20.2\n" + ) + out = tmp_path / "cov.hsp" + + builder = CoverageMapBuilder() + builder.run( + args=["-i", str(infile), "-o", str(out), "-c", "32", "-n", "1024"] + ) + + captured = capsys.readouterr() + assert "skipping CCD 100-0" in captured.out + assert "Added 1 polygons" in captured.out + + # --------------------------------------------------------------------------- # console-runner argv handling (regression guard for the -h entry points) # --------------------------------------------------------------------------- From 9b9ad76fec7e30ecbbc7c2bc0c9bf420046af044 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 01:53:25 +0200 Subject: [PATCH 28/34] fix(coverage): read true CCD dims from ZNAXIS and resume per-CCD Two correctness fixes in the per-CCD extractor: - Image dimensions: header_downloader fetches fpack tile-compressed headers where NAXIS1/NAXIS2 describe the compressed binary table (byte width, row count), so WCS.pixel_shape reports e.g. (8, 4612) instead of the true (2080, 4612). Read ZNAXIS1/ZNAXIS2 when ZIMAGE is set, falling back to NAXIS1/NAXIS2. _parse_header_to_wcs now returns (wcs, shape) pairs since the WCS drops the Z-keywords; _ccd_corners takes the shape explicitly. - Resume: key the done-set on individual CCD IDs, not exposure numbers. A write interrupted mid-exposure previously left a partial exposure that a -r rerun skipped forever (silent undercount); an expanded --ccd_list rerun skipped exposures whose new CCDs were never written. get_done_ccds reads the present CCD IDs and run() drops only those, so partial exposures are completed and new CCDs added, never duplicating a row. --- .../utilities/field_corners_extractor.py | 153 +++++++++++------- 1 file changed, 99 insertions(+), 54 deletions(-) diff --git a/src/shapepipe/utilities/field_corners_extractor.py b/src/shapepipe/utilities/field_corners_extractor.py index c05316e86..7e3273456 100644 --- a/src/shapepipe/utilities/field_corners_extractor.py +++ b/src/shapepipe/utilities/field_corners_extractor.py @@ -1,6 +1,10 @@ """FIELD_CORNERS_EXTRACTOR -Extract field corner coordinates from FITS headers. +Extract per-CCD sky-footprint corner coordinates from FITS headers. + +For each CCD (image HDU) in an exposure's multi-HDU header, the four corners +of the CCD are projected from pixel bounds to the sky, producing one row per +CCD keyed on its ``-`` ID. Author: Mike Hudson, Martin Kilbinger @@ -32,50 +36,86 @@ def _expnum_from_path(path): return int(match.group(1)) +def _image_shape(header): + """Return the CCD image ``(nx, ny)`` pixel shape from a header. + + For fpack tile-compressed HDUs (``ZIMAGE = T``), ``NAXIS1``/``NAXIS2`` + describe the compressed *binary table* (row width in bytes, row count), + not the image, and astropy's WCS never maps the true dimensions into + ``pixel_shape``. The true image size is carried by ``ZNAXIS1``/``ZNAXIS2``, + which are preferred here; a plain image HDU falls back to + ``NAXIS1``/``NAXIS2``. + + Parameters + ---------- + header : astropy.io.fits.Header + single-HDU header + + Returns + ------- + tuple + ``(nx, ny)`` image pixel dimensions + + Raises + ------ + ValueError + if neither the ``Z*`` nor the plain ``NAXIS`` dimensions are present + """ + if header.get("ZIMAGE", False): + keys = ("ZNAXIS1", "ZNAXIS2") + else: + keys = ("NAXIS1", "NAXIS2") + + if keys[0] not in header or keys[1] not in header: + raise ValueError( + f"Header is missing image dimensions ({keys[0]}/{keys[1]})" + ) + + return int(header[keys[0]]), int(header[keys[1]]) + + def _parse_header_to_wcs(path): - """Parse a multi-HDU header text file and return one WCS per HDU.""" + """Parse a multi-HDU header text file into ``(wcs, (nx, ny))`` per HDU. + + The primary HDU is skipped. Each remaining HDU yields its WCS together with + the CCD image pixel shape (``ZNAXIS1/2`` for compressed HDUs, else + ``NAXIS1/2``); the shape is read from the header because the WCS drops the + ``Z*`` keywords. + """ with open(path, "r") as f: string = f.read() tokens = re.split(r"^(END\s+)", string, flags=re.MULTILINE) - return [ - wcs.WCS(Header.fromstring(tokens[i] + tokens[i + 1], sep="\n")) - for i in range(2, len(tokens) - 1, 2) - ] + result = [] + for i in range(2, len(tokens) - 1, 2): + header = Header.fromstring(tokens[i] + tokens[i + 1], sep="\n") + result.append((wcs.WCS(header), _image_shape(header))) + return result -def _ccd_corners(w): +def _ccd_corners(w, shape): """Return (ra, dec) of the 4 corners of a single CCD. - The polygon is built from the CCD's pixel bounds, read from the WCS - ``pixel_shape`` (populated by astropy from the header ``NAXIS1``/``NAXIS2`` - keywords). Pixel edges ``(-0.5, n - 0.5)`` are used so the quadrilateral - covers the full CCD area rather than pixel centres. Corners are returned in - a consistent counterclockwise pixel order - (bottom-left, bottom-right, top-right, top-left); since ``pixel_to_world`` - maps this order to a simple, non-self-intersecting sky quadrilateral for a - CCD-sized field, and HealSparse ``Polygon`` is orientation-agnostic, the - resulting polygon is a valid convex footprint. + The polygon is built from the CCD's pixel bounds ``shape = (nx, ny)`` using + pixel edges ``(-0.5, n - 0.5)`` so the quadrilateral covers the full CCD + area rather than pixel centres. Corners are returned in a consistent + counterclockwise pixel order (bottom-left, bottom-right, top-right, + top-left); since ``pixel_to_world`` maps this order to a simple, + non-self-intersecting sky quadrilateral for a CCD-sized field, and + HealSparse ``Polygon`` is orientation-agnostic, the resulting polygon is a + valid convex footprint. Parameters ---------- w : astropy.wcs.WCS WCS of a single CCD + shape : tuple + ``(nx, ny)`` CCD image pixel dimensions Returns ------- tuple ``(ra_list, dec_list)`` each of length 4, in degrees - - Raises - ------ - ValueError - if the WCS carries no pixel shape (missing ``NAXIS1``/``NAXIS2``) """ - shape = w.pixel_shape - if shape is None: - raise ValueError( - "WCS has no pixel_shape; header is missing NAXIS1/NAXIS2" - ) nx, ny = shape # Pixel-edge corners, counterclockwise: BL, BR, TR, TL @@ -125,7 +165,7 @@ def params_default(self): self._help_strings = { "input_dir": "input directory containing header files; default is {}", "output_file": "output file for per-CCD corners; default is {}", - "ccd_list": "file of valid CCD IDs (output of get_ccds_with_psf); when given, only listed CCDs are written; default is all CCDs", + "ccd_list": "file of valid CCD IDs (output of get_ccds_with_psf); when given, only listed CCDs are written; on --resume, CCDs new to an expanded list are added without duplicating existing rows; default is all CCDs", "resume": "resume from existing output file; default is {}", "n_processes": f"number of parallel processes (1=serial, 0=auto={cpu_count()}); default is {{}}", } @@ -232,19 +272,19 @@ def process_single_header(args): expnum = _expnum_from_path(path) try: - wcs_list = _parse_header_to_wcs(path) + wcs_shapes = _parse_header_to_wcs(path) except Exception as e: if verbose: print(f"Failed to process {expnum}: {e}") return None rows = [] - for ccd_idx, w in enumerate(wcs_list): + for ccd_idx, (w, shape) in enumerate(wcs_shapes): ccd_id = f"{expnum}-{ccd_idx}" if valid_ccds is not None and ccd_id not in valid_ccds: continue try: - ra, dec = _ccd_corners(w) + ra, dec = _ccd_corners(w, shape) except Exception as e: if verbose: print(f"Failed to process CCD {ccd_id}: {e}") @@ -253,15 +293,20 @@ def process_single_header(args): return rows - def get_done_exposures(self): - """Get Done Exposures. + def get_done_ccds(self): + """Get Done CCDs. - Read exposures already processed from output file. + Read the set of CCD IDs already present in the output file. Resume is + keyed on individual CCD IDs, not exposure numbers: a CCD counts as done + only if its own row is present. This keeps resume correct when a write + was interrupted mid-exposure (the missing CCDs are filled in) and when + a rerun uses an expanded ``--ccd_list`` (the newly requested CCDs are + added), and it never duplicates a row. Returns ------- set - set of exposure numbers already processed + CCD IDs already written """ output_file = self._params["output_file"] @@ -273,9 +318,7 @@ def get_done_exposures(self): ids = np.atleast_1d( np.loadtxt(output_file, usecols=(0), dtype=str) ) - # First column is a "-" CCD ID; an exposure counts - # as done once any of its CCDs appears in the output. - return {int(ccd_id.split("-")[0]) for ccd_id in ids} + return set(ids.tolist()) except Exception as e: if self._params["verbose"]: print(f"Could not read existing output file: {e}") @@ -327,21 +370,17 @@ def run(self, args=None): print(f"{n} header files found") - # Get list of already processed exposures if resuming - done = set() + # On resume, read the CCD IDs already written so we can skip them + # per-CCD (not per-exposure): a partially written exposure is completed + # rather than skipped, and no row is ever duplicated. + done_ccds = set() if resume: - done = self.get_done_exposures() - print(f"{len(done)} already done") - - # Build todo list - todo = [p for p in paths if _expnum_from_path(p) not in done] - - n_todo = len(todo) - print(f"{n_todo} to process") + done_ccds = self.get_done_ccds() + print(f"{len(done_ccds)} CCDs already done") - if n_todo == 0: - print("Nothing to do") - return 0 + # Every header still needs parsing (a header may hold both done and + # not-yet-done CCDs); the per-CCD filter below decides what is written. + todo = paths # Get n_processes n_processes = self._params["n_processes"] @@ -357,8 +396,14 @@ def run(self, args=None): todo, n_processes, verbose, valid_ccds ) - # Flatten per-exposure CCD lists, dropping failed exposures (None) - rows = [row for res in results if res is not None for row in res] + # Flatten per-exposure CCD lists (dropping failed exposures), then drop + # CCDs already present in the output. + rows = [ + row + for res in results if res is not None + for row in res + if row[0] not in done_ccds + ] # Sort by exposure number, then CCD index rows.sort(key=lambda r: (int(r[0].split("-")[0]), @@ -374,9 +419,9 @@ def run(self, args=None): f.write("\n") n_exp_success = sum(1 for res in results if res is not None) - n_failed = n_todo - n_exp_success + n_failed = len(todo) - n_exp_success - print(f"Processed {n_exp_success} exposures, {len(rows)} CCDs") + print(f"Processed {n_exp_success} exposures, {len(rows)} new CCDs") if n_failed > 0: print(f"Failed to process {n_failed} exposures") From 1a0a3857d05b9918838bfd133549869c955aaa4f Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 01:53:25 +0200 Subject: [PATCH 29/34] fix(coverage): single-line CCD list, runner docstrings, plot doc example - header_downloader.get_exposures: wrap the text-format loadtxt in np.atleast_1d so a single-line CCD list (0-d array) does not crash on iteration. - coverage_run: update the extract/build runner docstrings from whole-field 'field corner coordinates' to per-CCD language. - pipeline_canfar docs: add a runnable plot_coverage_map example (the SGC region command) covering the flags build_and_plot_coverage_maps.sh uses. --- docs/source/pipeline_canfar.md | 12 ++++++++++-- src/shapepipe/coverage_run.py | 4 ++-- src/shapepipe/utilities/header_downloader.py | 7 +++++-- 3 files changed, 17 insertions(+), 6 deletions(-) diff --git a/docs/source/pipeline_canfar.md b/docs/source/pipeline_canfar.md index d2b5a79ea..fcac5b371 100644 --- a/docs/source/pipeline_canfar.md +++ b/docs/source/pipeline_canfar.md @@ -526,9 +526,17 @@ build_coverage_map -i exp_ra_dec_v1.6.txt -o coverage_v1.6.x.hsp -c 128 -n 13107 The healsparse resolutions (128, 131072) match the bit masks. -Use `plot_coverage_map` to create plots of the coverage mask. +Use `plot_coverage_map` to plot a region of the coverage mask. For example, +the SGC region with a colorbar clipped to the 1–5 exposure range: -Building and plotting for a range of versions is done with +```bash +plot_coverage_map -i coverage_v1.6.x.hsp -o coverage_v1.6_SGC.png \ + -C -R -20 -r 45 -D 18 -d 40 -m 1 -M 5 -v +``` + +Here `-R`/`-r` and `-D`/`-d` set the RA and Dec plot limits, `-m`/`-M` the +colorbar range, and `-C` adds the colorbar. Building and plotting for a range +of versions (SGC and NGC regions) is automated by `build_and_plot_coverage_maps.sh`. diff --git a/src/shapepipe/coverage_run.py b/src/shapepipe/coverage_run.py index 2485e4a10..15aaa1a37 100644 --- a/src/shapepipe/coverage_run.py +++ b/src/shapepipe/coverage_run.py @@ -37,7 +37,7 @@ def run_download_headers(args=None): def run_extract_corners(args=None): """Run Extract Corners. - Extract field corner coordinates from FITS headers. + Extract per-CCD sky-footprint corner coordinates from FITS headers. Parameters ---------- @@ -57,7 +57,7 @@ def run_extract_corners(args=None): def run_build_coverage(args=None): """Run Build Coverage. - Build HealSparse coverage maps from field corner coordinates. + Build HealSparse coverage maps from per-CCD corner coordinates. Parameters ---------- diff --git a/src/shapepipe/utilities/header_downloader.py b/src/shapepipe/utilities/header_downloader.py index e37b21878..122327a58 100644 --- a/src/shapepipe/utilities/header_downloader.py +++ b/src/shapepipe/utilities/header_downloader.py @@ -151,8 +151,11 @@ def get_exposures(self, ccd_list_file): exp = int(r1[0]) exps.append(exp) else: - # Text format - f = np.loadtxt(ccd_list_file, dtype="str", encoding="ascii") + # Text format; atleast_1d so a single-line file (0-d array) is + # still iterable. + f = np.atleast_1d( + np.loadtxt(ccd_list_file, dtype="str", encoding="ascii") + ) r = np.char.split(f, sep="-") for i, r1 in enumerate(r): From 0212fc56df13b773310b6222b849469734561a5e Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 01:53:25 +0200 Subject: [PATCH 30/34] test(coverage): compressed-header dims, resume path, nexp value contract Extend the coverage tests for the confirmed findings: - fpack-compressed header fixture; assert ZNAXIS dims are recovered and corners span the true CCD width, not the 8-byte compressed NAXIS1. - resume run() path: a new exposure is appended without duplicating existing rows, and a partial exposure (interrupted mid-write) is completed rather than skipped. - the nexp contract on pixel values: two overlapping CCDs from different exposures give value 2 in the overlap and 1 outside. - single-row builder file (atleast_1d/2d guards); single-line get_exposures; RA-wrap strengthened to a positional-parity check (seam CCD vs the same CCD shifted +10 deg agree in pixel count); get_ccds_with_psf subtraction pinned against the real get_all_shdus. --- tests/module/test_coverage.py | 329 ++++++++++++++++++++++++++-------- 1 file changed, 253 insertions(+), 76 deletions(-) diff --git a/tests/module/test_coverage.py b/tests/module/test_coverage.py index 96b104097..0c850ce4a 100644 --- a/tests/module/test_coverage.py +++ b/tests/module/test_coverage.py @@ -2,9 +2,11 @@ Covers the pure and lightly-fixtured logic behind the per-CCD coverage nexp masks: exposure-number parsing, the CCD-list -> unique-exposure reduction, the -handler's missing-CCD subtraction, per-CCD corner extraction from a multi-HDU -header, ``--ccd_list`` filtering, the builder's per-CCD row parsing and RA-wrap -guard, the power-of-two ``nside`` validation, the atomic-download rename, and +handler's missing-CCD subtraction (pinned against the real ID format), per-CCD +corner extraction from plain and fpack-compressed multi-HDU headers, +``--ccd_list`` filtering, the per-CCD resume path, the builder's per-CCD row +parsing, the accumulated exposure-count (nexp) contract, the RA-wrap and pole +guards, the power-of-two ``nside`` validation, the atomic-download rename, and the ``-h``/argv handling of the console runners. The heavy paths (real VOSpace download, production-resolution map building, plotting, multiprocessing) are exercised end to end by the pipeline, not here. @@ -12,17 +14,15 @@ import sys from types import SimpleNamespace -from unittest import mock import numpy as np import numpy.testing as npt import pytest from astropy import wcs -from astropy.io.fits import Header from hypothesis import given from hypothesis import strategies as st -from shapepipe.utilities import ccd_psf_handler as ccd_psf_handler_mod +from shapepipe.utilities import summary from shapepipe.utilities.ccd_psf_handler import CcdPsfHandler from shapepipe.utilities.coverage_map_builder import ( CoverageMapBuilder, @@ -32,6 +32,7 @@ FieldCornersExtractor, _ccd_corners, _expnum_from_path, + _image_shape, _parse_header_to_wcs, ) from shapepipe.utilities.header_downloader import HeaderDownloader @@ -40,8 +41,8 @@ def _tan_wcs(crval_ra, crval_dec=0.0, nx=2080, ny=4612): """Build a minimal 2-D TAN WCS with a populated pixel shape. - ``nx``/``ny`` set ``pixel_shape`` (populated by astropy from - ``NAXIS1``/``NAXIS2``), so ``_ccd_corners`` can read the CCD pixel bounds. + ``nx``/``ny`` set ``pixel_shape`` so ``_header_text`` can emit matching + ``NAXIS1``/``NAXIS2`` keywords. """ w = wcs.WCS(naxis=2) w.wcs.ctype = ["RA---TAN", "DEC--TAN"] @@ -53,10 +54,8 @@ def _tan_wcs(crval_ra, crval_dec=0.0, nx=2080, ny=4612): def _header_text(wcs_list): - """Render a multi-HDU header text file: a primary HDU then one per WCS. - - Each CCD header carries ``NAXIS1``/``NAXIS2`` so the reconstructed WCS has a - pixel shape. + """Render a multi-HDU header text file: a primary HDU then one plain image + HDU per WCS, each carrying ``NAXIS1``/``NAXIS2``. """ blocks = ["SIMPLE = T"] for w in wcs_list: @@ -69,6 +68,34 @@ def _header_text(wcs_list): return "".join(f"{block}\nEND \n" for block in blocks) +def _compressed_header_text(w): + """Render a one-CCD header mimicking an fpack tile-compressed HDU. + + ``NAXIS1``/``NAXIS2`` describe the compressed binary table (byte width, row + count); the true image dimensions live in ``ZNAXIS1``/``ZNAXIS2``. This is + the shape of headers fetched from 'p.fits.fz' with ``head=True``. + """ + nx, ny = w.pixel_shape + wcs_cards = "\n".join(str(card) for card in w.to_header().cards) + ccd = ( + "XTENSION= 'BINTABLE'\n" + "BITPIX = 8\n" + "NAXIS = 2\n" + "NAXIS1 = 8\n" + f"NAXIS2 = {ny}\n" + "PCOUNT = 1000000\n" + "GCOUNT = 1\n" + "TFIELDS = 1\n" + "ZIMAGE = T\n" + "ZBITPIX = -32\n" + "ZNAXIS = 2\n" + f"ZNAXIS1 = {nx}\n" + f"ZNAXIS2 = {ny}\n" + f"{wcs_cards}" + ) + return f"SIMPLE = T\nEND \n{ccd}\nEND \n" + + # --------------------------------------------------------------------------- # _expnum_from_path # --------------------------------------------------------------------------- @@ -113,6 +140,16 @@ def test_get_exposures_reduces_to_unique_exposures(tmp_path): npt.assert_array_equal(exps, np.array([2143523, 2143524])) +def test_get_exposures_single_line(tmp_path): + """A one-line CCD list (0-d loadtxt array) still yields one exposure.""" + ccd_list = tmp_path / "ccds.txt" + ccd_list.write_text("2143523-0\n") + + exps = HeaderDownloader().get_exposures(str(ccd_list)) + + npt.assert_array_equal(exps, np.array([2143523])) + + def test_get_exposures_csv_matches_txt(tmp_path): """The CSV and text code paths yield the same unique exposures.""" txt = tmp_path / "ccds.txt" @@ -177,29 +214,81 @@ def failing_copy(source, target, head=True): # --------------------------------------------------------------------------- def test_get_ccds_with_psf_subtracts_missing(monkeypatch): - """Valid CCDs are all exposure single-HDUs minus the missing set.""" + """Valid CCDs are all exposure single-HDUs minus the missing set. + + The real ``summary.get_all_shdus`` is used so the cross-component + ``-`` ID format is pinned end to end. + """ handler = CcdPsfHandler() # Two exposures, 3 CCDs each -> 6 candidate CCDs; two are missing. - monkeypatch.setattr( - handler, "get_exp", lambda patches: {"100", "200"} - ) + monkeypatch.setattr(handler, "get_exp", lambda patches: {"100", "200"}) monkeypatch.setattr( handler, "get_exp_shdu_missing", lambda patches: {"100-1", "200-2"}, ) - monkeypatch.setattr( - ccd_psf_handler_mod.summary, - "get_all_shdus", - lambda exposures, n_CCD: [ - f"{e}-{c}" for e in exposures for c in range(n_CCD) - ], - ) result = handler.get_ccds_with_psf(["P1"], n_CCD=3) + # get_all_shdus yields "-" for ccd in range(n_CCD). assert result == {"100-0", "100-2", "200-0", "200-1"} + # Guard the assumption that the missing IDs share the produced format. + assert set(summary.get_all_shdus({"100"}, 3)) == {"100-0", "100-1", "100-2"} + + +# --------------------------------------------------------------------------- +# image-shape resolution (fpack ZNAXIS vs plain NAXIS) +# --------------------------------------------------------------------------- + +def test_image_shape_prefers_znaxis_for_compressed_header(tmp_path): + """A compressed HDU reports ZNAXIS dims, not the binary-table NAXIS.""" + w = _tan_wcs(100.0, 20.0, nx=2080, ny=4612) + path = tmp_path / "1234567.txt" + path.write_text(_compressed_header_text(w)) + + (parsed_w, shape), = _parse_header_to_wcs(str(path)) + + # WCS pixel_shape would wrongly report the compressed byte width (8). + assert parsed_w.pixel_shape == (8, 4612) + # _image_shape recovers the true image dimensions. + assert shape == (2080, 4612) + + +def test_image_shape_falls_back_to_naxis(tmp_path): + """A plain image HDU (no ZIMAGE) uses NAXIS1/NAXIS2.""" + w = _tan_wcs(100.0, 20.0, nx=2080, ny=4612) + path = tmp_path / "1234567.txt" + path.write_text(_header_text([w])) + + (_, shape), = _parse_header_to_wcs(str(path)) + + assert shape == (2080, 4612) + + +def test_image_shape_raises_without_dimensions(): + """A header with no image dimensions raises a clear ``ValueError``.""" + from astropy.io.fits import Header + + header = Header() + header["CTYPE1"] = "RA---TAN" + + with pytest.raises(ValueError, match="image dimensions"): + _image_shape(header) + + +def test_compressed_header_corners_are_full_width(tmp_path): + """Corners from a compressed header span the true CCD width, not 8 px.""" + w = _tan_wcs(100.0, 20.0, nx=2080, ny=4612) + path = tmp_path / "1234567.txt" + path.write_text(_compressed_header_text(w)) + + (parsed_w, shape), = _parse_header_to_wcs(str(path)) + ra, dec = _ccd_corners(parsed_w, shape) + + # RA extent must reflect ~2080 px * 1e-5 deg/px * cos(dec), not 8 px. + ra_extent = max(ra) - min(ra) + assert ra_extent > 0.01 # --------------------------------------------------------------------------- @@ -211,7 +300,7 @@ def test_ccd_corners_returns_four_corners_around_centre(): nx, ny = 2080, 4612 w = _tan_wcs(100.0, 20.0, nx=nx, ny=ny) - ra, dec = _ccd_corners(w) + ra, dec = _ccd_corners(w, (nx, ny)) assert len(ra) == 4 assert len(dec) == 4 @@ -219,22 +308,11 @@ def test_ccd_corners_returns_four_corners_around_centre(): assert min(ra) < 100.0 < max(ra) assert min(dec) < 20.0 < max(dec) # Extent matches the CD-scaled pixel span (with edge padding). - npt.assert_allclose(max(dec) - min(dec), (ny) * 1e-5, rtol=1e-3) - - -def test_ccd_corners_raises_without_pixel_shape(): - """A WCS with no pixel shape raises a clear ``ValueError``.""" - w = wcs.WCS(naxis=2) - w.wcs.ctype = ["RA---TAN", "DEC--TAN"] - w.wcs.crval = [100.0, 20.0] - assert w.pixel_shape is None - - with pytest.raises(ValueError, match="pixel_shape"): - _ccd_corners(w) + npt.assert_allclose(max(dec) - min(dec), ny * 1e-5, rtol=1e-3) def test_parse_header_to_wcs_returns_one_wcs_per_extension(tmp_path): - """One WCS is returned per CCD HDU; the primary HDU is skipped.""" + """One (wcs, shape) pair is returned per CCD HDU; the primary is skipped.""" ccd_wcs = [_tan_wcs(10.0), _tan_wcs(20.0), _tan_wcs(30.0)] path = tmp_path / "1234567.txt" @@ -243,8 +321,7 @@ def test_parse_header_to_wcs_returns_one_wcs_per_extension(tmp_path): result = _parse_header_to_wcs(str(path)) assert len(result) == len(ccd_wcs) - # Reconstructed WCS carry the pixel shape from NAXIS1/NAXIS2. - assert all(w.pixel_shape == (2080, 4612) for w in result) + assert all(shape == (2080, 4612) for _, shape in result) # --------------------------------------------------------------------------- @@ -323,19 +400,81 @@ def test_run_extract_writes_per_ccd_rows(tmp_path): assert all(len(line.split()) == 9 for line in lines) -def test_get_done_exposures_keys_on_expnum(tmp_path): - """Resume treats an exposure as done if any of its CCD rows are present.""" +# --------------------------------------------------------------------------- +# FieldCornersExtractor resume path (per-CCD done-set) +# --------------------------------------------------------------------------- + +def _write_headers(header_dir, expnums, n_ccd=3): + """Write one plain multi-HDU header per exposure into ``header_dir``.""" + header_dir.mkdir(exist_ok=True) + for j, expnum in enumerate(expnums): + ccd_wcs = [_tan_wcs(10.0 + j + i) for i in range(n_ccd)] + (header_dir / f"{expnum}.txt").write_text(_header_text(ccd_wcs)) + + +def test_get_done_ccds_reads_present_ids(tmp_path): + """The done-set is the exact set of CCD IDs already in the output.""" out = tmp_path / "corners.txt" out.write_text( "1234567-0 10 10 10 10 20 20 20 20\n" - "1234567-5 10 10 10 10 20 20 20 20\n" - "7654321-3 30 30 30 30 40 40 40 40\n" + "1234567-2 10 10 10 10 20 20 20 20\n" ) extractor = FieldCornersExtractor() extractor._params["output_file"] = str(out) - assert extractor.get_done_exposures() == {1234567, 7654321} + assert extractor.get_done_ccds() == {"1234567-0", "1234567-2"} + + +def test_resume_adds_new_exposure_without_duplicating(tmp_path): + """Resume appends a new exposure and leaves existing rows untouched.""" + header_dir = tmp_path / "headers" + _write_headers(header_dir, [1000001, 1000002]) + + out = tmp_path / "corners.txt" + # Pre-populate with exposure 1's three CCD rows. + pre = FieldCornersExtractor().process_single_header( + (str(header_dir / "1000001.txt"), False, None) + ) + with open(out, "w") as f: + for ccd_id, ra, dec in pre: + f.write(f"{ccd_id} " + " ".join(f"{v:.5f}" for v in ra + dec) + "\n") + + FieldCornersExtractor().run( + args=["-i", str(header_dir), "-o", str(out), "-r"] + ) + + ids = [line.split()[0] for line in out.read_text().splitlines()] + # No duplicate exposure-1 rows; exposure 2's three CCDs added. + assert ids.count("1000001-0") == 1 + assert sorted(ids) == [ + "1000001-0", "1000001-1", "1000001-2", + "1000002-0", "1000002-1", "1000002-2", + ] + + +def test_resume_completes_partial_exposure(tmp_path): + """An exposure interrupted mid-write is completed, not skipped.""" + header_dir = tmp_path / "headers" + _write_headers(header_dir, [1000001]) + + out = tmp_path / "corners.txt" + # Simulate an interrupt: only the first of exposure 1's CCDs was written. + pre = FieldCornersExtractor().process_single_header( + (str(header_dir / "1000001.txt"), False, None) + ) + ccd_id, ra, dec = pre[0] + with open(out, "w") as f: + f.write(f"{ccd_id} " + " ".join(f"{v:.5f}" for v in ra + dec) + "\n") + + FieldCornersExtractor().run( + args=["-i", str(header_dir), "-o", str(out), "-r"] + ) + + ids = [line.split()[0] for line in out.read_text().splitlines()] + # The missing CCDs are filled in; the present one is not duplicated. + assert sorted(ids) == ["1000001-0", "1000001-1", "1000001-2"] + assert ids.count("1000001-0") == 1 # --------------------------------------------------------------------------- @@ -384,9 +523,54 @@ def test_check_params_missing_input_file_raises(tmp_path): # --------------------------------------------------------------------------- -# CoverageMapBuilder — RA-wrap guard and per-CCD row build +# CoverageMapBuilder — parsing, nexp contract, RA-wrap and pole guards # --------------------------------------------------------------------------- +def _read_map(path): + import healsparse as hsp + + return hsp.HealSparseMap.read(str(path)) + + +def test_build_map_single_row(tmp_path): + """A one-row corners file exercises the atleast_1d/2d parsing guards.""" + infile = tmp_path / "corners.txt" + infile.write_text("100-0 9.9 10.1 10.1 9.9 0.0 0.0 0.2 0.2\n") + out = tmp_path / "cov.hsp" + + CoverageMapBuilder().run( + args=["-i", str(infile), "-o", str(out), "-c", "32", "-n", "1024"] + ) + + m = _read_map(out) + assert 0 < len(m.valid_pixels) < 200 + assert m[m.valid_pixels].max() == 1 + + +def test_build_map_nexp_counts_overlapping_exposures(tmp_path): + """Two overlapping CCDs from different exposures give value 2 in overlap.""" + # Two 0.4x0.4 deg CCDs offset by 0.2 deg in RA -> a central overlap strip. + infile = tmp_path / "corners.txt" + infile.write_text( + "100-0 9.8 10.2 10.2 9.8 19.8 19.8 20.2 20.2\n" + "200-0 10.0 10.4 10.4 10.0 19.8 19.8 20.2 20.2\n" + ) + out = tmp_path / "cov.hsp" + + CoverageMapBuilder().run( + args=["-i", str(infile), "-o", str(out), "-c", "32", "-n", "1024"] + ) + + m = _read_map(out) + values = m[m.valid_pixels] + # The overlap is covered by both exposures (value 2); the union edges by + # one (value 1). Both must be present; nothing exceeds 2. + assert values.max() == 2 + assert values.min() == 1 + assert (values == 2).sum() > 0 + assert (values == 1).sum() > 0 + + def test_unwrap_ra_shifts_straddling_seam(): """A polygon straddling RA=0 is put on a common (negative) branch.""" out = unwrap_ra([359.9, 0.1, 0.1, 359.9]) @@ -401,41 +585,35 @@ def test_unwrap_ra_leaves_normal_polygon_unchanged(): ) -def test_build_map_wrapped_ccd_fills_small_patch(tmp_path): - """A CCD straddling RA=0 fills a small footprint, not the ~360° complement. +def test_build_map_invariant_under_ra_shift(tmp_path): + """A seam CCD's footprint matches an identical CCD shifted +10 deg in RA. - Verified empirically against a same-size CCD away from the seam: with the - unwrap guard the two footprints have the same pixel count to within a few - boundary pixels. (Without unwrapping, HealSparse happens to fill the small - region too for a convex CCD-sized polygon, but the guard makes the intent - explicit and is robust to degenerate orderings.) + Rotating both the seam CCD (via +360/unwrap) and a reference CCD onto the + same RA and comparing pixel counts fails if the seam polygon were filling + the ~360 deg complement. This is the real RA-wrap regression guard. """ - # CCD straddling RA=0 (0.2 x 0.2 deg) and a reference CCD at RA=10. - infile = tmp_path / "corners.txt" - infile.write_text( - "100-0 359.9 0.1 0.1 359.9 0.0 0.0 0.2 0.2\n" - "200-0 9.9 10.1 10.1 9.9 0.0 0.0 0.2 0.2\n" + # Seam CCD straddling RA=0, and the same CCD translated to RA~10. + seam = tmp_path / "seam.txt" + seam.write_text("100-0 359.9 0.1 0.1 359.9 20.0 20.0 20.2 20.2\n") + ref = tmp_path / "ref.txt" + ref.write_text("200-0 9.9 10.1 10.1 9.9 20.0 20.0 20.2 20.2\n") + + m_seam = tmp_path / "seam.hsp" + m_ref = tmp_path / "ref.hsp" + CoverageMapBuilder().run( + args=["-i", str(seam), "-o", str(m_seam), "-c", "32", "-n", "1024"] ) - out = tmp_path / "cov.hsp" - - builder = CoverageMapBuilder() - builder.run( - args=[ - "-i", str(infile), - "-o", str(out), - "-c", "32", - "-n", "1024", - ] + CoverageMapBuilder().run( + args=["-i", str(ref), "-o", str(m_ref), "-c", "32", "-n", "1024"] ) - import healsparse as hsp - - m = hsp.HealSparseMap.read(str(out)) - n_valid = len(m.valid_pixels) + n_seam = len(_read_map(m_seam).valid_pixels) + n_ref = len(_read_map(m_ref).valid_pixels) - # Two ~0.04 deg^2 patches at nside=1024 (~3.4 arcmin pixels): a small, - # bounded number of pixels — nowhere near a hemisphere. - assert 0 < n_valid < 200 + # Same-size footprints at the same declination: pixel counts agree to + # within a few boundary pixels, and are nowhere near a hemisphere. + assert abs(n_seam - n_ref) <= max(3, int(0.1 * n_ref)) + assert 0 < n_seam < 200 def test_build_map_pole_guard_skips_polygon(tmp_path, capsys): @@ -447,8 +625,7 @@ def test_build_map_pole_guard_skips_polygon(tmp_path, capsys): ) out = tmp_path / "cov.hsp" - builder = CoverageMapBuilder() - builder.run( + CoverageMapBuilder().run( args=["-i", str(infile), "-o", str(out), "-c", "32", "-n", "1024"] ) From 15401d964cf36b2b27be2074727dabfe570dff30 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 02:21:06 +0200 Subject: [PATCH 31/34] feat(coverage): accept catalogue version v2.0 (9-patch footprint, as v1.6) Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01MqH9wjbJicVmaQdhZ4kBYr --- log_-i | 1 + src/shapepipe/utilities/ccd_psf_handler.py | 5 +++-- tests/module/test_coverage.py | 21 +++++++++++++++++++++ 3 files changed, 25 insertions(+), 2 deletions(-) create mode 100644 log_-i diff --git a/log_-i b/log_-i new file mode 100644 index 000000000..35d62bcc5 --- /dev/null +++ b/log_-i @@ -0,0 +1 @@ +-i /tmp/pytest-of-cdaley/pytest-5/test_build_map_pole_guard_skip0/corners.txt -o /tmp/pytest-of-cdaley/pytest-5/test_build_map_pole_guard_skip0/cov.hsp -c 32 -n 1024 diff --git a/src/shapepipe/utilities/ccd_psf_handler.py b/src/shapepipe/utilities/ccd_psf_handler.py index 2d0dc90eb..2e361df28 100644 --- a/src/shapepipe/utilities/ccd_psf_handler.py +++ b/src/shapepipe/utilities/ccd_psf_handler.py @@ -51,7 +51,7 @@ def params_default(self): } self._help_strings = { - "version_cat": "catalogue major version, allowed are v1.3, v1.4, v1.5, v1.6; default is {}", + "version_cat": "catalogue major version, allowed are v1.3, v1.4, v1.5, v1.6, v2.0; default is {}", "n_CCD": "number of CCDs per exposure; default is {}", "output": "output file path; default is ccds_with_psf_.txt", @@ -93,7 +93,8 @@ def update_params(self): n_patch = 7 elif version == "v1.5": n_patch = 8 - elif version == "v1.6": + elif version in ("v1.6", "v2.0"): + # v2.0 (ngmix) covers the same 9-patch footprint as v1.6 n_patch = 9 else: raise ValueError(f"Invalid version {version}") diff --git a/tests/module/test_coverage.py b/tests/module/test_coverage.py index 0c850ce4a..8e79d427c 100644 --- a/tests/module/test_coverage.py +++ b/tests/module/test_coverage.py @@ -237,6 +237,27 @@ def test_get_ccds_with_psf_subtracts_missing(monkeypatch): assert set(summary.get_all_shdus({"100"}, 3)) == {"100-0", "100-1", "100-2"} +@pytest.mark.parametrize( + ("version", "n_patch"), + [("v1.3", 7), ("v1.4", 7), ("v1.5", 8), ("v1.6", 9), ("v2.0", 9)], +) +def test_version_to_patch_count(version, n_patch): + """Each catalogue version maps to its patch count; v2.0 mirrors v1.6.""" + handler = CcdPsfHandler() + handler._params["version_cat"] = version + handler.update_params() + assert handler._params["n_patch"] == n_patch + assert len(handler._params["patches"]) == n_patch + + +def test_invalid_version_raises(): + """An unknown catalogue version fails loudly.""" + handler = CcdPsfHandler() + handler._params["version_cat"] = "v9.9" + with pytest.raises(ValueError, match="v9.9"): + handler.update_params() + + # --------------------------------------------------------------------------- # image-shape resolution (fpack ZNAXIS vs plain NAXIS) # --------------------------------------------------------------------------- From a3aa4ec5b59f4c94a38a4b3722f156ed1ab08852 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 03:05:08 +0200 Subject: [PATCH 32/34] feat(coverage): v2.0 defaults to 10 patches; -p overrides the version map The v2 footprint is still growing (>5000 deg2 per the 2026-07 Paris meeting, new Dec<30 region); the newest patch definitions (auxdir/CFIS/tiles_202510) include P10. Default v2.0 to 10 patches and add an explicit -p/--n_patch override so the handler tracks evolving patch layouts without code changes. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01MqH9wjbJicVmaQdhZ4kBYr --- log_-i | 2 +- src/shapepipe/utilities/ccd_psf_handler.py | 20 ++++++++++++++++---- tests/module/test_coverage.py | 11 ++++++++++- 3 files changed, 27 insertions(+), 6 deletions(-) diff --git a/log_-i b/log_-i index 35d62bcc5..2d633c17b 100644 --- a/log_-i +++ b/log_-i @@ -1 +1 @@ --i /tmp/pytest-of-cdaley/pytest-5/test_build_map_pole_guard_skip0/corners.txt -o /tmp/pytest-of-cdaley/pytest-5/test_build_map_pole_guard_skip0/cov.hsp -c 32 -n 1024 +-i /tmp/pytest-of-cdaley/pytest-7/test_build_map_pole_guard_skip0/corners.txt -o /tmp/pytest-of-cdaley/pytest-7/test_build_map_pole_guard_skip0/cov.hsp -c 32 -n 1024 diff --git a/src/shapepipe/utilities/ccd_psf_handler.py b/src/shapepipe/utilities/ccd_psf_handler.py index 2e361df28..4a48e1a52 100644 --- a/src/shapepipe/utilities/ccd_psf_handler.py +++ b/src/shapepipe/utilities/ccd_psf_handler.py @@ -37,23 +37,30 @@ def params_default(self): self._params = { "version_cat": "v1.6", "n_CCD": 40, + "n_patch": 0, "output": None, } self._short_options = { "version_cat": "-V", "n_CCD": "-n", + "n_patch": "-p", "output": "-o", } self._types = { "n_CCD": "int", + "n_patch": "int", } self._help_strings = { "version_cat": "catalogue major version, allowed are v1.3, v1.4, v1.5, v1.6, v2.0; default is {}", "n_CCD": "number of CCDs per exposure; default is {}", + "n_patch": ( + "number of patches (P1..P); overrides the version" + + " default, 0 = use version default; default is {}" + ), "output": "output file path; default is ccds_with_psf_.txt", } @@ -87,15 +94,20 @@ def update_params(self): Set derived parameters based on input parameters. """ - # Determine number of patches based on version + # Determine number of patches: explicit -p wins, else version default version = self._params["version_cat"] - if version in ("v1.3", "v1.4"): + if self._params["n_patch"]: + n_patch = self._params["n_patch"] + elif version in ("v1.3", "v1.4"): n_patch = 7 elif version == "v1.5": n_patch = 8 - elif version in ("v1.6", "v2.0"): - # v2.0 (ngmix) covers the same 9-patch footprint as v1.6 + elif version == "v1.6": n_patch = 9 + elif version == "v2.0": + # Newest patch definitions (auxdir/CFIS/tiles_202510) include + # P10; override with -p if the v2 footprint settles differently. + n_patch = 10 else: raise ValueError(f"Invalid version {version}") diff --git a/tests/module/test_coverage.py b/tests/module/test_coverage.py index 8e79d427c..e70d31de6 100644 --- a/tests/module/test_coverage.py +++ b/tests/module/test_coverage.py @@ -239,7 +239,7 @@ def test_get_ccds_with_psf_subtracts_missing(monkeypatch): @pytest.mark.parametrize( ("version", "n_patch"), - [("v1.3", 7), ("v1.4", 7), ("v1.5", 8), ("v1.6", 9), ("v2.0", 9)], + [("v1.3", 7), ("v1.4", 7), ("v1.5", 8), ("v1.6", 9), ("v2.0", 10)], ) def test_version_to_patch_count(version, n_patch): """Each catalogue version maps to its patch count; v2.0 mirrors v1.6.""" @@ -250,6 +250,15 @@ def test_version_to_patch_count(version, n_patch): assert len(handler._params["patches"]) == n_patch +def test_n_patch_option_overrides_version_default(): + """An explicit -p patch count wins over the version mapping.""" + handler = CcdPsfHandler() + handler._params["version_cat"] = "v2.0" + handler._params["n_patch"] = 4 + handler.update_params() + assert handler._params["patches"] == ["P1", "P2", "P3", "P4"] + + def test_invalid_version_raises(): """An unknown catalogue version fails loudly.""" handler = CcdPsfHandler() From ae1017658dc7e18263f609e58f2075e4a7d9afee Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 03:10:47 +0200 Subject: [PATCH 33/34] =?UTF-8?q?fix(coverage):=20v2.0=20is=20patch-less?= =?UTF-8?q?=20=E2=80=94=20read=20summaries=20from=20the=20run=20root?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Martin (PR #812 review): the patch concept is removed in v2.0. For -V v2.0 the handler reads summary/missing_job_32_all.txt and exp_numbers.txt relative to the run root instead of iterating P1..P9; -p restores an explicit P1..P layout when needed. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01MqH9wjbJicVmaQdhZ4kBYr --- log_-i | 2 +- src/shapepipe/utilities/ccd_psf_handler.py | 14 ++++++++------ tests/module/test_coverage.py | 12 ++++++++++-- 3 files changed, 19 insertions(+), 9 deletions(-) diff --git a/log_-i b/log_-i index 2d633c17b..e90567d12 100644 --- a/log_-i +++ b/log_-i @@ -1 +1 @@ --i /tmp/pytest-of-cdaley/pytest-7/test_build_map_pole_guard_skip0/corners.txt -o /tmp/pytest-of-cdaley/pytest-7/test_build_map_pole_guard_skip0/cov.hsp -c 32 -n 1024 +-i /tmp/pytest-of-cdaley/pytest-8/test_build_map_pole_guard_skip0/corners.txt -o /tmp/pytest-of-cdaley/pytest-8/test_build_map_pole_guard_skip0/cov.hsp -c 32 -n 1024 diff --git a/src/shapepipe/utilities/ccd_psf_handler.py b/src/shapepipe/utilities/ccd_psf_handler.py index 4a48e1a52..495125624 100644 --- a/src/shapepipe/utilities/ccd_psf_handler.py +++ b/src/shapepipe/utilities/ccd_psf_handler.py @@ -59,7 +59,8 @@ def params_default(self): "n_CCD": "number of CCDs per exposure; default is {}", "n_patch": ( "number of patches (P1..P); overrides the version" - + " default, 0 = use version default; default is {}" + + " default, 0 = use version default (v2.0 is patch-less:" + + " summary files are read from the run root); default is {}" ), "output": "output file path; default is ccds_with_psf_.txt", } @@ -94,7 +95,8 @@ def update_params(self): Set derived parameters based on input parameters. """ - # Determine number of patches: explicit -p wins, else version default + # Determine patch directories: explicit -p wins, else version default. + # v2.0 removes the patch concept: summary files live at the run root. version = self._params["version_cat"] if self._params["n_patch"]: n_patch = self._params["n_patch"] @@ -105,14 +107,14 @@ def update_params(self): elif version == "v1.6": n_patch = 9 elif version == "v2.0": - # Newest patch definitions (auxdir/CFIS/tiles_202510) include - # P10; override with -p if the v2 footprint settles differently. - n_patch = 10 + n_patch = 0 else: raise ValueError(f"Invalid version {version}") self._params["n_patch"] = n_patch - self._params["patches"] = [f"P{x}" for x in np.arange(n_patch) + 1] + self._params["patches"] = ( + [f"P{x}" for x in np.arange(n_patch) + 1] if n_patch else ["."] + ) # Set output file if not specified if self._params["output"] is None: diff --git a/tests/module/test_coverage.py b/tests/module/test_coverage.py index e70d31de6..7f308b928 100644 --- a/tests/module/test_coverage.py +++ b/tests/module/test_coverage.py @@ -239,10 +239,10 @@ def test_get_ccds_with_psf_subtracts_missing(monkeypatch): @pytest.mark.parametrize( ("version", "n_patch"), - [("v1.3", 7), ("v1.4", 7), ("v1.5", 8), ("v1.6", 9), ("v2.0", 10)], + [("v1.3", 7), ("v1.4", 7), ("v1.5", 8), ("v1.6", 9)], ) def test_version_to_patch_count(version, n_patch): - """Each catalogue version maps to its patch count; v2.0 mirrors v1.6.""" + """Each v1.x catalogue version maps to its patch count.""" handler = CcdPsfHandler() handler._params["version_cat"] = version handler.update_params() @@ -250,6 +250,14 @@ def test_version_to_patch_count(version, n_patch): assert len(handler._params["patches"]) == n_patch +def test_v2_is_patchless(): + """v2.0 removes the patch concept: summaries read from the run root.""" + handler = CcdPsfHandler() + handler._params["version_cat"] = "v2.0" + handler.update_params() + assert handler._params["patches"] == ["."] + + def test_n_patch_option_overrides_version_default(): """An explicit -p patch count wins over the version mapping.""" handler = CcdPsfHandler() From 46214c729c69f8543d97b6702f7199fbecc9f811 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 16 Jul 2026 03:16:14 +0200 Subject: [PATCH 34/34] docs: state the patch convention per catalogue version (v2.0 is patch-less) Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01MqH9wjbJicVmaQdhZ4kBYr --- docs/source/pipeline_canfar.md | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/docs/source/pipeline_canfar.md b/docs/source/pipeline_canfar.md index fcac5b371..c548979bc 100644 --- a/docs/source/pipeline_canfar.md +++ b/docs/source/pipeline_canfar.md @@ -29,7 +29,15 @@ canfar auth switch default ### Set variables (optional) -Set the current patch in the shell as +> **Patch convention by catalogue version.** Runs up to `v1.6` are organised +> in sky patches `P1`–`P9` (7 patches for `v1.3`/`v1.4`, 8 for `v1.5`), each +> patch a run directory of its own with its `summary/` subdirectory. **The +> patch concept is removed in `v2.0`**: there is a single run root, and +> summary files (`summary/missing_job_32_all.txt`, `exp_numbers.txt`) live +> directly under it. Tools that take a catalogue version (e.g. +> `get_ccds_with_psf -V`) follow this convention automatically. + +For a `v1.x` run, set the current patch in the shell as ```bash patch=P[1-9] @@ -503,6 +511,12 @@ sky pixel, the number of exposures with a valid PSF. The list is written to get_ccds_with_psf -v -V v1.6 -o ccds_with_psf_v1.6.txt ``` +For a (patch-less) `v2.0` run, call this from the run root instead: + +```bash +get_ccds_with_psf -v -V v2.0 -o ccds_with_psf_v2.0.txt +``` + Next, download the exposure headers; indicate (with `-d`) a directory of already downloaded headers, which are linked so duplicate downloads are skipped: