diff --git a/src/mynteye/api/camodocal/src/gpl/gpl.cc b/src/mynteye/api/camodocal/src/gpl/gpl.cc index bfeacb7..591b00d 100644 --- a/src/mynteye/api/camodocal/src/gpl/gpl.cc +++ b/src/mynteye/api/camodocal/src/gpl/gpl.cc @@ -521,168 +521,8 @@ char UTMLetterDesignator(double latitude) { else letterDesignator = 'Z'; // This is here as an error flag to show that the // Latitude is outside the UTM limits - return letterDesignator; } -void LLtoUTM( - double latitude, double longitude, double &utmNorthing, double &utmEasting, // NOLINT - std::string &utmZone) { // NOLINT - // 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 decimal degrees - // Written by Chuck Gantz- chuck.gantz@globalstar.com +} // namespace camodocal - double k0 = 0.9996; - - double LongOrigin; - double eccPrimeSquared; - double N, T, C, A, M; - - double LatRad = latitude * M_PI / 180.0; - double LongRad = longitude * M_PI / 180.0; - double LongOriginRad; - - int ZoneNumber = static_cast((longitude + 180.0) / 6.0) + 1; - - if (latitude >= 56.0 && latitude < 64.0 && longitude >= 3.0 && - longitude < 12.0) { - ZoneNumber = 32; - } - - // Special zones for Svalbard - if (latitude >= 72.0 && latitude < 84.0) { - if (longitude >= 0.0 && longitude < 9.0) - ZoneNumber = 31; - else if (longitude >= 9.0 && longitude < 21.0) - ZoneNumber = 33; - else if (longitude >= 21.0 && longitude < 33.0) - ZoneNumber = 35; - else if (longitude >= 33.0 && longitude < 42.0) - ZoneNumber = 37; - } - LongOrigin = static_cast( - (ZoneNumber - 1) * 6 - 180 + 3); // +3 puts origin in middle of zone - LongOriginRad = LongOrigin * M_PI / 180.0; - - // compute the UTM Zone from the latitude and longitude - std::ostringstream oss; - oss << ZoneNumber << UTMLetterDesignator(latitude); - utmZone = oss.str(); - - eccPrimeSquared = WGS84_ECCSQ / (1.0 - WGS84_ECCSQ); - - N = WGS84_A / sqrt(1.0 - WGS84_ECCSQ * sin(LatRad) * sin(LatRad)); - T = tan(LatRad) * tan(LatRad); - C = eccPrimeSquared * cos(LatRad) * cos(LatRad); - A = cos(LatRad) * (LongRad - LongOriginRad); - - M = WGS84_A * - ((1.0 - WGS84_ECCSQ / 4.0 - 3.0 * WGS84_ECCSQ * WGS84_ECCSQ / 64.0 - - 5.0 * WGS84_ECCSQ * WGS84_ECCSQ * WGS84_ECCSQ / 256.0) * - LatRad - - (3.0 * WGS84_ECCSQ / 8.0 + 3.0 * WGS84_ECCSQ * WGS84_ECCSQ / 32.0 + - 45.0 * WGS84_ECCSQ * WGS84_ECCSQ * WGS84_ECCSQ / 1024.0) * - sin(2.0 * LatRad) + - (15.0 * WGS84_ECCSQ * WGS84_ECCSQ / 256.0 + - 45.0 * WGS84_ECCSQ * WGS84_ECCSQ * WGS84_ECCSQ / 1024.0) * - sin(4.0 * LatRad) - - (35.0 * WGS84_ECCSQ * WGS84_ECCSQ * WGS84_ECCSQ / 3072.0) * - sin(6.0 * LatRad)); - - utmEasting = - k0 * N * (A + (1.0 - T + C) * A * A * A / 6.0 + - (5.0 - 18.0 * T + T * T + 72.0 * C - 58.0 * eccPrimeSquared) * - A * A * A * A * A / 120.0) + - 500000.0; - - utmNorthing = - k0 * - (M + - N * tan(LatRad) * - (A * A / 2.0 + - (5.0 - T + 9.0 * C + 4.0 * C * C) * A * A * A * A / 24.0 + - (61.0 - 58.0 * T + T * T + 600.0 * C - 330.0 * eccPrimeSquared) * - A * A * A * A * A * A / 720.0)); - if (latitude < 0.0) { - utmNorthing += 10000000.0; // 10000000 meter offset for southern hemisphere - } -} - -void UTMtoLL( - double utmNorthing, double utmEasting, const std::string &utmZone, - double &latitude, double &longitude) { // NOLINT - // 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 decimal degrees. - // Written by Chuck Gantz- chuck.gantz@globalstar.com - - double k0 = 0.9996; - double eccPrimeSquared; - double e1 = (1.0 - sqrt(1.0 - WGS84_ECCSQ)) / (1.0 + sqrt(1.0 - WGS84_ECCSQ)); - double N1, T1, C1, R1, D, M; - double LongOrigin; - double mu, phi1, phi1Rad; - double x, y; - int ZoneNumber; - char ZoneLetter; - bool NorthernHemisphere; - - x = utmEasting - 500000.0; // remove 500,000 meter offset for longitude - y = utmNorthing; - - std::istringstream iss(utmZone); - iss >> ZoneNumber >> ZoneLetter; - if ((static_cast(ZoneLetter) - static_cast('N')) >= 0) { - NorthernHemisphere = true; // point is in northern hemisphere - } else { - NorthernHemisphere = false; // point is in southern hemisphere - y -= 10000000.0; // remove 10,000,000 meter offset used for southern - // hemisphere - } - - LongOrigin = (ZoneNumber - 1.0) * 6.0 - 180.0 + - 3.0; // +3 puts origin in middle of zone - - eccPrimeSquared = WGS84_ECCSQ / (1.0 - WGS84_ECCSQ); - - M = y / k0; - mu = M / (WGS84_A * - (1.0 - WGS84_ECCSQ / 4.0 - 3.0 * WGS84_ECCSQ * WGS84_ECCSQ / 64.0 - - 5.0 * WGS84_ECCSQ * WGS84_ECCSQ * WGS84_ECCSQ / 256.0)); - - phi1Rad = mu + (3.0 * e1 / 2.0 - 27.0 * e1 * e1 * e1 / 32.0) * sin(2.0 * mu) + - (21.0 * e1 * e1 / 16.0 - 55.0 * e1 * e1 * e1 * e1 / 32.0) * - sin(4.0 * mu) + - (151.0 * e1 * e1 * e1 / 96.0) * sin(6.0 * mu); - phi1 = phi1Rad / M_PI * 180.0; - - N1 = WGS84_A / sqrt(1.0 - WGS84_ECCSQ * sin(phi1Rad) * sin(phi1Rad)); - T1 = tan(phi1Rad) * tan(phi1Rad); - C1 = eccPrimeSquared * cos(phi1Rad) * cos(phi1Rad); - R1 = WGS84_A * (1.0 - WGS84_ECCSQ) / - pow(1.0 - WGS84_ECCSQ * sin(phi1Rad) * sin(phi1Rad), 1.5); - D = x / (N1 * k0); - - latitude = phi1Rad - - (N1 * tan(phi1Rad) / R1) * - (D * D / 2.0 - - (5.0 + 3.0 * T1 + 10.0 * C1 - 4.0 * C1 * C1 - - 9.0 * eccPrimeSquared) * - D * D * D * D / 24.0 + - (61.0 + 90.0 * T1 + 298.0 * C1 + 45.0 * T1 * T1 - - 252.0 * eccPrimeSquared - 3.0 * C1 * C1) * - D * D * D * D * D * D / 720.0); - latitude *= 180.0 / M_PI; - - longitude = (D - (1.0 + 2.0 * T1 + C1) * D * D * D / 6.0 + - (5.0 - 2.0 * C1 + 28.0 * T1 - 3.0 * C1 * C1 + - 8.0 * eccPrimeSquared + 24.0 * T1 * T1) * - D * D * D * D * D / 120.0) / - cos(phi1Rad); - longitude = LongOrigin + longitude / M_PI * 180.0; -} - -}