LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
healSparseMappingProperties.py
Go to the documentation of this file.
1# This file is part of pipe_tasks.
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
22__all__ = ["BasePropertyMapConfig", "PropertyMapRegistry", "register_property_map",
23 "PropertyMapMap", "BasePropertyMap", "ExposureTimePropertyMap",
24 "PsfSizePropertyMap", "PsfE1PropertyMap", "PsfE2PropertyMap",
25 "NExposurePropertyMap", "PsfMaglimPropertyMapConfig",
26 "PsfMaglimPropertyMap", "SkyBackgroundPropertyMap", "SkyNoisePropertyMap",
27 "DcrDraPropertyMap", "DcrDdecPropertyMap", "DcrE1PropertyMap",
28 "DcrE2PropertyMap", "EpochPropertyMap", "compute_approx_psf_size_and_shape"]
29
30import numpy as np
31import healsparse as hsp
32import copy
33
34import lsst.pex.config as pexConfig
35import lsst.geom
36from lsst.afw.math import ChebyshevBoundedField, ChebyshevBoundedFieldControl
37
38
39class 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
52class PropertyMapRegistry(pexConfig.Registry):
53 """Class for property map registry.
54
55 Notes
56 -----
57 This code is based on `lsst.meas.base.PluginRegistry`.
58 """
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.name = name
70 self.PropertyMapClass = PropertyMapClass
71
72 @property
73 def ConfigClass(self):
74 return self.PropertyMapClass.ConfigClass
75
76 def __call__(self, config):
77 return (self.name, config, self.PropertyMapClass)
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.Configurable(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
99def 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 # Protect against nans which can come in at the edges and masked regions.
170 good = np.isfinite(psf_size)
171 x = x[good]
172 y = y[good]
173
174 cheb_size = ChebyshevBoundedField.fit(lsst.geom.Box2I(bbox), x, y, psf_size[good], ctrl)
175 psf_array["psf_size"] = cheb_size.evaluate(pixel_x, pixel_y)
176 cheb_e1 = ChebyshevBoundedField.fit(lsst.geom.Box2I(bbox), x, y, psf_e1[good], ctrl)
177 psf_array["psf_e1"] = cheb_e1.evaluate(pixel_x, pixel_y)
178 cheb_e2 = ChebyshevBoundedField.fit(lsst.geom.Box2I(bbox), x, y, psf_e2[good], ctrl)
179 psf_array["psf_e2"] = cheb_e2.evaluate(pixel_x, pixel_y)
180 cheb_area = ChebyshevBoundedField.fit(lsst.geom.Box2I(bbox), x, y, psf_area[good], ctrl)
181 psf_array["psf_area"] = cheb_area.evaluate(pixel_x, pixel_y)
182
183 return psf_array
184
185
186class PropertyMapMap(dict):
187 """Map of property maps to be run for a given task.
188
189 Notes
190 -----
191 Property maps are classes derived from `BasePropertyMap`
192 """
193 def __iter__(self):
194 for property_map in self.values():
195 if (property_map.config.do_min or property_map.config.do_max or property_map.config.do_mean
196 or property_map.config.do_weighted_mean or property_map.config.do_sum):
197 yield property_map
198
199
201 """Base class for property maps.
202
203 Parameters
204 ----------
205 config : `BasePropertyMapConfig`
206 Property map configuration.
207 name : `str`
208 Property map name.
209 """
210 dtype = np.float64
211 requires_psf = False
212
213 description = ""
214 unit = ""
215
216 ConfigClass = BasePropertyMapConfig
217
218 registry = PropertyMapRegistry(BasePropertyMapConfig)
219
220 def __init__(self, config, name):
221 object.__init__(self)
222 self.config = config
223 self.name = name
224 self.zeropoint = 0.0
225
226 def initialize_tract_maps(self, nside_coverage, nside):
227 """Initialize the tract maps.
228
229 Parameters
230 ----------
231 nside_coverage : `int`
232 Healpix nside of the healsparse coverage map.
233 nside : `int`
234 Healpix nside of the property map.
235 """
236 base_metadata = {
237 "DESCRIPTION": self.descriptiondescription,
238 "UNIT": self.unitunit,
239 }
240 if self.config.do_min:
241 metadata = copy.copy(base_metadata)
242 metadata["OPERATION"] = "minimum"
243 self.min_map = hsp.HealSparseMap.make_empty(nside_coverage,
244 nside,
245 self.dtype,
246 metadata=metadata)
247 if self.config.do_max:
248 metadata = copy.copy(base_metadata)
249 metadata["OPERATION"] = "maximum"
250 self.max_map = hsp.HealSparseMap.make_empty(nside_coverage,
251 nside,
252 self.dtype,
253 metadata=metadata)
254 if self.config.do_mean:
255 metadata = copy.copy(base_metadata)
256 metadata["OPERATION"] = "mean"
257 self.mean_map = hsp.HealSparseMap.make_empty(nside_coverage,
258 nside,
259 self.dtype,
260 metadata=metadata)
261 if self.config.do_weighted_mean:
262 metadata = copy.copy(base_metadata)
263 metadata["OPERATION"] = "weighted mean"
264 self.weighted_mean_map = hsp.HealSparseMap.make_empty(nside_coverage,
265 nside,
266 self.dtype,
267 metadata=metadata)
268 if self.config.do_sum:
269 metadata = copy.copy(base_metadata)
270 metadata["OPERATION"] = "sum"
271 self.sum_map = hsp.HealSparseMap.make_empty(nside_coverage,
272 nside,
273 self.dtype,
274 metadata=metadata)
275
276 def initialize_values(self, n_pixels):
277 """Initialize the value arrays for accumulation.
278
279 Parameters
280 ----------
281 n_pixels : `int`
282 Number of pixels in the map.
283 """
284 if self.config.do_min:
285 self.min_values = np.zeros(n_pixels, dtype=self.dtype)
286 # This works for float types, need check for integers...
287 self.min_values[:] = np.nan
288 if self.config.do_max:
289 self.max_values = np.zeros(n_pixels, dtype=self.dtype)
290 self.max_values[:] = np.nan
291 if self.config.do_mean:
292 self.mean_values = np.zeros(n_pixels, dtype=self.dtype)
293 if self.config.do_weighted_mean:
294 self.weighted_mean_values = np.zeros(n_pixels, dtype=self.dtype)
295 if self.config.do_sum:
296 self.sum_values = np.zeros(n_pixels, dtype=self.dtype)
297
298 def accumulate_values(self, indices, ra, dec, weights, scalings, row,
299 psf_array=None):
300 """Accumulate values from a row of a visitSummary table.
301
302 Parameters
303 ----------
304 indices : `np.ndarray`
305 Indices of values that should be accumulated.
306 ra : `np.ndarray`
307 Array of right ascension for indices
308 dec : `np.ndarray`
309 Array of declination for indices
310 weights : `float` or `np.ndarray`
311 Weight(s) for indices to be accumulated.
312 scalings : `float` or `np.ndarray`
313 Scaling values to coadd zeropoint.
314 row : `lsst.afw.table.ExposureRecord`
315 Row of a visitSummary ExposureCatalog.
316 psf_array : `np.ndarray`, optional
317 Array of approximate psf values matched to ra/dec.
318
319 Raises
320 ------
321 ValueError : Raised if requires_psf is True and psf_array is None.
322 """
323 if self.requires_psf and psf_array is None:
324 name = self.__class__.__name__
325 raise ValueError(f"Cannot compute {name} without psf_array.")
326
327 values = self._compute(row, ra, dec, scalings, psf_array=psf_array)
328 if self.config.do_min:
329 self.min_values[indices] = np.fmin(self.min_values[indices], values)
330 if self.config.do_max:
331 self.max_values[indices] = np.fmax(self.max_values[indices], values)
332 if self.config.do_mean:
333 self.mean_values[indices] += values
334 if self.config.do_weighted_mean:
335 self.weighted_mean_values[indices] += weights*values
336 if self.config.do_sum:
337 self.sum_values[indices] += values
338
339 def finalize_mean_values(self, total_weights, total_inputs):
340 """Finalize the accumulation of the mean and weighted mean.
341
342 Parameters
343 ----------
344 total_weights : `np.ndarray`
345 Total accumulated weights, for each value index.
346 total_inputs : `np.ndarray`
347 Total number of inputs, for each value index.
348 """
349 if self.config.do_mean:
350 use, = np.where(total_inputs > 0)
351 self.mean_values[use] /= total_inputs[use]
352 if self.config.do_weighted_mean:
353 use, = np.where(total_weights > 0.0)
354 self.weighted_mean_values[use] /= total_weights[use]
355
356 # And perform any necessary post-processing
357 self._post_process(total_weights, total_inputs)
358
359 def set_map_values(self, pixels):
360 """Assign accumulated values to the maps.
361
362 Parameters
363 ----------
364 pixels : `np.ndarray`
365 Array of healpix pixels (nest scheme) to set in the map.
366 """
367 if self.config.do_min:
368 self.min_map[pixels] = self.min_values
369 if self.config.do_max:
370 self.max_map[pixels] = self.max_values
371 if self.config.do_mean:
372 self.mean_map[pixels] = self.mean_values
373 if self.config.do_weighted_mean:
374 self.weighted_mean_map[pixels] = self.weighted_mean_values
375 if self.config.do_sum:
376 self.sum_map[pixels] = self.sum_values
377
378 def _compute(self, row, ra, dec, scalings, psf_array=None):
379 """Compute map value from a row in the visitSummary catalog.
380
381 Parameters
382 ----------
383 row : `lsst.afw.table.ExposureRecord`
384 Row of a visitSummary ExposureCatalog.
385 ra : `np.ndarray`
386 Array of right ascensions
387 dec : `np.ndarray`
388 Array of declinations
389 scalings : `float` or `np.ndarray`
390 Scaling values to coadd zeropoint.
391 psf_array : `np.ndarray`, optional
392 Array of approximate psf values matched to ra/dec.
393 """
394 raise NotImplementedError("All property maps must implement _compute()")
395
396 def _post_process(self, total_weights, total_inputs):
397 """Perform post-processing on values.
398
399 Parameters
400 ----------
401 total_weights : `np.ndarray`
402 Total accumulated weights, for each value index.
403 total_inputs : `np.ndarray`
404 Total number of inputs, for each value index.
405 """
406 # Override of this method is not required.
407 pass
408
409
410@register_property_map("exposure_time")
412 """Exposure time property map."""
413 description = "Exposure time"
414 unit = "s"
415
416 def _compute(self, row, ra, dec, scalings, psf_array=None):
417 return row.getVisitInfo().getExposureTime()
418
419
420@register_property_map("psf_size")
422 """PSF size property map."""
423 description = "PSF size"
424 unit = "pixel"
425 requires_psf = True
426
427 def _compute(self, row, ra, dec, scalings, psf_array=None):
428 return psf_array["psf_size"]
429
430
431@register_property_map("psf_e1")
433 """PSF shape e1 property map."""
434 description = "PSF e1"
435 requires_psf = True
436
437 def _compute(self, row, ra, dec, scalings, psf_array=None):
438 return psf_array["psf_e1"]
439
440
441@register_property_map("psf_e2")
443 """PSF shape e2 property map."""
444 description = "PSF e2"
445 requires_psf = True
446
447 def _compute(self, row, ra, dec, scalings, psf_array=None):
448 return psf_array["psf_e2"]
449
450
451@register_property_map("n_exposure")
453 """Number of exposures property map."""
454 description = "Number of exposures"
455 dtype = np.int32
456
457 def _compute(self, row, ra, dec, scalings, psf_array=None):
458 return 1
459
460
462 """Configuration for the PsfMaglim property map."""
463 maglim_nsigma = pexConfig.Field(dtype=float, default=5.0,
464 doc="Number of sigma to compute magnitude limit.")
465
466 def validate(self):
467 super().validate()
468 if self.do_min or self.do_max or self.do_mean or self.do_sum:
469 raise ValueError("Can only use do_weighted_mean with PsfMaglimPropertyMap")
470
471
472@register_property_map("psf_maglim")
474 """PSF magnitude limit property map."""
475 description = "PSF magnitude limit"
476 unit = "mag(AB)"
477 requires_psf = True
478
479 ConfigClass = PsfMaglimPropertyMapConfig
480
481 def _compute(self, row, ra, dec, scalings, psf_array=None):
482 # Our values are the weighted mean of the psf area
483 return psf_array["psf_area"]
484
485 def _post_process(self, total_weights, total_inputs):
486 psf_area = self.weighted_mean_values.copy()
487 maglim = (self.zeropoint
488 - 2.5*np.log10(self.config.maglim_nsigma*np.sqrt(psf_area/total_weights)))
489 self.weighted_mean_values[:] = maglim
490
491
492@register_property_map("sky_background")
494 """Sky background property map."""
495 description = "Sky background"
496 unit = "nJy"
497
498 def _compute(self, row, ra, dec, scalings, psf_array=None):
499 return scalings*row["skyBg"]
500
501
502@register_property_map("sky_noise")
504 """Sky noise property map."""
505 description = "Sky noise"
506 unit = "nJy"
507
508 def _compute(self, row, ra, dec, scalings, psf_array=None):
509 return scalings*row["skyNoise"]
510
511
512@register_property_map("dcr_dra")
514 """Effect of DCR on delta-RA property map."""
515 description = "Effect of DCR on delta-RA"
516
517 def _compute(self, row, ra, dec, scalings, psf_array=None):
518 par_angle = row.getVisitInfo().getBoresightParAngle().asRadians()
519 return np.tan(np.deg2rad(row["zenithDistance"]))*np.sin(par_angle)
520
521
522@register_property_map("dcr_ddec")
524 """Effect of DCR on delta-Dec property map."""
525 description = "Effect of DCR on delta-Dec"
526
527 def _compute(self, row, ra, dec, scalings, psf_array=None):
528 par_angle = row.getVisitInfo().getBoresightParAngle().asRadians()
529 return np.tan(np.deg2rad(row["zenithDistance"]))*np.cos(par_angle)
530
531
532@register_property_map("dcr_e1")
534 """Effect of DCR on psf shape e1 property map."""
535 description = "Effect of DCR on PSF shape e1"
536
537 def _compute(self, row, ra, dec, scalings, psf_array=None):
538 par_angle = row.getVisitInfo().getBoresightParAngle().asRadians()
539 return (np.tan(np.deg2rad(row["zenithDistance"]))**2.)*np.cos(2.*par_angle)
540
541
542@register_property_map("dcr_e2")
544 """Effect of DCR on psf shape e2 property map."""
545 description = "Effect of DCR on PSF shape e2"
546
547 def _compute(self, row, ra, dec, scalings, psf_array=None):
548 par_angle = row.getVisitInfo().getBoresightParAngle().asRadians()
549 return (np.tan(np.deg2rad(row["zenithDistance"]))**2.)*np.sin(2.*par_angle)
550
551
552@register_property_map("epoch")
554 """Observation epoch (mjd) property map."""
555 description = "Observation epoch"
556 unit = "day"
557
558 def _compute(self, row, ra, dec, scalings, psf_array=None):
559 date = row.getVisitInfo().getDate()
560 return date.get(date.MJD)
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
accumulate_values(self, indices, ra, dec, weights, scalings, row, psf_array=None)
compute_approx_psf_size_and_shape(ccd_row, ra, dec, nx=20, ny=20, orderx=2, ordery=2)