PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ spheroid_distance()

double spheroid_distance ( const GEOGRAPHIC_POINT a,
const GEOGRAPHIC_POINT b,
const SPHEROID spheroid 
)

Computes the shortest distance along the surface of the spheroid between two points.

Based on Vincenty's formula for the geodetic inverse problem as described in "Geocentric Datum of Australia Technical Manual", Chapter 4. Tested against: http://mascot.gdbc.gov.bc.ca/mascot/util1a.html and http://www.ga.gov.au/nmd/geodesy/datums/vincenty_inverse.jsp

Parameters
a- location of first point.
b- location of second point.
s- spheroid to calculate on
Returns
spheroidal distance between a and b in spheroid units.

Definition at line 191 of file lwspheroid.c.

192 {
193  double lambda = (b->lon - a->lon);
194  double f = spheroid->f;
195  double omf = 1 - spheroid->f;
196  double u1, u2;
197  double cos_u1, cos_u2;
198  double sin_u1, sin_u2;
199  double big_a, big_b, delta_sigma;
200  double alpha, sin_alpha, cos_alphasq, c;
201  double sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqrsin_sigma, last_lambda, omega;
202  double cos_lambda, sin_lambda;
203  double distance;
204  int i = 0;
205 
206  /* Same point => zero distance */
207  if ( geographic_point_equals(a, b) )
208  {
209  return 0.0;
210  }
211 
212  u1 = atan(omf * tan(a->lat));
213  cos_u1 = cos(u1);
214  sin_u1 = sin(u1);
215  u2 = atan(omf * tan(b->lat));
216  cos_u2 = cos(u2);
217  sin_u2 = sin(u2);
218 
219  omega = lambda;
220  do
221  {
222  cos_lambda = cos(lambda);
223  sin_lambda = sin(lambda);
224  sqrsin_sigma = POW2(cos_u2 * sin_lambda) +
225  POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda));
226  sin_sigma = sqrt(sqrsin_sigma);
227  cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_lambda;
228  sigma = atan2(sin_sigma, cos_sigma);
229  sin_alpha = cos_u1 * cos_u2 * sin_lambda / sin(sigma);
230 
231  /* Numerical stability issue, ensure asin is not NaN */
232  if ( sin_alpha > 1.0 )
233  alpha = M_PI_2;
234  else if ( sin_alpha < -1.0 )
235  alpha = -1.0 * M_PI_2;
236  else
237  alpha = asin(sin_alpha);
238 
239  cos_alphasq = POW2(cos(alpha));
240  cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
241 
242  /* Numerical stability issue, cos2 is in range */
243  if ( cos2_sigma_m > 1.0 )
244  cos2_sigma_m = 1.0;
245  if ( cos2_sigma_m < -1.0 )
246  cos2_sigma_m = -1.0;
247 
248  c = (f / 16.0) * cos_alphasq * (4.0 + f * (4.0 - 3.0 * cos_alphasq));
249  last_lambda = lambda;
250  lambda = omega + (1.0 - c) * f * sin(alpha) * (sigma + c * sin(sigma) *
251  (cos2_sigma_m + c * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
252  i++;
253  }
254  while ( (i < 999) && (lambda != 0.0) && (fabs((last_lambda - lambda)/lambda) > 1.0e-9) );
255 
256  u2 = spheroid_mu2(alpha, spheroid);
257  big_a = spheroid_big_a(u2);
258  big_b = spheroid_big_b(u2);
259  delta_sigma = big_b * sin_sigma * (cos2_sigma_m + (big_b / 4.0) * (cos_sigma * (-1.0 + 2.0 * POW2(cos2_sigma_m)) -
260  (big_b / 6.0) * cos2_sigma_m * (-3.0 + 4.0 * sqrsin_sigma) * (-3.0 + 4.0 * POW2(cos2_sigma_m))));
261 
262  distance = spheroid->b * big_a * (sigma - delta_sigma);
263 
264  /* Algorithm failure, distance == NaN, fallback to sphere */
265  if ( distance != distance )
266  {
267  lwerror("spheroid_distance returned NaN: (%.20g %.20g) (%.20g %.20g) a = %.20g b = %.20g",a->lat, a->lon, b->lat, b->lon, spheroid->a, spheroid->b);
268  return spheroid->radius * sphere_distance(a, b);
269  }
270 
271  return distance;
272 }
double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Given two points on a unit sphere, calculate their distance apart in radians.
Definition: lwgeodetic.c:948
int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *g2)
Definition: lwgeodetic.c:170
#define POW2(x)
Definition: lwgeodetic.h:47
Datum distance(PG_FUNCTION_ARGS)
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
static double spheroid_big_b(double u2)
Definition: lwspheroid.c:60
static double spheroid_big_a(double u2)
Definition: lwspheroid.c:55
static double spheroid_mu2(double alpha, const SPHEROID *s)
Definition: lwspheroid.c:49
double radius
Definition: liblwgeom.h:321
double a
Definition: liblwgeom.h:316
double b
Definition: liblwgeom.h:317
double f
Definition: liblwgeom.h:318

References SPHEROID::a, SPHEROID::b, distance(), SPHEROID::f, geographic_point_equals(), GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, lwerror(), POW2, SPHEROID::radius, sphere_distance(), spheroid_big_a(), spheroid_big_b(), and spheroid_mu2().

Referenced by circ_tree_distance_tree(), ptarray_area_spheroid(), ptarray_distance_spheroid(), ptarray_length_spheroid(), and test_spheroid_distance().

Here is the call graph for this function:
Here is the caller graph for this function: