Add PROJ projection

CHANGELOG.md

@@ -56,6 +56,7 @@ release.
   override the timestamp style naming convention of the output cube with their
   own name; if not specified retains existing behavior [#5125](https://github.com/USGS-Astrogeology/ISIS3/issues/5162)
 - Added new parameters <b>ONERROR</b>, <b>ERRORLOG</b>, and <b>ERRORLIST</b> to <i>mosrange</i> to provide better control over error behavior and provide diagnostics when problems are encountered processing the input file list.[#3606](https://github.com/DOI-USGS/ISIS3/issues/3606)
+- Adds PROJ into ISIS, and exposes the capability with a new class called IProj. [#5317](https://github.com/DOI-USGS/ISIS3/pull/5317)
 
 ### Deprecated
 

isis/CMakeLists.txt

@@ -280,6 +280,7 @@ find_package(NN                REQUIRED)
 find_package(OpenCV    3.1.0   REQUIRED)
 find_package(PCL               REQUIRED)
 find_package(Protobuf  2.6.1   REQUIRED)
+find_package(PROJ              REQUIRED)
 find_package(Qwt       6       REQUIRED)
 find_package(SuperLU   4.3     REQUIRED)
 find_package(TIFF      4.0.0   REQUIRED)

isis/src/base/apps/cam2map/cam2map.cpp

@@ -30,7 +30,15 @@ namespace Isis {
 
     // Get the map projection file provided by the user
     Pvl userMap;
-    userMap.read(ui.GetFileName("MAP"));
+    if (ui.GetBoolean("USEPROJ")) {
+      PvlGroup mappingGroup("Mapping");
+      mappingGroup.addKeyword(PvlKeyword("ProjectionName", "IProj"));
+      mappingGroup.addKeyword(PvlKeyword("ProjStr", ui.GetAsString("PROJString")));
+      userMap.addGroup(mappingGroup);
+    }
+    else {
+      userMap.read(ui.GetFileName("MAP"));
+    }
     PvlGroup &userGrp = userMap.findGroup("Mapping", Pvl::Traverse);
 
     cam2map(&icube, userMap, userGrp, ui, log);

isis/src/base/apps/cam2map/cam2map.xml

@@ -381,6 +381,28 @@
           <item>LONSEAM</item>
         </exclusions>
       </parameter>
+
+      <parameter name="USEPROJ">
+        <type>boolean</type>
+        <default><item>false</item></default>
+        <exclusions>
+          <item>MAP</item>
+        </exclusions>
+        <inclusions>
+          <item>PROJString</item>
+        </inclusions>
+      </parameter>
+    </group>
+
+    <group name="PROJ">
+      <parameter name="PROJString">
+        <type>string</type>
+        <brief>A PROJ4 string that defines the projection</brief>
+        <description>
+          This parameters is used instead of the map file to define the projection to
+          transform into. For help with PROJ strings visit https://proj.org/en/9.3/index.html
+        </description>
+      </parameter>
     </group>
 
     <group name="Output Map Resolution">

isis/src/base/objs/IProj/IProj.cpp

@@ -0,0 +1,221 @@
+/** This is free and unencumbered software released into the public domain.
+The authors of ISIS do not claim copyright on the contents of this file.
+For more details about the LICENSE terms and the AUTHORS, you will
+find files of those names at the top level of this repository. **/
+
+/* SPDX-License-Identifier: CC0-1.0 */
+#include "IProj.h"
+
+#include <proj.h>
+
+#include "IException.h"
+#include "TProjection.h"
+#include "Pvl.h"
+#include "PvlGroup.h"
+#include "PvlKeyword.h"
+
+using namespace std;
+
+namespace Isis {
+  IProj::IProj(Pvl &label, bool allowDefaults) :
+      TProjection::TProjection(label) {
+    PvlGroup &mapGroup = label.findGroup("Mapping", Pvl::Traverse);
+    if (!mapGroup.hasKeyword("ProjStr")) {
+      QString message = "No ProjStr keyword in mapping group, either add a ProjStr or select a different projection method";
+      throw IException(IException::User, message, _FILEINFO_);
+    }
+    m_C = proj_context_create();
+
+    m_userOutputProjStr = new QString(mapGroup["ProjStr"]);
+
+    /* Create a projection. */
+    std::string projString = m_userOutputProjStr->toStdString();
+    m_outputProj = proj_create(m_C, projString.c_str());
+
+    if (0 == m_outputProj) {
+      QString msg = "Unable to create projection from [" + QString(projString.c_str()) + "]";
+      throw IException(IException::User, msg, _FILEINFO_);
+    }
+
+    PJ *outputEllipsoid = proj_get_ellipsoid(m_C, m_outputProj);
+
+    if (outputEllipsoid == 0) {
+      QString msg = "Unable to get ellipsoid from [" + *m_userOutputProjStr + "]. " + 
+      "Please add a radii definition to your proj string.";
+      throw IException(IException::User, msg, _FILEINFO_);
+    }
+
+    int res = proj_ellipsoid_get_parameters(m_C, outputEllipsoid, 
+                                            &m_equatorialRadius,
+                                            &m_polarRadius,
+                                            nullptr,
+                                            nullptr);
+    proj_destroy(outputEllipsoid);
+
+    if (res == 0) {
+      QString msg = "Unable to get ellipsoid information from [" + *m_userOutputProjStr + "]";
+      throw IException(IException::User, msg, _FILEINFO_);
+    }
+
+    m_llaProj = proj_crs_get_geodetic_crs(m_C, m_outputProj);
+
+    if (0 == m_llaProj) {
+      QString msg = "Unable to create projection from []";
+      throw IException(IException::Programmer, msg, _FILEINFO_);
+    }
+
+    m_llaProj2outputProj = proj_create_crs_to_crs_from_pj(m_C, m_llaProj, m_outputProj, 0, 0);
+    if (0 == m_llaProj2outputProj) {
+      QString msg = "Unable to create transform from [" + QString(projString.c_str()) + "]";
+      throw IException(IException::User, msg, _FILEINFO_);
+    }
+  }
+
+  //! Destroys the IProj object
+  IProj::~IProj() {
+    delete m_userOutputProjStr;
+    proj_destroy(m_llaProj);
+    proj_destroy(m_outputProj);
+    proj_destroy(m_llaProj2outputProj);
+    proj_context_destroy(m_C);
+  }
+
+  /**
+   * Returns the name of the map projection, "Proj"
+   *
+   * @return QString Name of projection, "Proj"
+   */
+  QString IProj::Name() const {
+    return "Proj";
+  }
+
+  /**
+   * This function returns the keywords that this projection uses.
+   * 
+   * We also include the generated PROJ string
+   *
+   * @return PvlGroup The keywords that this projection uses
+   */
+  PvlGroup IProj::Mapping()  {
+    PvlGroup mapping = TProjection::Mapping();
+
+    mapping.addKeyword(PvlKeyword("EquatorialRadius", toString(m_equatorialRadius, 15), "meters"), PvlContainer::InsertMode::Replace);
+    mapping.addKeyword(PvlKeyword("PolarRadius", toString(m_polarRadius, 15), "meters"), PvlContainer::InsertMode::Replace);
+    mapping.addKeyword(PvlKeyword("LatitudeType", "Planetographic"), PvlContainer::InsertMode::Replace);
+    mapping.addKeyword(PvlKeyword("LongitudeDomain", "180"), PvlContainer::InsertMode::Replace);
+
+    mapping += PvlKeyword("ProjStr", *m_userOutputProjStr);
+
+    return mapping;
+  }
+
+  /**
+   * Returns the version of the map projection
+   *
+   *
+   * @return QString Version number
+   */
+  QString IProj::Version() const {
+    return "1.0";
+  }
+
+  bool IProj::SetGround(const double lat, const double lon) {
+    m_longitude = lon;
+    m_latitude = lat;
+
+    PJ_COORD c_in;
+    c_in.lpz.lam = m_longitude;
+    // Convert to ographic as proj defaults to operating in ographic latitudes
+    c_in.lpz.phi = ToPlanetographic(m_latitude);
+
+    PJ_COORD c_out = proj_trans(m_llaProj2outputProj, PJ_FWD, c_in);
+    SetComputedXY(c_out.xy.x, c_out.xy.y);
+    m_good = true;
+    return m_good;
+  }
+
+  bool IProj::SetCoordinate(const double x, const double y) {
+    SetXY(x, y);
+
+    PJ_COORD c_in;
+    c_in.xy.x = x;
+    c_in.xy.y = y;
+
+    PJ_COORD c_out = proj_trans(m_llaProj2outputProj, PJ_INV, c_in);
+
+    m_longitude = c_out.lpz.lam;
+    // Convert back to ocentric as ISIS defaults to operating in ocentric latitudes
+    m_latitude = ToPlanetocentric(c_out.lpz.phi);
+    m_good = true;
+    return m_good;
+  }
+
+    /**
+   * This method is used to determine the x/y range which completely covers the
+   * area of interest specified by the lat/lon range. The latitude/longitude
+   * range may be obtained from the labels. The purpose of this method is to
+   * return the x/y range so it can be used to compute how large a map may need
+   * to be. For example, how big a piece of paper is needed or how large of an
+   * image needs to be created. The method may fail as indicated by its return
+   * value. This currently mimics the sinusoidal projection's XYRange check
+   * and should be made more robust for use in proj. This will likely be a 
+   * method that walks the boundary of the projection
+   *
+   * @param minX Minimum x projection coordinate which covers the latitude
+   *             longitude range specified in the labels.
+   *
+   * @param maxX Maximum x projection coordinate which covers the latitude
+   *             longitude range specified in the labels.
+   *
+   * @param minY Minimum y projection coordinate which covers the latitude
+   *             longitude range specified in the labels.
+   *
+   * @param maxY Maximum y projection coordinate which covers the latitude
+   *             longitude range specified in the labels.
+   *
+   * @return bool Indicates whether the method was successful.
+   */
+  bool IProj::XYRange(double &minX, double &maxX,
+                      double &minY, double &maxY) {
+    // Check the corners of the lat/lon range
+    XYRangeCheck(m_minimumLatitude, m_minimumLongitude);
+    XYRangeCheck(m_maximumLatitude, m_minimumLongitude);
+    XYRangeCheck(m_minimumLatitude, m_maximumLongitude);
+    XYRangeCheck(m_maximumLatitude, m_maximumLongitude);
+
+    // If the latitude crosses the equator check there
+    if ((m_minimumLatitude < 0.0) && (m_maximumLatitude > 0.0)) {
+      XYRangeCheck(0.0, m_minimumLongitude);
+      XYRangeCheck(0.0, m_maximumLongitude);
+    }
+
+    // Make sure everything is ordered
+    if (m_minimumX >= m_maximumX) return false;
+    if (m_minimumY >= m_maximumY) return false;
+
+    // Return X/Y min/maxs
+    minX = m_minimumX;
+    maxX = m_maximumX;
+    minY = m_minimumY;
+    maxY = m_maximumY;
+    return true;
+  }
+
+  /**
+   * This is the function that is called in order to instantiate a
+   * Sinusoidal object.
+   *
+   * @param lab Cube labels with appropriate Mapping information.
+   *
+   * @param allowDefaults Indicates whether CenterLongitude are allowed to
+   *                      be computed using the middle of the longitude
+   *                      range specified in the labels.
+   *
+   * @return @b Isis::Projection* Pointer to a Sinusoidal projection object.
+   */
+  extern "C" Isis::TProjection *IProjPlugin(Isis::Pvl &lab,
+                                                 bool allowDefaults)
+  {
+    return new Isis::IProj(lab, allowDefaults);
+  }
+}

isis/src/base/objs/IProj/IProj.h

@@ -0,0 +1,60 @@
+#ifndef IProj_h
+#define IProj_h
+/** This is free and unencumbered software released into the public domain.
+The authors of ISIS do not claim copyright on the contents of this file.
+For more details about the LICENSE terms and the AUTHORS, you will
+find files of those names at the top level of this repository. **/
+
+/* SPDX-License-Identifier: CC0-1.0 */
+#include <proj.h>
+
+#include "TProjection.h"
+
+namespace Isis {
+  class Pvl;
+  class PvlGroup;
+  /**
+   * @brief Proj Map Projection
+   *
+   * This class provides methods for the forward and inverse projection for
+   * any ISIS map file through PROJ. We take a map file and convert it into
+   * a PROJ string, that string is then fed into the PROJ projection engine.
+   *
+   *
+   * Please see the TProjection/Projection class for a full accounting of all the methods
+   * available.
+   *
+   * @ingroup MapProjection
+   *
+   * @author 2023-10-17 Adam Paquette
+   *
+   * @internal
+   *   @history 2023-10-17 Adam Paquette - Built Class
+   */
+
+  class IProj : public TProjection {
+    public:
+      IProj(Pvl &label, bool allowDefaults = false);
+      ~IProj();
+
+      QString Name() const;
+      QString Version() const;
+
+      PvlGroup Mapping();
+
+      bool SetGround(const double lat, const double lon);
+      bool SetCoordinate(const double x, const double y);
+
+      bool XYRange(double &minX, double &maxX,
+                   double &minY, double &maxY);
+
+    private:
+      QString *m_userOutputProjStr;
+      PJ_CONTEXT *m_C;
+      PJ *m_llaProj;
+      PJ *m_outputProj;
+      PJ *m_llaProj2outputProj;
+  };
+};
+
+#endif

isis/src/base/objs/IProj/Projection.plugin

@@ -0,0 +1,5 @@
+Group = IProj
+  Library = IProj
+  Routine = IProjPlugin
+  Keyword = ProjStr
+EndGroup

isis/src/qisis/objs/AdvancedTrackTool/AdvancedTrackTool.cpp

@@ -581,7 +581,7 @@ namespace Isis {
                                  setText(QString::number(TProjection::To180Domain(lon), 'f', 15));
             p_tableWin->table()->item(row, getIndex("360 Positive West Longitude"))->
                                  setText(QString::number(wlon, 'f', 15));
-            p_tableWin->table()->item(row, getIndex("180 Positive East Longitude"))->
+            p_tableWin->table()->item(row, getIndex("180 Positive West Longitude"))->
                                  setText(QString::number(TProjection::To180Domain(wlon), 'f', 15));
             p_tableWin->table()->item(row, getIndex("Local Radius"))->setText(QString::number(radius, 'f', 15));
           }