UTMPoint.cc

Go to the documentation of this file.

#include <utl/UTMPoint.h>
#include <utl/Point.h>
#include <utl/CoordinateSystem.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>

#include <boost/tuple/tuple.hpp>
#include <boost/tuple/tuple_io.hpp>

#include <sstream>
#include <cmath>

using namespace utl;


// some private constants for the implementation
namespace {
  // utm zones
  const double utmZoneWidth  = 6*deg;
  const double utmZoneOffset = 180*deg; // added to lattitude to move
                                        // angles from range
                                        // -180..180 to 0..360.

  const char utmBandLetter[] = "CDEFGHJKLMNPQRSTUVWX";
  const int numberOfUTMBands = sizeof(utmBandLetter) / sizeof(char) - 1;
  const double widthOfUTMBand = 8*deg;
  [[maybe_unused]] const double southernmostUTMBand = -80*deg;
  const char firstNorthernBand = 'N';  // 11th zone

  // utm definitions from transformation
  double k0 = 0.9996;
}


Point
UTMPoint::GetPoint(const CoordinateSystemPtr& theCS)
  const
{
  double latitude;
  double longitude;
  double height;
  boost::tie(latitude, longitude, height) = GetGeodeticCoordinates();
  const Point pt =
    fEllipsoid->LatitudeLongitudeHeightToPoint(latitude, longitude, height);

  if (theCS)
    pt.TransformTo(theCS);

  return pt;
}


boost::tuple<double, double, double>
UTMPoint::GetGeodeticCoordinates()
  const
{
  //converts UTM coords to lat/long.  Equations from USGS Bulletin
  //1532 East Longitudes are positive, West longitudes are negative.
  //North latitudes are positive, South latitudes are negative. Lat and
  //Long are in Auger units.

  const double a = fEllipsoid->GetEquatorialRadius();
  const double eccSquared = fEllipsoid->GetEccentricity2();

  const double e1 = (1 - sqrt(1-eccSquared)) / (1 + sqrt(1-eccSquared));

  const double x = fEasting - 500000;  //remove 500,000 meter offset for longitude
  double y = fNorthing;

  const bool inNorthernHemisphere = fBand >= firstNorthernBand;
  if (!inNorthernHemisphere)
    y -= 10000000;  // remove 10,000,000 meter offset used
                    // for southern hemisphere

  const double LongitudeOrigin = (fZone - 0.5) * utmZoneWidth - utmZoneOffset;

  const double eccPrimeSquared = eccSquared / (1-eccSquared);
  const double e4 = eccSquared*eccSquared;
  const double e6 = e4*eccSquared;
  const double e8 = e4*e4;
  const double e1_2 = e1*e1;

  const double M = y / k0;
  const double mu = M / (a * (1 - eccSquared/4 - 3*e4/64 - 5*e6/256 - 175*e8/16384));
  const double mur = mu/rad;
  const double phi1 = mur
    + (3*e1/2 - 27*e1_2*e1/32)*sin(2*mur)
    + (21*e1_2/16 - 55*e1_2*e1_2/32)*sin(4*mur)
    + (151*e1*e1_2/96)*sin(6*mur)
    + (1097*e1_2*e1_2/512)*sin(8*mur);

  const double phi1r = phi1/rad;
  const double cosPhi1 = cos(phi1r);
  const double sinPhi1 = sin(phi1r);
  const double tanPhi1 = sinPhi1/cosPhi1;

  const double N1 = a / sqrt(1 - eccSquared * sinPhi1*sinPhi1);
  const double T1 = tanPhi1*tanPhi1;
  const double C1 = eccPrimeSquared * cosPhi1*cosPhi1;
  const double R1 = a * (1-eccSquared) / pow(1 - eccSquared * sinPhi1*sinPhi1, 1.5);
  const double D = x / (N1*k0);
  const double D2 = D*D;
  const double D4 = D2*D2;
  const double C1_2 = C1*C1;
  const double T1_2 = T1*T1;

  const double latitude = phi1
    - N1*tanPhi1/R1 * (
        D2/2
        - (5 + 3*T1 + 10*C1 - 4*C1_2 - 9*eccPrimeSquared) * D4 / 24
        + (61 + 90*T1 + 298*C1 + 45*T1_2 - 252*eccPrimeSquared - 3*C1_2) * D2*D4 / 720
    );

  const double longitude =
    (D
     - (1 + 2*T1 + C1) * D*D2 / 6
     + (5 - 2*C1 + 28*T1 - 3*C1_2 + 8*eccPrimeSquared + 24*T1_2) * D*D4 / 120
    ) / cosPhi1 + LongitudeOrigin;

  return boost::make_tuple(latitude, longitude, fHeight);
}


/***************************************************************************
 * Private, internal functions
 ***************************************************************************/

void
UTMPoint::SetFromOtherDatum(const UTMPoint& theUTMPoint)
{
  if (*fEllipsoid != *theUTMPoint.fEllipsoid) {
    std::ostringstream msg;
    msg << "DATUM conversions not implemented in UTMPoint constructor";
    FATAL(msg);
    throw UTMException(msg.str());
  }

  fNorthing = theUTMPoint.fNorthing;
  fEasting  = theUTMPoint.fEasting;
  fHeight   = theUTMPoint.fHeight;
  fZone     = theUTMPoint.fZone;
  fBand     = theUTMPoint.fBand;
}


void
UTMPoint::SetFromCartesian(const Point& thePoint)
{
  double latitude;
  double longitude;
  double height;

  boost::tie(latitude, longitude, height) =
    fEllipsoid->PointToLatitudeLongitudeHeight(thePoint);

  fHeight = height;
  SetFromLatitudeLongitude(latitude, longitude);
}


void
UTMPoint::SetFromLatitudeLongitude(const double theLatitude,
                                   const double theLongitude)
{
  fZone = ZoneFromLatitudeLongitude(theLatitude, theLongitude);
  fBand = BandFromLatitude(theLatitude);

  const double LongitudeOrigin = (fZone-0.5)*utmZoneWidth - utmZoneOffset;

  //converts lat/long to UTM coords.  Equations from USGS Bulletin
  //1532 East Longitudes are positive, West longitudes are negative.
  //North latitudes are positive, South latitudes are negative. Lat and
  //Long are in Auger Units

  const double a = fEllipsoid->GetEquatorialRadius();
  const double eccSquared = fEllipsoid->GetEccentricity2();

  const double lat = theLatitude / rad;
  const double cosLatitude = cos(lat);
  const double sinLatitude = sin(lat);
  const double tanLatitude = sinLatitude / cosLatitude;

  const double eccPrimeSquared = eccSquared / (1 - eccSquared);

  const double N = a / sqrt(1 - eccSquared * sinLatitude*sinLatitude);
  const double T = tanLatitude*tanLatitude;
  const double C = eccPrimeSquared * cosLatitude*cosLatitude;
  const double A = cosLatitude * (theLongitude/rad - LongitudeOrigin/rad);

  const double e4 = eccSquared*eccSquared;
  const double e6 = e4 * eccSquared;
  const double e8 = e4*e4;

  const double c1 =
    1
    - eccSquared / 4
    - 3 * e4 / 64
    - 5 * e6 / 256
    - 175 * e8 / 16384;
  const double c2 =
    - 3 * eccSquared / 8
    - 3 * e4 / 32
    - 45 * e6 / 1024
    - 105 * e8 / 4096;
  const double c3 =
    15 * e4 / 256
    + 45 * e6 / 1024
    + 525 * e8 / 16384;
  const double c4 =
    - 35 * e6 / 3072
    - 175 * e8 / 12288;
  const double c5 = 315 * e8 / 131072;

  const double M = a *
    (c1 * lat
     + c2 * sin(2*lat)
     + c3 * sin(4*lat)
     + c4 * sin(6*lat)
     + c5 * sin(8*lat)
    );

  const double A2 = A*A;
  const double A3 = A2*A;
  const double T2 = T*T;

  fEasting = k0 * N * (A + (1-T+C) * A3/6
                       + (5 - 18*T + T2 + 72*C - 58*eccPrimeSquared) * A2*A3 / 120)
             + 500000;

  fNorthing = k0 * (M + N * tanLatitude *
                   (A2/2 + (5 - T + 9*C + 4*C*C) * A*A3 / 24
                    + (61 - 58*T + T2 + 600*C - 330*eccPrimeSquared) * A3*A3 / 720));

  if (theLatitude < 0)
    fNorthing += 10000000; //10 000 000 meter offset for southern hemisphere
}


void
UTMPoint::VerifyZone()
{
  if (fZone < eZoneMin || fZone > eZoneMax) {
    std::ostringstream msg;
    msg << "Zone out of bounds: " << fZone
        << " not in " << eZoneMin << ".." << eZoneMax << ".";
    FATAL(msg);
    throw UTMZoneException(msg.str());
  }
}


/***************************************************************************
 * static, private, internal functions
 ***************************************************************************/

int
UTMPoint::ZoneFromLatitudeLongitude(const double theLatitude,
                                    const double theLongitude)
{
  using std::floor;
  int zone = int(floor(theLongitude / utmZoneWidth)) + eZoneMax/2 + 1;

  // southern Norway zone exception
  if (56*deg <= theLatitude && theLatitude < 64*deg &&
      3*deg <= theLongitude && theLatitude < 12*deg) {
    zone = 32;
  }
  else if (72*deg <= theLatitude && theLatitude < 84*deg) {
    // Svalbard (Norway) zone exceptions
    if (0 <= theLongitude && theLongitude < 9*deg)
      zone = 31;
    else if (9*deg <= theLongitude && theLongitude < 21*deg)
      zone = 33;
    else if (21*deg <= theLongitude && theLongitude < 33*deg)
      zone = 35;
    else if (33*deg <= theLongitude && theLongitude < 42*deg)
      zone = 37;
  }

  return zone;
}


char
UTMPoint::BandFromLatitude(const double theLatitude)
{
  using std::floor;
  int index = int(floor(theLatitude / widthOfUTMBand)) + numberOfUTMBands/2;
  if (0 <= index && index < numberOfUTMBands)
    return utmBandLetter[index];
  else if (72*deg <= theLatitude && theLatitude < 84*deg)
    return 'X';
  else {
    std::ostringstream msg;
    msg << "Latitude " << theLatitude/deg
        << " deg out of range [-80 .. 84]. "
           "Cannot calculate band letter." << std::endl;
    ERROR(msg);
    throw UTMZoneException(msg.str());
    return ' '; // to keep the compiler happy
  }
}


/***************************************************************************
 * friends and related functions
 ***************************************************************************/

std::ostream&
utl::operator<<(std::ostream& s, const UTMPoint& p)
{
  s << p.GetCoordinates();
  return s;
}