00001 /** 00002 * \file LocalCartesian.hpp 00003 * \brief Header for GeographicLib::LocalCartesian class 00004 * 00005 * Copyright (c) Charles Karney (2008-2011) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_LOCALCARTESIAN_HPP) 00011 #define GEOGRAPHICLIB_LOCALCARTESIAN_HPP 1 00012 00013 #include <GeographicLib/Geocentric.hpp> 00014 #include <GeographicLib/Constants.hpp> 00015 00016 namespace GeographicLib { 00017 00018 /** 00019 * \brief Local cartesian coordinates 00020 * 00021 * Convert between geodetic coordinates latitude = \e lat, longitude = \e 00022 * lon, height = \e h (measured vertically from the surface of the ellipsoid) 00023 * to local cartesian coordinates (\e x, \e y, \e z). The origin of local 00024 * cartesian coordinate system is at \e lat = \e lat0, \e lon = \e lon0, \e h 00025 * = \e h0. The \e z axis is normal to the ellipsoid; the \e y axis points 00026 * due north. The plane \e z = - \e h0 is tangent to the ellipsoid. 00027 * 00028 * The conversions all take place via geocentric coordinates using a 00029 * Geocentric object (by default Geocentric::WGS84()). 00030 * 00031 * Example of use: 00032 * \include example-LocalCartesian.cpp 00033 * 00034 * <a href="CartConvert.1.html">CartConvert</a> is a command-line utility 00035 * providing access to the functionality of Geocentric and LocalCartesian. 00036 **********************************************************************/ 00037 00038 class GEOGRAPHICLIB_EXPORT LocalCartesian { 00039 private: 00040 typedef Math::real real; 00041 static const size_t dim_ = 3; 00042 static const size_t dim2_ = dim_ * dim_; 00043 Geocentric _earth; 00044 real _lat0, _lon0, _h0; 00045 real _x0, _y0, _z0, _r[dim2_]; 00046 void IntForward(real lat, real lon, real h, real& x, real& y, real& z, 00047 real M[dim2_]) const; 00048 void IntReverse(real x, real y, real z, real& lat, real& lon, real& h, 00049 real M[dim2_]) const; 00050 void MatrixMultiply(real M[dim2_]) const; 00051 public: 00052 00053 /** 00054 * Constructor setting the origin. 00055 * 00056 * @param[in] lat0 latitude at origin (degrees). 00057 * @param[in] lon0 longitude at origin (degrees). 00058 * @param[in] h0 height above ellipsoid at origin (meters); default 0. 00059 * @param[in] earth Geocentric object for the transformation; default 00060 * Geocentric::WGS84(). 00061 * 00062 * \e lat0 should be in the range [−90°, 90°]; \e 00063 * lon0 should be in the range [−540°, 540°). 00064 **********************************************************************/ 00065 LocalCartesian(real lat0, real lon0, real h0 = 0, 00066 const Geocentric& earth = Geocentric::WGS84()) 00067 : _earth(earth) 00068 { Reset(lat0, lon0, h0); } 00069 00070 /** 00071 * Default constructor. 00072 * 00073 * @param[in] earth Geocentric object for the transformation; default 00074 * Geocentric::WGS84(). 00075 * 00076 * Sets \e lat0 = 0, \e lon0 = 0, \e h0 = 0. 00077 **********************************************************************/ 00078 explicit LocalCartesian(const Geocentric& earth = Geocentric::WGS84()) 00079 : _earth(earth) 00080 { Reset(real(0), real(0), real(0)); } 00081 00082 /** 00083 * Reset the origin. 00084 * 00085 * @param[in] lat0 latitude at origin (degrees). 00086 * @param[in] lon0 longitude at origin (degrees). 00087 * @param[in] h0 height above ellipsoid at origin (meters); default 0. 00088 * 00089 * \e lat0 should be in the range [−90°, 90°]; \e 00090 * lon0 should be in the range [−540°, 540°). 00091 **********************************************************************/ 00092 void Reset(real lat0, real lon0, real h0 = 0); 00093 00094 /** 00095 * Convert from geodetic to local cartesian coordinates. 00096 * 00097 * @param[in] lat latitude of point (degrees). 00098 * @param[in] lon longitude of point (degrees). 00099 * @param[in] h height of point above the ellipsoid (meters). 00100 * @param[out] x local cartesian coordinate (meters). 00101 * @param[out] y local cartesian coordinate (meters). 00102 * @param[out] z local cartesian coordinate (meters). 00103 * 00104 * \e lat should be in the range [−90°, 90°]; \e lon 00105 * should be in the range [−540°, 540°). 00106 **********************************************************************/ 00107 void Forward(real lat, real lon, real h, real& x, real& y, real& z) 00108 const { 00109 IntForward(lat, lon, h, x, y, z, NULL); 00110 } 00111 00112 /** 00113 * Convert from geodetic to local cartesian coordinates and return rotation 00114 * matrix. 00115 * 00116 * @param[in] lat latitude of point (degrees). 00117 * @param[in] lon longitude of point (degrees). 00118 * @param[in] h height of point above the ellipsoid (meters). 00119 * @param[out] x local cartesian coordinate (meters). 00120 * @param[out] y local cartesian coordinate (meters). 00121 * @param[out] z local cartesian coordinate (meters). 00122 * @param[out] M if the length of the vector is 9, fill with the rotation 00123 * matrix in row-major order. 00124 * 00125 * \e lat should be in the range [−90°, 90°]; \e lon 00126 * should be in the range [−540°, 540°). 00127 * 00128 * Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can 00129 * express \e v as \e column vectors in one of two ways 00130 * - in east, north, up coordinates (where the components are relative to a 00131 * local coordinate system at (\e lat, \e lon, \e h)); call this 00132 * representation \e v1. 00133 * - in \e x, \e y, \e z coordinates (where the components are relative to 00134 * the local coordinate system at (\e lat0, \e lon0, \e h0)); call this 00135 * representation \e v0. 00136 * . 00137 * Then we have \e v0 = \e M ⋅ \e v1. 00138 **********************************************************************/ 00139 void Forward(real lat, real lon, real h, real& x, real& y, real& z, 00140 std::vector<real>& M) 00141 const { 00142 if (M.end() == M.begin() + dim2_) { 00143 real t[dim2_]; 00144 IntForward(lat, lon, h, x, y, z, t); 00145 std::copy(t, t + dim2_, M.begin()); 00146 } else 00147 IntForward(lat, lon, h, x, y, z, NULL); 00148 } 00149 00150 /** 00151 * Convert from local cartesian to geodetic coordinates. 00152 * 00153 * @param[in] x local cartesian coordinate (meters). 00154 * @param[in] y local cartesian coordinate (meters). 00155 * @param[in] z local cartesian coordinate (meters). 00156 * @param[out] lat latitude of point (degrees). 00157 * @param[out] lon longitude of point (degrees). 00158 * @param[out] h height of point above the ellipsoid (meters). 00159 * 00160 * The value of \e lon returned is in the range [−180°, 00161 * 180°). 00162 **********************************************************************/ 00163 void Reverse(real x, real y, real z, real& lat, real& lon, real& h) 00164 const { 00165 IntReverse(x, y, z, lat, lon, h, NULL); 00166 } 00167 00168 /** 00169 * Convert from local cartesian to geodetic coordinates and return rotation 00170 * matrix. 00171 * 00172 * @param[in] x local cartesian coordinate (meters). 00173 * @param[in] y local cartesian coordinate (meters). 00174 * @param[in] z local cartesian coordinate (meters). 00175 * @param[out] lat latitude of point (degrees). 00176 * @param[out] lon longitude of point (degrees). 00177 * @param[out] h height of point above the ellipsoid (meters). 00178 * @param[out] M if the length of the vector is 9, fill with the rotation 00179 * matrix in row-major order. 00180 * 00181 * Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can 00182 * express \e v as \e column vectors in one of two ways 00183 * - in east, north, up coordinates (where the components are relative to a 00184 * local coordinate system at (\e lat, \e lon, \e h)); call this 00185 * representation \e v1. 00186 * - in \e x, \e y, \e z coordinates (where the components are relative to 00187 * the local coordinate system at (\e lat0, \e lon0, \e h0)); call this 00188 * representation \e v0. 00189 * . 00190 * Then we have \e v1 = <i>M</i><sup>T</sup> ⋅ \e v0, where 00191 * <i>M</i><sup>T</sup> is the transpose of \e M. 00192 **********************************************************************/ 00193 void Reverse(real x, real y, real z, real& lat, real& lon, real& h, 00194 std::vector<real>& M) 00195 const { 00196 if (M.end() == M.begin() + dim2_) { 00197 real t[dim2_]; 00198 IntReverse(x, y, z, lat, lon, h, t); 00199 std::copy(t, t + dim2_, M.begin()); 00200 } else 00201 IntReverse(x, y, z, lat, lon, h, NULL); 00202 } 00203 00204 /** \name Inspector functions 00205 **********************************************************************/ 00206 ///@{ 00207 /** 00208 * @return latitude of the origin (degrees). 00209 **********************************************************************/ 00210 Math::real LatitudeOrigin() const { return _lat0; } 00211 00212 /** 00213 * @return longitude of the origin (degrees). 00214 **********************************************************************/ 00215 Math::real LongitudeOrigin() const { return _lon0; } 00216 00217 /** 00218 * @return height of the origin (meters). 00219 **********************************************************************/ 00220 Math::real HeightOrigin() const { return _h0; } 00221 00222 /** 00223 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00224 * the value of \e a inherited from the Geocentric object used in the 00225 * constructor. 00226 **********************************************************************/ 00227 Math::real MajorRadius() const { return _earth.MajorRadius(); } 00228 00229 /** 00230 * @return \e f the flattening of the ellipsoid. This is the value 00231 * inherited from the Geocentric object used in the constructor. 00232 **********************************************************************/ 00233 Math::real Flattening() const { return _earth.Flattening(); } 00234 ///@} 00235 00236 /// \cond SKIP 00237 /** 00238 * <b>DEPRECATED</b> 00239 * @return \e r the inverse flattening of the ellipsoid. 00240 **********************************************************************/ 00241 Math::real InverseFlattening() const 00242 { return _earth.InverseFlattening(); } 00243 /// \endcond 00244 }; 00245 00246 } // namespace GeographicLib 00247 00248 #endif // GEOGRAPHICLIB_LOCALCARTESIAN_HPP