Merge branch 'develop' of http://gitlab.mynt.com/mynteye/mynt-eye-s-sdk into develop

This commit is contained in:
Osenberg 2019-03-27 17:23:35 +08:00
commit 719c0b32ae

View File

@ -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<int>((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<double>(
(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<int>(ZoneLetter) - static_cast<int>('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;
}
}