Coverage for python/lsst/ip/isr/intrinsicZernikes.py: 95%
127 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-09 09:26 +0000
1# This file is part of ip_isr.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21"""
22Intrinsic Zernikes storage class.
23"""
25__all__ = ["IntrinsicZernikes"]
27import numpy as np
28from astropy import units as u
29from astropy.table import Table
30from scipy.interpolate import LinearNDInterpolator
32from lsst.ip.isr import IsrCalib
35class IntrinsicZernikes(IsrCalib):
36 """Intrinsic Zernike coefficients.
38 Stores Zernike wavefront-error coefficients sampled at a set of
39 focal-plane field angles. At query time the coefficients are
40 interpolated to an arbitrary field position provided in CCS.
42 The coefficients are stored in two coordinate systems, each as an
43 independent set of sample points and values:
45 - the Camera Coordinate System (CCS), which corresponds to the
46 focal plane heights and any other CCS contribution and
47 - the Optical Coordinate System (OCS), corresponding to
48 the intrinsics, which are by nature defined in the optical
49 coordinates.
51 See `LSE-349 <https://ls.st/LSE-349>`_ for the definitions.
53 The CCS sample points and values are stored on the un-suffixed
54 ``field_x``, ``field_y`` and ``values`` attributes (and serialized
55 under those same keys). These names are kept unchanged from version
56 1, which only stored the CCS system, so that version 1 calibrations
57 round-trip unchanged. The OCS system is stored alongside on the
58 ``*_ocs`` attributes/keys.
60 Parameters
61 ----------
62 table : `astropy.table.Table`, optional
63 Source table in the CCS. Must contain columns:
65 ``"x"``
66 Field x positions with angular units (e.g. ``u.deg``).
67 ``"y"``
68 Field y positions with angular units (e.g. ``u.deg``).
69 ``"Z{j}"``
70 One column per Noll index *j*, with length units
71 (e.g. ``u.um``).
72 table_ocs : `astropy.table.Table`, optional
73 Source table in the OCS, with the same column layout as
74 ``table``.
76 Attributes
77 ----------
78 field_x, field_y : `numpy.ndarray`
79 CCS x/y field positions in degrees for all sample points,
80 shape ``(n_points_ccs,)``.
81 field_x_ocs, field_y_ocs : `numpy.ndarray`
82 OCS x/y field positions in degrees for all sample points,
83 shape ``(n_points_ocs,)``.
84 noll_indices : `numpy.ndarray`
85 Noll indices of the stored Zernike terms, shape ``(n_zernikes,)``.
86 values, values_ocs : `numpy.ndarray`
87 Zernike coefficients in microns for the CCS and OCS sample
88 points, shape ``(n_points, n_zernikes)``.
89 interpolator, interpolator_ocs : `scipy.interpolate.LinearNDInterpolator`
90 or `None`
91 Interpolators built from the CCS and OCS sample points and
92 values. ``None`` until the corresponding system is populated.
94 Version 1.1 adds the OCS coordinate system alongside the CCS system
95 stored by version 1.
96 """
98 _OBSTYPE = "INTRINSIC_ZERNIKES"
99 _SCHEMA = "Intrinsic Zernikes"
100 _VERSION = 1.1
102 def __init__(self, table=None, table_ocs=None, **kwargs):
103 # CCS uses the un-suffixed names (field_x/field_y/values/
104 # interpolator) for backwards compatibility with version 1, which
105 # only stored the CCS system under those names.
106 self.field_x = np.array([])
107 self.field_y = np.array([])
108 self.values = np.array([])
109 self.field_x_ocs = np.array([])
110 self.field_y_ocs = np.array([])
111 self.values_ocs = np.array([])
112 self.noll_indices = np.array([])
113 self.interpolator = None
114 self.interpolator_ocs = None
116 super().__init__(**kwargs)
118 if table is not None:
119 (self.field_x, self.field_y, self.values, self.noll_indices) = (
120 self._unpackTable(table)
121 )
122 self.interpolator = self._makeInterpolator(
123 self.field_x, self.field_y, self.values
124 )
125 if table_ocs is not None:
126 (self.field_x_ocs, self.field_y_ocs, self.values_ocs, noll_indices_ocs) = (
127 self._unpackTable(table_ocs)
128 )
129 # The CCS and OCS systems are summed term-by-term at query time, so
130 # they must describe the same Noll indices. Store a single shared
131 # ``noll_indices`` and reject tables that disagree.
132 if table is not None and not np.array_equal(noll_indices_ocs, self.noll_indices): 132 ↛ 133line 132 didn't jump to line 133 because the condition on line 132 was never true
133 raise ValueError(
134 "CCS and OCS tables must share the same Noll indices; got "
135 f"{self.noll_indices.tolist()} (CCS) and "
136 f"{noll_indices_ocs.tolist()} (OCS)."
137 )
138 self.interpolator_ocs = self._makeInterpolator(
139 self.field_x_ocs, self.field_y_ocs, self.values_ocs
140 )
142 self.requiredAttributes.update(
143 [
144 "field_x",
145 "field_y",
146 "values",
147 "field_x_ocs",
148 "field_y_ocs",
149 "values_ocs",
150 "noll_indices",
151 ]
152 )
154 @staticmethod
155 def _unpackTable(table):
156 """Unpack a source table into field positions, values, and Noll
157 indices.
159 Parameters
160 ----------
161 table : `astropy.table.Table`
162 Source table with ``"x"``, ``"y"``, and ``"Z{j}"`` columns.
164 Returns
165 -------
166 field_x, field_y : `numpy.ndarray`
167 Field positions in degrees.
168 values : `numpy.ndarray`
169 Zernike coefficients in microns, shape
170 ``(n_points, n_zernikes)`` ordered by ascending Noll index.
171 noll_indices : `numpy.ndarray`
172 Sorted Noll indices.
173 """
174 field_x = table["x"].to("deg").value
175 field_y = table["y"].to("deg").value
176 zcols = [col for col in table.colnames if col.startswith("Z")]
177 noll_indices = np.array(sorted(int(col[1:]) for col in zcols))
178 values = np.column_stack([table[f"Z{j}"].to("um").value for j in noll_indices])
179 return field_x, field_y, values, noll_indices
181 @staticmethod
182 def _makeInterpolator(field_x, field_y, values):
183 """Build a field-position interpolator, or `None` if there are no
184 sample points."""
185 if np.asarray(field_x).size == 0: 185 ↛ 186line 185 didn't jump to line 186 because the condition on line 185 was never true
186 return None
187 return LinearNDInterpolator(np.column_stack((field_x, field_y)), values)
189 @classmethod
190 def fromDict(cls, dictionary):
191 """Construct an IntrinsicZernikes from dictionary of properties.
193 Parameters
194 ----------
195 dictionary : `dict`
196 Dictionary of properties.
198 Returns
199 -------
200 calib : `lsst.ip.isr.IntrinsicZernikes`
201 Constructed calibration.
203 Raises
204 ------
205 RuntimeError
206 Raised if the supplied dictionary is for a different
207 calibration type.
208 """
209 calib = cls()
211 if calib._OBSTYPE != dictionary["metadata"]["OBSTYPE"]:
212 raise RuntimeError(
213 f"Incorrect intrinsic zernikes supplied. "
214 f"Expected {calib._OBSTYPE}, found {dictionary['metadata']['OBSTYPE']}"
215 )
217 calib.setMetadata(dictionary["metadata"])
218 # CCS keys (field_x/field_y/values) are always present, including
219 # in version 1 dictionaries. The OCS keys are optional, so version
220 # 1 dictionaries (CCS only) load with an empty OCS system.
221 calib.field_x = np.array(dictionary["field_x"])
222 calib.field_y = np.array(dictionary["field_y"])
223 calib.values = np.array(dictionary["values"])
224 calib.field_x_ocs = np.array(dictionary.get("field_x_ocs", []))
225 calib.field_y_ocs = np.array(dictionary.get("field_y_ocs", []))
226 calib.values_ocs = np.array(dictionary.get("values_ocs", []))
227 calib.noll_indices = np.array(dictionary["noll_indices"])
228 calib.interpolator = cls._makeInterpolator(
229 calib.field_x, calib.field_y, calib.values
230 )
231 calib.interpolator_ocs = cls._makeInterpolator(
232 calib.field_x_ocs, calib.field_y_ocs, calib.values_ocs
233 )
235 calib.updateMetadata()
236 return calib
238 def toDict(self):
239 """Return a dictionary containing the calibration properties.
241 The dictionary should be able to be round-tripped through
242 `fromDict`.
244 Returns
245 -------
246 dictionary : `dict`
247 Dictionary of properties.
248 """
249 self.updateMetadata()
251 outDict = {}
252 outDict["metadata"] = self.getMetadata()
253 outDict["field_x"] = self.field_x.tolist()
254 outDict["field_y"] = self.field_y.tolist()
255 outDict["values"] = self.values.tolist()
256 outDict["field_x_ocs"] = self.field_x_ocs.tolist()
257 outDict["field_y_ocs"] = self.field_y_ocs.tolist()
258 outDict["values_ocs"] = self.values_ocs.tolist()
259 outDict["noll_indices"] = self.noll_indices.tolist()
261 return outDict
263 @classmethod
264 def fromTable(cls, tableList):
265 """Construct calibration from a list of tables.
267 Parameters
268 ----------
269 tableList : `list` [`astropy.table.Table`]
270 List of tables to use to construct the intrinsic zernikes
271 calibration. Each table is dispatched to the CCS or OCS
272 coordinate system according to its ``coord_sys`` metadata
273 entry (defaulting to ``"CCS"``). A version 1 single-table
274 calibration, which has no ``coord_sys`` entry, is therefore
275 read as the CCS system.
277 Returns
278 -------
279 calib : `lsst.ip.isr.IntrinsicZernikes`
280 The calibration defined in the tables.
281 """
282 tables = {}
283 for table in tableList:
284 coord_sys = table.meta.get("coord_sys", "CCS")
285 if coord_sys not in ("CCS", "OCS"): 285 ↛ 286line 285 didn't jump to line 286 because the condition on line 285 was never true
286 raise RuntimeError(
287 f"Invalid coordinate system {coord_sys} in table metadata; "
288 f"expected 'CCS' or 'OCS'"
289 )
290 tables[coord_sys] = table
292 calib = cls(table=tables.get("CCS"), table_ocs=tables.get("OCS", None))
293 # ``coord_sys`` is a per-table annotation used only to dispatch each
294 # table above; drop it so it does not leak into the calibration
295 # metadata (which must match across a toTable/fromTable round-trip).
296 meta = dict(tableList[0].meta)
297 meta.pop("coord_sys", None)
298 calib.setMetadata(meta)
299 calib.updateMetadata()
300 return calib
302 def toTable(self):
303 """Construct a list of tables containing the information in this
304 calibration.
306 One table is produced per populated coordinate system. The CCS
307 table is always emitted; the OCS table is only emitted when the
308 OCS system holds sample points, so a CCS-only calibration (e.g.
309 one read from a version 1 file) round-trips to a single table,
310 exactly as in version 1. The list of tables should be able to be
311 round-tripped through `fromTable`.
313 Returns
314 -------
315 tableList : `list` [`astropy.table.Table`]
316 List of tables containing the intrinsic zernikes calibration
317 information, one per populated coordinate system.
318 """
319 self.updateMetadata()
321 inMeta = self.getMetadata().toDict()
322 baseMeta = {k: v for k, v in inMeta.items() if v is not None}
323 baseMeta.update({k: "" for k, v in inMeta.items() if v is None})
325 systems = [("CCS", self.field_x, self.field_y, self.values)]
326 if np.asarray(self.field_x_ocs).size > 0: 326 ↛ 329line 326 didn't jump to line 329 because the condition on line 326 was always true
327 systems.append(("OCS", self.field_x_ocs, self.field_y_ocs, self.values_ocs))
329 tableList = []
330 for coord_sys, field_x, field_y, values in systems:
331 data = {
332 "x": field_x * u.deg,
333 "y": field_y * u.deg,
334 }
335 for i, j in enumerate(self.noll_indices):
336 column = values[:, i] if values.ndim == 2 else np.array([])
337 data[f"Z{j}"] = column * u.um
339 table = Table(data)
340 meta = dict(baseMeta)
341 meta["coord_sys"] = coord_sys
342 table.meta = meta
343 tableList.append(table)
345 return tableList
347 def writeText(self, filename, format="auto"):
348 raise NotImplementedError("Text output not implemented for IntrinsicZernikes")
350 def readText(self, filename, format="auto"):
351 raise NotImplementedError("Text input not implemented for IntrinsicZernikes")
353 def getIntrinsicZernikes(
354 self, field_x, field_y, rotTelPos=0.0, noll_indices=None
355 ):
356 """
357 Get the intrinsic Zernike coefficients at a given field position.
359 The returned coefficients are the sum of the CCS contribution
360 (heights_ccs), interpolated at the requested field position,
361 and the OCS contribution (measured_intrinsics),
362 interpolated at the field position rotated by
363 ``rotTelPos``. For calibrations that only store one
364 coordinate system (e.g. version 1 files, which only carry CCS),
365 the missing OCS contribution is simply omitted from the sum.
367 Parameters
368 ----------
369 field_x : `array-like`
370 x-field positions in degrees (CCS).
371 field_y : `array-like`
372 y-field positions in degrees (CCS).
373 rotTelPos : `float`, optional
374 Rotation angle in degrees applied to the query point before
375 interpolating the OCS contribution. Defaults to 0.
376 noll_indices : `list` [`int`], optional
377 List of Noll indices to return. If None, return all.
379 Returns
380 -------
381 zernikes : `array-like`
382 Array of Zernike coefficient values in microns corresponding to the
383 requested Noll indices and field positions.
384 """
385 if noll_indices is None:
386 # Default to all stored terms. The CCS and OCS systems share these
387 # indices, so the output shape is the same whether or not an OCS
388 # system is present (an absent OCS simply adds nothing).
389 noll_indices = self.noll_indices
390 noll_indices = np.array(noll_indices)
391 noll_mask = np.isin(self.noll_indices, noll_indices)
393 field_x = np.asarray(field_x)
394 field_y = np.asarray(field_y)
396 point = np.array([field_x, field_y]).T
397 total = self.interpolator(point)
399 if self.interpolator_ocs is not None: 399 ↛ 411line 399 didn't jump to line 411 because the condition on line 399 was always true
400 # Rotate the query point into the OCS frame before interpolating.
401 theta = np.deg2rad(rotTelPos)
402 cos_a, sin_a = np.cos(theta), np.sin(theta)
403 x_ocs = cos_a * field_x - sin_a * field_y
404 y_ocs = sin_a * field_x + cos_a * field_y
405 point_ocs = np.array([x_ocs, y_ocs]).T
406 # CCS and OCS share the same Noll indices (enforced in __init__),
407 # so the two interpolator outputs line up column-for-column and add
408 # element-wise, preserving the CCS output shape.
409 total = total + self.interpolator_ocs(point_ocs)
411 return total[..., noll_mask]