|
| 1 | +import numpy as np |
| 2 | + |
| 3 | +""" |
| 4 | +Coordinate transformation module. All methods accept arrays as input |
| 5 | +with each row as a position. |
| 6 | +""" |
| 7 | + |
| 8 | +a = 6378137 |
| 9 | +b = 6356752.3142 |
| 10 | +esq = 6.69437999014 * 0.001 |
| 11 | +e1sq = 6.73949674228 * 0.001 |
| 12 | + |
| 13 | + |
| 14 | +def geodetic2ecef(geodetic, radians=False): |
| 15 | + geodetic = np.array(geodetic) |
| 16 | + input_shape = geodetic.shape |
| 17 | + geodetic = np.atleast_2d(geodetic) |
| 18 | + |
| 19 | + ratio = 1.0 if radians else (np.pi / 180.0) |
| 20 | + lat = ratio * geodetic[:, 0] |
| 21 | + lon = ratio * geodetic[:, 1] |
| 22 | + alt = geodetic[:, 2] |
| 23 | + |
| 24 | + xi = np.sqrt(1 - esq * np.sin(lat) ** 2) |
| 25 | + x = (a / xi + alt) * np.cos(lat) * np.cos(lon) |
| 26 | + y = (a / xi + alt) * np.cos(lat) * np.sin(lon) |
| 27 | + z = (a / xi * (1 - esq) + alt) * np.sin(lat) |
| 28 | + ecef = np.array([x, y, z]).T |
| 29 | + return ecef.reshape(input_shape) |
| 30 | + |
| 31 | + |
| 32 | +def ecef2geodetic(ecef, radians=False): |
| 33 | + """ |
| 34 | + Convert ECEF coordinates to geodetic using ferrari's method |
| 35 | + """ |
| 36 | + # Save shape and export column |
| 37 | + ecef = np.atleast_1d(ecef) |
| 38 | + input_shape = ecef.shape |
| 39 | + ecef = np.atleast_2d(ecef) |
| 40 | + x, y, z = ecef[:, 0], ecef[:, 1], ecef[:, 2] |
| 41 | + |
| 42 | + ratio = 1.0 if radians else (180.0 / np.pi) |
| 43 | + |
| 44 | + # Conver from ECEF to geodetic using Ferrari's methods |
| 45 | + # https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#Ferrari.27s_solution |
| 46 | + r = np.sqrt(x * x + y * y) |
| 47 | + Esq = a * a - b * b |
| 48 | + F = 54 * b * b * z * z |
| 49 | + G = r * r + (1 - esq) * z * z - esq * Esq |
| 50 | + C = (esq * esq * F * r * r) / (pow(G, 3)) |
| 51 | + S = np.cbrt(1 + C + np.sqrt(C * C + 2 * C)) |
| 52 | + P = F / (3 * pow((S + 1 / S + 1), 2) * G * G) |
| 53 | + Q = np.sqrt(1 + 2 * esq * esq * P) |
| 54 | + r_0 = -(P * esq * r) / (1 + Q) + np.sqrt(0.5 * a * a * (1 + 1.0 / Q) - |
| 55 | + P * (1 - esq) * z * z / (Q * (1 + Q)) - 0.5 * P * r * r) |
| 56 | + U = np.sqrt(pow((r - esq * r_0), 2) + z * z) |
| 57 | + V = np.sqrt(pow((r - esq * r_0), 2) + (1 - esq) * z * z) |
| 58 | + Z_0 = b * b * z / (a * V) |
| 59 | + h = U * (1 - b * b / (a * V)) |
| 60 | + lat = ratio * np.arctan((z + e1sq * Z_0) / r) |
| 61 | + lon = ratio * np.arctan2(y, x) |
| 62 | + |
| 63 | + # stack the new columns and return to the original shape |
| 64 | + geodetic = np.column_stack((lat, lon, h)) |
| 65 | + return geodetic.reshape(input_shape) |
| 66 | + |
| 67 | + |
| 68 | +class LocalCoord: |
| 69 | + """ |
| 70 | + Allows conversions to local frames. In this case NED. |
| 71 | + That is: North East Down from the start position in |
| 72 | + meters. |
| 73 | + """ |
| 74 | + |
| 75 | + def __init__(self, init_geodetic, init_ecef): |
| 76 | + self.init_ecef = init_ecef |
| 77 | + lat, lon, _ = (np.pi / 180) * np.array(init_geodetic) |
| 78 | + self.ned2ecef_matrix = np.array([[-np.sin(lat) * np.cos(lon), -np.sin(lon), -np.cos(lat) * np.cos(lon)], |
| 79 | + [-np.sin(lat) * np.sin(lon), np.cos(lon), -np.cos(lat) * np.sin(lon)], |
| 80 | + [np.cos(lat), 0, -np.sin(lat)]]) |
| 81 | + self.ecef2ned_matrix = self.ned2ecef_matrix.T |
| 82 | + |
| 83 | + @classmethod |
| 84 | + def from_geodetic(cls, init_geodetic): |
| 85 | + init_ecef = geodetic2ecef(init_geodetic) |
| 86 | + return LocalCoord(init_geodetic, init_ecef) |
| 87 | + |
| 88 | + @classmethod |
| 89 | + def from_ecef(cls, init_ecef): |
| 90 | + init_geodetic = ecef2geodetic(init_ecef) |
| 91 | + return LocalCoord(init_geodetic, init_ecef) |
| 92 | + |
| 93 | + def ecef2ned(self, ecef): |
| 94 | + ecef = np.array(ecef) |
| 95 | + return np.dot(self.ecef2ned_matrix, (ecef - self.init_ecef).T).T |
| 96 | + |
| 97 | + def ned2ecef(self, ned): |
| 98 | + ned = np.array(ned) |
| 99 | + # Transpose so that init_ecef will broadcast correctly for 1d or 2d ned. |
| 100 | + return np.dot(self.ned2ecef_matrix, ned.T).T + self.init_ecef |
| 101 | + |
| 102 | + def geodetic2ned(self, geodetic): |
| 103 | + ecef = geodetic2ecef(geodetic) |
| 104 | + return self.ecef2ned(ecef) |
| 105 | + |
| 106 | + def ned2geodetic(self, ned): |
| 107 | + ecef = self.ned2ecef(ned) |
| 108 | + return ecef2geodetic(ecef) |
0 commit comments