LocalCartesian.hpp

Go to the documentation of this file.
00001 /**
00002  * \file LocalCartesian.hpp
00003  * \brief Header for GeographicLib::LocalCartesian class
00004  *
00005  * Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com>
00006  * and licensed under the LGPL.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 #if !defined(GEOGRAPHICLIB_LOCALCARTESIAN_HPP)
00011 #define GEOGRAPHICLIB_LOCALCARTESIAN_HPP "$Id: LocalCartesian.hpp 6952 2011-02-14 20:26:44Z karney $"
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 
00032   class LocalCartesian {
00033   private:
00034     typedef Math::real real;
00035     static const size_t dim = 3, dim2 = dim * dim;
00036     const Geocentric _earth;
00037     real _lat0, _lon0, _h0;
00038     real _x0, _y0, _z0, _r[dim2];
00039     void IntForward(real lat, real lon, real h, real& x, real& y, real& z,
00040                     real M[dim2]) const throw();
00041     void IntReverse(real x, real y, real z, real& lat, real& lon, real& h,
00042                     real M[dim2]) const throw();
00043     void MatrixMultiply(real M[dim2]) const throw();
00044   public:
00045 
00046     /**
00047      * Constructor setting the origin.
00048      *
00049      * @param[in] lat0 latitude at origin (degrees).
00050      * @param[in] lon0 longitude at origin (degrees).
00051      * @param[in] h0 height above ellipsoid at origin (meters); default 0.
00052      * @param[in] earth Geocentric object for the transformation; default
00053      *   Geocentric::WGS84.
00054      **********************************************************************/
00055     LocalCartesian(real lat0, real lon0, real h0 = 0,
00056                    const Geocentric& earth = Geocentric::WGS84) throw()
00057       : _earth(earth)
00058     { Reset(lat0, lon0, h0); }
00059 
00060     /**
00061      * Default constructor.
00062      *
00063      * @param[in] earth Geocentric object for the transformation; default
00064      *   Geocentric::WGS84.
00065      *
00066      * Sets \e lat0 = 0, \e lon0 = 0, \e h0 = 0.
00067      **********************************************************************/
00068     explicit LocalCartesian(const Geocentric& earth = Geocentric::WGS84)
00069       throw()
00070       : _earth(earth)
00071     { Reset(real(0), real(0), real(0)); }
00072 
00073     /**
00074      * Reset the origin.
00075      *
00076      * @param[in] lat0 latitude at origin (degrees).
00077      * @param[in] lon0 longitude at origin (degrees).
00078      * @param[in] h0 height above ellipsoid at origin (meters); default 0.
00079      **********************************************************************/
00080     void Reset(real lat0, real lon0, real h0 = 0)
00081       throw();
00082 
00083     /**
00084      * Convert from geodetic to local cartesian coordinates.
00085      *
00086      * @param[in] lat latitude of point (degrees).
00087      * @param[in] lon longitude of point (degrees).
00088      * @param[in] h height of point above the ellipsoid (meters).
00089      * @param[out] x local cartesian coordinate (meters).
00090      * @param[out] y local cartesian coordinate (meters).
00091      * @param[out] z local cartesian coordinate (meters).
00092      *
00093      * \e lat should be in the range [-90, 90]; \e lon and \e lon0 should be in
00094      * the range [-180, 360].
00095      **********************************************************************/
00096     void Forward(real lat, real lon, real h, real& x, real& y, real& z)
00097       const throw() {
00098       IntForward(lat, lon, h, x, y, z, NULL);
00099     }
00100 
00101     /**
00102      * Convert from geodetic to local cartesian coordinates and return rotation
00103      * matrix.
00104      *
00105      * @param[in] lat latitude of point (degrees).
00106      * @param[in] lon longitude of point (degrees).
00107      * @param[in] h height of point above the ellipsoid (meters).
00108      * @param[out] x local cartesian coordinate (meters).
00109      * @param[out] y local cartesian coordinate (meters).
00110      * @param[out] z local cartesian coordinate (meters).
00111      * @param[out] M if the length of the vector is 9, fill with the rotation
00112      *   matrix in row-major order.
00113      *
00114      * Pre-multiplying a unit vector in local cartesian coordinates at (lat,
00115      * lon, h) by \e M transforms the vector to local cartesian coordinates at
00116      * (lat0, lon0, h0).
00117      **********************************************************************/
00118     void Forward(real lat, real lon, real h, real& x, real& y, real& z,
00119                  std::vector<real>& M)
00120       const throw()  {
00121       real t[dim2];
00122       IntForward(lat, lon, h, x, y, z, t);
00123       if (M.end() == M.begin() + dim2)
00124         copy(t, t + dim2, M.begin());
00125     }
00126 
00127     /**
00128      * Convert from local cartesian to geodetic coordinates.
00129      *
00130      * @param[in] x local cartesian coordinate (meters).
00131      * @param[in] y local cartesian coordinate (meters).
00132      * @param[in] z local cartesian coordinate (meters).
00133      * @param[out] lat latitude of point (degrees).
00134      * @param[out] lon longitude of point (degrees).
00135      * @param[out] h height of point above the ellipsoid (meters).
00136      *
00137      * The value of \e lon returned is in the range [-180, 180).
00138      **********************************************************************/
00139     void Reverse(real x, real y, real z, real& lat, real& lon, real& h)
00140       const throw() {
00141       IntReverse(x, y, z, lat, lon, h, NULL);
00142     }
00143 
00144     /**
00145      * Convert from local cartesian to geodetic coordinates and return rotation
00146      * matrix.
00147      *
00148      * @param[in] x local cartesian coordinate (meters).
00149      * @param[in] y local cartesian coordinate (meters).
00150      * @param[in] z local cartesian coordinate (meters).
00151      * @param[out] lat latitude of point (degrees).
00152      * @param[out] lon longitude of point (degrees).
00153      * @param[out] h height of point above the ellipsoid (meters).
00154      * @param[out] M if the length of the vector is 9, fill with the rotation
00155      *   matrix in row-major order.
00156      *
00157      * Pre-multiplying a unit vector in local cartesian coordinates at (lat0,
00158      * lon0, h0) by the transpose of \e M transforms the vector to local
00159      * cartesian coordinates at (lat, lon, h).
00160      **********************************************************************/
00161     void Reverse(real x, real y, real z, real& lat, real& lon, real& h,
00162                  std::vector<real>& M)
00163       const throw() {
00164       real t[dim2];
00165       IntReverse(x, y, z, lat, lon, h, t);
00166       if (M.end() == M.begin() + dim2)
00167         copy(t, t + dim2, M.begin());
00168     }
00169 
00170     /** \name Inspector functions
00171      **********************************************************************/
00172     ///@{
00173     /**
00174      * @return latitude of the origin (degrees).
00175      **********************************************************************/
00176     Math::real LatitudeOrigin() const throw() { return _lat0; }
00177 
00178     /**
00179      * @return longitude of the origin (degrees).
00180      **********************************************************************/
00181     Math::real LongitudeOrigin() const throw() { return _lon0; }
00182 
00183     /**
00184      * @return height of the origin (meters).
00185      **********************************************************************/
00186     Math::real HeightOrigin() const throw() { return _h0; }
00187 
00188     /**
00189      * @return \e a the equatorial radius of the ellipsoid (meters).  This is
00190      *   the value of \e a inherited from the Geocentric object used in the
00191      *   constructor.
00192      **********************************************************************/
00193     Math::real MajorRadius() const throw() { return _earth.MajorRadius(); }
00194 
00195     /**
00196      * @return \e r the inverse flattening of the ellipsoid.  This is the
00197      *   value of \e r inherited from the Geocentric object used in the
00198      *   constructor.  A value of 0 is returned for a sphere (infinite inverse
00199      *   flattening).
00200      **********************************************************************/
00201     Math::real InverseFlattening() const throw()
00202     { return _earth.InverseFlattening(); }
00203     ///@}
00204   };
00205 
00206 } // namespace GeographicLib
00207 
00208 #endif