This is a C++ library meant to provide implementations of the most commonly used
geodetic calculations. The whole library is wrapped around the dso
namespace.
Since December 2021, the build system has been changed from GNU Autotools to scons.
First off, clone the project.
To buid the project just use scons:
scons [OPTIONS]
in the top-level directory. The optional [OPTIONS]
argument, can be any/multiple of:
debug=1
to trigger aDEBUG
build,boost=1
to trigger a build including boost geometry comparisson/test programseigen=1
to use eigen library for vector/matrix operations. Note that because the library uses extensively template code, if you includeggeodesy
header files in another project you should also include the-DUSE_EIGEN
compilation flag when building the dependent project
Installation is trivial, use:
sudo scons install
The folder boost includes source code for testing the algorithms in ggeodesy against the ones implemented in (boost/geometry)[https://www.boost.org/doc/libs/1_74_0/libs/geometry/doc/html/index.html]. To compile these, you will need the respective developement files. Installing boost is usually pretty trivial; relevant documentation can be found on the (boost webpage)[https://www.boost.org/].
Once you have downloaded boost/geometry
you can include the boost
source code in the build process via scons boost=1
TODO
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).
This comparisson function is defined in the test_help.hpp file and used throughout the test/check programs.
Here is a list of the provided utilities:
In general, when refering to coordinate arrays/vector, the elements are in the following order:
- for cartesian vector
[x, y, z]
- for ellipsoidal/geodetic (longitude, latitude, height)
[longitude, latitude, height]
- topocentric
[east, north, up]
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 -π to π).
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);
// ....
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 [φ, λ, height])
- Ellipsoidal to Cartesian (aka [φ, λ, height] to [x, y, z])
- Cartesian to Topocentric (aka [δx, δy, δz] to [north, east, up])
- Topocentric to Cartesian (aka [north, east, up] to [δx, δy, δz])
As is well known in geodesy, the meridian arc length S(\f$\phi\f$) on the earth ellipsoid from the equator to the geographic latitude \f$\phi\f$ 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 range [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 dso
namespace
- static
- dynamic
- build dox with doxygen (or link to dox)
Xanthos, xanthos@mail.ntua.gr Mitsos, danast@mail.ntua.gr
[1] 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
[2] Wiki page on Meridian Arc
[3] H. Moritz, GEODETIC REFERENCE SYSTEM 1980, available here
[4] Li, X.; Li, H.; Liu, G.; Bian, S.; Jiao, C. Simplified Expansions of Common Latitudes with Geodetic Latitude and Geocentric Latitude as Variables. Appl. Sci. 2022, 12, 7818. https://doi.org/10.3390/app12157818
[5] Snyder, J. P. Map projections - A working manual. Professional Paper 1395, U.S. Geological Survey, 1987. doi:10.3133/pp1395.
[6] https://proj.org/operations/conversions/geoc.html
[7] https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf
[8] W. Torge, Geodesy, De Gruyter; 3rd ed. edition (May 30, 2001)