Using GDAL/C++ to calculate distance in meters

Question

I'm using GDAL in C++, with WGS84. I want to calculate distance between two points in meters. The code snippet below tries to find the distance between 25N/75E and 25N/76E. The code, as is, returns 1.0000, which is just a Cartesian number. I need the answer in meters or feet or miles. How do I do that?

OGRSpatialReference  wgs84;
OGRPoint   point1,
           point2;

assert(wgs84.SetWellKnownGeogCS("WGS84")==OGRERR_NONE);
point1.assignSpatialReference(&wgs84);
point2.assignSpatialReference(&wgs84);
point1.setY(25.0000); point1.setX(75.0000);
printf("Distance %f\n",point1.Distance(&point2));

point2.setY(25.0000); point2.setX(76.0000);
printf("Distance %f\n",point1.Distance(&point2)); 

Answer 1

GDAL depends on Proj.4 and the latest version of Proj.4 (4.9.1) includes an API for geodesic calculations (and the geod command uses this API). This is a translation of the C++ geodesic routines in GeographicLib to C. The documentation on the C interface is here. Alternatively, you can use the C++ interface of GeographicLib directly (GeographicLib and GDAL won't interfere with one another). I recommend against implementing Vincenty's algorithm, because it doesn't always converge.

An example of using the C interface is inverse.c.

An example of using the C++ interface is given here (search for "Example of use").

Answer 2

Before you use this, please check out this link for a solid Vincenty example. Luckily, Javascript is very C-like so its not too hard to make it work for C++.

Please also forgive any coding no-nos as this was thrown together rather quickly. My Datum is WGS84. I ran this through g++ and got a matching answer as their online calculator.

I have never been able to find a good Vincenty algorithm in C++ in either GDAL or Proj. Boost has one but it only seems to be in the newer versions of Boost Geometry.

// C++ Standard Libraries
#include <cmath>
#include <iostream>

/**
 * Vincenty Distance
*/
double Vincenty_Distance( const double& latitude_01, const double& longitude_01,
                          const double& latitude_02, const double& longitude_02,
                          const double& a,
                          const double& b )
{
    // Flattening
    const double f = (a-b)/a;

    // tan U1
    const double tan_U1 = (1-f) * std::tan(latitude_01);
    const double tan_U2 = (1-f) * std::tan(latitude_02);

    // Longitudinal Distance
    const double cos_U1 = 1 / std::sqrt(1 + tan_U1 * tan_U1);
    const double cos_U2 = 1 / std::sqrt(1 + tan_U2 * tan_U2);
    const double sin_U1 = tan_U1 * cos_U1;
    const double sin_U2 = tan_U2 * cos_U2;

    // Iterate until complete
    const double L = longitude_02 - longitude_01;
    double lambda = L;
    double diff, sigma;
    double cos_alpha_sq, cos_2sigma_m;
    double cos_sigma, sin_sigma;

    while( true ){

        // 
        double sin_lambda = std::sin( lambda );
        double cos_lambda = std::cos( lambda );

        double c1 = (cos_U2 * sin_lambda)*(cos_U2 * sin_lambda);
        double c2 = (cos_U1 * sin_U2);
        double c3 = (sin_U1 * cos_U2 * cos_lambda);


        //  sin sigma
        sin_sigma = std::sqrt( c1 + ( c2 - c3 )*( c2 - c3 ) );

        // cos sigma
        cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda;

        // sigma
        sigma = std::atan2( sin_sigma, cos_sigma );

        // sin alpha
        double sin_alpha = (cos_U1 * cos_U2 * sin_lambda)/(sin_sigma);

        // cos^2 alpha
        cos_alpha_sq = 1 - (sin_alpha*sin_alpha);

        // cos^2 2sigmam
        cos_2sigma_m = cos_sigma - (2 * sin_U1 * sin_U2)/(cos_alpha_sq);

        // C
        double C = (f/16.0) * cos_alpha_sq * (4 + f * (4 - 3 * cos_alpha_sq));

        // Update Lambda
        diff = lambda;
        lambda = L + (1-C) * f * sin_alpha * (sigma + C * sin_sigma * ( cos_2sigma_m + C * cos_sigma * (-1 + 2 * cos_2sigma_m*cos_2sigma_m)));
        diff = lambda - diff;
        if( std::fabs(diff) < 0.00001 ){ break; }
    }

    // U2
    double u_sq = cos_alpha_sq  * (a*a - b*b)/(b*b);

    // Compute A, B
    double A = 1 + (u_sq/16384) * (4096 + u_sq * (-768 + u_sq * (320 - 175 * u_sq)));

    double B = (u_sq / 1024) * (256 + u_sq * (-128 + u_sq * (-128 + u_sq * (74 - 47 * u_sq))));

    // Sigma
    double cos_2sigma_m_sq = cos_2sigma_m * cos_2sigma_m;
    double delta_sigma = B * sin_sigma * ( cos_2sigma_m + (B/4.0) * (cos_sigma * (-1 * 2 * cos_2sigma_m_sq ) - (B/6.0) * cos_2sigma_m * (-3 + 4 * sin_sigma*sin_sigma) * (-3 + 4 * cos_2sigma_m_sq)));

    // Distance
    double s = b * A * (sigma - delta_sigma);
    return s;
}

/**
 * @brief Main Function
 */
int main( int argc, char* argv[] )
{

    // Set our coordinates
    const double latitude_01  =   39.5 * M_PI / 180.0;
    const double longitude_01 = -120.5 * M_PI / 180.0;

    const double latitude_02 =    40 * M_PI / 180.0;
    const double longitude_02 = -119 * M_PI / 180.0;

    // Set the datum components
    const double a = 6378137.0;
    const double b = 6356752.314245;

    // Vincenty Distance
    double distance = Vincenty_Distance( latitude_01, longitude_01,
                                         latitude_02, longitude_02,
                                         a, b);

    // Distance
    std::cout << std::fixed << "Distance: " << distance << std::endl;

    // Exit
    return 0;
}

Answer 3

You could try using Proj4. The geod command will allow you to calculate the geodesic distance between to points. https://trac.osgeo.org/proj/wiki/man_geod

I have not used the C++ API, so perhaps someone else can provide an example.

Answer 4

msmith81886's example works well. Vincenty's method, however is flawed and old, and can go out to lunch in certain circumstances... when the two points are almost opposite each other on the globe (antipodal).

Karney improved this method in 2013 by three orders of magnitude, and fixed the bounds case issues. See this wiki article for info.

https://en.wikipedia.org/wiki/Geographical_distance

Karney has also, awesomely, built a very nice library available in ports to most languages you are likely to use, including C++, called geographiclib, available on sourceforge.

https://geographiclib.sourceforge.io/

enjoy!