23 import healsparse
as hsp
27 from lsst.afw.math import ChebyshevBoundedField, ChebyshevBoundedFieldControl
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"]
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.")
53 """Class for property map registry.
57 This code is based on `lsst.meas.base.PluginRegistry`.
60 """Class used as the element in the property map registry.
65 Name under which the property map is registered.
66 PropertyMapClass : subclass of `BasePropertyMap`
80 """Register a property map class with the given name.
85 The name of the property map.
86 PropertyMapClass : subclass of `BasePropertyMap`
88 pexConfig.Registry.register(self, name, self.
ConfigurableConfigurable(name, PropertyMapClass))
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
100 """Compute the approximate psf size and shape.
102 This routine fits how the psf size and shape varies over a field by approximating
103 with a Chebyshev bounded field.
107 ccd_row : `lsst.afw.table.ExposureRecord`
108 Exposure metadata for a given detector exposure.
110 Right ascension of points to compute size and shape (degrees).
112 Declination of points to compute size and shape (degrees).
114 Number of sampling points in the x direction.
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.
124 psf_array : `np.ndarray`
125 Record array with "psf_size", "psf_e1", "psf_e2".
128 r, d
in zip(ra, dec)]
129 pixels = ccd_row.getWcs().skyToPixel(pts)
134 ctrl.triangular =
False
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)
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)
147 psf = ccd_row.getPsf()
148 for i
in range(x.size):
150 psf_size[i] = shape.getDeterminantRadius()
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.)
159 psf_area[i] = np.sum(im.array)/np.sum(im.array**2.)
161 pixel_x = np.array([pix.getX()
for pix
in pixels])
162 pixel_y = np.array([pix.getY()
for pix
in pixels])
164 psf_array = np.zeros(pixel_x.size, dtype=[(
"psf_size",
"f8"),
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)
182 """Map of property maps to be run for a given task.
186 Property maps are classes derived from `BasePropertyMap`
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):
196 """Base class for property maps.
200 config : `BasePropertyMapConfig`
201 Property map configuration.
208 ConfigClass = BasePropertyMapConfig
213 object.__init__(self)
219 """Initialize the tract maps.
223 nside_coverage : `int`
224 Healpix nside of the healsparse coverage map.
226 Healpix nside of the property map.
228 if self.
configconfig.do_min:
229 self.
min_mapmin_map = hsp.HealSparseMap.make_empty(nside_coverage,
232 if self.
configconfig.do_max:
233 self.
max_mapmax_map = hsp.HealSparseMap.make_empty(nside_coverage,
236 if self.
configconfig.do_mean:
237 self.
mean_mapmean_map = hsp.HealSparseMap.make_empty(nside_coverage,
240 if self.
configconfig.do_weighted_mean:
244 if self.
configconfig.do_sum:
245 self.
sum_mapsum_map = hsp.HealSparseMap.make_empty(nside_coverage,
250 """Initialize the value arrays for accumulation.
255 Number of pixels in the map.
257 if self.
configconfig.do_min:
261 if self.
configconfig.do_max:
264 if self.
configconfig.do_mean:
266 if self.
configconfig.do_weighted_mean:
268 if self.
configconfig.do_sum:
273 """Accumulate values from a row of a visitSummary table.
277 indices : `np.ndarray`
278 Indices of values that should be accumulated.
280 Array of right ascension for indices
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.
294 ValueError : Raised if requires_psf is True and psf_array is None.
297 name = self.__class__.__name__
298 raise ValueError(f
"Cannot compute {name} without psf_array.")
300 values = self.
_compute_compute(row, ra, dec, scalings, psf_array=psf_array)
301 if self.
configconfig.do_min:
303 if self.
configconfig.do_max:
305 if self.
configconfig.do_mean:
307 if self.
configconfig.do_weighted_mean:
309 if self.
configconfig.do_sum:
313 """Finalize the accumulation of the mean and weighted mean.
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.
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)
333 """Assign accumulated values to the maps.
337 pixels : `np.ndarray`
338 Array of healpix pixels (nest scheme) to set in the map.
340 if self.
configconfig.do_min:
342 if self.
configconfig.do_max:
344 if self.
configconfig.do_mean:
346 if self.
configconfig.do_weighted_mean:
348 if self.
configconfig.do_sum:
351 def _compute(self, row, ra, dec, scalings, psf_array=None):
352 """Compute map value from a row in the visitSummary catalog.
356 row : `lsst.afw.table.ExposureRecord`
357 Row of a visitSummary ExposureCatalog.
359 Array of right ascensions
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.
367 raise NotImplementedError(
"All property maps must implement _compute()")
369 def _post_process(self, total_weights, total_inputs):
370 """Perform post-processing on values.
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.
385 """Exposure time property map."""
387 def _compute(self, row, ra, dec, scalings, psf_array=None):
388 return row.getVisitInfo().getExposureTime()
391 @register_property_map("psf_size")
393 """PSF size property map."""
396 def _compute(self, row, ra, dec, scalings, psf_array=None):
397 return psf_array[
"psf_size"]
400 @register_property_map("psf_e1")
402 """PSF shape e1 property map."""
405 def _compute(self, row, ra, dec, scalings, psf_array=None):
406 return psf_array[
"psf_e1"]
409 @register_property_map("psf_e2")
411 """PSF shape e2 property map."""
414 def _compute(self, row, ra, dec, scalings, psf_array=None):
415 return psf_array[
"psf_e2"]
418 @register_property_map("n_exposure")
420 """Number of exposures property map."""
423 def _compute(self, row, ra, dec, scalings, psf_array=None):
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.")
435 raise ValueError(
"Can only use do_weighted_mean with PsfMaglimPropertyMap")
438 @register_property_map("psf_maglim")
440 """PSF magnitude limit property map."""
443 ConfigClass = PsfMaglimPropertyMapConfig
445 def _compute(self, row, ra, dec, scalings, psf_array=None):
447 return psf_array[
"psf_area"]
449 def _post_process(self, total_weights, total_inputs):
452 - 2.5*np.log10(self.
configconfig.maglim_nsigma*np.sqrt(psf_area/total_weights)))
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"]
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"]
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)
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)
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)
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.
Point in an unspecified spherical coordinate system.
def _compute(self, row, ra, dec, scalings, psf_array=None)
def accumulate_values(self, indices, ra, dec, weights, scalings, row, psf_array=None)
def set_map_values(self, pixels)
def __init__(self, config, name)
def initialize_values(self, n_pixels)
def initialize_tract_maps(self, nside_coverage, nside)
def finalize_mean_values(self, total_weights, total_inputs)
def _post_process(self, total_weights, total_inputs)
def __call__(self, config)
def __init__(self, name, PropertyMapClass)
def register(self, name, PropertyMapClass)
def register_property_map(name)
def compute_approx_psf_size_and_shape(ccd_row, ra, dec, nx=20, ny=20, orderx=2, ordery=2)