From e51bd84ba2002d0c7485df582d476c7dc797d2ed Mon Sep 17 00:00:00 2001 From: Clara Brasseur Date: Thu, 7 May 2020 18:45:53 -0400 Subject: [PATCH] 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