LSSTApplications  18.0.0+54,19.0.0,19.0.0+1,19.0.0+16,19.0.0+17,19.0.0+19,19.0.0+24,19.0.0+3,19.0.0-1-g20d9b18+10,19.0.0-1-g425ff20,19.0.0-1-g5549ca4,19.0.0-1-g580fafe+10,19.0.0-1-g6fe20d0+2,19.0.0-1-g8c57eb9+10,19.0.0-1-gbfe0924+1,19.0.0-1-gdc0e4a7+15,19.0.0-1-ge272bc4+10,19.0.0-1-ge3aa853+1,19.0.0-14-gbb28fe44+1,19.0.0-16-g8258e2a+1,19.0.0-2-g0d9f9cd+17,19.0.0-2-g260436e+1,19.0.0-2-g9b11441+4,19.0.0-2-gd955cfd+23,19.0.0-3-g6513920+1,19.0.0-3-gc4f6e04,19.0.0-4-g41ffa1d+3,19.0.0-4-g725f80e+19,19.0.0-4-ga8eba22,19.0.0-5-g0745e3f+1,19.0.0-5-gd943061d,19.0.0-6-gb6b8b0a+1,19.0.0-7-ge358e0fc4,19.0.0-7-gea0a0fe+6,w.2020.03
LSSTDataManagementBasePackage
testUtils.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 #
4 # Copyright 2008-2017 AURA/LSST.
5 #
6 # This product includes software developed by the
7 # LSST Project (http://www.lsst.org/).
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 LSST License Statement and
20 # the GNU General Public License along with this program. If not,
21 # see <https://www.lsstcorp.org/LegalNotices/>.
22 #
23 
24 import numpy as np
25 
26 import lsst.geom
27 import lsst.afw.image as afwImage
28 from . import SingleGaussianPsf
29 
30 
31 def plantSources(bbox, kwid, sky, coordList, addPoissonNoise=True):
32  """Make an exposure with stars (modelled as Gaussians)
33 
34  @param bbox: parent bbox of exposure
35  @param kwid: kernel width (and height; kernel is square)
36  @param sky: amount of sky background (counts)
37  @param coordList: a list of [x, y, counts, sigma], where:
38  * x,y are relative to exposure origin
39  * counts is the integrated counts for the star
40  * sigma is the Gaussian sigma in pixels
41  @param addPoissonNoise: add Poisson noise to the exposure?
42  """
43  # make an image with sources
44  img = afwImage.ImageD(bbox)
45  meanSigma = 0.0
46  for coord in coordList:
47  x, y, counts, sigma = coord
48  meanSigma += sigma
49 
50  # make a single gaussian psf
51  psf = SingleGaussianPsf(kwid, kwid, sigma)
52 
53  # make an image of it and scale to the desired number of counts
54  thisPsfImg = psf.computeImage(lsst.geom.PointD(x, y))
55  thisPsfImg *= counts
56 
57  # bbox a window in our image and add the fake star image
58  psfBox = thisPsfImg.getBBox()
59  psfBox.clip(bbox)
60  if psfBox != thisPsfImg.getBBox():
61  thisPsfImg = thisPsfImg[psfBox, afwImage.PARENT]
62  imgSeg = img[psfBox, afwImage.PARENT]
63  imgSeg += thisPsfImg
64  meanSigma /= len(coordList)
65 
66  img += sky
67 
68  # add Poisson noise
69  if (addPoissonNoise):
70  np.random.seed(seed=1) # make results reproducible
71  imgArr = img.getArray()
72  imgArr[:] = np.random.poisson(imgArr)
73 
74  # bundle into a maskedimage and an exposure
75  mask = afwImage.Mask(bbox)
76  var = img.convertFloat()
77  img -= sky
78  mimg = afwImage.MaskedImageF(img.convertFloat(), mask, var)
79  exposure = afwImage.makeExposure(mimg)
80 
81  # insert an approximate psf
82  psf = SingleGaussianPsf(kwid, kwid, meanSigma)
83  exposure.setPsf(psf)
84 
85  return exposure
86 
87 
88 def makeRandomTransmissionCurve(rng, minWavelength=4000.0, maxWavelength=7000.0, nWavelengths=200,
89  maxRadius=80.0, nRadii=30, perturb=0.05):
90  """Create a random TransmissionCurve with nontrivial spatial and
91  wavelength variation.
92 
93  Parameters
94  ----------
95  rng : numpy.random.RandomState
96  Random number generator.
97  minWavelength : float
98  Average minimum wavelength for generated TransmissionCurves (will be
99  randomly perturbed).
100  maxWavelength : float
101  Average maximum wavelength for generated TransmissionCurves (will be
102  randomly perturbed).
103  nWavelengths : int
104  Number of samples in the wavelength dimension.
105  maxRadius : float
106  Average maximum radius for spatial variation (will be perturbed).
107  nRadii : int
108  Number of samples in the radial dimension.
109  perturb: float
110  Fraction by which wavelength and radius bounds should be randomly
111  perturbed.
112  """
113  dWavelength = maxWavelength - minWavelength
114 
115  def perturbed(x, s=perturb*dWavelength):
116  return x + 2.0*s*(rng.rand() - 0.5)
117 
118  wavelengths = np.linspace(perturbed(minWavelength), perturbed(maxWavelength), nWavelengths)
119  radii = np.linspace(0.0, perturbed(maxRadius, perturb*maxRadius), nRadii)
120  throughput = np.zeros(wavelengths.shape + radii.shape, dtype=float)
121  # throughput will be a rectangle in wavelength, shifting to higher wavelengths and shrinking
122  # in height with radius, going to zero at all bounds.
123  peak0 = perturbed(0.9, 0.05)
124  start0 = perturbed(minWavelength + 0.25*dWavelength)
125  stop0 = perturbed(minWavelength + 0.75*dWavelength)
126  for i, r in enumerate(radii):
127  mask = np.logical_and(wavelengths >= start0 + r, wavelengths <= stop0 + r)
128  throughput[mask, i] = peak0*(1.0 - r/1000.0)
129  return afwImage.TransmissionCurve.makeRadial(throughput, wavelengths, radii)
def plantSources(bbox, kwid, sky, coordList, addPoissonNoise=True)
Definition: testUtils.py:31
def makeRandomTransmissionCurve(rng, minWavelength=4000.0, maxWavelength=7000.0, nWavelengths=200, maxRadius=80.0, nRadii=30, perturb=0.05)
Definition: testUtils.py:89
std::shared_ptr< Exposure< ImagePixelT, MaskPixelT, VariancePixelT > > makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, std::shared_ptr< geom::SkyWcs const > wcs=std::shared_ptr< geom::SkyWcs const >())
A function to return an Exposure of the correct type (cf.
Definition: Exposure.h:442
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...