# Using GDAL/C++ to calculate distance in meters

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));
``````
• I forgot to include the last line of the code snippet:
– Rex
Apr 26, 2015 at 18:21
• printf("Distance %f\n",point1.Distance(&point2));
– Rex
Apr 26, 2015 at 18:21
• Please edit the question to make corrections/changes. Apr 26, 2015 at 18:37
• You want to lookup vincenty formula. Sadly this is not in ogr as far as I know. Apr 26, 2015 at 20:01
• You need to project your points to a spatial reference system with the horizontal units in metres. See the GDAL OSR projections tutorial gdal.org/osr_tutorial.html Apr 26, 2015 at 21:07

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").

• I like the idea of using GeographicLib. It's quite simple and does exactly what I need.
– Rex
Apr 28, 2015 at 23:16

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;
}
``````

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.

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!