This is a C++ library meant to provide implementations of the most commonly used
geodetic calculations. The whole library is wrapped around the ngpt
namespace.
Source code is ISO C++17. Compilation should be trivial using any C++ compiler
supporting the c++17
standard (option -std=c++17
in gcc and clang).
This software is meant to be implemented on Unix-type OS's. No effort will be undertaken for compatibility with other OS types.
To compile the library, just follow the basic steps: (note that the library is still at development phase so users need to configure the project before compiling)
For the following, ROOTDIR
will be the root directory of this repository,
aka the directory under which /src
, /test
and /doc
folders live.
If you do not need the DEBUG version (which most probably you don't), create the Makefile.am
templates. This means that you
should rename Makefile.am.production and Makefile.am.production to
src/Makefile.am
and test/Makefile.am
respectively, that is:
mv src/Makefile.am.production src/Makefile.am
mv test/Makefile.am.production test/Makefile.am
Then run Autotools and compile:
autoreconf -if
./configure
make
sudo make install
In the ggeodesy/test
folder you should be able to see a list of executables; run ggeodesy/test/testGeodesy
to validate the library.
After a succesefull installation, users should have:
- all library header files in
/usr/local/include/ggeodesy/
- the library (both a static and shared) in
/usr/local/lib/
To run a validity check, just run: make check
at the root directory. Hopefully,
you 'll see all checks passing!
Link, include and have fun!
For testing purposes, we need to compare (floating point) results, aka check if two floating point numbers are equal. But what does equal mean? Well, in this case for most of the test programs we define the floating point comparisson (when comparing numbers other than zero) as:
diff = abs(a - b)
diff < max(abs(a), abs(b)) * tolerance or diff < tolerance
where tolerance is defined as: std::numeric_limits<double>::epsilon()
. That is,
we use a tolerance value proportional to max(a,b) when comparing two (floating point) numbers.
(see also this SO thread,
and this also).
Here is a list of the provided utilities:
The library includes functions for transforming between radians, decimal degrees and hexicondal degrees. Namely:
ngpt::deg2rad
transforms radians to decimal degreesngpt::decd2hexd
transforms decimal degrees to hexicondal degreesngpt::rad2deg
transforms decimal degrees to radiansngpt::rad2hexd
transforms radians to hexicondal degreesngpt::hexd2decd
transforms hexicondal degrees to decimal degreesngpt::hexd2rad
transforms hexicondal degrees to radiansngpt::rad2sec
convert radians to seconds (of degree)
Additionaly the function ngpt::normalize_angle
can normalize a given angle
to a specified range (e.g. in range -\f$\pi\f$ to \f$\pi\f$).
All of the above functions are defined in the header file units.hpp. For usage examples, see test_units.hpp and test_angle_normalization.hpp.
Currently there are implementations for the (reference) ellipsoids:
Reference ellipsoids can be used in either on of two ways:
- via using the
enum
classngpt::ellipsoid
(e.g.ngpt::ellipsoid::grs80
), or - via the class
ngpt::Ellipsoid
Users can easily add more reference ellipsoids if they need to, via constructing
an Ellipsoid
instance with the wanted parameters (aka semi-major axis and flattening), e.g.
using namespace ngpt;
Ellipsoid myell = Ellipsoid(6378136.0e0/*semi-major axis*/, 1/298.25784/*flattening factor*/);
// use the created ellipsoid in some way ....
double semi_major = myell.semi_major();
double N = myell.N(some_latitude);\f$\pi\f$
// ....
Users can also expand the source code to add a reference ellipsoid in the enum
(class)
ellipsoid
. In this case, you will need one entry in the ellipsoid
enum and a
respective specialization of the template <ellipsoid E> struct ellipsoid_traits {};
class.
Note that most of the constructors and function (for the Ellipsoid
class and the
ellipsoid
enum are constexpr
. Hence, the following code is computed at compile-time:
using namespace ngpt;
constexpr auto wgs84 = Ellipsoid(ellipsoid::wgs84);
constexpr auto grs80 = Ellipsoid(ellipsoid_traits<ellipsoid::grs80>::a,
ellipsoid_traits<ellipsoid::grs80>::f);
constexpr auto pz90 = Ellipsoid(ellipsoid::pz90);
static_assert(wgs84.eccentricity_squared() ==
eccentricity_squared<ellipsoid::wgs84>());
static_assert(grs80.semi_minor() == semi_minor<ellipsoid::grs80>());
static_assert(std::abs(pz90.eccentricity_squared()-0.0066943662)<1e-9);
For more information on how to use the reference ellipsoids, see e.g. test_ellipsoid.hpp.
The following coordinate transformations are provided (for points on some reference ellipsoid):
- Cartesian to Ellipsoidal (aka [x, y, z] to [\f$\phi\f$, \f$\lambda\f$, height])
- Ellipsoidal to Cartesian (aka [\f$\phi\f$, \f$\lambda\f$, height] to [x, y, z])
- Cartesian to Topocentric (aka [\f$delta\f$x, \f$delta\f$y, \f$delta\f$z] to [north, east, up])
- Topocentric to Cartesian (aka [north, east, up] to [\f$delta\f$x, \f$delta\f$y, \f$delta\f$z])
As is well known in geodesy, the meridian arc length S( \phi ) on the earth ellipsoid from the equator to the geographic latitude \phi includes an elliptic integral and cannot be expressed explicitly using a combination of elementary functions. The library provides three implementations for computing the meridian arc length using approximations.
-
ngpt::core::meridian_arc_length_impl1
: implementation using a binomial series with respect to e^2 obtained by truncating the expansion at order e^10 -
ngpt::core::meridian_arc_length_impl2
: the Bessel’s formula implementation -
ngpt::core::meridian_arc_length_impl3
: the Helmert's formula implementation
The above implementations agree within the following precision limits:
- impl1 vs impl2 : 3e-5 meters
- impl2 vs impl3 : 1e-6 meters
- impl2 vs impl3 : 3e-5 meters
Users can select one of the algorithms via the (optional) intput parameter alg
(using values in rabge [0,2]) when calling the function
template<ellipsoid E> double meridian_arc_length(double lat, int alg=0)
(aka
by default the function will use the
More information and implementation details can be found in Kawase, 2011. Usage examples of the algorithm(s) can be found in the file test_meridian_arc.hpp
Note that if we only want the arc length of an infinitesimal element of the meridian the
computation is way more straight-forward Meridian Arc; for this computation users may use the
function (template) infinitesimal_meridian_arc
.
The whole of the library is wrapped around the ngpt
namespace
- static
- dynamic
- build dox with doxygen (or link to dox)
Xanthos, [email protected] Mitsos, [email protected]
K. Kawase, A General Formula for Calculating Meridian Arc Length and its Application to Coordinate Conversion in the Gauss-Krüger Projection, Bulletin of the Geospatial Information Authority of Japan, Vol.59 December, 2011
Wiki page on Meridian Arc