11 #include "boost/regex.hpp" 12 #include "boost/multi_index_container.hpp" 13 #include "boost/multi_index/sequenced_index.hpp" 14 #include "boost/multi_index/ordered_index.hpp" 15 #include "boost/multi_index/hashed_index.hpp" 16 #include "boost/multi_index/member.hpp" 35 template <std::
string FitsSchemaItem::*Member>
36 struct SetFitsSchemaString {
37 void operator()(FitsSchemaItem &item) { item.*Member = _v; }
38 explicit SetFitsSchemaString(
std::string const &v) : _v(v) {}
52 typedef boost::multi_index_container<
54 boost::multi_index::indexed_by<
55 boost::multi_index::ordered_non_unique<
56 boost::multi_index::member<FitsSchemaItem, int, &FitsSchemaItem::column>>,
57 boost::multi_index::ordered_non_unique<
58 boost::multi_index::member<FitsSchemaItem, int, &FitsSchemaItem::bit>>,
59 boost::multi_index::hashed_unique<
60 boost::multi_index::member<FitsSchemaItem, std::string, &FitsSchemaItem::ttype>>,
61 boost::multi_index::sequenced<>>>
65 typedef SetFitsSchemaString<&FitsSchemaItem::ttype>
SetTTYPE;
66 typedef SetFitsSchemaString<&FitsSchemaItem::tform>
SetTFORM;
67 typedef SetFitsSchemaString<&FitsSchemaItem::tccls>
SetTCCLS;
68 typedef SetFitsSchemaString<&FitsSchemaItem::tunit>
SetTUNIT;
69 typedef SetFitsSchemaString<&FitsSchemaItem::doc>
SetDoc;
78 ByColumn &
byColumn() {
return inputs.get<0>(); }
79 ByBit &
byBit() {
return inputs.get<1>(); }
80 ByName &
byName() {
return inputs.get<2>(); }
81 AsList &
asList() {
return inputs.get<3>(); }
101 : _impl(
std::make_shared<
Impl>()) {
104 if (!metadata.
exists(
"AFW_TYPE")) {
107 _impl->version = metadata.
get(
"AFW_TABLE_VERSION", _impl->version);
108 _impl->type = metadata.
get(
"AFW_TYPE", _impl->type);
110 metadata.
remove(
"AFW_TABLE_VERSION");
113 metadata.
remove(
"AFW_TYPE");
117 _impl->archiveHdu = metadata.
get(
"AR_HDU", -1);
118 if (_impl->archiveHdu > 0) {
121 metadata.
remove(
"AR_HDU");
130 if (pos == std::string::npos) {
132 (
boost::format(
"Malformed alias definition: '%s'") % (*i)).str());
134 _impl->schema.getAliasMap()->set(i->substr(0, pos), i->substr(pos + 1, std::string::npos));
143 if (_impl->version == 0) {
154 if (!target.
empty()) {
155 _impl->schema.getAliasMap()->set(oldSlotKeys[i].
second, target);
158 metadata.
remove(oldSlotKeys[i].first +
"_ERR_SLOT");
159 metadata.
remove(oldSlotKeys[i].first +
"_FLAG_SLOT");
168 if (
key->compare(0, 5,
"TTYPE") == 0) {
170 auto iter = _impl->byColumn().lower_bound(column);
171 if (iter == _impl->byColumn().end() || iter->column != column) {
172 iter = _impl->byColumn().insert(iter,
FitsSchemaItem(column, -1));
176 if (iter->doc.empty()) {
182 }
else if (
key->compare(0, 5,
"TFLAG") == 0) {
184 auto iter = _impl->byBit().lower_bound(bit);
185 if (iter == _impl->byBit().end() || iter->bit != bit) {
190 if (iter->doc.empty()) {
196 }
else if (
key->compare(0, 4,
"TDOC") == 0) {
198 auto iter = _impl->byColumn().lower_bound(column);
199 if (iter == _impl->byColumn().end() || iter->column != column) {
200 iter = _impl->byColumn().insert(iter,
FitsSchemaItem(column, -1));
206 }
else if (
key->compare(0, 5,
"TFDOC") == 0) {
208 auto iter = _impl->byBit().lower_bound(bit);
209 if (iter == _impl->byBit().end() || iter->bit != bit) {
216 }
else if (
key->compare(0, 5,
"TUNIT") == 0) {
218 auto iter = _impl->byColumn().lower_bound(column);
219 if (iter == _impl->byColumn().end() || iter->column != column) {
220 iter = _impl->byColumn().insert(iter,
FitsSchemaItem(column, -1));
226 }
else if (
key->compare(0, 5,
"TCCLS") == 0) {
228 auto iter = _impl->byColumn().lower_bound(column);
229 if (iter == _impl->byColumn().end() || iter->column != column) {
230 iter = _impl->byColumn().insert(iter,
FitsSchemaItem(column, -1));
236 }
else if (
key->compare(0, 5,
"TFORM") == 0) {
238 auto iter = _impl->byColumn().lower_bound(column);
239 if (iter == _impl->byColumn().end() || iter->column != column) {
240 iter = _impl->byColumn().insert(iter,
FitsSchemaItem(column, -1));
246 }
else if (
key->compare(0, 5,
"TZERO") == 0) {
250 }
else if (
key->compare(0, 5,
"TSCAL") == 0) {
254 }
else if (
key->compare(0, 5,
"TNULL") == 0) {
258 }
else if (
key->compare(0, 5,
"TDISP") == 0) {
266 _impl->flagColumn = metadata.
get(
"FLAGCOL", 0);
267 if (_impl->flagColumn > 0) {
269 metadata.
remove(
"FLAGCOL");
272 auto iter = _impl->byColumn().find(_impl->flagColumn);
273 if (
iter == _impl->byColumn().end()) {
276 (
boost::format(
"Column for flag data not found; FLAGCOL=%d") % _impl->flagColumn).str());
281 static boost::regex
const regex(
"(\\d+)?X\\(?(\\d)*\\)?", boost::regex::perl);
283 if (!boost::regex_match(
iter->tform, m, regex)) {
286 (
boost::format(
"Invalid TFORM key for flags column: '%s'") %
iter->tform).str());
292 _impl->flagKeys.resize(nFlags);
293 _impl->flagWorkspace.reset(
new bool[nFlags]);
296 _impl->byColumn().erase(
iter);
309 int oldHdu = fits.
getHdu();
310 if (_impl->archiveHdu < 0) _impl->archiveHdu = oldHdu + 1;
312 fits.
setHdu(_impl->archiveHdu);
319 _impl->archiveHdu = -1;
327 auto iter = _impl->byName().find(ttype);
328 if (
iter == _impl->byName().end()) {
335 auto iter = _impl->byColumn().lower_bound(column);
336 if (
iter == _impl->byColumn().end() ||
iter->column != column) {
343 auto iter = _impl->byColumn().lower_bound(item->
column);
344 assert(
iter != _impl->byColumn().end() &&
iter->column == item->
column);
345 _impl->byColumn().erase(
iter);
349 auto iter = _impl->byName().find(ttype);
350 if (
iter != _impl->byName().end() &&
iter->ttype == ttype) {
351 _impl->byName().erase(
iter);
356 auto iter = _impl->byColumn().lower_bound(column);
357 if (
iter != _impl->byColumn().end() &&
iter->column == column) {
358 _impl->byColumn().erase(
iter);
362 void erase(
int column);
365 _impl->readers.push_back(
std::move(reader));
370 template <
typename T>
380 _cache(), _cacheFirstRow(0)
387 if (_key.getElementCount() == 1u) {
388 std::size_t nElements = nRows*_key.getElementCount();
389 _cache.resize(nElements);
390 _cacheFirstRow = firstRow;
391 fits.
readTableArray(firstRow, _column, nElements, &_cache.front());
397 if (_cache.empty()) {
400 assert(row >= _cacheFirstRow);
402 assert(offset < _cache.size());
431 "Angle fields must be persisted in radians (TUNIT='rad').");
436 assert(_key.getElementCount() == 1u);
437 _cache.resize(nRows);
438 _cacheFirstRow = firstRow;
444 if (_cache.empty()) {
449 assert(row >= _cacheFirstRow);
451 assert(offset < _cache.size());
472 _isVariableLength(size == 0) {}
484 bool _isVariableLength;
487 template <
typename T>
500 ndarray::Array<T, 1, 1> array = ndarray::allocate(size);
502 record.
set(_key, array);
512 template <
typename T>
584 template <
typename T,
int N>
588 static boost::regex
const regex(
"(.*)(\\^(\\d+))?", boost::regex::perl);
590 if (!boost::regex_match(oldUnits, m, regex)) {
613 for (
int i = 0; i <
_size; ++i) {
614 for (
int j = i; j <
_size; ++j) {
631 static boost::regex
const regex(
"(\\d+)?([PQ])?(\\u)\\(?(\\d)*\\)?", boost::regex::perl);
634 if (!boost::regex_match(item.
tform, m, regex)) {
641 char code = m[3].str()[0];
651 if (item.
tccls ==
"Array") {
652 return StandardReader<Array<std::uint8_t>>::make(schema, item, size);
654 return StandardReader<std::uint8_t>::make(schema, item);
657 return VariableLengthArrayReader<std::uint8_t>::make(schema, item);
659 return StandardReader<Array<std::uint8_t>>::make(schema, item, size);
664 if (item.
tccls ==
"Array") {
665 return StandardReader<Array<std::uint16_t>>::make(schema, item, size);
667 return StandardReader<std::uint16_t>::make(schema, item);
670 return VariableLengthArrayReader<std::uint16_t>::make(schema, item);
672 return StandardReader<Array<std::uint16_t>>::make(schema, item, size);
675 return VariableLengthArrayReader<std::int32_t>::make(schema, item);
677 if (item.
tccls ==
"Point") {
678 return PointConversionReader<std::int32_t>::make(schema, item);
680 if (size > 1 || item.
tccls ==
"Array") {
681 return StandardReader<Array<std::int32_t>>::make(schema, item, size);
683 return StandardReader<std::int32_t>::make(schema, item);
686 return StandardReader<std::int64_t>::make(schema, item);
690 return VariableLengthArrayReader<float>::make(schema, item);
693 if (item.
tccls ==
"Array") {
694 return StandardReader<Array<float>>::make(schema, item, 1);
698 return StandardReader<float>::make(schema, item);
700 if (size == 3 && item.
tccls ==
"Covariance(Point)") {
702 return CovarianceConversionReader<float, 2>::make(schema, item, names);
704 if (size == 6 && item.
tccls ==
"Covariance(Moments)") {
706 return CovarianceConversionReader<float, 3>::make(schema, item, names);
708 if (item.
tccls ==
"Covariance") {
709 double v = 0.5 * (
std::sqrt(1 + 8 * size) - 1);
711 if (n * (n + 1) != size * 2) {
715 for (
int i = 0; i < n; ++i) {
718 return CovarianceConversionReader<float, Eigen::Dynamic>::make(schema, item, names);
720 return StandardReader<Array<float>>::make(schema, item, size);
723 return VariableLengthArrayReader<double>::make(schema, item);
726 if (item.
tccls ==
"Angle") {
727 return AngleReader::make(schema, item);
729 if (item.
tccls ==
"Array") {
730 return StandardReader<Array<double>>::make(schema, item, 1);
732 return StandardReader<double>::make(schema, item);
735 if (item.
tccls ==
"Point") {
736 return PointConversionReader<double>::make(schema, item);
738 if (item.
tccls ==
"Coord") {
739 return CoordConversionReader::make(schema, item);
742 if (size == 3 && item.
tccls ==
"Moments") {
743 return MomentsConversionReader::make(schema, item);
745 return StandardReader<Array<double>>::make(schema, item, size);
748 return StringReader::make(schema, item, size);
763 if (!includes(item.
ttype,
"flux"))
return false;
772 return includes(units,
"count") || includes(units,
"dn") || includes (units,
"adu");
783 if (_impl->version == 0) {
784 AliasMap &aliases = *_impl->schema.getAliasMap();
785 for (
auto iter = _impl->asList().begin();
iter != _impl->asList().end(); ++
iter) {
787 if (flagPos != std::string::npos) {
796 ttype.
replace(flagPos, 5,
"flag");
802 aliases.
set(prefix +
"flags", prefix +
"flag");
803 }
else if (isInstFlux(*
iter)) {
805 if (endswith(
iter->ttype,
"_err")) {
806 aliases.
set(replace(
iter->ttype,
"_err",
"_instFluxErr"),
iter->ttype);
808 aliases.
set(
iter->ttype +
"_instFlux",
iter->ttype);
810 }
else if (endswith(
iter->ttype,
"_err")) {
816 if (
iter->tccls ==
"Covariance(Point)") {
817 aliases.
set(replace(
iter->ttype,
"_err",
"_yErr"),
iter->ttype +
"_yErr");
818 aliases.
set(replace(
iter->ttype,
"_err",
"_xErr"),
iter->ttype +
"_xErr");
819 aliases.
set(replace(
iter->ttype,
"_err",
"_x_y_Cov"),
iter->ttype +
"_x_y_Cov");
820 }
else if (
iter->tccls ==
"Covariance(Moments)") {
821 aliases.
set(replace(
iter->ttype,
"_err",
"_xxErr"),
iter->ttype +
"_xxErr");
822 aliases.
set(replace(
iter->ttype,
"_err",
"_yyErr"),
iter->ttype +
"_yyErr");
823 aliases.
set(replace(
iter->ttype,
"_err",
"_xyErr"),
iter->ttype +
"_xyErr");
824 aliases.
set(replace(
iter->ttype,
"_err",
"_xx_yy_Cov"),
iter->ttype +
"_xx_yy_Cov");
825 aliases.
set(replace(
iter->ttype,
"_err",
"_xx_xy_Cov"),
iter->ttype +
"_xx_xy_Cov");
826 aliases.
set(replace(
iter->ttype,
"_err",
"_yy_xy_Cov"),
iter->ttype +
"_yy_xy_Cov");
830 }
else if (_impl->version < 3) {
834 AliasMap &aliases = *_impl->schema.getAliasMap();
835 for (
auto iter = _impl->asList().begin();
iter != _impl->asList().end(); ++
iter) {
837 if (_impl->version < 2 && endswith(name,
"Sigma")) {
838 name = replace(
std::move(name),
"Sigma",
"Err");
840 if (_impl->version < 3 && isInstFlux(*
iter)) {
841 name = replace(
std::move(name),
"flux",
"instFlux");
843 if (name !=
iter->ttype) {
844 aliases.
set(name,
iter->ttype);
848 for (
auto iter = _impl->asList().begin();
iter != _impl->asList().end(); ++
iter) {
852 _impl->readers.push_back(
std::move(reader));
854 LOGLS_WARN(
"afw.FitsSchemaInputMapper",
"Format " <<
iter->tform <<
" for column " 856 <<
" not supported; skipping.");
859 if (static_cast<std::size_t>(
iter->bit) >= _impl->flagKeys.size()) {
861 (
boost::format(
"Flag field '%s' is is in bit %d (0-indexed) of only %d") %
862 iter->ttype %
iter->bit % _impl->flagKeys.size())
865 _impl->flagKeys[
iter->bit] = _impl->schema.addField<Flag>(
iter->ttype,
iter->doc);
868 _impl->asList().clear();
869 if (_impl->schema.getRecordSize() <= 0) {
872 (
boost::format(
"Non-positive record size: %d; file is corrupt or invalid.") %
873 _impl->schema.getRecordSize()).str()
877 return _impl->schema;
881 if (!_impl->flagKeys.empty()) {
882 fits.
readTableArray<
bool>(
row, _impl->flagColumn, _impl->flagKeys.size(), _impl->flagWorkspace.get());
883 for (
std::size_t bit = 0; bit < _impl->flagKeys.size(); ++bit) {
884 record.
set(_impl->flagKeys[bit], _impl->flagWorkspace[bit]);
887 if (_impl->nRowsToPrep != 1 && row % _impl->nRowsToPrep == 0) {
891 for (
auto & reader : _impl->readers) {
892 reader->prepRead(row, size, fits);
895 for (
auto const & reader : _impl->readers) {
896 reader->readCell(record, row, fits, _impl->archive);
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
An ellipse core with quadrupole moments as parameters.
Defines the fields and offsets for a table.
std::size_t countRows()
Return the number of row in a table.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Field base class default implementation (used for numeric scalars and lsst::geom::Angle).
table::Key< table::Array< int > > _size
A coordinate class intended to represent absolute positions (2-d specialization). ...
Key< Flag > const & target
Class for storing ordered metadata with comments.
Reports attempts to exceed implementation-defined length limits for some classes. ...
std::string const & getComment(std::string const &name) const
Get the comment for a string property name (possibly hierarchical).
std::vector< T > getArray(std::string const &name) const
Get the vector of values for a property name (possibly hierarchical).
void set(std::string const &alias, std::string const &target)
Add an alias to the schema or replace an existing one.
void readTableScalar(std::size_t row, int col, T &value)
Read an array scalar from a binary table.
Mapping class that holds aliases for a Schema.
A simple struct that combines the two arguments that must be passed to most cfitsio routines and cont...
int computeCovariancePackedSize(int size)
Defines the packed size of a covariance matrices.
long getTableArraySize(int col)
Return the size of an array column.
AngleUnit constexpr radians
constant with units of radians
A class representing an angle.
A structure that describes a field as a collection of related strings read from the FITS header...
LSST DM logging module built on log4cxx.
virtual void remove(std::string const &name)
Remove all values for a property name (possibly hierarchical).
An exception thrown when problems are found when reading or writing FITS files.
Reports attempts to access elements using an invalid key.
A base class for image defects.
A FunctorKey used to get or set a lsst::geom::Point from an (x,y) pair of int or double Keys...
static QuadrupoleKey addFields(Schema &schema, std::string const &name, std::string const &doc, CoordinateType coordType=CoordinateType::PIXEL)
Add a set of quadrupole subfields to a schema and return a QuadrupoleKey that points to them...
int indexCovariance(int i, int j)
Defines the ordering of packed covariance matrices.
bool exists(std::string const &name) const
Determine if a name (possibly hierarchical) exists.
Tag types used to declare specialized field types.
void setHdu(int hdu, bool relative=false)
Set the current HDU.
Field< T >::Element * getElement(Key< T > const &key)
Return a pointer to the underlying elements of a field (non-const).
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
T get(std::string const &name) const
Get the last value for a property name (possibly hierarchical).
Base class for all records.
Point in an unspecified spherical coordinate system.
A class used as a handle to a particular field in a table.
static CoordKey addFields(afw::table::Schema &schema, std::string const &name, std::string const &doc)
Add a pair of _ra, _dec fields to a Schema, and return a CoordKey that points to them.
void readTableArray(std::size_t row, int col, int nElements, T *value)
Read an array value from a binary table.
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Polymorphic reader interface used to read different kinds of objects from one or more FITS binary tab...
int getHdu()
Return the current HDU (0-indexed; 0 is the Primary HDU).
std::vector< std::string > getOrderedNames() const
Get the list of property names, in the order they were added.
A FunctorKey used to get or set a geom::ellipses::Quadrupole from a tuple of constituent Keys...
Key< T > addField(Field< T > const &field, bool doReplace=false)
Add a new field to the Schema, and return the associated Key.
A FunctorKey used to get or set celestial coordinates from a pair of lsst::geom::Angle keys...