From 668e05d539ba8556a25ccf03d0d557ee729d5f5d Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Fri, 9 Jan 2026 14:18:07 +0000 Subject: [PATCH 1/9] 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 2/9] 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 4d7ebe3b..9a4153db 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 5f917ad0..e284e264 --- 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 00000000..fb18b0a7 --- /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 00000000..2707cc62 --- /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 3/9] 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 00000000..742acf38 --- /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 4/9] 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 9a4153db..b20220ec 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 00000000..3b30ce04 --- /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 00000000..3f91f763 --- /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 00000000..cd8d2c2f --- /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 00000000..323dd5ac --- /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 5/9] 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 2707cc62..78bbcca7 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 3f91f763..152598e1 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 cd8d2c2f..14f84fc6 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 323dd5ac..f0765a04 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 742acf38..e7cb70e6 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 6/9] 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 8b3af63d..ee4c6199 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 7/9] 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 78bbcca7..3ee29303 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 337b69ffec889ff79a4d3508b31a823a7ad7c8ba Mon Sep 17 00:00:00 2001 From: "martin.kilbinger" Date: Mon, 2 Feb 2026 11:46:13 +0100 Subject: [PATCH 8/9] Updated contributors --- docs/source/contributors.md | 7 +++++-- src/shapepipe/info.py | 10 +++++----- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/docs/source/contributors.md b/docs/source/contributors.md index a9ab324e..4762aa95 100644 --- a/docs/source/contributors.md +++ b/docs/source/contributors.md @@ -41,12 +41,15 @@ Martin Kilbinger - [Axel Guinot](https://github.com/aguinot) : Primary module developer, validation, performance and debugging - [Martin Kilbinger](https://github.com/martinkilbinger) : Scientific lead, module developer and debugging -## Principal Contributors +## Contributors - [Lucie Baumont](https://github.com/lbaumo) : Testing - [Jérôme Bonnin](https://github.com/jerome-bonnin) : PSF module developer +- [Cail Daley](https://www.cosmostat.org/people/cail-daley) : Container development +- [Fabian Hervas Peters](https://www.cosmostat.org/people/former-members/fabian-hervas-peters) : Simulations, validation - Marc Gentile : File IO, initial architecture -- Xavier Jimenez : External catalogue matching module developer +- [Mike Hudson](https://uwaterloo.ca/physics-astronomy/profile/mjhudson) : Further contributors +- Xavier Jimenez : External catalogue matching - [François Lanusse](https://github.com/eiffl) : Further contributors - [Tobias Liaudat](https://github.com/tobias-liaudat) : PSF module developer and validation - [Austin Peel](https://github.com/austinpeel) : Validation diff --git a/src/shapepipe/info.py b/src/shapepipe/info.py index 7ce7f1d7..de3f8894 100644 --- a/src/shapepipe/info.py +++ b/src/shapepipe/info.py @@ -69,11 +69,11 @@ def shapepipe_logo(colour=False): Martin Kilbinger Main Contributors: - Tobias Liaudat - Morgan Schmitz - Andre Zamorano Vitorelli - Francois Lanusse - Xavier Jimenez + Lucie Baumont + Cail Daley + Fabian Hervas Peters + Tobias Liaudat + Morgan Schmitz Version: {} From 44dfb7c3ce9549708d996bc5511dd9e971a65519 Mon Sep 17 00:00:00 2001 From: "martin.kilbinger" Date: Wed, 8 Apr 2026 10:26:51 +0200 Subject: [PATCH 9/9] Updated contributors --- src/shapepipe/info.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/shapepipe/info.py b/src/shapepipe/info.py index de3f8894..4748c2bc 100644 --- a/src/shapepipe/info.py +++ b/src/shapepipe/info.py @@ -69,11 +69,8 @@ def shapepipe_logo(colour=False): Martin Kilbinger Main Contributors: - Lucie Baumont Cail Daley Fabian Hervas Peters - Tobias Liaudat - Morgan Schmitz Version: {}