LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
mockObservation.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010, 2011, 2012 LSST Corporation.
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 
23 import numpy
24 
25 import lsst.pex.config
26 import lsst.afw.table
27 import lsst.afw.geom
28 from lsst.afw.cameraGeom import PIXELS, FOCAL_PLANE
29 import lsst.afw.image
30 import lsst.afw.math
31 import lsst.afw.detection
32 import lsst.pipe.base
33 
34 class MockObservationConfig(lsst.pex.config.Config):
35  pixelScale = lsst.pex.config.Field(
36  dtype=float, default=0.2, optional=False,
37  doc="Pixel scale for mock WCSs in arcseconds/pixel"
38  )
39  doRotate = lsst.pex.config.Field(
40  dtype=bool, default=True, optional=False,
41  doc="Whether to randomly rotate observations relative to the tract Wcs"
42  )
43  fluxMag0 = lsst.pex.config.Field(
44  dtype=float, default=1E11, optional=False,
45  doc="Flux at zero magnitude used to define Calibs."
46  )
47  fluxMag0Sigma = lsst.pex.config.Field(
48  dtype=float, default=100.0, optional=False,
49  doc="Error on flux at zero magnitude used to define Calibs; used to add scatter as well."
50  )
51  expTime = lsst.pex.config.Field(
52  dtype=float, default=60.0, optional=False,
53  doc="Exposure time set in generated Calibs (does not affect flux or noise level)"
54  )
55  psfImageSize = lsst.pex.config.Field(
56  dtype=int, default=21, optional=False,
57  doc="Image width and height of generated Psfs."
58  )
59  psfMinSigma = lsst.pex.config.Field(
60  dtype=float, default=1.5, optional=False,
61  doc="Minimum radius for generated Psfs."
62  )
63  psfMaxSigma = lsst.pex.config.Field(
64  dtype=float, default=3.0, optional=False,
65  doc="Maximum radius for generated Psfs."
66  )
67  seed = lsst.pex.config.Field(dtype=int, default=1, doc="Seed for numpy random number generator")
68 
69 class MockObservationTask(lsst.pipe.base.Task):
70  """Task to generate mock Exposure parameters (Wcs, Psf, Calib), intended for use as a subtask
71  of MockCoaddTask.
72 
73  @todo:
74  - document "pa" in detail; angle of what to what?
75  - document the catalog parameter of the run method
76  """
77 
78  ConfigClass = MockObservationConfig
79 
80  def __init__(self, **kwds):
81  lsst.pipe.base.Task.__init__(self, **kwds)
83  self.ccdKey = self.schema.addField("ccd", type=int, doc="CCD number")
84  self.visitKey = self.schema.addField("visit", type=int, doc="visit number")
85  self.pointingKey = lsst.afw.table.CoordKey.addFields(self.schema, "pointing", "center of visit")
86  self.rng = numpy.random.RandomState(self.config.seed)
87 
88  def run(self, butler, n, tractInfo, camera, catalog=None):
89  """Driver that generates an ExposureCatalog of mock observations.
90 
91  @param[in] butler: a data butler
92  @param[in] n: number of pointings
93  @param[in] camera: camera geometry (an lsst.afw.cameraGeom.Camera)
94  @param[in] catalog: catalog to which to add observations (an ExposureCatalog);
95  if None then a new catalog is created.
96  """
97  if catalog is None:
99  else:
100  if not catalog.getSchema().contains(self.schema):
101  raise ValueError("Catalog schema does not match Task schema")
102  visit = 1
103  for position, pa in self.makePointings(n, tractInfo):
104  for detector in camera:
105  calib = self.buildCalib()
106  record = catalog.addNew()
107  record.setI(self.ccdKey, detector.getId())
108  record.setI(self.visitKey, visit)
109  record.set(self.pointingKey, position)
110  record.setWcs(self.buildWcs(position, pa, detector))
111  record.setCalib(calib)
112  record.setPsf(self.buildPsf(detector))
113  record.setBBox(detector.getBBox())
114  record.setId(butler.get("ccdExposureId", visit=visit, ccd=detector.getId(), immediate=True))
115  visit += 1
116  return catalog
117 
118  def makePointings(self, n, tractInfo):
119  """Generate (celestial) positions and rotation angles that define field locations.
120 
121  Default implementation draws random pointings that are uniform in the tract's image
122  coordinate system.
123 
124  @param[in] n: number of pointings
125  @param[in] tractInfo: skymap tract (a lsst.skymap.TractInfo)
126  @return a Python iterable over (coord, angle) pairs:
127  - coord is an object position (an lsst.afw.coord.Coord)
128  - angle is a position angle (???) (an lsst.afw.geom.Angle)
129 
130  The default implementation returns an iterator (i.e. the function is a "generator"),
131  but derived-class overrides may return any iterable.
132  """
133  wcs = tractInfo.getWcs()
134  bbox = lsst.afw.geom.Box2D(tractInfo.getBBox())
135  bbox.grow(lsst.afw.geom.Extent2D(-0.1 * bbox.getWidth(), -0.1 * bbox.getHeight()))
136  for i in xrange(n):
137  x = self.rng.rand() * bbox.getWidth() + bbox.getMinX()
138  y = self.rng.rand() * bbox.getHeight() + bbox.getMinY()
139  pa = 0.0 * lsst.afw.geom.radians
140  if self.config.doRotate:
141  pa = self.rng.rand() * 2.0 * numpy.pi * lsst.afw.geom.radians
142  yield wcs.pixelToSky(x, y), pa
143 
144  def buildWcs(self, position, pa, detector):
145  """Build a simple TAN Wcs with no distortion and exactly-aligned CCDs.
146 
147  @param[in] position: object position on sky (an lsst.afw.coord.Coord)
148  @param[in] pa: position angle (an lsst.afw.geom.Angle)
149  @param[in] detector: detector information (an lsst.afw.cameraGeom.Detector)
150  """
151  crval = position
152  pixelScale = (self.config.pixelScale * lsst.afw.geom.arcseconds).asDegrees()
155  fpCtr = detector.makeCameraPoint(lsst.afw.geom.Point2D(0, 0), FOCAL_PLANE)
156  crpix = detector.transform(fpCtr, PIXELS).getPoint()
157 
158  wcs = lsst.afw.image.makeWcs(crval, crpix, *cd.getMatrix().flatten())
159  return wcs
160 
161  def buildCalib(self):
162  """Build a simple Calib object with exposure time fixed by config, fluxMag0 drawn from
163  a Gaussian defined by config, and mid-time set to DateTime.now().
164  """
165  calib = lsst.afw.image.Calib()
166  calib.setMidTime(lsst.daf.base.DateTime.now())
167  calib.setExptime(self.config.expTime)
168  calib.setFluxMag0(
169  self.rng.randn() * self.config.fluxMag0Sigma + self.config.fluxMag0,
170  self.config.fluxMag0Sigma
171  )
172  return calib
173 
174  def buildPsf(self, detector):
175  """Build a simple Gaussian Psf with linearly-varying ellipticity and size.
176 
177  The Psf pattern increases sigma_x linearly along the x direction, and sigma_y
178  linearly along the y direction.
179 
180  @param[in] detector: detector information (an lsst.afw.cameraGeom.Detector)
181  @return a psf (an instance of lsst.meas.algorithms.KernelPsf)
182  """
183  bbox = detector.getBBox()
184  dx = (self.config.psfMaxSigma - self.config.psfMinSigma) / bbox.getWidth()
185  dy = (self.config.psfMaxSigma - self.config.psfMinSigma) / bbox.getHeight()
186  sigmaXFunc = lsst.afw.math.PolynomialFunction2D(1)
187  sigmaXFunc.setParameter(0, self.config.psfMinSigma - dx * bbox.getMinX() - dy * bbox.getMinY())
188  sigmaXFunc.setParameter(1, dx)
189  sigmaXFunc.setParameter(2, 0.0)
190  sigmaYFunc = lsst.afw.math.PolynomialFunction2D(1)
191  sigmaYFunc.setParameter(0, self.config.psfMinSigma)
192  sigmaYFunc.setParameter(1, 0.0)
193  sigmaYFunc.setParameter(2, dy)
194  angleFunc = lsst.afw.math.PolynomialFunction2D(0)
195  spatialFuncList = lsst.afw.math.Function2DList()
196  spatialFuncList.append(sigmaXFunc)
197  spatialFuncList.append(sigmaYFunc)
198  spatialFuncList.append(angleFunc)
200  self.config.psfImageSize, self.config.psfImageSize,
201  lsst.afw.math.GaussianFunction2D(self.config.psfMinSigma, self.config.psfMinSigma),
202  spatialFuncList
203  )
204  return lsst.meas.algorithms.KernelPsf(kernel)
static Schema makeMinimalSchema()
Return a minimal schema for Exposure tables and records.
Definition: Exposure.h:164
static DateTime now(void)
Definition: DateTime.cc:554
A Psf defined by a Kernel.
Definition: KernelPsf.h:33
static LinearTransform makeRotation(Angle t)
static CoordKey addFields(afw::table::Schema &schema, std::string const &name, std::string const &doc)
static LinearTransform makeScaling(double s)
Wcs::Ptr makeWcs(coord::Coord const &crval, geom::Point2D const &crpix, double CD11, double CD12, double CD21, double CD22)
Create a Wcs object from crval, crpix, CD, using CD elements (useful from python) ...
Definition: makeWcs.cc:141
boost::enable_if_c< ((C+Nf-N)>=1), ArrayRef< T, Nf,(C+Nf-N)> >::type flatten(Array< T, N, C > const &input)
Create a view into an array with trailing contiguous dimensions merged.
Definition: casts.h:150
Custom catalog class for ExposureRecord/Table.
Definition: Exposure.h:55
A floating-point coordinate rectangle geometry.
Definition: Box.h:271
A kernel described by a function.
Definition: Kernel.h:628