Skip to content

Commit

Permalink
Merge 5aa7938 into ec7d7bd
Browse files Browse the repository at this point in the history
  • Loading branch information
C. E. Brasseur committed May 7, 2020
2 parents ec7d7bd + 5aa7938 commit 7334cda
Show file tree
Hide file tree
Showing 2 changed files with 90 additions and 43 deletions.
1 change: 1 addition & 0 deletions CHANGES.rst
View file
Expand Up @@ -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)
Expand Down
132 changes: 89 additions & 43 deletions astrocut/cutouts.py
View file
Expand Up @@ -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

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -126,9 +125,38 @@ 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.
Parameters
----------
hdu_header : ~astropy.io.fits.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:
del hdu_header["{}_ORDER".format(lets)]
except KeyError:
pass

try:
del hdu_header["{}_DMAX".format(lets)]
except KeyError:
pass

for i, j in product([0, 1, 2, 3], [0, 1, 2, 3]):
try:
del hdu_header["{}_{}_{}".format(lets, i, j)]
except KeyError:
pass

def _hducut(img_hdu, center_coord, cutout_size, correct_wcs=False, drop_after=None, verbose=False):

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,
Expand All @@ -146,13 +174,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.
Expand All @@ -161,18 +183,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:
Expand Down Expand Up @@ -223,7 +266,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"
Expand Down Expand Up @@ -260,7 +306,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):
Expand Down Expand Up @@ -330,8 +380,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.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):
"""
Takes one or more fits files with the same WCS/pointing, makes the same cutout in each file,
Expand All @@ -357,12 +407,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
Expand Down Expand Up @@ -402,11 +446,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)
Expand Down Expand Up @@ -538,8 +583,9 @@ 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,
cutout_prefix="cutout", output_dir='.', drop_after=None, verbose=False):
Expand Down Expand Up @@ -632,7 +678,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
Expand Down

0 comments on commit 7334cda

Please sign in to comment.