LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
healSparseMappingProperties.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008-2021 AURA/LSST.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 import numpy as np
23 import healsparse as hsp
24 
25 import lsst.pex.config as pexConfig
26 import lsst.geom
27 from lsst.afw.math import ChebyshevBoundedField, ChebyshevBoundedFieldControl
28 
29 
30 __all__ = ["BasePropertyMapConfig", "PropertyMapRegistry", "register_property_map",
31  "PropertyMapMap", "BasePropertyMap", "ExposureTimePropertyMap",
32  "PsfSizePropertyMap", "PsfE1PropertyMap", "PsfE2PropertyMap",
33  "NExposurePropertyMap", "PsfMaglimPropertyMapConfig",
34  "PsfMaglimPropertyMap", "SkyBackgroundPropertyMap", "SkyNoisePropertyMap",
35  "DcrDraPropertyMap", "DcrDdecPropertyMap", "DcrE1PropertyMap",
36  "DcrE2PropertyMap", "compute_approx_psf_size_and_shape"]
37 
38 
39 class BasePropertyMapConfig(pexConfig.Config):
40  do_min = pexConfig.Field(dtype=bool, default=False,
41  doc="Compute map of property minima.")
42  do_max = pexConfig.Field(dtype=bool, default=False,
43  doc="Compute map of property maxima.")
44  do_mean = pexConfig.Field(dtype=bool, default=False,
45  doc="Compute map of property means.")
46  do_weighted_mean = pexConfig.Field(dtype=bool, default=False,
47  doc="Compute map of weighted property means.")
48  do_sum = pexConfig.Field(dtype=bool, default=False,
49  doc="Compute map of property sums.")
50 
51 
52 class PropertyMapRegistry(pexConfig.Registry):
53  """Class for property map registry.
54 
55  Notes
56  -----
57  This code is based on `lsst.meas.base.PluginRegistry`.
58  """
59  class Configurable:
60  """Class used as the element in the property map registry.
61 
62  Parameters
63  ----------
64  name : `str`
65  Name under which the property map is registered.
66  PropertyMapClass : subclass of `BasePropertyMap`
67  """
68  def __init__(self, name, PropertyMapClass):
69  self.namename = name
70  self.PropertyMapClassPropertyMapClass = PropertyMapClass
71 
72  @property
73  def ConfigClass(self):
74  return self.PropertyMapClassPropertyMapClass.ConfigClass
75 
76  def __call__(self, config):
77  return (self.namename, config, self.PropertyMapClassPropertyMapClass)
78 
79  def register(self, name, PropertyMapClass):
80  """Register a property map class with the given name.
81 
82  Parameters
83  ----------
84  name : `str`
85  The name of the property map.
86  PropertyMapClass : subclass of `BasePropertyMap`
87  """
88  pexConfig.Registry.register(self, name, self.ConfigurableConfigurable(name, PropertyMapClass))
89 
90 
92  """A decorator to register a property map class in its base class's registry."""
93  def decorate(PropertyMapClass):
94  PropertyMapClass.registry.register(name, PropertyMapClass)
95  return PropertyMapClass
96  return decorate
97 
98 
99 def compute_approx_psf_size_and_shape(ccd_row, ra, dec, nx=20, ny=20, orderx=2, ordery=2):
100  """Compute the approximate psf size and shape.
101 
102  This routine fits how the psf size and shape varies over a field by approximating
103  with a Chebyshev bounded field.
104 
105  Parameters
106  ----------
107  ccd_row : `lsst.afw.table.ExposureRecord`
108  Exposure metadata for a given detector exposure.
109  ra : `np.ndarray`
110  Right ascension of points to compute size and shape (degrees).
111  dec : `np.ndarray`
112  Declination of points to compute size and shape (degrees).
113  nx : `int`, optional
114  Number of sampling points in the x direction.
115  ny : `int`, optional
116  Number of sampling points in the y direction.
117  orderx : `int`, optional
118  Chebyshev polynomial order for fit in x direction.
119  ordery : `int`, optional
120  Chebyshev polynomial order for fit in y direction.
121 
122  Returns
123  -------
124  psf_array : `np.ndarray`
125  Record array with "psf_size", "psf_e1", "psf_e2".
126  """
127  pts = [lsst.geom.SpherePoint(r*lsst.geom.degrees, d*lsst.geom.degrees) for
128  r, d in zip(ra, dec)]
129  pixels = ccd_row.getWcs().skyToPixel(pts)
130 
132  ctrl.orderX = orderx
133  ctrl.orderY = ordery
134  ctrl.triangular = False
135 
136  bbox = ccd_row.getBBox()
137  xSteps = np.linspace(bbox.getMinX(), bbox.getMaxX(), nx)
138  ySteps = np.linspace(bbox.getMinY(), bbox.getMaxY(), ny)
139  x = np.tile(xSteps, nx)
140  y = np.repeat(ySteps, ny)
141 
142  psf_size = np.zeros(x.size)
143  psf_e1 = np.zeros(x.size)
144  psf_e2 = np.zeros(x.size)
145  psf_area = np.zeros(x.size)
146 
147  psf = ccd_row.getPsf()
148  for i in range(x.size):
149  shape = psf.computeShape(lsst.geom.Point2D(x[i], y[i]))
150  psf_size[i] = shape.getDeterminantRadius()
151  ixx = shape.getIxx()
152  iyy = shape.getIyy()
153  ixy = shape.getIxy()
154 
155  psf_e1[i] = (ixx - iyy)/(ixx + iyy + 2.*psf_size[i]**2.)
156  psf_e2[i] = (2.*ixy)/(ixx + iyy + 2.*psf_size[i]**2.)
157 
158  im = psf.computeKernelImage(lsst.geom.Point2D(x[i], y[i]))
159  psf_area[i] = np.sum(im.array)/np.sum(im.array**2.)
160 
161  pixel_x = np.array([pix.getX() for pix in pixels])
162  pixel_y = np.array([pix.getY() for pix in pixels])
163 
164  psf_array = np.zeros(pixel_x.size, dtype=[("psf_size", "f8"),
165  ("psf_e1", "f8"),
166  ("psf_e2", "f8"),
167  ("psf_area", "f8")])
168 
169  cheb_size = ChebyshevBoundedField.fit(lsst.geom.Box2I(bbox), x, y, psf_size, ctrl)
170  psf_array["psf_size"] = cheb_size.evaluate(pixel_x, pixel_y)
171  cheb_e1 = ChebyshevBoundedField.fit(lsst.geom.Box2I(bbox), x, y, psf_e1, ctrl)
172  psf_array["psf_e1"] = cheb_e1.evaluate(pixel_x, pixel_y)
173  cheb_e2 = ChebyshevBoundedField.fit(lsst.geom.Box2I(bbox), x, y, psf_e2, ctrl)
174  psf_array["psf_e2"] = cheb_e2.evaluate(pixel_x, pixel_y)
175  cheb_area = ChebyshevBoundedField.fit(lsst.geom.Box2I(bbox), x, y, psf_area, ctrl)
176  psf_array["psf_area"] = cheb_area.evaluate(pixel_x, pixel_y)
177 
178  return psf_array
179 
180 
181 class PropertyMapMap(dict):
182  """Map of property maps to be run for a given task.
183 
184  Notes
185  -----
186  Property maps are classes derived from `BasePropertyMap`
187  """
188  def __iter__(self):
189  for property_map in self.values():
190  if (property_map.config.do_min or property_map.config.do_max or property_map.config.do_mean
191  or property_map.config.do_weighted_mean or property_map.config.do_sum):
192  yield property_map
193 
194 
196  """Base class for property maps.
197 
198  Parameters
199  ----------
200  config : `BasePropertyMapConfig`
201  Property map configuration.
202  name : `str`
203  Property map name.
204  """
205  dtype = np.float64
206  requires_psf = False
207 
208  ConfigClass = BasePropertyMapConfig
209 
210  registry = PropertyMapRegistry(BasePropertyMapConfig)
211 
212  def __init__(self, config, name):
213  object.__init__(self)
214  self.configconfig = config
215  self.namename = name
216  self.zeropointzeropoint = 0.0
217 
218  def initialize_tract_maps(self, nside_coverage, nside):
219  """Initialize the tract maps.
220 
221  Parameters
222  ----------
223  nside_coverage : `int`
224  Healpix nside of the healsparse coverage map.
225  nside : `int`
226  Healpix nside of the property map.
227  """
228  if self.configconfig.do_min:
229  self.min_mapmin_map = hsp.HealSparseMap.make_empty(nside_coverage,
230  nside,
231  self.dtypedtype)
232  if self.configconfig.do_max:
233  self.max_mapmax_map = hsp.HealSparseMap.make_empty(nside_coverage,
234  nside,
235  self.dtypedtype)
236  if self.configconfig.do_mean:
237  self.mean_mapmean_map = hsp.HealSparseMap.make_empty(nside_coverage,
238  nside,
239  self.dtypedtype)
240  if self.configconfig.do_weighted_mean:
241  self.weighted_mean_mapweighted_mean_map = hsp.HealSparseMap.make_empty(nside_coverage,
242  nside,
243  self.dtypedtype)
244  if self.configconfig.do_sum:
245  self.sum_mapsum_map = hsp.HealSparseMap.make_empty(nside_coverage,
246  nside,
247  self.dtypedtype)
248 
249  def initialize_values(self, n_pixels):
250  """Initialize the value arrays for accumulation.
251 
252  Parameters
253  ----------
254  n_pixels : `int`
255  Number of pixels in the map.
256  """
257  if self.configconfig.do_min:
258  self.min_valuesmin_values = np.zeros(n_pixels, dtype=self.dtypedtype)
259  # This works for float types, need check for integers...
260  self.min_valuesmin_values[:] = np.nan
261  if self.configconfig.do_max:
262  self.max_valuesmax_values = np.zeros(n_pixels, dtype=self.dtypedtype)
263  self.max_valuesmax_values[:] = np.nan
264  if self.configconfig.do_mean:
265  self.mean_valuesmean_values = np.zeros(n_pixels, dtype=self.dtypedtype)
266  if self.configconfig.do_weighted_mean:
267  self.weighted_mean_valuesweighted_mean_values = np.zeros(n_pixels, dtype=self.dtypedtype)
268  if self.configconfig.do_sum:
269  self.sum_valuessum_values = np.zeros(n_pixels, dtype=self.dtypedtype)
270 
271  def accumulate_values(self, indices, ra, dec, weights, scalings, row,
272  psf_array=None):
273  """Accumulate values from a row of a visitSummary table.
274 
275  Parameters
276  ----------
277  indices : `np.ndarray`
278  Indices of values that should be accumulated.
279  ra : `np.ndarray`
280  Array of right ascension for indices
281  dec : `np.ndarray`
282  Array of declination for indices
283  weights : `float` or `np.ndarray`
284  Weight(s) for indices to be accumulated.
285  scalings : `float` or `np.ndarray`
286  Scaling values to coadd zeropoint.
287  row : `lsst.afw.table.ExposureRecord`
288  Row of a visitSummary ExposureCatalog.
289  psf_array : `np.ndarray`, optional
290  Array of approximate psf values matched to ra/dec.
291 
292  Raises
293  ------
294  ValueError : Raised if requires_psf is True and psf_array is None.
295  """
296  if self.requires_psfrequires_psf and psf_array is None:
297  name = self.__class__.__name__
298  raise ValueError(f"Cannot compute {name} without psf_array.")
299 
300  values = self._compute_compute(row, ra, dec, scalings, psf_array=psf_array)
301  if self.configconfig.do_min:
302  self.min_valuesmin_values[indices] = np.fmin(self.min_valuesmin_values[indices], values)
303  if self.configconfig.do_max:
304  self.max_valuesmax_values[indices] = np.fmax(self.max_valuesmax_values[indices], values)
305  if self.configconfig.do_mean:
306  self.mean_valuesmean_values[indices] += values
307  if self.configconfig.do_weighted_mean:
308  self.weighted_mean_valuesweighted_mean_values[indices] += weights*values
309  if self.configconfig.do_sum:
310  self.sum_valuessum_values[indices] += values
311 
312  def finalize_mean_values(self, total_weights, total_inputs):
313  """Finalize the accumulation of the mean and weighted mean.
314 
315  Parameters
316  ----------
317  total_weights : `np.ndarray`
318  Total accumulated weights, for each value index.
319  total_inputs : `np.ndarray`
320  Total number of inputs, for each value index.
321  """
322  if self.configconfig.do_mean:
323  use, = np.where(total_inputs > 0)
324  self.mean_valuesmean_values[use] /= total_inputs[use]
325  if self.configconfig.do_weighted_mean:
326  use, = np.where(total_weights > 0.0)
327  self.weighted_mean_valuesweighted_mean_values[use] /= total_weights[use]
328 
329  # And perform any necessary post-processing
330  self._post_process_post_process(total_weights, total_inputs)
331 
332  def set_map_values(self, pixels):
333  """Assign accumulated values to the maps.
334 
335  Parameters
336  ----------
337  pixels : `np.ndarray`
338  Array of healpix pixels (nest scheme) to set in the map.
339  """
340  if self.configconfig.do_min:
341  self.min_mapmin_map[pixels] = self.min_valuesmin_values
342  if self.configconfig.do_max:
343  self.max_mapmax_map[pixels] = self.max_valuesmax_values
344  if self.configconfig.do_mean:
345  self.mean_mapmean_map[pixels] = self.mean_valuesmean_values
346  if self.configconfig.do_weighted_mean:
347  self.weighted_mean_mapweighted_mean_map[pixels] = self.weighted_mean_valuesweighted_mean_values
348  if self.configconfig.do_sum:
349  self.sum_mapsum_map[pixels] = self.sum_valuessum_values
350 
351  def _compute(self, row, ra, dec, scalings, psf_array=None):
352  """Compute map value from a row in the visitSummary catalog.
353 
354  Parameters
355  ----------
356  row : `lsst.afw.table.ExposureRecord`
357  Row of a visitSummary ExposureCatalog.
358  ra : `np.ndarray`
359  Array of right ascensions
360  dec : `np.ndarray`
361  Array of declinations
362  scalings : `float` or `np.ndarray`
363  Scaling values to coadd zeropoint.
364  psf_array : `np.ndarray`, optional
365  Array of approximate psf values matched to ra/dec.
366  """
367  raise NotImplementedError("All property maps must implement _compute()")
368 
369  def _post_process(self, total_weights, total_inputs):
370  """Perform post-processing on values.
371 
372  Parameters
373  ----------
374  total_weights : `np.ndarray`
375  Total accumulated weights, for each value index.
376  total_inputs : `np.ndarray`
377  Total number of inputs, for each value index.
378  """
379  # Override of this method is not required.
380  pass
381 
382 
383 @register_property_map("exposure_time")
385  """Exposure time property map."""
386 
387  def _compute(self, row, ra, dec, scalings, psf_array=None):
388  return row.getVisitInfo().getExposureTime()
389 
390 
391 @register_property_map("psf_size")
393  """PSF size property map."""
394  requires_psf = True
395 
396  def _compute(self, row, ra, dec, scalings, psf_array=None):
397  return psf_array["psf_size"]
398 
399 
400 @register_property_map("psf_e1")
402  """PSF shape e1 property map."""
403  requires_psf = True
404 
405  def _compute(self, row, ra, dec, scalings, psf_array=None):
406  return psf_array["psf_e1"]
407 
408 
409 @register_property_map("psf_e2")
411  """PSF shape e2 property map."""
412  requires_psf = True
413 
414  def _compute(self, row, ra, dec, scalings, psf_array=None):
415  return psf_array["psf_e2"]
416 
417 
418 @register_property_map("n_exposure")
420  """Number of exposures property map."""
421  dtype = np.int32
422 
423  def _compute(self, row, ra, dec, scalings, psf_array=None):
424  return 1
425 
426 
428  """Configuration for the PsfMaglim property map."""
429  maglim_nsigma = pexConfig.Field(dtype=float, default=5.0,
430  doc="Number of sigma to compute magnitude limit.")
431 
432  def validate(self):
433  super().validate()
434  if self.do_mindo_min or self.do_maxdo_max or self.do_meando_mean or self.do_sumdo_sum:
435  raise ValueError("Can only use do_weighted_mean with PsfMaglimPropertyMap")
436 
437 
438 @register_property_map("psf_maglim")
440  """PSF magnitude limit property map."""
441  requires_psf = True
442 
443  ConfigClass = PsfMaglimPropertyMapConfig
444 
445  def _compute(self, row, ra, dec, scalings, psf_array=None):
446  # Our values are the weighted mean of the psf area
447  return psf_array["psf_area"]
448 
449  def _post_process(self, total_weights, total_inputs):
450  psf_area = self.weighted_mean_valuesweighted_mean_values.copy()
451  maglim = (self.zeropointzeropoint
452  - 2.5*np.log10(self.configconfig.maglim_nsigma*np.sqrt(psf_area/total_weights)))
453  self.weighted_mean_valuesweighted_mean_values[:] = maglim
454 
455 
456 @register_property_map("sky_background")
458  """Sky background property map."""
459  def _compute(self, row, ra, dec, scalings, psf_array=None):
460  return scalings*row["skyBg"]
461 
462 
463 @register_property_map("sky_noise")
465  """Sky noise property map."""
466  def _compute(self, row, ra, dec, scalings, psf_array=None):
467  return scalings*row["skyNoise"]
468 
469 
470 @register_property_map("dcr_dra")
472  """Effect of DCR on delta-RA property map."""
473  def _compute(self, row, ra, dec, scalings, psf_array=None):
474  par_angle = row.getVisitInfo().getBoresightParAngle().asRadians()
475  return np.tan(np.deg2rad(row["zenithDistance"]))*np.sin(par_angle)
476 
477 
478 @register_property_map("dcr_ddec")
480  """Effect of DCR on delta-Dec property map."""
481  def _compute(self, row, ra, dec, scalings, psf_array=None):
482  par_angle = row.getVisitInfo().getBoresightParAngle().asRadians()
483  return np.tan(np.deg2rad(row["zenithDistance"]))*np.cos(par_angle)
484 
485 
486 @register_property_map("dcr_e1")
488  """Effect of DCR on psf shape e1 property map."""
489  def _compute(self, row, ra, dec, scalings, psf_array=None):
490  par_angle = row.getVisitInfo().getBoresightParAngle().asRadians()
491  return (np.tan(np.deg2rad(row["zenithDistance"]))**2.)*np.cos(2.*par_angle)
492 
493 
494 @register_property_map("dcr_e2")
496  """Effect of DCR on psf shape e2 property map."""
497  def _compute(self, row, ra, dec, scalings, psf_array=None):
498  par_angle = row.getVisitInfo().getBoresightParAngle().asRadians()
499  return (np.tan(np.deg2rad(row["zenithDistance"]))**2.)*np.sin(2.*par_angle)
A control object used when fitting ChebyshevBoundedField to data (see ChebyshevBoundedField::fit)
An integer coordinate rectangle.
Definition: Box.h:55
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
def _compute(self, row, ra, dec, scalings, psf_array=None)
def accumulate_values(self, indices, ra, dec, weights, scalings, row, psf_array=None)
def compute_approx_psf_size_and_shape(ccd_row, ra, dec, nx=20, ny=20, orderx=2, ordery=2)