Skip to content

Commit

Permalink
standardizeWCS() uses P and QOFFSETs rather than RA and DEC because o…
Browse files Browse the repository at this point in the history 
…f problems. checkWCS() calls standardizeWCS()
  • Loading branch information
@chris-simpson
chris-simpson committed Dec 17, 2024
1 parent 8a68dd5 commit d48e740
Show file tree
Hide file tree
Showing 2 changed files with 66 additions and 64 deletions.
2 changes: 1 addition & 1 deletion geminidr/gemini/parameters_gemini.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,4 +29,4 @@ class standardizeWCSConfig(parameters_standardize.standardizeWCSConfig):


class checkWCSConfig(config.Config):
tolerance = config.RangeField("Positional tolerance (arcsec)", float, 2, min=0.1)
pass
128 changes: 65 additions & 63 deletions geminidr/gemini/primitives_gemini.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
# primitives_gemini.py
# ------------------------------------------------------------------------------
import datetime
from contextlib import suppress
from copy import deepcopy
import numpy as np

Expand All @@ -12,7 +13,7 @@
from astropy.modeling import models

from gempy.gemini import gemini_tools as gt
from gempy.library import astrotools as at
from gempy.library import astromodels as am, astrotools as at

from geminidr.core import Bookkeeping, CalibDB, Preprocess
from geminidr.core import Visualize, Standardize, Stack
Expand Down Expand Up @@ -223,7 +224,7 @@ def standardizeWCS(self, adinputs=None, **params):
limit = params["debug_consistency_limit"]
max_deadtime = params["debug_max_deadtime"]

want_to_fix = bad_wcs not in ('exit', 'ignore')
want_to_fix = bad_wcs not in ('check', 'exit', 'ignore')
bad_wcs_list = []
base_pointing = None
last_pointing = None
Expand Down Expand Up @@ -260,9 +261,12 @@ def standardizeWCS(self, adinputs=None, **params):
base_pointing = None
last_pointing = None

p = Pointing(ad)
needs_fixing = (bad_wcs == 'new' or
not p.self_consistent(limit=limit))
p = Pointing(ad, log_function=self.log.debug)
selfcon = p.self_consistent(limit=limit)
if not selfcon:
log.stdinfo(f"{ad.filename}: WCS inconsistent with target "
"coordinates")
needs_fixing = (bad_wcs == 'new' or not selfcon)
if base_pointing is not None:
needs_fixing |= not base_pointing.consistent_with(p)

Expand All @@ -273,6 +277,7 @@ def standardizeWCS(self, adinputs=None, **params):
for ext, wcs in zip(ad, p.wcs):
ext.wcs = wcs
base_pointing = p
log.stdinfo(f"Using {ad.filename} as base pointing")
elif not want_to_fix or base_pointing is None:
# Do not want to, or cannot yet, fix
bad_wcs_list.append(ad)
Expand All @@ -293,6 +298,7 @@ def standardizeWCS(self, adinputs=None, **params):
if base_pointing is None:
# Found a reliable base WCS
base_pointing = p
log.stdinfo(f"Using {ad.filename} as base pointing")
# Fix all backed up ADs if we want to
if want_to_fix:
while bad_wcs_list:
Expand All @@ -316,15 +322,18 @@ def standardizeWCS(self, adinputs=None, **params):
last_obsid = this_obsid

if not want_to_fix and bad_wcs_list:
log.stdinfo("The following files were identified as having bad "
"WCS information:")
log.stdinfo("\nThe following files were identified as having bad"
" WCS information:")
for ad in bad_wcs_list:
log.stdinfo(f" {ad.filename}")
if base_pointing is None:
log.stdinfo("\nNo valid base pointing was identified, so a "
"'new' WCS is required, 'fix' will not work.")
if bad_wcs == 'exit':
raise ValueError("Some files have bad WCS information and "
"user has requested an exit")
else:
log.stdinfo("This is being ignored as requested.")
log.stdinfo("No changes are being made to the WCS information.")

for ad in adinputs:
# Timestamp and update filename
Expand All @@ -333,15 +342,10 @@ def standardizeWCS(self, adinputs=None, **params):

return adinputs

def checkWCS(self, adinputs=None, tolerance=None):
def checkWCS(self, adinputs=None):
"""
This primitive checks for consistency within the WCS by comparing the
header offsets with the WCS coordinates
Parameters
----------
tolerance: float
positional tolerance in arcseconds
This primitive checks for consistency within the WCS by calling
standardizeWCS().
"""
log = self.log
log.debug(gt.log_message("primitive", self.myself(), "starting"))
Expand All @@ -351,32 +355,12 @@ def checkWCS(self, adinputs=None, tolerance=None):
if any('PREPARED' in ad.tags for ad in adinputs):
raise ValueError(f"'{self.myself()}' requires unprepared data")

adref = adinputs[0]
refindex = (len(adref) - 1) // 2
refpixel = tuple(length // 2 for length in adref[refindex].shape[::-1])

bad = []
for i, ad in enumerate(adinputs):
offsets = ad.detector_x_offset(), ad.detector_y_offset()
pixel = tuple(r + o for r, o in zip(refpixel, offsets)) + refpixel[2:]
# Be wary of third F2 axis
coord = SkyCoord(*ad[refindex].wcs(*pixel)[:2], unit='deg')
if i == 0:
log.stdinfo(f"Using {ad.filename} as the reference")
ref_coord = coord
else:
sep = ref_coord.separation(coord).arcsec
if sep > tolerance:
bad.append((ad.filename, sep))

for (fname, sep) in bad:
log.stdinfo(f"{fname} has a discrepancy of {sep:.2f} arcsec")

if len(bad) > len(adinputs) // 2:
log.stdinfo(f"The first frame ({adref.filename}) is likely to "
"be the one in error")
elif not bad:
log.stdinfo("All files appear to be OK")
# We don't want to modify the filename or PHU in case such changes
# cause the outputs to be written to disk
Gemini([]).standardizeWCS(adinputs, bad_wcs='ignore', suffix=None)
for ad in adinputs:
with suppress(KeyError):
del ad.phu[self.timestamp_keys['standardizeWCS']]

return adinputs

Expand All @@ -391,20 +375,39 @@ class Pointing:
# (x, y), 0-indexed
center_of_rotation_dict = {'GNIRS': (629., 519.)}

def __init__(self, ad):
def __init__(self, ad, log_function=None):
self.phu = ad.phu.copy()
self.instrument = ad.instrument()
self.wcs = [ext.wcs for ext in ad]
self.target_coords = SkyCoord(ad.target_ra(), ad.target_dec(),
unit=u.deg)
self.coords = SkyCoord(ad.wcs_ra(), ad.wcs_dec(), unit=u.deg)
self.expected_coords = self.target_coords.spherical_offsets_by(
self.phu['RAOFFSET']*u.arcsec, self.phu['DECOFFSE']*u.arcsec)
self.pa = ad.position_angle()
self.xoffset = ad.detector_x_offset()
self.yoffset = ad.detector_y_offset()
self.pixel_scale = ad.pixel_scale()
self.filename = ad.filename
self.wcs = [ext.wcs for ext in ad]
self.logit = log_function or print

self.logit(f"{self.filename} {self.xoffset} {self.yoffset} {self.pa}")
self.target_coords = SkyCoord(ad.target_ra(), ad.target_dec(),
unit=u.deg)
# This returns the CRVALi keywords, which are expected to change
# because the CRPIXj keywords are expected to be constant.
self.coords = SkyCoord(ad.wcs_ra(), ad.wcs_dec(), unit=u.deg)

# We want to know whether East is to the left or right when North
# is up, which may depend on the port the instrument is on.
# Unfortunately, we can't trust the CD matrix because values close
# to zero (when the instrument is oriented along one of the cardinal
# axes) can have random signs. So we need to do this empirically.
self.flipped = (self.instrument == "GNIRS" or
self.phu.get("INPORT") == 1 and (
"GMOS" in self.instrument or
self.instrument == "NIRI" and not ad.is_ao()))
# If an offset is positive, it means that an object's pixel coords
# will *decrease*, so the telescope has moved either up or right.
self.expected_coords = self.target_coords.directional_offset_by(
self.pa * u.deg, -self.yoffset * self.pixel_scale * u.arcsec).directional_offset_by(
(self.pa + 90 * (-1, 1)[self.flipped]) * u.deg, -self.xoffset * self.pixel_scale * u.arcsec
)

def __repr__(self):
return f"Pointing object from {self.filename}"
Expand All @@ -426,7 +429,14 @@ def self_consistent(self, limit=10):
-------
bool: is the pointing self-consistent?
"""
return self.coords.separation(self.expected_coords).arcsec <= limit
sep = self.coords.separation(self.expected_coords).arcsec
self.logit(f"Self-consistency report on {self.filename} "
f"(flipped={self.flipped})")
self.logit(f"Base coords: {self.coords.to_string(precision=5)}")
self.logit("Expected coords: "
f"{self.expected_coords.to_string(precision=5)}")
self.logit(f"Separation: {sep:.2f} arcsec")
return sep <= limit

def consistent_with(self, other):
"""
Expand Down Expand Up @@ -455,8 +465,8 @@ def consistent_with(self, other):
ra, dec = at.get_center_of_projection(wcs1)
except TypeError: # if this returns None
return False
x, y = wcs1.invert(ra, dec)
x2, y2 = wcs2.invert(ra, dec)
x, y = am.get_named_submodel(wcs1.forward_transform, 'SKY').inverse(ra, dec)
x2, y2 = am.get_named_submodel(wcs2.forward_transform, 'SKY').inverse(ra, dec)
dx = other.xoffset - self.xoffset
dy = other.yoffset - self.yoffset
distsq = dx * dx + dy * dy
Expand Down Expand Up @@ -486,19 +496,11 @@ def fix_pointing(self):
# Update projection center with coordinates
nat2cel.lon = self.expected_coords.ra.value
nat2cel.lat = self.expected_coords.dec.value
# Try to construct a new CD matrix. Whether East is to the left or
# right when North is up may depend on the port the instrument is
# on, but we can assume that the matrix in the header is correct
# in this regard since the instrument won't have switched ports.
# Not true for GNIRS where there was a sign error, but GNIRS is
# always flipped (East is to the right if North is up).
flipped = (aftran.matrix[0, 0] * aftran.matrix[1, 1] > 0 or
aftran.matrix[0, 1] * aftran.matrix[1, 0] < 0 or
self.instrument == "GNIRS")
# Try to construct a new CD matrix.
new_matrix = self.pixel_scale / 3600 * np.identity(2)
pa = -self.pa if flipped else self.pa
pa = -self.pa if self.flipped else self.pa
new_matrix = np.asarray(models.Rotation2D(pa)(*new_matrix))
if not flipped:
if not self.flipped:
new_matrix[0] *= -1
aftran.matrix = new_matrix

Expand Down

0 comments on commit d48e740

Please sign in to comment.