LSSTApplications  20.0.0
LSSTDataManagementBasePackage
positionGalSimFakes.py
Go to the documentation of this file.
1 import os
2 import fcntl
3 
4 import numpy as np
5 from astropy.io import fits
6 
7 import lsst.afw.image
8 import lsst.afw.geom
9 import lsst.afw.math
11 import lsst.pex.config as lsstConfig
12 
13 from lsst.pipe.tasks.fakes import BaseFakeSourcesConfig, BaseFakeSourcesTask
14 
15 from . import makeFakeGalaxy as makeFake
16 
17 
19  galList = lsstConfig.Field(dtype=str,
20  doc="catalog of galaxies to add")
21  maxMargin = lsstConfig.Field(dtype=int, default=600,
22  optional=True,
23  doc="Size of margin")
24  seed = lsstConfig.Field(dtype=int, default=1,
25  doc="Seed for random number generator")
26  addShear = lsstConfig.Field(dtype=bool, default=False,
27  doc='include shear in the galaxies')
28  addMask = lsstConfig.Field(dtype=bool, default=False,
29  doc='add FAKE mask plane')
30  sersic_prec = lsstConfig.Field(dtype=float, default=0.0,
31  doc='The desired precision for n')
32  cosStr = 'Use Galsim.COSMOSlog()'
33  serStr = 'Single sersic galaxies added'
34  dSerStr = 'Double sersic galaxies added'
35  realStr = 'Real HST galaxy images added'
36  galType = lsstConfig.ChoiceField(dtype=str, default='sersic',
37  allowed={'dsersic': dSerStr,
38  'sersic': serStr,
39  'real': realStr,
40  'cosmos': cosStr},
41  doc='type of GalSim galaxies to add')
42  exclusionLevel = lsstConfig.ChoiceField(dtype=str, default='none',
43  allowed={'none': "None",
44  'marginal': "Marginal",
45  'bad_stamp': "Bad Stamp",
46  'bad_fits': "Bad Fits"},
47  doc='Exclusion level')
48 
49 
51  ConfigClass = PositionGalSimFakesConfig
52 
53  def __init__(self, **kwargs):
54  BaseFakeSourcesTask.__init__(self, **kwargs)
55  print("RNG seed:", self.config.seed)
56  self.rng = lsst.afw.math.Random(seed=self.config.seed)
57  self.npRand = np.random.RandomState(self.config.seed)
58  self.galData = fits.open(self.config.galList)[1].data
59 
60  def run(self, exposure, background):
61  self.log.info("Adding fake galaxies at real positions")
62  PARENT = lsst.afw.image.PARENT
63  psf = exposure.getPsf()
64  md = exposure.getMetadata()
65  calib = exposure.getCalib()
66  expBBox = exposure.getBBox(PARENT)
67  wcs = exposure.getWcs()
68 
69  """Deal with the skipped ones."""
70  skipLog = 'runAddFake.skipped'
71  if not os.path.isfile(skipLog):
72  os.system('touch ' + skipLog)
73 
74  if self.config.galType == 'cosmos':
75  import galsim
76  exLevel = self.config.exclusionLevel
77  cosmosCat = galsim.COSMOSCatalog(exclusion_level=exLevel)
78  else:
79  cosmosCat = None
80 
81  for igal, gal in enumerate(self.galData):
82  try:
83  galident = gal["ID"]
84  except KeyError:
85  galident = igal + 1
86 
87  try:
88  flux = calib.getFlux(float(gal['mag']))
89  except KeyError:
90  raise KeyError("No mag column in %s" % self.config.galList)
91 
92  try:
93  galCoord = lsst.afw.geom.SpherePoint(gal['RA'], gal['DEC'], lsst.afw.geom.degrees)
94  except KeyError:
95  raise KeyError("No RA/DEC column in {} table".format(self.config.galList))
96 
97  skyToPixelMatrix = wcs.linearizeSkyToPixel(galCoord, lsst.afw.geom.degrees)
98  skyToPixelMatrix = skyToPixelMatrix.getLinear().getParameterVector() / 3600.0
99 
100  galXY = wcs.skyToPixel(galCoord)
101  bboxI = exposure.getBBox(PARENT)
102  bboxI.grow(self.config.maxMargin)
103  if not bboxI.contains(lsst.afw.geom.Point2I(galXY)):
104  # Will just skip this object
105  continue
106 
107  # Check the magnitude
108  if gal['mag'] <= 0:
109  self.log.info("Mag <= 0: Skipping %d" % galident)
110  self.log.info(" mag: %7.3d" % gal['mag'])
111  with open(skipLog, "a") as slog:
112  try:
113  fcntl.flock(slog, fcntl.LOCK_EX)
114  slog.write("%8d , negMag\n" % galident)
115  fcntl.flock(slog, fcntl.LOCK_UN)
116  except IOError:
117  continue
118  continue
119 
120  # This is extrapolating for the PSF, probably not a good idea
121  # Return an Image of the PSF, in a form suitable for convolution.
122  # The returned image is normalized to sum to unity.
123  try:
124  psfImage = psf.computeKernelImage(galXY)
125  except Exception:
126  # There may not be any data at this point, and the object
127  # should be skipped
128  continue
129  try:
130  addShear = self.config.addShear
131  prec = self.config.sersic_prec
132  galType = self.config.galType
133  if self.config.galType != 'cosmos':
134  galArray = makeFake.makeGalaxy(flux, gal,
135  psfImage.getArray(),
136  galType=galType,
137  cosmosCat=None,
138  calib=None,
139  addShear=addShear,
140  transform=skyToPixelMatrix)
141  else:
142  galArray = makeFake.makeGalaxy(flux, gal,
143  psfImage.getArray(),
144  cosmosCat=cosmosCat,
145  calib=calib,
146  galType=galType,
147  addShear=addShear,
148  sersic_prec=prec,
149  transform=skyToPixelMatrix)
150  except IndexError as ierr:
151  self.log.info("GalSim Index Error: Skipping %d" % galident)
152  self.log.info(ierr.message)
153  with open(skipLog, "a") as slog:
154  try:
155  fcntl.flock(slog, fcntl.LOCK_EX)
156  slog.write("%8d , galsimI\n" % galident)
157  fcntl.flock(slog, fcntl.LOCK_UN)
158  except IOError:
159  continue
160  continue
161  except KeyError as kerr:
162  self.log.info("GalSim Key Error: Skipping %d" % galident)
163  self.log.info(kerr.message)
164  with open(skipLog, "a") as slog:
165  try:
166  fcntl.flock(slog, fcntl.LOCK_EX)
167  slog.write("%8d , galsimK\n" % galident)
168  fcntl.flock(slog, fcntl.LOCK_UN)
169  except IOError:
170  continue
171  continue
172  except ValueError as verr:
173  self.log.info("GalSim Value Error: Skipping %d" % galident)
174  self.log.info(verr.message)
175  with open(skipLog, "a") as slog:
176  try:
177  fcntl.flock(slog, fcntl.LOCK_EX)
178  slog.write("%8d , galsimV\n" % galident)
179  fcntl.flock(slog, fcntl.LOCK_UN)
180  except IOError:
181  continue
182  continue
183  except RuntimeError as rerr:
184  self.log.info("GalSim Runtime Error: Skipping %d" % galident)
185  self.log.info(rerr.message)
186  with open(skipLog, "a") as slog:
187  try:
188  fcntl.flock(slog, fcntl.LOCK_EX)
189  slog.write("%8d , galsimR\n" % galident)
190  fcntl.flock(slog, fcntl.LOCK_UN)
191  except IOError:
192  continue
193  continue
194  except Exception as uerr:
195  self.log.info("Unexpected Error: Skipping %d" % galident)
196  self.log.info(uerr.message)
197  with open(skipLog, "a") as slog:
198  try:
199  fcntl.flock(slog, fcntl.LOCK_EX)
200  slog.write("%8d , Unexpected\n" % galident)
201  fcntl.flock(slog, fcntl.LOCK_UN)
202  except IOError:
203  continue
204  continue
205 
206  galImage = lsst.afw.image.ImageF(galArray.astype(np.float32))
207  galBBox = galImage.getBBox(PARENT)
208  galX0 = (galXY.getX() - galBBox.getWidth()/2.0 + 0.5)
209  galY0 = (galXY.getY() - galBBox.getHeight()/2.0 + 0.5)
210  galImage = lsst.afw.math.offsetImage(galImage,
211  galX0, galY0,
212  'lanczos3')
213  galBBox = galImage.getBBox(PARENT)
214 
215  # Check that we're within the larger exposure, otherwise crop
216  parentBox = galImage.getBBox(PARENT)
217  if expBBox.contains(parentBox) is False:
218  newBBox = galImage.getBBox(PARENT)
219  newBBox.clip(expBBox)
220  if newBBox.getArea() <= 0:
221  self.log.info("BBoxEdge Error: Skipping %d" % galident)
222  with open(skipLog, "a") as slog:
223  try:
224  fcntl.flock(slog, fcntl.LOCK_EX)
225  slog.write("%8d , bboxEdge\n" % galident)
226  fcntl.flock(slog, fcntl.LOCK_UN)
227  except IOError:
228  continue
229  continue
230  self.log.info("Cropping FAKE%d from %s to %s" % (galident,
231  str(galBBox), str(newBBox)))
232  galImage = galImage.Factory(galImage, newBBox,
233  PARENT)
234  galBBox = newBBox
235 
236  galMaskedImage = lsst.afw.image.MaskedImageF(galImage)
237 
238  # Put information of the added fake galaxies into the header
239  md.set("FAKE%s" % str(galident), "%.3f, %.3f" % (galXY.getX(),
240  galXY.getY()))
241  self.log.info("Adding fake %s at: %.1f,%.1f" % (str(galident),
242  galXY.getX(),
243  galXY.getY()))
244 
245  galMaskedImage.getMask().set(self.bitmask)
246  if not self.config.addMask:
247  try:
248  galMaskedImage.getMask().removeAndClearMaskPlane('FAKE',
249  True)
250  except Exception:
251  pass
252  try:
253  galMaskedImage.getMask().removeAndClearMaskPlane('CROSSTALK',
254  True)
255  except Exception:
256  pass
257  try:
258  galMaskedImage.getMask().removeAndClearMaskPlane('UNMASKEDNAN',
259  True)
260  except Exception:
261  pass
262 
263  maskedImage = exposure.getMaskedImage()
264  try:
265  maskedImage.getMask().removeAndClearMaskPlane('CROSSTALK',
266  True)
267  except Exception:
268  pass
269  try:
270  maskedImage.getMask().removeAndClearMaskPlane('UNMASKEDNAN',
271  True)
272  except Exception:
273  pass
274  if not self.config.addMask:
275  try:
276  maskedImage.getMask().removeAndClearMaskPlane('FAKE', True)
277  except Exception:
278  pass
279 
280  BBox = galMaskedImage.getBBox(PARENT)
281  subMaskedImage = maskedImage.Factory(exposure.getMaskedImage(),
282  BBox,
283  PARENT)
284  subMaskedImage += galMaskedImage
285 
286  """
287  #
288  galMaskedImage.getMask().set(self.bitmask)
289 
290  maskedImage = exposure.getMaskedImage()
291  try:
292  maskedImage.getMask().removeAndClearMaskPlane('CROSSTALK',
293  True)
294  except Exception:
295  pass
296  try:
297  maskedImage.getMask().removeAndClearMaskPlane('UNMASKEDNAN',
298  True)
299  except Exception:
300  pass
301  try:
302  maskedImage.getMask().removeAndClearMaskPlane('FAKE', True)
303  except Exception:
304  pass
305  """
lsst.synpipe.positionGalSimFakes.PositionGalSimFakesTask.__init__
def __init__(self, **kwargs)
Definition: positionGalSimFakes.py:53
lsst::afw::image
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Definition: imageAlgorithm.dox:1
lsst::log.log.logContinued.info
def info(fmt, *args)
Definition: logContinued.py:198
pex.config.history.format
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
lsst.synpipe.positionGalSimFakes.PositionGalSimFakesTask.run
def run(self, exposure, background)
Definition: positionGalSimFakes.py:60
lsst.pipe.tasks.fakes.BaseFakeSourcesTask.bitmask
bitmask
Definition: fakes.py:61
lsst.pipe.tasks.fakes.BaseFakeSourcesConfig
Definition: fakes.py:29
lsst.synpipe.positionGalSimFakes.PositionGalSimFakesTask
Definition: positionGalSimFakes.py:50
lsst.pipe.tasks.fakes.BaseFakeSourcesTask
Definition: fakes.py:36
lsst.synpipe.positionGalSimFakes.PositionGalSimFakesConfig
Definition: positionGalSimFakes.py:18
lsst.synpipe.positionGalSimFakes.PositionGalSimFakesTask.npRand
npRand
Definition: positionGalSimFakes.py:57
lsst.pipe.base.task.Task.config
config
Definition: task.py:149
lsst.pipe.base.task.Task.log
log
Definition: task.py:148
lsst.pipe.tasks.fakes
Definition: fakes.py:1
lsst::afw::cameraGeom
Definition: Amplifier.h:33
lsst::afw::math::Random
A class that can be used to generate sequences of random numbers according to a number of different a...
Definition: Random.h:57
lsst::afw::math
Definition: statistics.dox:6
lsst::afw::math::offsetImage
std::shared_ptr< ImageT > offsetImage(ImageT const &image, float dx, float dy, std::string const &algorithmName="lanczos5", unsigned int buffer=0)
Return an image offset by (dx, dy) using the specified algorithm.
Definition: offsetImage.cc:41
set
daf::base::PropertySet * set
Definition: fits.cc:912
lsst.synpipe.positionGalSimFakes.PositionGalSimFakesTask.rng
rng
Definition: positionGalSimFakes.py:56
lsst.synpipe.positionGalSimFakes.PositionGalSimFakesTask.galData
galData
Definition: positionGalSimFakes.py:58
lsst::afw::geom
Definition: frameSetUtils.h:40