From e51bd84ba2002d0c7485df582d476c7dc797d2ed Mon Sep 17 00:00:00 2001 From: Clara Brasseur Date: Thu, 7 May 2020 18:45:53 -0400 Subject: [PATCH 1/4] replacing drop_after with something better --- astrocut/cutouts.py | 127 +++++++++++++++++++++++++++++--------------- 1 file changed, 84 insertions(+), 43 deletions(-) diff --git a/astrocut/cutouts.py b/astrocut/cutouts.py index 822664c9..be91c8c3 100644 --- a/astrocut/cutouts.py +++ b/astrocut/cutouts.py @@ -2,24 +2,25 @@ """This module implements cutout functionality similar to fitscut.""" - +import os +import warnings import numpy as np -import astropy.units as u +from time import time +from datetime import date +from itertools import product + +from astropy import log +from astropy import units as u from astropy.io import fits from astropy.coordinates import SkyCoord from astropy import wcs +from astropy.utils.decorators import deprecated_renamed_argument from astropy.visualization import (SqrtStretch, LogStretch, AsinhStretch, SinhStretch, LinearStretch, MinMaxInterval, ManualInterval, AsymmetricPercentileInterval) from PIL import Image -from time import time -from datetime import date - -import os -import warnings - from . import __version__ from .exceptions import InputWarning, DataWarning, InvalidQueryError, InvalidInputError @@ -54,8 +55,6 @@ def _get_cutout_limits(img_wcs, center_coord, cutout_size): center_pixel = center_coord.to_pixel(img_wcs, 1) except wcs.NoConvergence: # If wcs can't converge, center coordinate is far from the footprint raise InvalidQueryError("Cutout location is not in image footprint!") - - # print("Center pixel: {center_pixel}".format(center_pixel)) lims = np.zeros((2, 2), dtype=int) @@ -126,9 +125,34 @@ def _get_cutout_wcs(img_wcs, cutout_lims): #### FUNCTIONS FOR UTILS #### +def remove_sip_coefficients(hdu_header): + """ + Remove standard sip coefficient keywords for a fits header. + + Everything is wrapped in try catches because we don't care is a keyword is not actually there. + """ + for lets in product(["A","B"], ["","P"]): + lets = ''.join(lets) + + try: + del hdu_header[f"{lets}_ORDER"] + except KeyError: + pass -def _hducut(img_hdu, center_coord, cutout_size, correct_wcs=False, drop_after=None, verbose=False): + try: + del hdu_header[f"{lets}_DMAX"] + except KeyError: + pass + + for i,j in product([0,1,2,3],[0,1,2,3]): + try: + del hdu_header[f"{lets}_{i}_{j}"] + except KeyError: + pass + + +def _hducut(img_hdu, center_coord, cutout_size, correct_wcs=False, verbose=False): """ Takes an ImageHDU (image and associated metatdata in the fits format), as well as a center coordinate and size and make a cutout of that image, which is returned as another ImageHDU, @@ -146,13 +170,7 @@ def _hducut(img_hdu, center_coord, cutout_size, correct_wcs=False, drop_after=No or `~astropy.Quantity` values, either pixels or angular quantities. correct_wcs : bool Default False. If true a new WCS will be created for the cutout that is tangent projected - and does not include distortions. - drop_after : str or None - Default None. When creating the header for the cutout (and crucially, before - building the WCS object) drop all header keywords starting with the one given. This is - useful particularly for drizzle files that contain a multitude of extranious keywords - and sometimes leftover WCS keywords that astropy will try to parse even thought they should be - ignored. + and does not include distortions. verbose : bool Default False. If true intermediate information is printed. @@ -161,18 +179,39 @@ def _hducut(img_hdu, center_coord, cutout_size, correct_wcs=False, drop_after=No response : `~astropy.io.fits.hdu.image.ImageHDU` The cutout image and associated metadata. """ - - # Pulling out the header - max_ind = len(img_hdu.header) - if drop_after is not None: - try: - max_ind = img_hdu.header.index(drop_after) - except ValueError: - warnings.warn("Last desired keyword not found in image header, using the entire header.", - DataWarning) - hdu_header = fits.Header(img_hdu.header[:max_ind], copy=True) - img_wcs = wcs.WCS(hdu_header) + hdu_header = fits.Header(img_hdu.header, copy=True) + + # We are going to reroute the logging to a string stream temporarily so we can + # intercept any message from astropy, chiefly the "Inconsistent SIP distortion information" + # INFO message which will indicate that we need to remove existing SIP keywords + # from a WCS whose CTYPE does not include SIP. In this we are taking the CTYPE to be + # correct and adjusting the header keywords to match. + with log.log_to_list() as log_list: + hdlr = log.handlers[0] + log.removeHandler(log.handlers[0]) + + img_wcs = wcs.WCS(hdu_header, relax=True) + + log.addHandler(hdlr) + + no_sip = False + if (len(log_list) > 0): + if ("Inconsistent SIP distortion information" in log_list[0].msg): + # Delete standard sip keywords + remove_sip_coefficients(hdu_header) + + # load wcs ignoring any nonstandard keywords + img_wcs = wcs.WCS(hdu_header, relax=False) + + # As an extra precaution make sure the img wcs has no sip coeefficients + img_wcs.sip = None + no_sip = True + + else: # Message(s) we didn't prepare for we want to go ahead and display + for log_rec in log_list: + log.log(log_rec.levelno, log_rec.msg, extra={"origin": log_rec.name}) + img_data = img_hdu.data if verbose: @@ -223,7 +262,10 @@ def _hducut(img_hdu, center_coord, cutout_size, correct_wcs=False, drop_after=No cutout_wcs = _get_cutout_wcs(img_wcs, cutout_lims) # Updating the header with the new wcs info - hdu_header.update(cutout_wcs.to_header(relax=True)) # relax arg is for sip distortions if they exist + if no_sip: + hdu_header.update(cutout_wcs.to_header(relax=False)) + else: + hdu_header.update(cutout_wcs.to_header(relax=True)) # relax arg is for sip distortions if they exist # Naming the extension hdu_header["EXTNAME"] = "CUTOUT" @@ -260,7 +302,11 @@ def _save_single_fits(cutout_hdus, output_path, center_coord): ('DEC_OBJ', center_coord.dec.deg, '[deg] declination')]) cutout_hdulist = fits.HDUList([primary_hdu] + cutout_hdus) - cutout_hdulist.writeto(output_path, overwrite=True, checksum=True) + + # Writing out the hdu often causes a warning as the ORIG_FLE card description is truncated + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + cutout_hdulist.writeto(output_path, overwrite=True, checksum=True) def _save_multiple_fits(cutout_hdus, output_paths, center_coord): @@ -330,8 +376,8 @@ def _parse_size_input(cutout_size): return cutout_size - -def fits_cut(input_files, coordinates, cutout_size, correct_wcs=False, drop_after=None, +@deprecated_renamed_argument('drop_after', None, '0.5') +def fits_cut(input_files, coordinates, cutout_size, correct_wcs=False, drop_after=None, single_outfile=True, cutout_prefix="cutout", output_dir='.', verbose=False): """ Takes one or more fits files with the same WCS/pointing, makes the same cutout in each file, @@ -357,12 +403,6 @@ def fits_cut(input_files, coordinates, cutout_size, correct_wcs=False, drop_afte correct_wcs : bool Default False. If true a new WCS will be created for the cutout that is tangent projected and does not include distortions. - drop_after : str or None - Default None. When creating the header for the cutout (and crucially, before - building the WCS object) drop all header keywords starting with the one given. This is - useful particularly for drizzle files that contain a multitude of extranious keywords - and sometimes leftover WCS keywords that astropy will try to parse even thought they should be - ignored. single_outfile : bool Default True. If true return all cutouts in a single fits file with one cutout per extension, if False return cutouts in individual fits files. If returing a single file the filename will @@ -402,11 +442,12 @@ def fits_cut(input_files, coordinates, cutout_size, correct_wcs=False, drop_afte for in_fle in input_files: if verbose: print("\n{}".format(in_fle)) - + + warnings.filterwarnings("ignore", category=wcs.FITSFixedWarning) with fits.open(in_fle) as hdulist: try: cutout = _hducut(hdulist[0], coordinates, cutout_size, - correct_wcs=correct_wcs, drop_after=drop_after, verbose=verbose) + correct_wcs=correct_wcs, verbose=verbose) except OSError as err: warnings.warn("Error {} encountered when performing cutout on {}, skipping...".format(err, in_fle), DataWarning) @@ -539,7 +580,7 @@ def normalize_img(img_arr, stretch='asinh', minmax_percent=None, minmax_value=No return norm_img - +@deprecated_renamed_argument('drop_after', None, '0.5') def img_cut(input_files, coordinates, cutout_size, stretch='asinh', minmax_percent=None, minmax_value=None, invert=False, img_format='jpg', colorize=False, cutout_prefix="cutout", output_dir='.', drop_after=None, verbose=False): @@ -632,7 +673,7 @@ def img_cut(input_files, coordinates, cutout_size, stretch='asinh', minmax_perce print("\n{}".format(in_fle)) hdulist = fits.open(in_fle) cutout = _hducut(hdulist[0], coordinates, cutout_size, - correct_wcs=False, drop_after=drop_after, verbose=verbose) + correct_wcs=False, verbose=verbose) hdulist.close() # We just want the data array From 5aa79385023120c90a57ddeb17c3e11671b99bb9 Mon Sep 17 00:00:00 2001 From: Clara Brasseur Date: Thu, 7 May 2020 19:05:02 -0400 Subject: [PATCH 2/4] pep8 and changelog --- CHANGES.rst | 1 + astrocut/cutouts.py | 23 ++++++++++++++--------- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 90fbf533..3e1e48fd 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -2,6 +2,7 @@ ---------------- - Update wcs fitting to match Astropy (and use Astropy when available) [#29] - Limit the number of pixels used for WCS fitting to 100 [#30] +- Deprecate drop_after and handle inconsistant wcs keywords automatically [#31] 0.5 (2020-01-13) diff --git a/astrocut/cutouts.py b/astrocut/cutouts.py index be91c8c3..0e57344a 100644 --- a/astrocut/cutouts.py +++ b/astrocut/cutouts.py @@ -129,25 +129,29 @@ def remove_sip_coefficients(hdu_header): """ Remove standard sip coefficient keywords for a fits header. - Everything is wrapped in try catches because we don't care is a keyword is not actually there. + Parameters + ---------- + hdu_header : ~astropy.io.fits.Header + The header from which SIP keywords will be removed. This is done in place. """ - for lets in product(["A","B"], ["","P"]): + # Everything is wrapped in try catches because if a keyword is missing we just move on. + for lets in product(["A", "B"], ["", "P"]): lets = ''.join(lets) try: - del hdu_header[f"{lets}_ORDER"] + del hdu_header["{}_ORDER".format(lets)] except KeyError: pass try: - del hdu_header[f"{lets}_DMAX"] + del hdu_header["{}_DMAX".format(lets)] except KeyError: pass - for i,j in product([0,1,2,3],[0,1,2,3]): + for i, j in product([0, 1, 2, 3], [0, 1, 2, 3]): try: - del hdu_header[f"{lets}_{i}_{j}"] + del hdu_header["{}_{}_{}".format(lets, i, j)] except KeyError: pass @@ -208,7 +212,7 @@ def _hducut(img_hdu, center_coord, cutout_size, correct_wcs=False, verbose=False img_wcs.sip = None no_sip = True - else: # Message(s) we didn't prepare for we want to go ahead and display + else: # Message(s) we didn't prepare for we want to go ahead and display for log_rec in log_list: log.log(log_rec.levelno, log_rec.msg, extra={"origin": log_rec.name}) @@ -376,7 +380,7 @@ def _parse_size_input(cutout_size): return cutout_size -@deprecated_renamed_argument('drop_after', None, '0.5') +@deprecated_renamed_argument('drop_after', None, '0.6') def fits_cut(input_files, coordinates, cutout_size, correct_wcs=False, drop_after=None, single_outfile=True, cutout_prefix="cutout", output_dir='.', verbose=False): """ @@ -579,7 +583,8 @@ def normalize_img(img_arr, stretch='asinh', minmax_percent=None, minmax_value=No norm_img = 255 - norm_img return norm_img - + + @deprecated_renamed_argument('drop_after', None, '0.5') def img_cut(input_files, coordinates, cutout_size, stretch='asinh', minmax_percent=None, minmax_value=None, invert=False, img_format='jpg', colorize=False, From 3c1ce0cc6b5c8212c0e71b6ec27224812c07e16a Mon Sep 17 00:00:00 2001 From: Clara Brasseur Date: Thu, 7 May 2020 19:55:16 -0400 Subject: [PATCH 3/4] changing deprecation for backwards compatibility --- astrocut/cutouts.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/astrocut/cutouts.py b/astrocut/cutouts.py index 0e57344a..6c54a12b 100644 --- a/astrocut/cutouts.py +++ b/astrocut/cutouts.py @@ -15,7 +15,7 @@ from astropy.io import fits from astropy.coordinates import SkyCoord from astropy import wcs -from astropy.utils.decorators import deprecated_renamed_argument +from astropy.utils.exceptions import AstropyDeprecationWarning from astropy.visualization import (SqrtStretch, LogStretch, AsinhStretch, SinhStretch, LinearStretch, MinMaxInterval, ManualInterval, AsymmetricPercentileInterval) @@ -379,8 +379,7 @@ def _parse_size_input(cutout_size): return cutout_size - -@deprecated_renamed_argument('drop_after', None, '0.6') + def fits_cut(input_files, coordinates, cutout_size, correct_wcs=False, drop_after=None, single_outfile=True, cutout_prefix="cutout", output_dir='.', verbose=False): """ @@ -427,6 +426,11 @@ def fits_cut(input_files, coordinates, cutout_size, correct_wcs=False, drop_afte the output filepaths. """ + # Dealing with deprecation + if drop_after is not None: + warnings.warn("Argument 'drop_after' is deprecated and will be ignored", + AstropyDeprecationWarning) + if verbose: start_time = time() @@ -585,7 +589,6 @@ def normalize_img(img_arr, stretch='asinh', minmax_percent=None, minmax_value=No return norm_img -@deprecated_renamed_argument('drop_after', None, '0.5') def img_cut(input_files, coordinates, cutout_size, stretch='asinh', minmax_percent=None, minmax_value=None, invert=False, img_format='jpg', colorize=False, cutout_prefix="cutout", output_dir='.', drop_after=None, verbose=False): @@ -644,7 +647,11 @@ def img_cut(input_files, coordinates, cutout_size, stretch='asinh', minmax_perce the output filepaths. """ - + # Dealing with deprecation + if drop_after is not None: + warnings.warn("Argument 'drop_after' is deprecated and will be ignored", + AstropyDeprecationWarning) + if verbose: start_time = time() From de1f5ce9a0b5870b991d6138957b5bc8684334d2 Mon Sep 17 00:00:00 2001 From: Clara Brasseur Date: Fri, 8 May 2020 10:29:03 -0400 Subject: [PATCH 4/4] making log capture more robust --- astrocut/cutouts.py | 32 ++++++++++++++------------------ 1 file changed, 14 insertions(+), 18 deletions(-) diff --git a/astrocut/cutouts.py b/astrocut/cutouts.py index 6c54a12b..d2e1821c 100644 --- a/astrocut/cutouts.py +++ b/astrocut/cutouts.py @@ -124,7 +124,6 @@ def _get_cutout_wcs(img_wcs, cutout_lims): return wcs.WCS(wcs_header) -#### FUNCTIONS FOR UTILS #### def remove_sip_coefficients(hdu_header): """ Remove standard sip coefficient keywords for a fits header. @@ -135,26 +134,23 @@ def remove_sip_coefficients(hdu_header): The header from which SIP keywords will be removed. This is done in place. """ - # Everything is wrapped in try catches because if a keyword is missing we just move on. for lets in product(["A", "B"], ["", "P"]): lets = ''.join(lets) - - try: + + key = "{}_ORDER".format(lets) + if key in hdu_header.keys(): del hdu_header["{}_ORDER".format(lets)] - except KeyError: - pass - try: + key = "{}_DMAX".format(lets) + if key in hdu_header.keys(): del hdu_header["{}_DMAX".format(lets)] - except KeyError: - pass for i, j in product([0, 1, 2, 3], [0, 1, 2, 3]): - try: + key = "{}_{}_{}".format(lets, i, j) + if key in hdu_header.keys(): del hdu_header["{}_{}_{}".format(lets, i, j)] - except KeyError: - pass +#### FUNCTIONS FOR UTILS #### def _hducut(img_hdu, center_coord, cutout_size, correct_wcs=False, verbose=False): """ @@ -191,13 +187,13 @@ def _hducut(img_hdu, center_coord, cutout_size, correct_wcs=False, verbose=False # INFO message which will indicate that we need to remove existing SIP keywords # from a WCS whose CTYPE does not include SIP. In this we are taking the CTYPE to be # correct and adjusting the header keywords to match. - with log.log_to_list() as log_list: - hdlr = log.handlers[0] - log.removeHandler(log.handlers[0]) - + hdlrs = log.handlers + log.handlers = [] + with log.log_to_list() as log_list: img_wcs = wcs.WCS(hdu_header, relax=True) - - log.addHandler(hdlr) + + for hd in hdlrs: + log.addHandler(hd) no_sip = False if (len(log_list) > 0):