IProj.cpp

 
/* 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();
    proj_log_level(m_C, PJ_LOG_ERROR);
 
    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 (!m_outputProj) {
      QString msg = "Unable to create projection from [" + QString(projString.c_str()) + "]";
      throw IException(IException::User, msg, _FILEINFO_);
    }
 
    // Read the radii information from the PROJ projection
    PJ *ellipsoid = proj_get_ellipsoid(m_C, m_outputProj);
 
    if (ellipsoid == nullptr) {
      QString msg = "Unable to create ellipsoid from [" + *m_userOutputProjStr + "]";
      throw IException(IException::User, msg, _FILEINFO_);
    }
 
    int res = proj_ellipsoid_get_parameters(m_C, ellipsoid, 
                                            &m_equatorialRadius,
                                            &m_polarRadius,
                                            nullptr,
                                            nullptr);
 
    proj_destroy(ellipsoid);
 
    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 lla projection from [" + QString(projString.c_str()) + "]";
      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_);
    }
  }

  IProj::~IProj() {
    delete m_userOutputProjStr;
    proj_destroy(m_llaProj);
    proj_destroy(m_outputProj);
    proj_destroy(m_llaProj2outputProj);
    proj_context_destroy(m_C);
  }

  QString IProj::Name() const {
    return "Proj";
  }

  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 += PvlKeyword("ProjStr", *m_userOutputProjStr);
 
    return mapping;
  }

  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;
    c_in.lpz.phi = 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;
    m_latitude = c_out.lpz.phi;
    m_good = true;
    return m_good;
  }

  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);
 
    // 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;
  }

  extern "C" Isis::TProjection *IProjPlugin(Isis::Pvl &lab,
                                                 bool allowDefaults)
  {
    return new Isis::IProj(lab, allowDefaults);
  }
}