diff --git a/scripts/sh/find_exposures_theli.sh b/scripts/sh/find_exposures_theli.sh new file mode 100755 index 000000000..c42140704 --- /dev/null +++ b/scripts/sh/find_exposures_theli.sh @@ -0,0 +1,70 @@ +#!/usr/bin/env bash + +# Name: find_exposures.sh +# Description: Script to convert THELI images to shapepipe +# bookkeeping conventions +# Author: Lucie Baumont +# Date: v1.0 1/2023 + + +# Command line arguments + +## Help string +usage="Usage: $(basename "$0") THELI_DIR OUTPUT_DIR\n + THELI_DIR\n + \thome director of THELI images, eg. /n17data/baumont/THELI_W3/W3/CFISCOLLAB_V2.0.0A\n + OUTPUT_DIR\n + \troot directory where outputs will go \n +" + +## Help if no arguments +if [ -z $1 ]; then + echo -ne $usage + exit 1 +fi +theli_dir=$1 +output_dir=$2 + +#make necessary directories +mkdir -p ${output_dir}/tiles +mkdir -p ${output_dir}/tiles/exposure_lists +mkdir -p ${output_dir}/exposures +mkdir -p ${output_dir}/exposures/headers + +for dir in ${theli_dir}/CFIS*; + do + # create link for all tiles and rename them + coadd_dir=$dir/r.MP9602/coadd_V2.0.0A/ + single_dir=${coadd_dir/coadd/single} + echo $single_dir + header_dir=${coadd_dir/coadd/headers} + file=${coadd_dir}/*.cut.fits + base=$(basename ${file} _r.MP9602.V2.0.0A.swarp.cut.fits) + num_pattern=${base:5} + num=${num_pattern//[-._]/} + ln -s ${coadd_dir}/${base}_r.MP9602.V2.0.0A.swarp.cut.fits ${output_dir}/tiles/CFIS_$num.fits + ln -s ${coadd_dir}/${base}_r.MP9602.V2.0.0A.swarp.cut.weight.fits ${output_dir}/tiles/CFIS_$num.weight.fits + ln -s ${coadd_dir}/${base}_r.MP9602.V2.0.0A.swarp.cut.flag.fits ${output_dir}/tiles/CFIS_$num.flag.fits + # create exposure list, link files but exclude broken links + for file in $(find ${single_dir}/*.sub.fits -type l -not -xtype l); do + base=$(basename $file C.sub.fits) + weight_link=${single_dir}/${base}C.weight.fits.gz + if [ -L ${weight_link} ] ; then + if [ -e ${weight_link} ] ; then + echo $(basename $file C.sub.fits) + ln -s ${single_dir}/${base}C.sub.fits ${output_dir}/exposures + ln -s ${single_dir}/${base}C.weight.fits.gz ${output_dir}/exposures + # remove header line 2 which has french accents that cause problems with astropy, change wcs keys + awk 'NR!~/^(2)$/ {sub(/TAN/,"TPV"); print}' ${header_dir}/${base}.head > ${output_dir}/exposures/headers/${base}.head + + fi + fi + done >> ${output_dir}/tiles/exposure_lists/exp_numbers-$num.txt + + #list=$(find ${single_dir}/*.sub.fits | awk -F/ '{print substr($NF,1,8)}'|uniq) + #echo "${list}">exposure_lists/$num.txt + + done + + #uncompress step + gunzip -f ${output_dir}/exposures/*.gz diff --git a/shapepipe/modules/find_exposures_package/__init__.py b/shapepipe/modules/find_exposures_package/__init__.py index 0b0debcc7..579107c26 100644 --- a/shapepipe/modules/find_exposures_package/__init__.py +++ b/shapepipe/modules/find_exposures_package/__init__.py @@ -14,8 +14,8 @@ =========== Identify the exposure images that were co-added to produce the tiles -(stacked image). The image names are listed in the tile FITS header, -which is read by this module to extract the names. +(stacked image). The image names are listed in the megapipe tile FITS +header,which is read by this module to extract the names. The output ASCII file contains the image base names (without file extension). diff --git a/shapepipe/modules/get_images_package/get_images.py b/shapepipe/modules/get_images_package/get_images.py index f0a4220c6..e7fa6365f 100644 --- a/shapepipe/modules/get_images_package/get_images.py +++ b/shapepipe/modules/get_images_package/get_images.py @@ -59,10 +59,10 @@ def in2out_pattern(number): # replace dots ('.') with dashes ('-') to avoid confusion # with file extension delimiters number_final = re.sub(r'\.', '-', number) - + # replace underscore with dashes + number_final = re.sub(r'_', '-', number_final) # remove letters in number number_final = re.sub('[a-zA-Z]', '', number_final) - return number_final @@ -228,16 +228,17 @@ def get_file_list( list_files_per_type = [] for number in image_number_list: - + # find all files with exposure number and extension with glob + # loop over this list, make sure it is robust for a single exposure + # allow for prefix if use_output_file_pattern: # Transform input to output number patterns - number_final = in2out_pattern(number) - - # Keep initial dot in extension - x = in_ext[1:] - x2 = re.sub(r'\.', '', x) - ext_final = in_ext[0] + x2 + # remove leading suffix from suffix.fits files + x = re.sub(r'.+(?=\b.fits\b)','',in_ext) + # remove all but leading . + ext_final='.'+re.sub(r'\.','',x) + fbase = ( f'{self._output_file_pattern[idx]}{number_final}' ) diff --git a/shapepipe/modules/merge_headers_package/__init__.py b/shapepipe/modules/merge_headers_package/__init__.py index 48a972c14..c8a0eb962 100644 --- a/shapepipe/modules/merge_headers_package/__init__.py +++ b/shapepipe/modules/merge_headers_package/__init__.py @@ -2,7 +2,7 @@ This package contains the module for ``merge_headers``. -:Author: Axel Guinot +:Author: Axel Guinot, Lucie Baumont :Parent module: ``split_exp_runner`` diff --git a/shapepipe/modules/merge_headers_package/merge_headers.py b/shapepipe/modules/merge_headers_package/merge_headers.py index 95d3e9acc..c3f22f3a2 100644 --- a/shapepipe/modules/merge_headers_package/merge_headers.py +++ b/shapepipe/modules/merge_headers_package/merge_headers.py @@ -1,10 +1,10 @@ """MERGE HEADERS. This module merges the output *header* files of the ``split_exp`` -module. It creates a binnary file that contains the WCS of each CCD for each -exposure. +module or external scamp headers. It creates a binary file that +contains the WCS of each CCD for each exposure. -:Author: Axel Guinot +:Author: Axel Guinot, Lucie Baumont """ @@ -13,9 +13,10 @@ import numpy as np from sqlitedict import SqliteDict +from astropy.io.fits import Header +from astropy.wcs import WCS - -def merge_headers(input_file_list, output_dir): +def merge_headers(input_file_list, output_dir, pattern, ext_header, n_hdu): """Merge Headers. This function opens the files in the input file list and merges them into @@ -27,7 +28,12 @@ def merge_headers(input_file_list, output_dir): List of input files output_dir : str Output path - + pattern : str + File pattern + ext_header : bool + Use external scamp header if ``True`` + n_hdu : int + number of ccds Raises ------ TypeError @@ -42,24 +48,103 @@ def merge_headers(input_file_list, output_dir): # Open SqliteDict file final_file = SqliteDict(f'{output_dir}/log_exp_headers.sqlite') - # Set matching pattern - pattern = 'headers-' + # define file properties for bookkeeping + header_dir = os.path.split(input_file_list[0][0])[0] + ext = os.path.splitext(input_file_list[0][0])[1] + keys = unique_exposure_list(input_file_list, header_dir, pattern, ext, ext_header) + + for key in keys: + # Load header, numpy binary or external header + if ext_header: + final_file[key] = create_joint_header(n_hdu, key, header_dir, pattern, ext) + else: + file_path = datadir+pattern+"key"+ext + final_file[key] = np.load(file_path, allow_pickle=True) + # Close file + final_file.commit() + final_file.close() + +def unique_exposure_list(input_file_list, header_dir, file_pattern, ext, ext_header): + """unique_exposure_list. + + Extract unique exposure ids from file list. + + Parameters + ---------- + input_file_list : list + list of input files + header_dir : str + directory of headers + file_pattern : str + text in filename + ext : str + file extension, eg .npy, .head + ext_header : bool + Uses external scamp header if ``True`` - for file_path in input_file_list: - # Extract file path - file_path_scalar = file_path[0] - # Check file name against pattern + Returns + ------- + list + List of unique exposure numbers + """ + file_list = [] + if ext_header: + ext = '\-\d{1,2}' + ext + # extract keys + for file in input_file_list: + file_path_scalar = file[0] file_name = os.path.split(file_path_scalar)[1] - file_base_name = os.path.splitext(file_name)[0] - pattern_split = re.split(pattern, file_base_name) + root = re.split(ext, file_name)[0] + pattern_split = re.split(file_pattern, root) if len(pattern_split) < 2: raise IndexError( f'Regex "{pattern}" not found in base name "{file_base_name}".' ) key = pattern_split[1] - # Load Numpy binary file - final_file[key] = np.load(file_path_scalar, allow_pickle=True) + + file_list.append(key) + # only keep unique keys + unique_keys = list(set(file_list)) - # Close file - final_file.commit() - final_file.close() + return unique_keys + +def create_joint_header(n_hdu, key, header_dir, pattern, ext): + """create_joint_header. + + Packages scamp headers for ccds into a numpy array per exposure. + + Parameters + ---------- + n_hdu : int + total number of ccds + key : str + exposure number + header_dir : str + directory where headers are stored + pattern : str + file pattern, e.g. headers- + ext : str + header extension + Returns + ------- + list + compilation of headers corresponding to exposure number + """ + + header_file = np.zeros(n_hdu, dtype='O') + for idx in range(1, n_hdu + 1): + filepath= header_dir+'/'+pattern+key+'-'+str(idx)+ext + print(filepath) + # check for file + if os.path.exists(filepath): + h=Header.fromfile(filepath,sep='\n',endcard=False,padding=False) + try: + w = WCS(h) + except Exception: + print(f'WCS error for file {exp_path}') + raise + header_file[idx - 1] = {'WCS': w, 'header': h.tostring()} + else: + header_file[idx - 1] = 'empty' + + return header_file diff --git a/shapepipe/modules/merge_headers_runner.py b/shapepipe/modules/merge_headers_runner.py index 9e692bba7..dc41c98a5 100644 --- a/shapepipe/modules/merge_headers_runner.py +++ b/shapepipe/modules/merge_headers_runner.py @@ -31,12 +31,14 @@ def merge_headers_runner( output_dir = run_dirs['output'] if config.has_option(module_config_sec, 'OUTPUT_PATH'): output_dir = config.getexpanded(module_config_sec, 'OUTPUT_PATH') - + file_pattern = config.get(module_config_sec, 'FILE_PATTERN') + ext_header = config.get(module_config_sec, 'EXT_HEADER') + n_hdu = int(config.get(module_config_sec, 'N_HDU')) # Log output directory w_log.info(f'output_dir = {output_dir}') # Merge header files - merge_headers(input_file_list, output_dir) + merge_headers(input_file_list, output_dir, file_pattern, ext_header, n_hdu) # No return objects return None, None diff --git a/shapepipe/modules/ngmix_package/ngmix.py b/shapepipe/modules/ngmix_package/ngmix.py index 7cb614a42..d54b1c10a 100644 --- a/shapepipe/modules/ngmix_package/ngmix.py +++ b/shapepipe/modules/ngmix_package/ngmix.py @@ -813,7 +813,6 @@ def do_ngmix_metacal( weight_map = np.copy(weights[n_e]) weight_map[np.where(flags[n_e] != 0)] = 0. - weight_map[weight_map != 0] = 1 psf_guess = np.array([0., 0., 0., 0., psf_T, 1.]) try: diff --git a/shapepipe/modules/split_exp_package/split_exp.py b/shapepipe/modules/split_exp_package/split_exp.py index 2dd1a5875..06aedeba9 100644 --- a/shapepipe/modules/split_exp_package/split_exp.py +++ b/shapepipe/modules/split_exp_package/split_exp.py @@ -93,7 +93,7 @@ def create_hdus( output_suffix : str Suffix for the output file transf_coord : bool - Transform the WCS (``pv`` to ``sip``) if ``True`` + Transform the WCS (``TAN`` to ``TPV``) if ``True`` transf_int : bool Set data types to int if ``True`` save_header : bool @@ -103,11 +103,11 @@ def create_hdus( header_file = np.zeros(self._n_hdu, dtype='O') for idx in range(1, self._n_hdu + 1): - h = fits.getheader(exp_path, idx) if transf_coord: - stp.pv_to_sip(h) - + # correct WCS convention so astropy can read PVs + h['CTYPE1']='RA--TPV' + h['CTYPE2']='RA--TPV' d = fits.getdata(exp_path, idx) if transf_int: d = d.astype(np.int16) diff --git a/shapepipe/pipeline/file_handler.py b/shapepipe/pipeline/file_handler.py index b7479f6f3..a9124fa25 100644 --- a/shapepipe/pipeline/file_handler.py +++ b/shapepipe/pipeline/file_handler.py @@ -986,6 +986,7 @@ def _save_num_patterns( ) # Save file list + os.makedirs(os.path.dirname(output_file), exist_ok=True) np.save(output_file, np.array(final_file_list)) del true_file_list, final_file_list