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
Classes | Functions
lsst::meas::astrom::detail Namespace Reference

Classes

struct  MagColInfo
 
struct  IndexManager
 

Functions

lsst::afw::table::SimpleCatalog getCatalogImpl (std::vector< index_t * > inds, lsst::afw::coord::Coord const &ctrCoord, lsst::afw::geom::Angle const &radius, const char *idCol, std::vector< MagColInfo > const &magColInfoList, const char *starGalCol, const char *varCol, bool uniqueIds=true)
 

Function Documentation

afwTable::SimpleCatalog lsst::meas::astrom::detail::getCatalogImpl ( std::vector< index_t * >  inds,
lsst::afw::coord::Coord const &  ctrCoord,
lsst::afw::geom::Angle const &  radius,
const char *  idCol,
std::vector< MagColInfo > const &  magColInfoList,
const char *  starGalCol,
const char *  varCol,
bool  uniqueIds = true 
)

Implementation for index_t::getCatalog method

Parameters
[in]indsstar kd-trees from astrometry.net
[in]ctrCoordcenter of search region
[in]radiussearch radius
[in]idColname of ID column in astrometry.net data
[in]magColInfoListlist of information about magnitude columns in astrometry.net data
[in]starGalColname of "starGal" column (true if object is a star) in astrometry.net data
[in]varColname of "var" column (true if brightness is variable) in astrometry.net data
[in]uniqueIdsif true then only return unique IDs (the first of each seen)

Returned schema:

  • id
  • coord: sky position (an lsst::afw::coord::IcrsCoord)
  • centroid: centroid on some exposure, if relevant (an lsst::afw::geom::Point2D); returned value is not set
  • hasCentroid: if true then centroid has been set; returned value is false
  • filterName_flux: flux in the specified filter (double)
  • filterName_fluxSigma: flux uncertainty in the specified filter (double)
  • resolved (if starGalCol specified): true if object is not resolved
  • variable (if varCol specified): true if brightness is variable
  • photometric: true if not resolved (or starGalCol blank) and not variable (or varCol blank); note that if starGalCol and varCol both blank then all objects are claimed to be photometric

Definition at line 58 of file utils.cc.

66 {
67  /*
68  If uniqueIds == true: return only reference sources with unique IDs;
69  arbitrarily keep the first star found with each ID.
70  */
71 
72  size_t const nMag = magColInfoList.size(); /* number of magnitude[error] columns */
73 
74  for (auto mc = magColInfoList.cbegin(); mc != magColInfoList.cend(); ++mc) {
75  if (mc->filterName.empty()) {
76  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
77  "Magnitude names cannot be empty strings.");
78  }
79  // We enforce this condition because we convert the mags to fluxes, and
80  // we need a flux to compute a flux error!
81  if (mc->magCol.empty()) {
82  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
83  "Magnitude column names cannot be empty strings.");
84  }
85  //printf("mag col \"%s\", \"%s\", \"%s\"\n", mc->filterName.c_str(), mc->magCol.c_str(), mc->magErrCol.c_str());
86  }
87 
88  auto icrsCoord = ctrCoord.toIcrs();
89  double raDeg = icrsCoord.getLongitude().asDegrees();
90  double decDeg = icrsCoord.getLatitude().asDegrees();
91  double xyz[3];
92  radecdeg2xyzarr(raDeg, decDeg, xyz);
93  double r2 = deg2distsq(radius.asDegrees());
94 
96 
97  afw::table::PointKey<double>::addFields(schema, "centroid",
98  "centroid on some exposure; invalid unless \"hasCentroid\" is true)", "pixels");
99  auto hasCentroidKey = schema.addField<afwTable::Flag>("hasCentroid",
100  "true if centroid field has been set");
101 
102  std::vector<afwTable::Key<double> > fluxKey; // these are double for consistency with measured fluxes;
103  std::vector<afwTable::Key<double> > fluxErrKey; // double may be unnecessary, but less surprising.
104  fluxKey.reserve(nMag);
105  fluxErrKey.reserve(nMag);
106 
107  for (auto mc = magColInfoList.cbegin(); mc != magColInfoList.cend(); ++mc) {
108  // Add schema elements for each requested flux (and optionally flux error)
109  // avoid the comment "flux flux"
110  fluxKey.push_back(
111  schema.addField<double>(
112  mc->filterName + "_flux",
113  mc->filterName + " flux"));
114  if (mc->hasErr()) {
115  fluxErrKey.push_back(
116  schema.addField<double>(
117  mc->filterName + "_fluxSigma",
118  mc->filterName + " flux uncertainty (sigma)"));
119  }
120  }
121 
122  afwTable::Key<afwTable::Flag> resolvedKey;
123  if (isStarCol) {
124  resolvedKey = schema.addField<afwTable::Flag>(
125  "resolved",
126  "set if the reference object is resolved");
127  }
128  afwTable::Key<afwTable::Flag> variableKey;
129  if (isVarCol) {
130  variableKey = schema.addField<afwTable::Flag>(
131  "variable",
132  "set if the reference object is variable");
133  }
134  afwTable::Key<afwTable::Flag> photometricKey = schema.addField<afwTable::Flag>(
135  "photometric",
136  "set if the reference object can be used in photometric calibration");
137 
139  if (idCol) {
140  // make catalog with no IdFactory, since IDs are external
142  } else {
143  // let the catalog assign IDs
145  }
146 
147  // for uniqueIds: keep track of the IDs we have already added to the result set.
148  std::set<boost::int64_t> uids;
149 
150  for (std::vector<index_t*>::iterator pind = inds.begin(); pind != inds.end(); ++pind) {
151  index_t* ind = (*pind);
152  // Find nearby stars
153  double *radecs = NULL;
154  int *starinds = NULL;
155  int nstars = 0;
156  startree_search_for(ind->starkd, xyz, r2, NULL, &radecs, &starinds, &nstars);
157  //printf("found %i in \"%s\"\n", nstars, ind->indexname);
158  if (nstars == 0) {
159  continue;
160  }
161 
162  std::vector<float*> mag;
163  std::vector<float*> magErr;
164  mag.reserve(nMag);
165  magErr.reserve(nMag);
166  boost::int64_t* id = NULL;
167  bool* stargal = NULL;
168  bool* var = NULL;
169  if (idCol || nMag || isStarCol || isVarCol) {
170  fitstable_t* tag = startree_get_tagalong(ind->starkd);
171  tfits_type flt = fitscolumn_float_type();
172  tfits_type boo = fitscolumn_boolean_type();
173  tfits_type i64 = fitscolumn_i64_type();
174 
175  if (!tag) {
176  std::string msg = boost::str(boost::format(
177  "astrometry_net_data index file %s does not contain a tag-along table, "
178  "so can't retrieve extra columns. idCol=%s, isStarCol=%s, isVarCol=%s") %
179  ind->indexname % idCol % isStarCol % isVarCol);
180  msg += ", mag columns=[";
181  for (unsigned int i=0; i<nMag; i++) {
182  if (i) {
183  msg += ",";
184  }
185  msg += " name='" + magColInfoList[i].filterName +
186  "', mag='" + magColInfoList[i].magCol +
187  "', magErr='" + magColInfoList[i].magErrCol + "'";
188  }
189  msg += " ]. You may need to edit the $ASTROMETRY_NET_DATA_DIR/andConfig.py file to set idColumn=None, etc.";
190  throw LSST_EXCEPT(lsst::pex::exceptions::NotFoundError, msg);
191  }
192 
193  if (idCol) {
194  id = static_cast<int64_t*>(fitstable_read_column_inds(tag, idCol, i64, starinds, nstars));
195  if (!id) {
196  throw LSST_EXCEPT(lsst::pex::exceptions::NotFoundError,
197  str(boost::format("Unable to read data for %s from %s") % idCol % ind->indexname));
198  }
199  }
200  if (id && uniqueIds) {
201  // remove duplicate IDs.
202 
203  // FIXME -- this shouldn't be necessary once we get astrometry_net 0.40
204  // multi-index functionality in place.
205 
206  if (uids.empty()) {
207  uids = std::set<boost::int64_t>(id, id+nstars);
208  } else {
209  int nkeep = 0;
210  for (int i=0; i<nstars; i++) {
211  //std::pair<std::set<boost::int64_t>::iterator, bool>
212  if (uids.insert(id[i]).second) {
213  // inserted; keep this one.
214  if (nkeep != i) {
215  // compact the arrays.
216  starinds[nkeep] = starinds[i];
217  radecs[nkeep*2+0] = radecs[i*2+0];
218  radecs[nkeep*2+1] = radecs[i*2+1];
219  id[nkeep] = id[i];
220  }
221  nkeep++;
222  } else {
223  // did not insert (this id has already been found);
224  // drop this star.
225  }
226  }
227  nstars = nkeep;
228  // if they were all duplicate IDs...
229  if (nstars == 0) {
230  free(starinds);
231  free(radecs);
232  free(id);
233  continue;
234  }
235  }
236  }
237 
238  for (auto mc = magColInfoList.cbegin(); mc != magColInfoList.cend(); ++mc) {
239  char const* col = mc->magCol.c_str();
240  mag.push_back(read_column(tag, col, flt, starinds, nstars, ind->indexname));
241  if (mc->hasErr()) {
242  char const* col = mc->magErrCol.c_str();
243  magErr.push_back(read_column(tag, col, flt, starinds, nstars, ind->indexname));
244  }
245  }
246  if (isStarCol) {
247  /* There is something weird going on with handling of bools; maybe "T" vs "F"?
248  stargal = static_cast<bool*>(fitstable_read_column_inds(tag, isStarCol, boo, starinds, nstars));
249  for (int j=0; j<nstars; j++) {
250  printf(" sg %i = %i, vs %i\n", j, (int)sg[j], stargal[j] ? 1:0);
251  }
252  */
253  uint8_t* sg = static_cast<uint8_t*>(fitstable_read_column_inds(
254  tag, isStarCol, fitscolumn_u8_type(), starinds, nstars));
255  stargal = static_cast<bool*>(malloc(nstars));
256  if (!stargal) {
257  throw LSST_EXCEPT(lsst::pex::exceptions::NotFoundError,
258  str(boost::format("Unable to read data for %s from %s") % isStarCol % ind->indexname));
259  }
260  for (int j=0; j<nstars; j++) {
261  stargal[j] = (sg[j] > 0);
262  }
263  free(sg);
264  }
265  if (isVarCol) {
266  var = static_cast<bool*>(fitstable_read_column_inds(tag, isVarCol, boo, starinds, nstars));
267  if (!var) {
268  throw LSST_EXCEPT(lsst::pex::exceptions::NotFoundError,
269  str(boost::format("Unable to read data for %s from %s") % isVarCol % ind->indexname));
270  }
271  }
272  }
273 
274  for (int i=0; i<nstars; i++) {
275  PTR(afwTable::SimpleRecord) src = cat.addNew();
276 
277  // Note that all coords in afwTable catalogs are ICRS; hopefully that's what the
278  // reference catalogs are (and that's what the code assumed before JFB modified it).
279  src->setCoord(
280  lsst::afw::coord::IcrsCoord(
281  radecs[i * 2 + 0] * afwGeom::degrees,
282  radecs[i * 2 + 1] * afwGeom::degrees
283  )
284  );
285 
286  if (id) {
287  src->setId(id[i]);
288  }
289 
290  src->set(hasCentroidKey, false);
291 
292  assert(fluxKey.size() == nMag);
293  // only non-empty error columns are populated in these vectors.
294  assert(fluxErrKey.size() == magErr.size());
295  // index into non-empty error columns.
296  size_t ej = 0;
297  size_t j = 0;
298  for (auto mc = magColInfoList.cbegin(); mc != magColInfoList.cend(); ++mc, ++j) {
299  // There is some dispute about returning flux or magnitudes
300  // flux is not conventional, but is convenient for several reasons
301  // - it simplifies matching to sources for astrometric calibration
302  // - flux errors are easier to combine than magnitude errors
303  //printf("mag %s = %g\n", mc->filterName.c_str(), mag[j][i]);
304  double flux = lsst::afw::image::fluxFromABMag(mag[j][i]);
305  //printf("flux %g\n", flux);
306  src->set(fluxKey[j], flux);
307  if (mc->hasErr()) {
308  //printf("mag err = %g\n", magErr[ej][i]);
309  double fluxErr = lsst::afw::image::fluxErrFromABMagErr(magErr[ej][i], mag[j][i]);
310  //printf("flux err = %g\n", fluxErr);
311  src->set(fluxErrKey[ej], fluxErr);
312  ej++;
313  }
314  }
315  assert(ej == fluxErrKey.size());
316 
317  bool photometric = true;
318  if (stargal) {
319  src->set(resolvedKey, !stargal[i]);
320  photometric &= stargal[i];
321  }
322  if (var) {
323  src->set(variableKey, var[i]);
324  photometric &= (!var[i]);
325  }
326  src->set(photometricKey, photometric);
327  }
328 
329  free(id);
330  for (size_t j=0; j<mag.size(); ++j) {
331  free(mag[j]);
332  }
333  for (size_t j=0; j<magErr.size(); ++j) {
334  free(magErr[j]);
335  }
336  free(stargal);
337  free(var);
338  free(radecs);
339  free(starinds);
340  }
341  return cat;
342 }
Defines the fields and offsets for a table.
Definition: Schema.h:46
#define PTR(...)
Definition: base.h:41
static Schema makeMinimalSchema()
Return a minimal schema for Simple tables and records.
Definition: Simple.h:121
double fluxErrFromABMagErr(double magErr, double mag)
Compute flux error in Janskys from AB magnitude error and AB magnitude.
Definition: Calib.h:74
static boost::shared_ptr< SimpleTable > make(Schema const &schema, boost::shared_ptr< IdFactory > const &idFactory)
Construct a new table.
lsst::afw::coord::IcrsCoord IcrsCoord
Definition: misc.h:38
double fluxFromABMag(double mag)
Compute flux in Janskys from AB magnitude.
Definition: Calib.h:69
AngleUnit const degrees
Definition: Angle.h:92
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:55
SortedCatalogT< SimpleRecord > SimpleCatalog
Definition: fwd.h:69
if(width!=gim.getWidth()||height!=gim.getHeight()||x0!=gim.getX0()||y0!=gim.getY0())
Definition: saturated.cc:47
tbl::Schema schema
A polymorphic functor base class for generating record IDs for a table.
Definition: IdFactory.h:19
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
A class used as a handle to a particular field in a table.
Definition: fwd.h:44
int id
Definition: CR.cc:151
Record class that must contain a unique ID field and a celestial coordinate field.
Definition: Simple.h:46
int col
Definition: CR.cc:152