#!/usr/bin/env python
""" make_custom_mosaic.py - Module to control the generation of user-defined multi-skycell mosaics. This
script extends the capabilities of the HAP multi-visit mosaics (MVM) processing code (hapmultisequencer.py).
Based on user input, this script will produce drizzle-combined mosaics that will be as large as necessary to
incorporate all input images. Mosaic images can span multiple skycells, but NOT multiple projection cells. In
the exceedingly rare case that a custom mosaic contains observations that fall on multiple adjacent
projection cells, the bounds of the mosaic will be automatically clipped at the projection cell edge.
As this script is simply a wrapper around hapmultisequencer.py, the final output products are the exact same
as what one would expect from a pipeline MVM run. For a complete details about all hapmultisequencer.py
inputs and output products, please refer to the hapmultisequencer.py documentation.
hapmultisequencer.py was explicitly written to process data from only a single skycell. However, this script
gives users the ability to create mosaics spanning multiple skycells. As such, all output products produced
by hapmultisequencer.py normally include both the projection cell and skycell names in the all output
products. To avoid any confusion that might arise during the creation of mosaics that span multiple skycells,
all custom mosaic output product filenames will by default include projection cell name and the projection
cell reference right ascension and declination instead of the skycell name. Users can also optionally
specify an output product filename prefix of their choosing.
The output world cooordinate system (WCS) information is based on that of the projection cell in which the
observations reside. If the input observations happen to fall in a region of the sky where more than one
projection cells overlap, the WCS information of the output products will be based on the projection cell
whose center is closest to the geometric center of the input observations.
For all detectors, output products will be generated with a plate scale of 0.04 arcseconds per pixel. The
script produces an additional set up output products for WFC3/IR observations. This additional set of
products is generated using the native detector platescale of 0.12 arcseconds per pixel. These files end in
"_coarse_all_flt.fits".
The output world coordinate system (WCS) information is based on that of the projection cell in which the
observations reside. If the input observations happen to fall in a region of the sky where more than one
projection cells overlap, the WCS information of the output products will be based on the projection cell
whose center is closest to the geometric center of the input observations.
USAGE:
>>> python drizzlepac/make_custom_mosaic.py -[los] <search pattern enclosed in quotes>
or
>>> python drizzlepac/make_custom_mosaic.py -[los] <text file with list of input files>
- The '-l' option allows the user to control the level of verboseness in the log statements displayed on
the screen and written to the .log file. The level of verboseness from left to right, and includes all
log statements with a log_level left of the specified level. Specifying "critical" will only
record/display "critical" log statements, and specifying "error" will record/display both "error" and
"critical" log statements, and so on.
- The '-o' option allows the user to specify the text string that will be used as the filename prefix all
files created by hapmultisequencer.py during the MVM custom mosaic generation process. If not
explicitly specified, all output files will start with the following formatted text string:
"hst-skycell-p<pppp>-ra<##>d<####>-dec<n|s><##>d<####>", where p<pppp> is the projection cell ID,
ra<##>d<####> are the whole-number and decimal portions of the right ascention, respectively, and
dec<n|s><##>d<####> are the whole-number and decimal portions of the declination, respectively. Note
that the "<n|s>" denotes if the declination is north (positive) or south (negative). Example: For
skycell = 1974, ra = 201.9512, and dec = +26.0012, The filename prefix would be
"skycell-p1974-ra201d9512-decn26d0012".
- When turned on, the '-s' option allows users to skip alignment of all input images to known Gaia/HSC
sources in the input image footprint and instead use the existing input image alignment solution.
Python USAGE:
>>> python
>>> from drizzlepac import make_custom_mosaic
>>> make_custom_mosaic.perform(<list file or search pattern>, log_level='info', output_file_prefix=None, skip_gaia_alignment=False)
"""
import argparse
import glob
import math
import os
import sys
import tempfile
import traceback
from astropy.io import fits
import numpy as np
from drizzlepac import hapmultisequencer
from drizzlepac.haputils import cell_utils
from drizzlepac.haputils import make_poller_files
from stsci.tools import logutil
__taskname__ = 'make_custom_mosaic'
MSG_DATEFMT = '%Y%j%H%M%S'
SPLUNK_MSG_FORMAT = '%(asctime)s %(levelname)s src=%(name)s- %(message)s'
log = logutil.create_logger(__name__, level=logutil.logging.NOTSET, stream=sys.stdout,
format=SPLUNK_MSG_FORMAT, datefmt=MSG_DATEFMT)
# ------------------------------------------------------------------------------------------------------------
[docs]
def calc_skycell_dist(x, y, x_ref, y_ref):
"""Calculate distance from one skyframe to another
Parameters
----------
x : int
skyframe x index
y : int
skyframe y index
x_ref : int
reference x index
y_ref : int
reference y index
Returns
-------
dist : float
distance from (x, y) to (x_ref, y_ref)
"""
dist = math.sqrt(float((x - x_ref)**2) + float((y - y_ref)**2))
return(dist)
# ------------------------------------------------------------------------------------------------------------
# ------------------------------------------------------------------------------------------------------------
[docs]
def create_poller_file(img_list, proj_cell_dict):
"""Subroutine that executes make_poller_files.generate_poller_file() to generate custom MVM poller file
for execution of hapmultisequencer.run_mvm_processing().
Parameters
----------
img_list : list
list of images to process
proj_cell_dict : dict
Dictionary containing projection cell information
Returns
-------
poller_filename : str
Name of the newly created poller_filename
"""
# Locate bottom-left most skycell
pcell_filename = os.path.join(os.path.abspath(os.path.dirname(__file__)), 'pars', 'allsky_cells.fits')
sc_nxy = fits.getval(pcell_filename, "SC_NXY",
ext=0) # Projection cell side length (in skycells) in X or Y
closest_skycell = ""
min_dist = calc_skycell_dist(sc_nxy, sc_nxy, 1, 1) + 2.0
for sc in proj_cell_dict.keys():
x_index = proj_cell_dict[sc].x_index
y_index = proj_cell_dict[sc].y_index
dist = calc_skycell_dist(x_index, y_index, 1, 1)
if dist < min_dist:
min_dist = dist
closest_skycell = sc
# Create poller file
poller_filename = "temp_{}_mvm.out".format(closest_skycell)
rootname_list = []
for imgname in img_list:
rootname_list.append(imgname+"\n")
tf = tempfile.NamedTemporaryFile(mode='w+t', dir=os.getcwd())
with open(tf.name, 'w') as f:
f.writelines(rootname_list)
f.close()
make_poller_files.generate_poller_file(tf.name, poller_file_type="mvm",
output_poller_filename=poller_filename,
skycell_name=closest_skycell)
return poller_filename
# ------------------------------------------------------------------------------------------------------------
[docs]
def compute_mosaic_limits(proj_cell_dict):
"""Compute min and max limits of the rectangle that encloses the mosaic observations
Parameters
----------
proj_cell_dict : dictionary
Dictionary containing projection cell information
Returns
-------
mosaic_limits : list
4-element list containing the mosaic bounding rectangle X min and max and Y min and max values
"""
# set up storage arrays
array_size = len(proj_cell_dict.keys()) * 4
ra_values = np.empty(array_size)
dec_values = np.empty(array_size)
x_values = np.empty(array_size)
y_values = np.empty(array_size)
# Use wcs.calc_footprint() to get skycell corners, convert RA, Dec to X, Y in projection cell frame of reference
i = 0
for sc_name in proj_cell_dict.keys():
for ra_dec in proj_cell_dict[sc_name].wcs.calc_footprint().tolist():
x_y = proj_cell_dict[sc_name].projection_cell.wcs.all_world2pix([ra_dec], 0).tolist()[0]
ra_values[i] = ra_dec[0]
dec_values[i] = ra_dec[1]
x_values[i] = x_y[0]
y_values[i] = x_y[1]
i += 1
# determine min and max X and Y values, convert them from floating points to integers
x_min = int(np.rint(x_values.min()))
x_max = int(np.rint(x_values.max()))
y_min = int(np.rint(y_values.min()))
y_max = int(np.rint(y_values.max()))
# Make sure all min and max X and Y values are within the projection cell. If not, reset to the projection cell limit.
ps_x_max = proj_cell_dict[sc_name].projection_cell.wcs.pixel_shape[0]
ps_y_max = proj_cell_dict[sc_name].projection_cell.wcs.pixel_shape[1]
if x_min < 0:
log.warning("Custom mosaic minimum X value '{}' is beyond the projection cell minimum X value of '{}' Resetting custom mosaic minimum X value to projection cell minimum X value.".format(x_min, 0))
x_min = 0
if x_max > ps_x_max:
log.warning("Custom mosaic maximum X value '{}' is beyond the projection cell maximum X value of '{}' Resetting custom mosaic maximum X value to projection cell maximum X value.".format(x_max, ps_x_max))
x_max = ps_x_max
if y_min < 0:
log.warning("Custom mosaic minimum Y value '{}' is beyond the projection cell minimum Y value of '{}' Resetting custom mosaic minimum Y value to projection cell minimum Y value.".format(y_min, 0))
y_min = 0
if y_max > ps_y_max:
log.warning("Custom mosaic maximum Y value '{}' is beyond the projection cell maximum Y value of '{}' Resetting custom mosaic maximum Y value to projection cell maximum Y value.".format(y_max, ps_y_max))
y_max = ps_y_max
mosaic_limits = [x_min, x_max, y_min, y_max]
if log.level <= 10: # only display x and y limits if logging level is set to "debug".
for mosaic_limit, limit_label in zip(mosaic_limits, ["minimum X", "maximum X", "minimum Y", "maximum Y"]):
log.debug("Custom mosaic {} value (in projection cell pixels) : {}".format(limit_label, mosaic_limit))
return mosaic_limits
# ------------------------------------------------------------------------------------------------------------
[docs]
def determine_projection_cell(img_list):
"""Determine which projection cell should be used as the basis for the WCS of the output mosaic
product(s) based on which projection cell center is closest to the center of the observations.
Parameters
----------
img_list : list
A list of images to process
Returns
-------
best_pc_dict : dict
Dictionary of SkyCell objects for the input observations. These SkyCell objects contain projection
cell information from the projection cell whose center is closest to the geometric center of the
input observations.
"""
# Create dictionary containing information about the skycells that contain input images
skycell_dict = cell_utils.get_sky_cells(img_list)
# recast skycell_dict into nested dictionary with the top-most layer keyed by projection cell name
proj_cell_dict = {}
for key in skycell_dict.keys():
proj_cell = key[9:13]
if proj_cell not in proj_cell_dict.keys():
proj_cell_dict[proj_cell] = {}
proj_cell_dict[proj_cell][key] = skycell_dict[key]
# Determine which skycell's WCS information should be used as the basis for WCS of the output product(s)
if len(proj_cell_dict.keys()) == 1:
log.info("Observations are present in only a single projection cell.")
best_pc = list(proj_cell_dict)[0]
else:
log.info("Observations are present in multiple projection cells.")
log.info("Output WCS will be based on WCS from the projection cell whose center is closest to the "
"center of the input observations.")
# Determine which projection cell's center is closest to the center of the observations
pcell_filename = os.path.join(os.path.abspath(os.path.dirname(__file__)), 'pars', 'allsky_cells.fits')
sc_nxy = fits.getval(pcell_filename, "SC_NXY", ext=0) # Projection cell side length (in skycells) in X or Y
min_dist = calc_skycell_dist(sc_nxy, sc_nxy, 1, 1) + 2.0
best_pc = ""
for pc in proj_cell_dict.keys():
dist_ra = np.empty(len(proj_cell_dict[pc].keys()))
for i, skycell_name in zip(range(0, len(proj_cell_dict[pc].keys())), proj_cell_dict[pc].keys()):
pc_center_length = math.ceil(21/2.0)
x_index = proj_cell_dict[pc][skycell_name].x_index
y_index = proj_cell_dict[pc][skycell_name].y_index
dist_ra[i] = calc_skycell_dist(x_index, y_index, pc_center_length, pc_center_length)
if min_dist > dist_ra.mean():
min_dist = dist_ra.mean()
best_pc = pc
log.info("Output WCS will be based on WCS from projection cell {}".format(best_pc))
best_pc_dict = proj_cell_dict[best_pc]
return best_pc_dict
# ------------------------------------------------------------------------------------------------------------
# ------------------------------------------------------------------------------------------------------------
[docs]
def main():
"""Command-line interface
Parameters
----------
None
Returns
-------
return_value : int
return value from the run. 0 for successful run, something else otherwise.
"""
# Parse command-line inputs
parser = argparse.ArgumentParser(description='Create custom mosaic based on user-specified images and '
'world coordinate system (WCS) information')
parser.add_argument('input_image_source',
help='Search pattern to be used to identify images to process (NOTE: Pattern must be '
'enclosed in single or double quotes) or alternately, the '
'name of a text file containing a list of images to process')
parser.add_argument('-l', '--log_level', required=False, default='info',
choices=['critical', 'error', 'warning', 'info', 'debug'],
help='The desired level of verboseness in the log statements displayed on the screen '
'and written to the .log file. The level of verboseness from left to right, and '
'includes all log statements with a log_level left of the specified level. '
'Specifying "critical" will only record/display "critical" log statements, and '
'specifying "error" will record/display both "error" and "critical" log '
'statements, and so on.')
parser.add_argument('-o', '--output_file_prefix', required=False,
help='Text string that will be used as the filename prefix all files created by '
'hapmultisequencer.py during the MVM custom mosaic generation process. If not '
'explicitly specified, all output files will start with the following formatted '
'text string: "hst-skycell-p<pppp>-ra<##>d<####>-dec<n|s><##>d<####>", where '
'p<pppp> is the projection cell ID, ra<##>d<####> are the whole-number and '
'decimal portions of the right ascention, respectively, and dec<n|s><##>d<####> '
'are the whole-number and decimal portions of the declination, respectively. '
'Note that the "<n|s>" denotes if the declination is north (positive) or south '
'(negative). Example: For skycell = 1974, ra = 201.9512, and dec = +26.0012, '
'The filename prefix would be "skycell-p1974-ra201d9512-decn26d0012".')
parser.add_argument('-s', '--skip_gaia_alignment', required=False, action='store_true',
help='Skip alignment of all input images to known Gaia/HSC sources in the input '
'image footprint? If this option is turned on, the existing input image '
'alignment solution will be used instead. The default is False.')
user_args = parser.parse_args()
return_value = perform(user_args.input_image_source,
log_level=user_args.log_level,
output_file_prefix=user_args.output_file_prefix,
skip_gaia_alignment=user_args.skip_gaia_alignment)
# ------------------------------------------------------------------------------------------------------------
if __name__ == '__main__':
main()