|
| 1 | +#include <math.h> |
| 2 | +#include <tuple> |
| 3 | + |
| 4 | +std::tuple<double, double, int> gpsToUTM(double lat, double lon) { |
| 5 | + // Constants for UTM conversion |
| 6 | + constexpr double a = 6378137.0; // WGS-84 major axis |
| 7 | + constexpr double f = 1 / 298.257223563; // WGS-84 flattening |
| 8 | + constexpr double k0 = 0.9996; // UTM scale factor |
| 9 | + constexpr double e = sqrt(f * (2 - f)); // Eccentricity |
| 10 | + constexpr double e2 = e * e; |
| 11 | + |
| 12 | + // UTM Zone |
| 13 | + int zone = static_cast<int>((lon + 180) / 6) + 1; |
| 14 | + |
| 15 | + // Convert latitude and longitude to radians |
| 16 | + double lat_rad = lat * M_PI / 180.0; |
| 17 | + double lon_rad = lon * M_PI / 180.0; |
| 18 | + |
| 19 | + // Central meridian of the UTM zone |
| 20 | + double lon_origin = (zone - 1) * 6 - 180 + 3; |
| 21 | + double lon_origin_rad = lon_origin * M_PI / 180.0; |
| 22 | + |
| 23 | + // Calculations for UTM coordinates |
| 24 | + double N = a / sqrt(1 - e2 * sin(lat_rad) * sin(lat_rad)); |
| 25 | + double T = tan(lat_rad) * tan(lat_rad); |
| 26 | + double C = e2 / (1 - e2) * cos(lat_rad) * cos(lat_rad); |
| 27 | + double A = cos(lat_rad) * (lon_rad - lon_origin_rad); |
| 28 | + |
| 29 | + double M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * lat_rad |
| 30 | + - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * sin(2 * lat_rad) |
| 31 | + + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * sin(4 * lat_rad) |
| 32 | + - (35 * e2 * e2 * e2 / 3072) * sin(6 * lat_rad)); |
| 33 | + |
| 34 | + double easting = k0 * N * (A + (1 - T + C) * A * A * A / 6 |
| 35 | + + (5 - 18 * T + T * T + 72 * C - 58 * e2) * A * A * A * A * A / 120) + 500000.0; |
| 36 | + |
| 37 | + double northing = k0 * (M + N * tan(lat_rad) * (A * A / 2 |
| 38 | + + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 |
| 39 | + + (61 - 58 * T + T * T + 600 * C - 330 * e2) * A * A * A * A * A * A / 720)); |
| 40 | + |
| 41 | + // Correct for southern hemisphere |
| 42 | + if (lat < 0) { |
| 43 | + northing += 10000000.0; |
| 44 | + } |
| 45 | + |
| 46 | + return std::make_tuple(easting, northing, zone); |
| 47 | +} |
| 48 | + |
| 49 | +std::tuple<double, double> gpsToGlobalCoord(double lat0, double lon0, double lat1, double lon1) { |
| 50 | + // Convert latitude and longitude to utm coordinates |
| 51 | + const auto& [e0, n0, zone0] = gpsToUTM(lat0, lon0); |
| 52 | + const auto& [e1, n1, zone1] = gpsToUTM(lat1, lon1); |
| 53 | + return {e1 - e0, n1 - n0}; |
| 54 | +} |
0 commit comments