|
| 1 | +/**************************************************************************** |
| 2 | + * |
| 3 | + * Copyright (c) 2023 Jaeyoung Lim. All rights reserved. |
| 4 | + * |
| 5 | + * Redistribution and use in source and binary forms, with or without |
| 6 | + * modification, are permitted provided that the following conditions |
| 7 | + * are met: |
| 8 | + * |
| 9 | + * 1. Redistributions of source code must retain the above copyright |
| 10 | + * notice, this list of conditions and the following disclaimer. |
| 11 | + * 2. Redistributions in binary form must reproduce the above copyright |
| 12 | + * notice, this list of conditions and the following disclaimer in |
| 13 | + * the documentation and/or other materials provided with the |
| 14 | + * distribution. |
| 15 | + * 3. Neither the name terrain-navigation nor the names of its contributors may be |
| 16 | + * used to endorse or promote products derived from this software |
| 17 | + * without specific prior written permission. |
| 18 | + * |
| 19 | + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
| 20 | + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
| 21 | + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS |
| 22 | + * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE |
| 23 | + * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, |
| 24 | + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, |
| 25 | + * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS |
| 26 | + * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED |
| 27 | + * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT |
| 28 | + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN |
| 29 | + * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
| 30 | + * POSSIBILITY OF SUCH DAMAGE. |
| 31 | + * |
| 32 | + ****************************************************************************/ |
| 33 | + |
| 34 | +#ifndef GEOCONVERSIONS_H |
| 35 | +#define GEOCONVERSIONS_H |
| 36 | + |
| 37 | +#include <math.h> |
| 38 | + |
| 39 | +class GeoConversions { |
| 40 | + public: |
| 41 | + GeoConversions(); |
| 42 | + virtual ~GeoConversions(); |
| 43 | + |
| 44 | + /** |
| 45 | + * @brief Convert WGS84 (LLA) to LV03/CH1903 |
| 46 | + * |
| 47 | + * @param lat latitude (degrees) WGS84 |
| 48 | + * @param lon lontitude (degrees) WGS84 |
| 49 | + * @param alt Altitude WGS84 |
| 50 | + * @param x |
| 51 | + * @param y |
| 52 | + * @param h |
| 53 | + */ |
| 54 | + static void forward(const double lat, const double lon, const double alt, double &y, double &x, double &h) { |
| 55 | + // 1. Convert the ellipsoidal latitudes φ and longitudes λ into arcseconds ["] |
| 56 | + const double lat_arc = lat * 3600.0; |
| 57 | + const double lon_arc = lon * 3600.0; |
| 58 | + |
| 59 | + // 2. Calculate the auxiliary values (differences of latitude and longitude relative to Bern in the unit [10000"]): |
| 60 | + // φ' = (φ – 169028.66 ")/10000 |
| 61 | + // λ' = (λ – 26782.5 ")/10000 |
| 62 | + const double lat_aux = (lat_arc - 169028.66) / 10000.0; |
| 63 | + const double lon_aux = (lon_arc - 26782.5) / 10000.0; |
| 64 | + |
| 65 | + // 3. Calculate projection coordinates in LV95 (E, N, h) or in LV03 (y, x, h) |
| 66 | + // E [m] = 2600072.37 + 211455.93 * λ' - 10938.51 * λ' * φ' - 0.36 * λ' * φ'2 - 44.54 * λ'3 |
| 67 | + // y [m] = E – 2000000.00 N [m] = 1200147.07 + 308807.95 * φ' + 3745.25 * λ'2 + 76.63 * φ'2 - 194.56 * λ'2 * φ' + |
| 68 | + // 119.79 * φ'3 x [m] = N – 1000000.00 |
| 69 | + // hCH [m] =hWGS – 49.55 + 2.73 * λ' + 6.94 * φ' |
| 70 | + const double E = 2600072.37 + 211455.93 * lon_aux - 10938.51 * lon_aux * lat_aux - |
| 71 | + 0.36 * lon_aux * std::pow(lat_aux, 2) - 44.54 * std::pow(lon_aux, 3); |
| 72 | + y = E - 2000000.00; |
| 73 | + const double N = 1200147.07 + 308807.95 * lat_aux + 3745.25 * std::pow(lon_aux, 2) + 76.63 * std::pow(lat_aux, 2) - |
| 74 | + 194.56 * std::pow(lon_aux, 2) * lat_aux + 119.79 * std::pow(lat_aux, 3); |
| 75 | + x = N - 1000000.00; |
| 76 | + |
| 77 | + h = alt - 49.55 + 2.73 * lon_aux + 6.84 * lat_aux; |
| 78 | + }; |
| 79 | + |
| 80 | + /** |
| 81 | + * @brief LV03/CH1903 to Convert WGS84 (LLA) |
| 82 | + * |
| 83 | + * @param x |
| 84 | + * @param y |
| 85 | + * @param h |
| 86 | + * @param lat latitude |
| 87 | + * @param lon longitude |
| 88 | + * @param alt altitude |
| 89 | + */ |
| 90 | + static void reverse(const double y, const double x, const double h, double &lat, double &lon, double &alt) { |
| 91 | + // 1. Convert the projection coordinates E (easting) and N (northing) in LV95 (or y / x in LV03) into the civilian |
| 92 | + // system (Bern = 0 / 0) and express in the unit [1000 km]: E' = (E – 2600000 m)/1000000 = (y – 600000 m)/1000000 |
| 93 | + // N' = (N – 1200000 m)/1000000 = (x – 200000 m)/1000000 |
| 94 | + const double y_aux = (y - 600000.0) / 1000000.0; |
| 95 | + const double x_aux = (x - 200000.0) / 1000000.0; |
| 96 | + |
| 97 | + // 2. Calculate longitude λ and latitude φ in the unit [10000"]: |
| 98 | + // λ' = 2.6779094 + 4.728982 * y' + 0.791484* y' * x' + 0.1306 * y' * x'2 - 0.0436 * y'3 |
| 99 | + // φ' = 16.9023892 + 3.238272 * x' - 0.270978 * y'2 - 0.002528 * x'2 - 0.0447 * y'2 * x' - 0.0140 * x'3 |
| 100 | + // hWGS [m] = hCH + 49.55 - 12.60 * y' - 22.64 * x' |
| 101 | + const double lon_aux = 2.6779094 + 4.728982 * y_aux + 0.791484 * y_aux * x_aux + |
| 102 | + 0.1306 * y_aux * std::pow(x_aux, 2) - 0.0436 * std::pow(y_aux, 3); |
| 103 | + const double lat_aux = 16.9023892 + 3.238272 * x_aux - 0.270978 * std::pow(y_aux, 2) - |
| 104 | + 0.002528 * std::pow(x_aux, 2) - 0.0447 * std::pow(y_aux, 2) * x_aux - |
| 105 | + 0.0140 * std::pow(x_aux, 3); |
| 106 | + alt = h + 49.55 - 12.60 * y_aux - 22.64 * x_aux; |
| 107 | + |
| 108 | + lon = lon_aux * 100.0 / 36.0; |
| 109 | + lat = lat_aux * 100.0 / 36.0; |
| 110 | + }; |
| 111 | + |
| 112 | + private: |
| 113 | +}; |
| 114 | + |
| 115 | +#endif |
0 commit comments