Merge pull request #9 from lsst/tickets/DM-31226

DM-31226: Add functionality needed for WCSFitTask in drp_tasks
lsst · Dec 16, 2022 · ce05b32 · ce05b32
2 parents 00d849c + b11d828
commit ce05b32
Show file tree

Hide file tree

Showing 22 changed files with 955 additions and 949 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,10 +1,6 @@
 # Stuff I've stashed away:
 OLD
 
-# CMS additions:
-delete_prob
-*.fits
-
 # The pdfs for the paper:
 doc/WCSFit.pdf
 doc/WCSFoF.pdf
@@ -20,6 +16,8 @@ testbin/*
 
 # Standard stuff:
 build
+install
+build_cmake
 *.o
 .depend
 *.log

diff --git a/astrometry/include/PixelMap.h b/astrometry/include/PixelMap.h
@@ -118,6 +118,7 @@ namespace astrometry {
     IdentityMap(string name_="Identity"): PixelMap(name_) {}
     virtual PixelMap* duplicate() const {return new IdentityMap(*this);}
     static string type() {return "Identity";}
+    virtual string getType() const {return type();}
     void toWorld(double xpix, double ypix,
 		 double& xworld, double& yworld,
 		 double color=astrometry::NODATA) const;
@@ -256,6 +257,9 @@ namespace astrometry {
     virtual void write(YAML::Emitter& os) const;
 #endif
     virtual string getType() const {return type();}
+    Matrix33 getPixMatrix() const;
+    Matrix33 getWorldMatrix() const;
+    Vector2 getWorldPole() const;
   private:
     SphericalCoords* pix;
     SphericalCoords* world;

diff --git a/astrometry/include/PixelMapCollection.h b/astrometry/include/PixelMapCollection.h
@@ -204,6 +204,7 @@ namespace astrometry {
     // Set/get the master parameter vector for all PixelMaps
     void setParams(const DVector& p);
     DVector getParams() const;
+    std::map<std::string, astrometry::DVector> getParamDict() const;
     int nParams() const {return parameterCount;}
 
     // Set parameters of a member map by copying from
@@ -219,11 +220,14 @@ namespace astrometry {
     int nFreeMaps() const {return freeCount;}
     bool mapExists(string name) const {return mapElements.count(name);}
     vector<string> allMapNames() const;
+    string getMapType(string mapName) const;
 
     // And the WCS specified for this collection:
     int nWcs() const {return wcsElements.size();}
     bool wcsExists(string name) const {return wcsElements.count(name);}
     vector<string> allWcsNames() const;
+    // Get center orientation of a WCS element, default is in radians:
+    DVector getWcsNativeCoords(string wcsName, bool degrees=true) const;
 
     // Return whether the map given by first name depends on
     // a map of the second name.
@@ -233,6 +237,9 @@ namespace astrometry {
     // Includes self.  Assumes no dependence cycles.
     set<string> dependencies(string mapName) const;
 
+    // Produce a list giving the atomic transformation sequence needed to implement a
+    // specified map. Assumes no dependence cycles.
+    list<string> orderAtoms(string mapName) const;
 
     // This is a routine useful for debugging: return the name of the atomic
     // map that a certain parameter in the vector belongs to.
@@ -314,10 +321,6 @@ namespace astrometry {
     // Check that all referenced names exist, and that there are no circular dependences.
     void checkCompleteness() const;
 
-    // Produce a list giving the atomic transformation sequence needed to implement a
-    // specified map. Assumes no dependence cycles.
-    list<string> orderAtoms(string mapName) const;
-
     // **** Static structures / methods for serialization: ****
     // This routine will write a complete YAML key/value to emitter for this PhotoMap
     void writeSingleMap(YAML::Emitter& os, const MapElement& mel, string name) const;

diff --git a/astrometry/src/Astrometry.cpp b/astrometry/src/Astrometry.cpp
@@ -1106,16 +1106,12 @@ namespace astrometry {
   void
   Gnomonic::setLonLat(double xi, double eta) {
     x[2] = pow(1.+xi*xi+eta*eta , -0.5);
-    //if (x[2] < 0.5){
-    //   cerr << xi << "  " << eta << "  " << x[2] << endl;
-    //}
     x[0] = xi*x[2];
     x[1] = eta*x[2];
   }
   void
   Gnomonic::getLonLat(double& xi, double& eta) const {
     if (x[2]<=0.){ 
-      //cerr << x[0] << "  " << x[1] << "  " << x[2] << endl;
       throw AstrometryError("Gnomonic projection invalid for z<=0.");
     }
     xi = x[0]/x[2];

diff --git a/astrometry/src/PixelMap.cpp b/astrometry/src/PixelMap.cpp
@@ -210,6 +210,42 @@ ReprojectionMap::write(YAML::Emitter& os) const {
 }
 #endif
 
+Matrix33
+ReprojectionMap::getPixMatrix() const {
+  auto gnom = dynamic_cast<const astrometry::Gnomonic*>(pix);
+  if (gnom == nullptr) { 
+    if (dynamic_cast<const astrometry::SphericalICRS*>(pix)) {
+      Vector2 lonlat = pix->getLonLat();
+    }
+    Matrix33 tmpMatrix;
+    tmpMatrix.setToIdentity();
+    return tmpMatrix;
+  }
+  return gnom->getOrient()->m();
+}
+
+Matrix33
+ReprojectionMap::getWorldMatrix() const {
+  auto gnom = dynamic_cast<const astrometry::Gnomonic*>(world);
+  if (gnom == nullptr) {
+    Matrix33 tmpMatrix;
+    tmpMatrix.setToIdentity();
+    return tmpMatrix;
+  }
+  return gnom->getOrient()->m();
+}
+
+Vector2
+ReprojectionMap::getWorldPole() const {
+  auto gnom = dynamic_cast<const astrometry::Gnomonic*>(world);
+  if (gnom == nullptr) {
+    Vector2 lonlat;
+    return lonlat;
+  }
+  SphericalICRS pole = gnom->getOrient()->getPole();
+  Vector2 lonlat = pole.getLonLat();
+  return lonlat;
+}
 ///////////////////////////////////////////////////////////////////
 // Color term calculations
 ///////////////////////////////////////////////////////////////////

diff --git a/astrometry/src/PixelMapCollection.cpp b/astrometry/src/PixelMapCollection.cpp
@@ -135,6 +135,62 @@ PixelMapCollection::getParams() const {
   return p;
 }
 
+std::map<std::string, astrometry::DVector>
+PixelMapCollection::getParamDict() const {
+  map<string, linalg::Vector<double>> outMap;
+  for (auto& melpair : mapElements) {
+    const MapElement& map = melpair.second;
+    if (map.isFixed) continue;
+    int nSub = map.nParams;
+    if (nSub<=0) continue;
+    if (!map.atom) 
+      cerr << "mapElement is not atomic: " << melpair.first << " params: " << nSub << endl;
+    Assert(map.atom);
+    DVector p = map.atom->getParams().subVector(0,nSub);
+    string mapName = melpair.first;
+    outMap.insert(std::pair<string, astrometry::DVector>(mapName, p));
+    }
+  return outMap;
+}
+
+std::string
+PixelMapCollection::getMapType(string mapName) const {
+  std::string mapType;
+  auto melpair = mapElements.find(mapName);
+  if (melpair == mapElements.end())
+    throw AstrometryError("Could not find mapName " + mapName);
+
+  const MapElement& mel = melpair->second;
+  if (mel.atom) {
+    mapType = mel.atom->getType();
+  }
+  else {
+    mapType = SubMap::type();
+  }
+  return mapType;
+}
+
+DVector
+PixelMapCollection::getWcsNativeCoords(string wcsName, bool degrees) const {
+  auto wcspair = wcsElements.find(wcsName);
+  if (wcspair == wcsElements.end())
+    throw AstrometryError("Could not find wcsName " + wcsName);
+
+  const WcsElement& wcs = wcspair->second;
+  const SphericalCustomBase *coords = dynamic_cast<const SphericalCustomBase*>(wcs.nativeCoords);
+  if (coords == nullptr) {
+    throw AstrometryError("NativeCoords subclass does not have 'orient' property");
+  }
+  SphericalICRS pole = coords->getOrient()->getPole();
+  DVector lonLat(2);
+  pole.getRADec(lonLat[0], lonLat[1]);
+  if (degrees) {
+    lonLat[0] /= DEGREE;
+    lonLat[1] /= DEGREE;
+  }
+  return lonLat;
+}
+
 void
 PixelMapCollection::copyParamsFrom(const PixelMap& pm) {
   auto melpair = mapElements.find(pm.getName());
@@ -873,12 +929,8 @@ PixelMapCollection::read(istream& is, string namePrefix) {
       if (keyName == magicKey) continue;
 
       string name = namePrefix + keyName;
-      //cerr << i->second << endl;
+
       const YAML::Node& node=i->second;
-      //cerr << node["XMin"] << endl;
-      //cerr << node["XMin"].as<string>() << endl;
-      //string s = typeid(node["XMin"]).name();
-      //cerr << "XMin type: " << s << endl;
       // If there is a node named "WCS" it holds all of our WCS definitions
       if (keyName=="WCS") {
 	if (!node.IsMap()) 

diff --git a/astrometry/src/YAMLCollector.cpp b/astrometry/src/YAMLCollector.cpp
@@ -13,6 +13,10 @@ YAMLCollector::YAMLCollector(string specs, string magic_): magic(magic_) {
   // Split specs at commas
   auto speclist = split(specs, ',');
   for (auto spec : speclist) {
+    // Allow empty inputs:
+    if (spec.size() == 0) {
+      continue;
+    }
     // Split each spec at '@' sign
     auto regname = split(spec,'@');
     string expr = "";  // Default will be to match anything

diff --git a/gbfits/CMakeLists.txt b/gbfits/CMakeLists.txt
@@ -1,5 +1,5 @@
-add_library(gbfits SHARED FITS.cpp FTableExpression.cpp FitsTable.cpp Header.cpp Image.cpp 
+add_library(gbfits FITS.cpp FTableExpression.cpp FitsTable.cpp Header.cpp Image.cpp 
             FTable.cpp FitsImage.cpp Hdu.cpp HeaderFromStream.cpp)
 
 target_include_directories(gbfits PUBLIC ${PROJECT_SOURCE_DIR}/gbfits)
-target_link_libraries(gbfits gbutil cfitsio)
+target_link_libraries(gbfits PUBLIC gbutil cfitsio)
diff --git a/gbutil/include/Poly2d.h b/gbutil/include/Poly2d.h
@@ -52,6 +52,9 @@ namespace poly2d {
     void write(std::ostream& os, int precision=DEFAULT_PRECISION) const;  
     static Poly2d* create(std::istream& is);  // Build polynomial from serialized string
 
+    // Multiply by another polynomial:
+    Poly2d operator* (Poly2d polyB);
+
 #ifdef USE_YAML
     // Serialize to/from YAML
     void write(YAML::Emitter& os) const;

diff --git a/gbutil/src/Poly2d.cpp b/gbutil/src/Poly2d.cpp
@@ -319,3 +319,26 @@ Poly2d::fillFromVector(const DVector& v) {
   }
 }
 
+Poly2d
+Poly2d::operator*(Poly2d polyB) {
+  int orderX = this->getOrderX() + polyB.getOrderX();
+  int orderY = this->getOrderY() + polyB.getOrderY();
+  DMatrix product(orderX + 1, orderY + 1);
+  int iA, iB, jA, jB;
+  DVector AVector = getC();
+  DVector BVector = polyB.getC();
+  for (int m = 0; m < AVector.size(); m++) {
+    for (int n = 0; n < BVector.size(); n++) {
+      this->powersOfIndex(m, iA, jA);
+      polyB.powersOfIndex(n, iB, jB);
+      int outOrderX = iA + iB;
+      int outOrderY = jA + jB;
+      if ((outOrderX > orderX) || (outOrderY > orderY)) {
+        throw Poly2dError("Poly coefficients don't correspond to poly order");
+      }
+      product(outOrderX, outOrderY) = AVector[m] * BVector[n];
+    }
+  }
+  return Poly2d(product);
+}
+
diff --git a/include/FitSubroutines.h b/include/FitSubroutines.h
@@ -69,6 +69,66 @@ struct Astro {
     // saving results wants an array giving the field projection to use for each catalog.
     static void saveResults(const astrometry::MCat &matches, string outCatalog, string starCatalog,
                             vector<astrometry::SphericalCoords *> catalogProjections);
+    // writes out the final residuals catalog as in `saveResults`, but returns results as a map
+    struct outputCatalog {
+        vector<int> matchID;
+        vector<long> catalogNumber;
+        vector<long> objectNumber;
+        vector<bool> clip;
+        vector<bool> reserve;
+        vector<bool> hasPM;
+        vector<double> color;
+        vector<double> xresw;
+        vector<double> yresw;
+        vector<double> xpix;
+        vector<double> ypix;
+        vector<double> sigpix;
+        vector<double> xrespix;
+        vector<double> yrespix;
+        vector<double> xworld;
+        vector<double> yworld;
+        vector<double> covTotalW_00;
+        vector<double> covTotalW_11;
+        vector<double> covTotalW_01;
+        vector<double> chisq;
+        vector<double> chisqExpected;
+    };
+    static outputCatalog getOutputCatalog(const astrometry::MCat &matches);
+
+    struct PMCatalog {
+        vector<int> pmMatchID;
+        vector<long> pmCatalogNumber;
+        vector<long> pmObjectNumber;
+        vector<bool> pmClip;
+        vector<bool> pmReserve;
+        vector<vector<double>> pmMean;
+        vector<vector<double>> pmInvCov;
+        vector<double> pmChisq;
+        vector<double> pmChisqExpected;
+    };
+    // writes out the residuals for PMDetections as in `saveResults`, but returns results as a map
+    static PMCatalog getPMCatalog(const astrometry::MCat &matches);
+
+    struct StarCatalog {
+        vector<int> starMatchID;
+        vector<bool> starReserve;
+        vector<double> starColor;
+        vector<int> starPMCount;
+        vector<int> starDetCount;
+        vector<int> starClipCount;
+        vector<int> starDOF;
+        vector<double> starChisq;
+        vector<double> starX;
+        vector<double> starY;
+        vector<double> starPMx;
+        vector<double> starPMy;
+        vector<double> starParallax;
+        vector<vector<double>> starInvCov;
+    };
+    // writes out the catalog of fit star positions (and proper motions if fit) as in `saveResults`, but
+    // returns results as a map
+    static StarCatalog getStarCatalog(const astrometry::MCat &matches,
+                                      vector<astrometry::SphericalCoords *> catalogProjections);
 
     static void reportStatistics(const astrometry::MCat &matches,
                                  const vector<unique_ptr<Exposure>> &exposures,

diff --git a/include/Instrument.h b/include/Instrument.h
@@ -81,15 +81,16 @@ class ExposuresHelper {
 public:
     ExposuresHelper(vector<string> expNames_, vector<int> fieldNumbers_, vector<int> instrumentNumbers_,
                     vector<double> ras_, vector<double> decs_, vector<double> airmasses_,
-                    vector<double> exposureTimes_, vector<double> mjds_)
+                    vector<double> exposureTimes_, vector<double> mjds_, vector<vector<double>> observatories_)
             : expNames(expNames_),
               fieldNumbers(fieldNumbers_),
               instrumentNumbers(instrumentNumbers_),
               ras(ras_),
               decs(decs_),
               airmasses(airmasses_),
               exposureTimes(exposureTimes_),
-              mjds(mjds_){};
+              mjds(mjds_),
+              observatories(observatories_){};
     vector<string> expNames;
     // Index of the exposures' fields:
     vector<int> fieldNumbers;
@@ -100,8 +101,9 @@ class ExposuresHelper {
     vector<double> airmasses;
     vector<double> exposureTimes;
     vector<double> mjds;
+    vector<vector<double>> observatories;
 
-    std::vector<std::unique_ptr<Exposure>> getExposuresVector() {
+    std::vector<std::unique_ptr<Exposure>> getExposuresVector(const vector<double> &fieldEpochs = std::vector<double>()) {
         std::vector<std::unique_ptr<Exposure>> tmpExposures;
         tmpExposures.reserve(expNames.size());
         for (int i = 0; i < expNames.size(); i++) {
@@ -113,6 +115,14 @@ class ExposuresHelper {
             expo->airmass = airmasses[i];
             expo->exptime = exposureTimes[i];
             expo->mjd = mjds[i];
+            if (fieldEpochs.size() > 0) {
+                // Also calculate time in years after field's reference epoch
+                astrometry::UT ut;
+                ut.setMJD(expo->mjd);
+                // Note that getTTyr returns years since J2000.
+                expo->pmTDB = ut.getTTyr() - (fieldEpochs[expo->field] - 2000.);
+            }
+            for (int j = 0; j < 3; ++j) expo->observatory[j] = observatories[i][j];
             tmpExposures.emplace_back(std::move(expo));
         }
         return tmpExposures;