Commit 6d22c1fa authored by Carl Schaffer's avatar Carl Schaffer
Browse files

refactoring IFU wcs generator for maintainability

parent d6bf1d4e
......@@ -49,6 +49,7 @@ class GrisIFUTranslator(GrisTranslator):
return False
def to_EXTNAME(self):
"""Generate extension name (date_run_map)"""
iserie = int(self._header["ISERIE"])
value = self.to_observation_id()
value += f"_{iserie:03d}"
......
......@@ -9,3 +9,11 @@ class WCSGenerator(ABC):
@abstractmethod
def reordered_data(self):
"""Generate WCS cards for the underlying file"""
@abstractmethod
def data_shape(self):
"""Return the shape of the data array after translation"""
@abstractmethod
def rotate_reference_coords(self, X, Y, angle):
"""Rotate a set of reference coords by a given angle"""
......@@ -180,7 +180,7 @@ class GrisWCSGenerator(WCSGenerator):
the SOLARX axis. They the also contain a rotation matrix to orient the slit correctly within the
coordinate system. This is necessary to make the specification of the values more intuitive."""
shape = self.reordered_data().T.shape
shape = self.data_shape()
def get_angle(x, y):
v = np.array([x, y])
......@@ -193,7 +193,7 @@ class GrisWCSGenerator(WCSGenerator):
X, Y = self.coords
angle = get_angle(X[0] - X[-1], Y[0] - Y[-1])
X_rot, Y_rot = rotate_around_first_coord(X, Y, angle)
X_rot, Y_rot = self.rotate_reference_coords(X, Y, angle)
assert np.isclose(0, sum(X_rot - np.mean(X_rot)), 0.01)
cards = []
......@@ -252,7 +252,17 @@ class GrisWCSGenerator(WCSGenerator):
return cards
def rotate_reference_coords(self, X, Y, angle):
"""
Rotate the reference coords by a given angle in radians
for slit data we can directly rotate the x and y coords of the slit.
"""
X_rot, Y_rot = rotate_around_first_coord(X, Y, angle)
return X_rot, Y_rot
def data_shape(self):
shape = self.reordered_data().T.shape
return shape
class GrisIFUWCSGenerator(GrisWCSGenerator):
......@@ -265,13 +275,11 @@ class GrisIFUWCSGenerator(GrisWCSGenerator):
self.coords = self.get_coords_from_header()
def get_coords_from_header(self):
"""coordinates for slit, the Keywords 'SLITPOSX' and 'SLITPOSY' are assumed to be the center of the slit
while the keyword 'SLITORIE' describes the rotation of the slit w.r.t. the helioprojective axes. The
algorithm further assumes that each pixel represents a square area with side length of 'STEPSIZE'.
These assumptions are used to calculate the coordinates for the center of each pixel within the data array.
"""
Extract coordinates from the underlying data.
Returns:
(X,Y):ndarray(float) x and y coordinates for each pixel of the slit
(X,Y):ndarray(float) x and y coordinates in HPC/Arcsec for the leftmost column and topmost row respectively
"""
coords = self.infile.coords.value
......@@ -292,84 +300,21 @@ class GrisIFUWCSGenerator(GrisWCSGenerator):
data = self.infile.data
return data
def make_wcs_cards(self):
"""Generate WCS keywords. The spatial extent is represented by introducing a degenerate axis.
The values specified in the header contain axes of SOLARX and SOLARY rotated to be parallel with
the SOLARX axis. They the also contain a rotation matrix to orient the slit correctly within the
coordinate system. This is necessary to make the specification of the values more intuitive."""
shape = self.reordered_data().T.shape
def get_angle(x, y):
v = np.array([x, y])
u = np.array([0, 1])
temp = np.dot(v, u) / (math.sqrt(np.dot(v, v)) * math.sqrt(np.dot(u, u)))
angle = math.acos(temp)
angle = -1 * angle if x < 0 else angle
return angle
X, Y = self.coords
angle = get_angle(X[0] - X[-1], Y[0] - Y[-1])
def rotate_reference_coords(self, X, Y, angle):
"""
Rotate the reference coords by a given angle in radians
for IFU data we need to account for different sitzes in x and y direction.
Essentially the code below is meant to take a horizontal slice with varying X values and a vertical slice with
varying Y values, pad each with a constant value of the other coordinate and rotate them around the first
coordinate.
"""
X, Y_placeholder = [X, np.ones(X.shape) * Y.mean()]
X_placeholder, Y = [np.ones(Y.shape) * X.mean(), Y]
X_rot, _ = rotate_around_first_coord(X, Y_placeholder, angle)
_, Y_rot = rotate_around_first_coord(X_placeholder, Y, angle)
assert np.isclose(0, sum(X_rot - np.mean(X_rot)), 0.01)
cards = []
# COORDINATE TYPES
c_types = ["HPLN-TAN", "HPLT-TAN", "WAVE", "STOKES"]
# UNITS
c_units = ["arcsec", "arcsec", "Angstrom", ""]
# Reference Pixels
# CRPIX values are 1-indexed according to the fits standard
c_rpix = [1, 1, 1, 1]
return X_rot, Y_rot
# REFERENCE VALUES
# the coordinates are rotated in such a way that they are parallel to the y-axis
# which means that the y- value is the same along the entire axis. The information on the rotation is stored
# in the PCi_j keywords
c_rval = [X_rot[0], Y_rot[0], self.infile.wavelength_region[0], 1]
# COORDINATE UNCERTAINTIES
# set systematic error: Use estimated 2 arcsec if cross-correlation is successful, else use one solar disk
# with 1000" because the only thing we can be sure of with current GREGOR coordinates is that we are indeed
# pointing at the sun -.-
has_ccorr = bool(self.cross_correlation_file)
# Uncertainties estimated from difference between cross correlated coordinates and coordinates
# from the telescope. They showed a significantly larger uncertainty for the Y coordinate. See documentation of
# https://gitlab.leibniz-kis.de/sdc/sdc_importer/issues/199
spatial_uncertainty_x = 5 if has_ccorr else self.coord_uncertainties[0]
spatial_uncertainty_y = 5 if has_ccorr else self.coord_uncertainties[1]
c_syer = [
spatial_uncertainty_x,
spatial_uncertainty_y,
self.infile.wavelength_step_size,
0,
]
# STEP SIZES
step_size = np.diff(Y_rot).mean()
c_delt = [step_size, step_size, self.infile.wavelength_step_size, 1]
for i in range(len(shape)):
cards.append(self.make_card(f"NAXIS{i + 1}", shape[i]))
cards.append(self.make_card(f"CTYPE{i + 1}", c_types[i]))
cards.append(self.make_card(f"CUNIT{i + 1}", c_units[i]))
cards.append(self.make_card(f"CRPIX{i + 1}", c_rpix[i]))
cards.append(self.make_card(f"CRVAL{i + 1}", c_rval[i]))
cards.append(self.make_card(f"CDELT{i + 1}", c_delt[i]))
cards.append(self.make_card(f"CSYER{i + 1}", c_syer[i]))
# ROTATION MATRIX
# add rotation of spatial axes
c, s = np.cos(-angle), np.sin(-angle)
cards.append(self.make_card("PC1_1", c))
cards.append(self.make_card("PC1_2", -s))
cards.append(self.make_card("PC2_1", s))
cards.append(self.make_card("PC2_2", c))
return cards
def data_shape(self):
shape = self.reordered_data().T.shape
return shape
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment