Skip to content

Run THELI on shapepipe part 1 #617

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 11 commits into
base: develop
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 70 additions & 0 deletions scripts/sh/find_exposures_theli.sh
Original file line number Diff line number Diff line change
@@ -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 <[email protected]>
# 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
4 changes: 2 additions & 2 deletions shapepipe/modules/find_exposures_package/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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).

Expand Down
19 changes: 10 additions & 9 deletions shapepipe/modules/get_images_package/get_images.py
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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}'
)
Expand Down
2 changes: 1 addition & 1 deletion shapepipe/modules/merge_headers_package/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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``

Expand Down
123 changes: 104 additions & 19 deletions shapepipe/modules/merge_headers_package/merge_headers.py
Original file line number Diff line number Diff line change
@@ -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

"""

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
6 changes: 4 additions & 2 deletions shapepipe/modules/merge_headers_runner.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
1 change: 0 additions & 1 deletion shapepipe/modules/ngmix_package/ngmix.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will create an issue. Line 859 there is weight_map *= 1 / sig_noise ** 2 if the weight_map is not equal to 1 it will lead to weird behavior.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

for theli weights, there are relative weights between pixels, which then need to be rescaled by the sky background variance to become variance maps. Thomas Erban seemed to think this operation was okay. are you talking about the weirdness associated with some megaprime weight maps that are not always 1 or 0?


psf_guess = np.array([0., 0., 0., 0., psf_T, 1.])
try:
Expand Down
8 changes: 4 additions & 4 deletions shapepipe/modules/split_exp_package/split_exp.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand Down
1 change: 1 addition & 0 deletions shapepipe/pipeline/file_handler.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down