PostGIS  2.2.8dev-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 172 of file lwspheroid.c.

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(), spheroid_big_b(), and test_spheroid_distance().

173 {
174  double lambda = (b->lon - a->lon);
175  double f = spheroid->f;
176  double omf = 1 - spheroid->f;
177  double u1, u2;
178  double cos_u1, cos_u2;
179  double sin_u1, sin_u2;
180  double big_a, big_b, delta_sigma;
181  double alpha, sin_alpha, cos_alphasq, c;
182  double sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqrsin_sigma, last_lambda, omega;
183  double cos_lambda, sin_lambda;
184  double distance;
185  int i = 0;
186 
187  /* Same point => zero distance */
188  if ( geographic_point_equals(a, b) )
189  {
190  return 0.0;
191  }
192 
193  u1 = atan(omf * tan(a->lat));
194  cos_u1 = cos(u1);
195  sin_u1 = sin(u1);
196  u2 = atan(omf * tan(b->lat));
197  cos_u2 = cos(u2);
198  sin_u2 = sin(u2);
199 
200  omega = lambda;
201  do
202  {
203  cos_lambda = cos(lambda);
204  sin_lambda = sin(lambda);
205  sqrsin_sigma = POW2(cos_u2 * sin_lambda) +
206  POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda));
207  sin_sigma = sqrt(sqrsin_sigma);
208  cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_lambda;
209  sigma = atan2(sin_sigma, cos_sigma);
210  sin_alpha = cos_u1 * cos_u2 * sin_lambda / sin(sigma);
211 
212  /* Numerical stability issue, ensure asin is not NaN */
213  if ( sin_alpha > 1.0 )
214  alpha = M_PI_2;
215  else if ( sin_alpha < -1.0 )
216  alpha = -1.0 * M_PI_2;
217  else
218  alpha = asin(sin_alpha);
219 
220  cos_alphasq = POW2(cos(alpha));
221  cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
222 
223  /* Numerical stability issue, cos2 is in range */
224  if ( cos2_sigma_m > 1.0 )
225  cos2_sigma_m = 1.0;
226  if ( cos2_sigma_m < -1.0 )
227  cos2_sigma_m = -1.0;
228 
229  c = (f / 16.0) * cos_alphasq * (4.0 + f * (4.0 - 3.0 * cos_alphasq));
230  last_lambda = lambda;
231  lambda = omega + (1.0 - c) * f * sin(alpha) * (sigma + c * sin(sigma) *
232  (cos2_sigma_m + c * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
233  i++;
234  }
235  while ( (i < 999) && (lambda != 0.0) && (fabs((last_lambda - lambda)/lambda) > 1.0e-9) );
236 
237  u2 = spheroid_mu2(alpha, spheroid);
238  big_a = spheroid_big_a(u2);
239  big_b = spheroid_big_b(u2);
240  delta_sigma = big_b * sin_sigma * (cos2_sigma_m + (big_b / 4.0) * (cos_sigma * (-1.0 + 2.0 * POW2(cos2_sigma_m)) -
241  (big_b / 6.0) * cos2_sigma_m * (-3.0 + 4.0 * sqrsin_sigma) * (-3.0 + 4.0 * POW2(cos2_sigma_m))));
242 
243  distance = spheroid->b * big_a * (sigma - delta_sigma);
244 
245  /* Algorithm failure, distance == NaN, fallback to sphere */
246  if ( distance != distance )
247  {
248  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);
249  return spheroid->radius * sphere_distance(a, b);
250  }
251 
252  return distance;
253 }
static double spheroid_big_b(double u2)
Definition: lwspheroid.c:46
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:898
static double spheroid_big_a(double u2)
Definition: lwspheroid.c:41
double b
Definition: liblwgeom.h:298
double radius
Definition: liblwgeom.h:302
static double spheroid_mu2(double alpha, const SPHEROID *s)
Definition: lwspheroid.c:35
double f
Definition: liblwgeom.h:299
#define POW2(x)
Definition: lwgeodetic.h:27
Datum distance(PG_FUNCTION_ARGS)
int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *g2)
Definition: lwgeodetic.c:146
double a
Definition: liblwgeom.h:297
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:74
Here is the call graph for this function:
Here is the caller graph for this function: