LSSTApplications  20.0.0
LSSTDataManagementBasePackage
makeRaDecCat.py
Go to the documentation of this file.
1 """
2 Create N random selected (RA, Dec) pairs within a desired region.
3 The inputs can be:
4  1. RA, Dec range: (minRa, maxRa, minDec, maxDec)
5  2. DataId for single frame image: (visit, ccd)
6  3. DataId for coadd image: (tract, patch, filter)
7 """
8 
9 import os
10 import warnings
11 
12 import numpy as np
13 from numpy.random import uniform
14 
15 
16 def polyReadWkb(wkbName, load=True):
17  from shapely.wkt import loads
18 
19  wkbFile = open(wkbName, 'r')
20  polyWkb = wkbFile.read().decode('hex')
21  wkbFile.close()
22 
23  if load is True:
24  return loads(polyWkb)
25  else:
26  return polyWkb
27 
28 
29 def getImageRaDecRange(rootDir, dataId, dataType='calexp'):
30  """
31  Get the Ra,Dec range for certain single frame or coadded image
32  """
33 
34  import lsst.daf.persistence as dafPersist
35  import lsst.afw.image as afwImage
36 
37  butler = dafPersist.Butler(rootDir)
38  exposure = butler.get(dataType, dataId)
39 
40  bboxI = exposure.getBBox(afwImage.PARENT)
41  wcs = exposure.getWcs()
42 
43  minPix = wcs.pixelToSky(bboxI.getMinX(), bboxI.getMinY())
44  maxPix = wcs.pixelToSky(bboxI.getMaxX(), bboxI.getMaxY())
45 
46  ra1 = minPix.getLongitude().asDegrees()
47  ra2 = maxPix.getLongitude().asDegrees()
48  dec1 = minPix.getLatitude().asDegrees()
49  dec2 = maxPix.getLatitude().asDegrees()
50 
51  minRa, maxRa = min(ra1, ra2), max(ra1, ra2)
52  minDec, maxDec = min(dec1, dec2), max(dec1, dec2)
53 
54  return [minRa, maxRa, minDec, maxDec]
55 
56 
57 def getRandomRaDec(nRand, minRa, maxRa, minDec, maxDec, rad=None):
58  """
59  Randomly select Ra,Dec pairs from the input Ra,Dec range
60  """
61  if minRa > maxRa or minDec > maxDec:
62  raise Exception('Please provide appropriate Ra,Dec range !')
63 
64  if rad is None:
65  raArr = uniform(low=minRa, high=maxRa, size=nRand)
66  decArr = uniform(low=minDec, high=maxDec, size=nRand)
67  return list(zip(raArr, decArr))
68  else:
69  import lsst.afw.geom as afwGeom
70 
71  minSep = float(rad)
72  raArr = []
73  decArr = []
74  numTry = 0
75  while len(raArr) < nRand:
76  if numTry == 0:
77  raArr.append(uniform(low=minRa, high=maxRa))
78  decArr.append(uniform(low=minDec, high=maxDec))
79  numTry += 1
80  else:
81  raTry = uniform(low=minRa, high=maxRa)
82  decTry = uniform(low=minDec, high=maxDec)
83  coordTry = afwGeom.SpherePoint(raTry, decTry, afwGeom.degrees)
84  nExist = len(raArr)
85  sepGood = True
86  for ii in range(nExist):
87  coordTest = afwGeom.SpherePoint(raArr[ii], decArr[ii], afwGeom.degrees)
88  sep = coordTry.separation(coordTest).asArcseconds()
89  if sep <= minSep:
90  sepGood = False
91  break
92  if sepGood:
93  raArr.append(raTry)
94  decArr.append(decTry)
95  return list(zip(raArr, decArr))
96 
97 
98 def plotRandomRaDec(randomRaDec, rangeRaDec=None):
99  """
100  Plot the distribution of radom Ra, Dec for examination
101  """
102  import matplotlib.pyplot as plt
103 
104  plt.scatter(*list(zip(*randomRaDec)))
105  plt.xlabel(r'RA (J2000)', fontsize=20, labelpad=20)
106  plt.ylabel(r'DEC (J2000)', fontsize=20, labelpad=20)
107 
108  # TODO : To be finished
109  """
110  if rangeRaDec is not None:
111  if type(rangeRaDec) is dict:
112  raMin, raMax = rangeRaDec['raMin'], rangeRaDec['raMax']
113  decMin, decMax = rangeRaDec['raMin'], rangeRaDec['raMax']
114  else:
115  raMin, raMax = rangeRaDec[0], rangeRaDec[1]
116  decMin, decMax = rangeRaDec[0], rangeRaDec[1]
117  raDec0 = (raMin, decMin)
118  raRange = (raMax - raMin)
119  decRange = (decMax - decMin)
120  """
121 
122  plt.gcf().savefig('randomRaDec.png')
123 
124  return None
125 
126 
127 def makeRaDecCat(nRand, dataId=None, rangeRaDec=None, rad=None,
128  rootDir='/lustre/Subaru/SSP/rerun/yasuda/SSP3.8.5_20150725/',
129  inputCat=None, plot=False, acpMask=None, rejMask=None):
130  """
131  Generate nRand random RA,Dec pairs in a desired region of sky
132  The region can be defined by:
133  1) Single frame image
134  2) Coadd image
135  3) RA, Dec range
136  """
137 
138  if dataId is not None:
139  if 'visit' in list(dataId.keys()) and 'ccd' in list(dataId.keys()):
140  # Input should be a single frame image
141  rangeRaDec = getImageRaDecRange(rootDir, dataId)
142  randomRaDec = getRandomRaDec(nRand, rangeRaDec[0], rangeRaDec[1],
143  rangeRaDec[2], rangeRaDec[3], rad=rad)
144  elif ('tract' in list(dataId.keys()) and
145  'patch' in list(dataId.keys()) and
146  'filter' in dataId.keys):
147  # Input should be a coadd image
148  rangeRaDec = getImageRaDecRange(rootDir, dataId,
149  dataType='deepCoadd_calexp')
150  randomRaDec = getRandomRaDec(nRand, rangeRaDec[0], rangeRaDec[1],
151  rangeRaDec[2], rangeRaDec[3], rad=rad)
152  else:
153  raise KeyError('Please provide the correct dataId !')
154  elif rangeRaDec is not None:
155  if type(rangeRaDec) is dict:
156  # rKeys = rangeRaDec.keys()
157  randomRaDec = getRandomRaDec(nRand, rangeRaDec['minRa'],
158  rangeRaDec['maxRa'],
159  rangeRaDec['minDec'],
160  rangeRaDec['maxDec'], rad=rad)
161  elif (type(rangeRaDec) is list or type(rangeRaDec).__module__ == 'numpy'):
162  if len(rangeRaDec) >= 4:
163  randomRaDec = getRandomRaDec(nRand, rangeRaDec[0],
164  rangeRaDec[1], rangeRaDec[2],
165  rangeRaDec[3], rad=rad)
166  else:
167  raise Exception('randomRaDec should have at least 4 elements!')
168  else:
169  raise Exception('randomRaDec need to be Dict/List/Numpy.Array')
170  else:
171  raise Exception("Need to provide either dataId or rangeRaDec")
172 
173  """
174  Add by Song Huang 15-09-01
175  Filter the random catalog through two masks
176  """
177  if acpMask is not None or rejMask is not None:
178  try:
179  from shapely.geometry import Point
180  from shapely.prepared import prep
181 
182  raArr, decArr = np.array(list(zip(*randomRaDec)))
183  if os.path.isfile(acpMask):
184  acpRegs = polyReadWkb(acpMask)
185  acpPrep = prep(acpRegs)
186  inside = list(map(lambda x, y: acpPrep.contains(Point(x, y)),
187  raArr, decArr))
188  else:
189  inside = np.isfinite(raArr)
190  if os.path.isfile(rejMask):
191  rejRegs = polyReadWkb(rejMask)
192  rejPrep = prep(rejRegs)
193  masked = list(map(lambda x, y: rejPrep.contains(Point(x, y)),
194  raArr, decArr))
195  else:
196  masked = np.isnan(raArr)
197  useful = list(map(lambda x, y: x and (not y), inside, masked))
198  randomUse = list(zip(raArr[useful], decArr[useful]))
199  except ImportError:
200  warnings.warn('Please install the Shapely')
201  randomUse = randomRaDec
202  else:
203  randomUse = randomRaDec
204 
205  if plot:
206  plotRandomRaDec(randomUse, rangeRaDec=rangeRaDec)
207 
208  if inputCat is not None:
209 
210  import os
211  import astropy.table
212 
213  if os.path.exists(inputCat):
214 
215  galCat = astropy.table.Table.read(inputCat, format='fits')
216  outCat = inputCat.strip().replace('.fits', '_radec.fits')
217 
218  nGal = len(galCat)
219 
220  if nGal == nRand:
221  raArr, decArr = np.array(list(zip(*randomUse)))
222  raCol = astropy.table.Column(name='RA', data=raArr)
223  decCol = astropy.table.Column(name='Dec', data=decArr)
224  galCat.add_columns([raCol, decCol])
225  elif nGal < nRand:
226  import random
227  raArr, decArr = np.array(list(zip(*random.sample(randomUse, nGal))))
228  raCol = astropy.table.Column(name='RA', data=raArr)
229  decCol = astropy.table.Column(name='Dec', data=decArr)
230  galCat.add_columns([raCol, decCol])
231  else:
232  import random
233  indGal = np.arange(nGal)
234  galCatRand = galCat[random.sample(indGal, nRand)]
235  raArr, decArr = np.array(list(zip(*randomUse)))
236  raCol = astropy.table.Column(name='RA', data=raArr)
237  decCol = astropy.table.Column(name='Dec', data=decArr)
238  galCatRand.add_columns([raCol, decCol])
239  galCat = galCatRand
240 
241  galCat.write(outCat, format='fits', overwrite=True)
242  else:
243  raise Exception('Can not find input catalog %s!' % inputCat)
244 
245  return randomUse
lsst::afw::image
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Definition: imageAlgorithm.dox:1
lsst.synpipe.makeRaDecCat.polyReadWkb
def polyReadWkb(wkbName, load=True)
Definition: makeRaDecCat.py:16
lsst.synpipe.makeRaDecCat.makeRaDecCat
def makeRaDecCat(nRand, dataId=None, rangeRaDec=None, rad=None, rootDir='/lustre/Subaru/SSP/rerun/yasuda/SSP3.8.5_20150725/', inputCat=None, plot=False, acpMask=None, rejMask=None)
Definition: makeRaDecCat.py:127
lsst.synpipe.makeRaDecCat.getImageRaDecRange
def getImageRaDecRange(rootDir, dataId, dataType='calexp')
Definition: makeRaDecCat.py:29
lsst.synpipe.makeRaDecCat.getRandomRaDec
def getRandomRaDec(nRand, minRa, maxRa, minDec, maxDec, rad=None)
Definition: makeRaDecCat.py:57
lsst.synpipe.makeRaDecCat.plotRandomRaDec
def plotRandomRaDec(randomRaDec, rangeRaDec=None)
Definition: makeRaDecCat.py:98
lsst::daf::persistence.butler.Butler
Definition: butler.py:321
max
int max
Definition: BoundedField.cc:104
list
daf::base::PropertyList * list
Definition: fits.cc:913
min
int min
Definition: BoundedField.cc:103
type
table::Key< int > type
Definition: Detector.cc:163
lsst::daf::persistence
Definition: Utils.h:50
lsst::afw::geom
Definition: frameSetUtils.h:40