LSSTApplications  17.0+103,17.0+11,17.0+61,18.0.0+13,18.0.0+25,18.0.0+5,18.0.0+52,18.0.0-4-g68ffd23,18.1.0-1-g0001055+8,18.1.0-1-g03d53ef+1,18.1.0-1-g1349e88+28,18.1.0-1-g2505f39+22,18.1.0-1-g380d4d4+27,18.1.0-1-g5315e5e+1,18.1.0-1-g5e4b7ea+10,18.1.0-1-g7e8fceb+1,18.1.0-1-g85f8cd4+23,18.1.0-1-g9a6769a+13,18.1.0-1-ga1a4c1a+22,18.1.0-1-gd55f500+17,18.1.0-12-g42eabe8e+10,18.1.0-14-gd04256d+15,18.1.0-16-g430f6a53+1,18.1.0-17-gd2166b6e4,18.1.0-18-gb5d19ff+1,18.1.0-2-gfbf3545+7,18.1.0-2-gfefb8b5+16,18.1.0-3-g52aa583+13,18.1.0-3-g62b5e86+14,18.1.0-3-g8f4a2b1+17,18.1.0-3-g9bc06b8+7,18.1.0-3-gb69f684+9,18.1.0-4-g1ee41a7+1,18.1.0-5-g6dbcb01+13,18.1.0-5-gc286bb7+3,18.1.0-6-g48bdcd3+2,18.1.0-6-gd05e160+9,18.1.0-7-gc4d902b+2,18.1.0-7-gebc0338+8,18.1.0-9-gae7190a+10,w.2019.38
LSSTDataManagementBasePackage
Classes | Functions
lsst::meas::extensions::astrometryNet::detail Namespace Reference

Classes

struct  IndexManager
 RAII manager for astrometry.net indices. More...
 
struct  MagColInfo
 

Functions

lsst::afw::table::SimpleCatalog getCatalogImpl (std::vector< index_t *> inds, lsst::geom::SpherePoint const &ctrCoord, lsst::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. More...
 

Function Documentation

◆ getCatalogImpl()

afwTable::SimpleCatalog lsst::meas::extensions::astrometryNet::detail::getCatalogImpl ( std::vector< index_t *>  inds,
lsst::geom::SpherePoint const &  ctrCoord,
lsst::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: ICRS sky position (an lsst::geom::SpherePoint)
  • centroid: centroid on some exposure, if relevant (an lsst::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_fluxErr: 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 61 of file utils.cc.

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