Mercator.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include "Mercator.h"
 
#include <cmath>
#include <cfloat>
 
#include "IException.h"
#include "Constants.h"
#include "TProjection.h"
#include "Pvl.h"
#include "PvlGroup.h"
#include "PvlKeyword.h"
 
using namespace std;
namespace Isis {
  Mercator::Mercator(Pvl &label, bool allowDefaults) :
    TProjection::TProjection(label) {
    try {
      // Try to read the mapping group
      PvlGroup &mapGroup = label.findGroup("Mapping", Pvl::Traverse);
 
      // Compute and write the default center longitude if allowed and
      // necessary
      if ((allowDefaults) && (!mapGroup.hasKeyword("CenterLongitude"))) {
        double lon = (m_minimumLongitude + m_maximumLongitude) / 2.0;
        mapGroup += PvlKeyword("CenterLongitude", toString(lon));
      }
 
      // Compute and write the default center latitude if allowed and
      // necessary
      if ((allowDefaults) && (!mapGroup.hasKeyword("CenterLatitude"))) {
        double lat = (m_minimumLatitude + m_maximumLatitude) / 2.0;
        mapGroup += PvlKeyword("CenterLatitude", toString(lat));
      }
 
      // Get the center longitude  & latitude
      m_centerLongitude = mapGroup["CenterLongitude"];
      m_centerLatitude = mapGroup["CenterLatitude"];
      if (IsPlanetocentric()) {
        m_centerLatitude = ToPlanetographic(m_centerLatitude);
      }
 
      // convert to radians, adjust for longitude direction
      m_centerLongitude *= PI / 180.0;
      m_centerLatitude *= PI / 180.0;
      if (m_longitudeDirection == PositiveWest) m_centerLongitude *= -1.0;
 
      // Compute the scale factor
      double cos_clat = cos(m_centerLatitude);
      double sin_clat = sin(m_centerLatitude);
      double m_eccsq = Eccentricity() * Eccentricity();
      m_scalefactor = cos_clat / sqrt(1.0 - m_eccsq * sin_clat * sin_clat);
    }
    catch (IException &e) {
      QString message = "Invalid label group [Mapping]";
      throw IException(e, IException::Io, message, _FILEINFO_);
    }
  }
 
  Mercator::~Mercator() {
  }
 
  bool Mercator::operator== (const Projection &proj) {
    if (!Projection::operator==(proj)) return false;
    // don't do the below it is a recusive plunge
    //  if (Projection::operator!=(proj)) return false;
    Mercator *merc = (Mercator *) &proj;
    if ((merc->m_centerLongitude != m_centerLongitude) ||
        (merc->m_centerLatitude != m_centerLatitude)) return false;
    return true;
  }
 
  QString Mercator::Name() const {
    return "Mercator";
  }
 
  QString Mercator::Version() const {
    return "1.0";
  }
 
  double Mercator::TrueScaleLatitude() const {
    return m_centerLatitude * 180.0 / PI;
  }
 
  bool Mercator::IsEquatorialCylindrical() {
    return true;
  }
 
  bool Mercator::SetGround(const double lat, const double lon) {
    // Convert longitude to radians & clean up
    m_longitude = lon;
    double lonRadians = lon * PI / 180.0;
    if (m_longitudeDirection == PositiveWest) lonRadians *= -1.0;
 
    // Now convert latitude to radians & clean up ... it must be planetographic
    m_latitude = lat;
    double latRadians = lat;
    if (IsPlanetocentric()) latRadians = ToPlanetographic(latRadians);
    latRadians *= PI / 180.0;
 
    // Make sure latitude value is not too close to either pole
    if (fabs(fabs(m_latitude) - 90.0) <= DBL_EPSILON) {
      m_good = false;
      return m_good;
    }
 
    // Compute the coordinate
    double deltaLon = (lonRadians - m_centerLongitude);
    double x = m_equatorialRadius * deltaLon * m_scalefactor;
    double sinphi = sin(latRadians);
    double t = tCompute(latRadians, sinphi);
    double y = -m_equatorialRadius * m_scalefactor * log(t);
    SetComputedXY(x, y);
    m_good = true;
    return m_good;
  }
 
  bool Mercator::SetCoordinate(const double x, const double y) {
    // Save the coordinate
    SetXY(x, y);
 
    // Compute Snyder's t
    double snyders_t = exp(-GetY() / (m_equatorialRadius * m_scalefactor));
 
    // Compute latitude and make sure it is not above 90
    m_latitude = phi2Compute(snyders_t);
    if (fabs(m_latitude) > HALFPI) {
      if (fabs(HALFPI - fabs(m_latitude)) > DBL_EPSILON) {
        m_good = false;
        return m_good;
      }
      else if (m_latitude < 0.0) {
        m_latitude = -HALFPI;
      }
      else {
        m_latitude = HALFPI;
      }
    }
 
    // Compute longitude
    double coslat = cos(m_latitude);
    if (coslat <= DBL_EPSILON) {
      m_longitude = m_centerLongitude;
    }
    else {
      m_longitude = m_centerLongitude + GetX() /
                    (m_equatorialRadius * m_scalefactor);
    }
 
    // Convert to degrees
    m_latitude *= 180.0 / PI;
    m_longitude *= 180.0 / PI;
 
    // Cleanup the longitude
    if (m_longitudeDirection == PositiveWest) m_longitude *= -1.0;
    // These need to be done for circular type projections
    // m_longitude = To360Domain (m_longitude);
    // if (m_longitudeDomain == 180) m_longitude = To180Domain(m_longitude);
 
    // Cleanup the latitude
    if (IsPlanetocentric()) m_latitude = ToPlanetocentric(m_latitude);
 
    m_good = true;
    return m_good;
  }
 
  bool Mercator::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;
  }
 
  PvlGroup Mercator::Mapping() {
    PvlGroup mapping = TProjection::Mapping();
 
    mapping += m_mappingGrp["CenterLatitude"];
    mapping += m_mappingGrp["CenterLongitude"];
 
    return mapping;
  }
 
  PvlGroup Mercator::MappingLatitudes() {
    PvlGroup mapping = TProjection::MappingLatitudes();
 
    mapping += m_mappingGrp["CenterLatitude"];
 
    return mapping;
  }
 
  PvlGroup Mercator::MappingLongitudes() {
    PvlGroup mapping = TProjection::MappingLongitudes();
 
    mapping += m_mappingGrp["CenterLongitude"];
 
    return mapping;
  }
 
} // end namespace isis
 
extern "C" Isis::Projection *MercatorPlugin(Isis::Pvl &lab,
    bool allowDefaults) {
  return new Isis::Mercator(lab, allowDefaults);
}