# -*- coding: utf-8 -*-
#pylint: disable=C0103,W0311
'''
brutifus: a set of Python modules to process datacubes from integral field spectrographs.\n
Copyright (C) 2018-2019, F.P.A. Vogt
----------------------------------------------------------------------------------------------------
This file contains WCS-related function for the higher-leve brutifus routines.
Created August 2019, F.P.A. Vogt - frederic.vogt@alumni.anu.edu.au
'''
# --------------------------------------------------------------------------------------------------
import warnings
import os
import numpy as np
from scipy import signal
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
from astropy.io import fits
import astropy.units as u
from astropy.coordinates.sky_coordinate import SkyCoord
from astropy.stats import sigma_clipped_stats
from astropy.time import Time
from astropy.wcs import WCS
from astroquery.gaia import Gaia
from photutils import DAOStarFinder
from . import brutifus_metadata as bifus_m
from . import brutifus_plots as bifus_p
# --------------------------------------------------------------------------------------------------
[docs]def get_linear_WCS_corr(fn, obstime=None, verbose=True, suffix=None, target=None):
''' Identify the linear offset between Gaia and the WCS solution of a given 2D image.
:param fn: 2D FITS image filename (and location)
:type fn: string
:param obstime: Observing time and date
:type obstime: astropy.Time()
:param suffix: the tag of this step (for plot filename)
:type suffix: str
:param target: target name (for plot filename)
:type target: str
:param verbose: print some info to screen
:type verbose: bool
:return: (dx,dy), the linear offsets in pixels
:rtype: tuple
.. note:: Assumes the image is in the primary extension.
'''
# Get the data
hdu = fits.open(fn)
header = hdu[0].header
wl_im = hdu[0].data
hdu.close()
# First, let's get some of the statistics on this image
wl_mean, wl_median, wl_std = sigma_clipped_stats(wl_im, sigma=3.0)
if verbose:
print(' Looking for stars in the field ...')
# Source: https://photutils.readthedocs.io/en/stable/detection.html
# Get ready to look for stars in there
daofind = DAOStarFinder(fwhm=3.0, threshold=5. * wl_std)
# Start the search
sources = daofind(wl_im - wl_median)
# Create a list of sources in pixels
sources_list = np.array([sources['xcentroid'], sources['ycentroid']])
# Reshape this properly as Nx2
sources_list = sources_list.T
# Store the 2D WCS info for later
wcs2d = WCS(header)
# Find the coordinate at the center of the image.
# To be robust (E.g. OCam), should be CRVAl 1/2 + CDELT1/2*(0.5*NAXIS 1/2 - CRPIX1/2)
# Fixed v. 2019.08.1
# Credits to J. Suherli for stumbling upon and reporting the bug.
core_coord = SkyCoord(ra=header['CRVAL1'] +
header['CD1_1'] * (0.5*header['NAXIS1'] - header['CRPIX1']),
dec=header['CRVAL2'] +
header['CD2_2'] * (0.5*header['NAXIS2'] - header['CRPIX2']),
unit=(u.deg, u.deg), frame='icrs')
search_radius = 0.75 * max([np.abs(header['NAXIS1'] * header['CD1_1']*u.deg),
np.abs(header['NAXIS2'] * header['CD2_2']*u.deg)])
if verbose:
print(' Querying GAIA ...')
# Make it a sync search, because I don't think I need more than 2000 targets ...
j = Gaia.cone_search_async(core_coord, search_radius, verbose=False)
r = j.get_results()
if len(r) >= 2000:
warnings.warn(' Maximum number of GAIA entries returned. Some stars might be missing.')
# Some of these proper motions are NaN's ... need to deal with this before feeding them
# to SkyCoord
r['pmra'].fill_value = 0
r['pmdec'].fill_value = 0
# Turn this into a "catalogue" of SkyCoord
gaia_cat = SkyCoord(ra=r['ra'], dec=r['dec'], unit=(u.deg, u.deg), frame='icrs',
pm_ra_cosdec=r['pmra'].filled(),
pm_dec=r['pmdec'].filled(),
obstime=[Time(item, format='decimalyear') for item in r['ref_epoch']],
equinox='J2000',
distance=1.e3 * u.pc, # Assumes the default distance
)
# Propagate Gaia to the observation date!
if not obstime is None:
if verbose:
print(' Propagating proper motions to %s ...' % (obstime))
gaia_cat = gaia_cat.apply_space_motion(new_obstime=obstime)
# Clean this up, to only show the points within the image
gaia_cat = gaia_cat[wcs2d.footprint_contains(gaia_cat)]
# Transform this in pixel coords
gaia_cat_pix = wcs2d.wcs_world2pix(gaia_cat.ra, gaia_cat.dec, 0)
# Reshape it properly
gaia_cat_pix = np.array(gaia_cat_pix).T
# Ok, let's make two fake image from the source list, that I then cross-correlate
nx = header['NAXIS1']
ny = header['NAXIS2']
im_ref = np.zeros((ny, nx))
im_obs = np.zeros((ny, nx))
# Fill the two images with fake stars
# v2019.08.2: this gets slow with large numbers of stars
# (thanks Janette Suherli for reporting this)
# let's be smarter about this ... only create stars on small sub-arrays,
# and add these to the final image (rather than creating full images each time).
# Note: the FWHM of the stars in the fake images is fixed in pixels ... should I change this ?
fwhm_pix = 2
# Let's get ready to fill the images.
for (im, entries) in [(im_obs, sources_list), (im_ref, gaia_cat_pix)]:
# Let's create the fake observed image ...
for i in range(len(entries)):
# Pixel ID closest from the star centroid
xc = np.floor(entries[i][0])
yc = np.floor(entries[i][1])
# If we are too close from the edge, skip the star
# (otherwise I have array size issues)
# That's the easy way out of course ... !
if (xc < 2 * fwhm_pix + 3) or (xc > nx - (2 * fwhm_pix + 3)) or \
(yc < 2 * fwhm_pix + 3) or (yc > ny - (2 * fwhm_pix + 3)):
continue
# Find the x and y indices of the sub-array pixels
x_ind, y_ind = np.meshgrid(np.arange(xc-(2 * fwhm_pix + 1), xc+(2 * fwhm_pix + 1), 1),
np.arange(yc-(2 * fwhm_pix + 1), yc+(2 * fwhm_pix + 1), 1))
# For each of these, compute the distance to the star center
d = np.sqrt((x_ind - entries[i][0])**2 + (y_ind - entries[i][1])**2)
# And now add the star to the full image
im[y_ind.astype(int), x_ind.astype(int)] += np.exp(-((d)**2 / (2.0 * fwhm_pix**2)))
# Having done that, let's convolve the two images
# Credit to Nick Whyborn @ ALMA for the idea
# Don't forget to flip one of the image to get the correlation
im_corr = signal.convolve(im_ref, im_obs[::-1, ::-1], mode='same', method='auto')
# Now, let's find the brightest pixel
peak_loc = np.unravel_index(im_corr.argmax(), im_corr.shape)
# Let's also find the "fine" position of the peak with photutils
mean, median, std = sigma_clipped_stats(im_corr, sigma=10.0)
daofind = DAOStarFinder(fwhm=5.0, threshold=10.*std)
peaks = daofind(im_corr)
# Add a clean error, in case I find no peak.
if not peaks: # Neat trick ... empty sequences == False !
raise Exception('Ouch! No WCS correlation peak identified.')
# Ok, which one is the brightest peak ?
peak_id = np.argmax(peaks['peak'])
peak_loc_fine = (peaks['xcentroid'][peak_id], peaks['ycentroid'][peak_id])
# Issue a warning if the brightest peak is too far from the brightest pixel
if np.abs(peak_loc_fine[0] - peak_loc[1]) > 2 or \
np.abs(peak_loc_fine[1] - peak_loc[0]) > 2:
warnings.warn('WCS correlation peak does not coincide with the brightest pixel!')
# Ok, so what are the offsets I should apply ?
dx = peak_loc_fine[0]-np.floor(nx/2) # Added the floor... it is required by 'same' ?
dy = peak_loc_fine[1]-np.floor(ny/2) # v2019.07.1
if verbose:
print(' Best offset: %+.2f %+.2f pixels' % (dx, dy))
# Very well, now apply the corrections to the image
fits.setval(fn, 'CRPIX1', value=header['CRPIX1'] - dx, ext=0)
fits.setval(fn, 'CRPIX2', value=header['CRPIX2'] - dy, ext=0)
# Let's make a plot to show how good I'm doing with getting a reasonable offset
plt.close(99)
fig99 = plt.figure(99, figsize=(6.94, 7))
gs = gridspec.GridSpec(1, 1, height_ratios=[1], width_ratios=[1],
left=0.12, right=0.98, bottom=0.12, top=0.9,
wspace=0.1, hspace=0.1)
ax99 = plt.subplot(gs[0, 0])
ax99.imshow(im_corr, cmap='gray', origin='lower', interpolation='nearest')
ax99.plot(peak_loc_fine[0], peak_loc_fine[1], marker=bifus_p.crosshair(pa=0, inner_r=0.5),
c='darkorange', markersize=50)
ax99.set_aspect('equal')
ax99.set_xlabel('x [pix]')
ax99.set_ylabel('y [pix]')
ax99.set_title('Correlation peak: (%+.2f,%+.2f) spaxels' % (dx, dy))
fig99.savefig(os.path.join(bifus_m.plot_loc, suffix + '_' + target +
'_get-WCS-offsets.pdf'))
plt.close(99)
# Now plot the white-light image to show how good I did ...
(fig, ax, ofn) = bifus_p.make_2Dplot(fn, ext=0,
ofn=os.path.join(bifus_m.plot_loc, suffix + '_' + target +
'_check-WCS.pdf'),
plims=[10, 99],
cmap=None, cblabel=None)
# And add them to the plot
ax.scatter(gaia_cat.ra, gaia_cat.dec, marker='o', color='darkorange',
s=20, facecolor='none', transform=ax.get_transform('world'))
# Show the stars I found in the data ... these are in pix coordinates ...
# -> always "right" irrespective of the WCS!
ax.scatter(sources['xcentroid'], sources['ycentroid'],
color='crimson', s=50, marker='o', facecolor='none')
# Save the updated plot
fig.savefig(ofn)
return (dx, dy)
# --------------------------------------------------------------------------------------------------