PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ spheroid_project()

int spheroid_project ( const GEOGRAPHIC_POINT r,
const SPHEROID spheroid,
double  distance,
double  azimuth,
GEOGRAPHIC_POINT g 
)

Given a location, an azimuth and a distance, computes the location of the projected point.

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

Parameters
r- location of first point.
distance- distance in meters.
azimuth- azimuth in radians.
Returns
s - location of projected point.

Definition at line 364 of file lwspheroid.c.

365 {
366  double omf = 1 - spheroid->f;
367  double tan_u1 = omf * tan(r->lat);
368  double u1 = atan(tan_u1);
369  double sigma, last_sigma, delta_sigma, two_sigma_m;
370  double sigma1, sin_alpha, alpha, cos_alphasq;
371  double u2, A, B;
372  double lat2, lambda, lambda2, C, omega;
373  int i = 0;
374 
375  if (azimuth < 0.0)
376  {
377  azimuth = azimuth + M_PI * 2.0;
378  }
379  if (azimuth > (M_PI * 2.0))
380  {
381  azimuth = azimuth - M_PI * 2.0;
382  }
383 
384  sigma1 = atan2(tan_u1, cos(azimuth));
385  sin_alpha = cos(u1) * sin(azimuth);
386  alpha = asin(sin_alpha);
387  cos_alphasq = 1.0 - POW2(sin_alpha);
388 
389  u2 = spheroid_mu2(alpha, spheroid);
390  A = spheroid_big_a(u2);
391  B = spheroid_big_b(u2);
392 
393  sigma = (distance / (spheroid->b * A));
394  do
395  {
396  two_sigma_m = 2.0 * sigma1 + sigma;
397  delta_sigma = B * sin(sigma) * (cos(two_sigma_m) + (B / 4.0) * (cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)) - (B / 6.0) * cos(two_sigma_m) * (-3.0 + 4.0 * POW2(sin(sigma))) * (-3.0 + 4.0 * POW2(cos(two_sigma_m))))));
398  last_sigma = sigma;
399  sigma = (distance / (spheroid->b * A)) + delta_sigma;
400  i++;
401  }
402  while (i < 999 && fabs((last_sigma - sigma) / sigma) > 1.0e-9);
403 
404  lat2 = atan2((sin(u1) * cos(sigma) + cos(u1) * sin(sigma) *
405  cos(azimuth)), (omf * sqrt(POW2(sin_alpha) +
406  POW2(sin(u1) * sin(sigma) - cos(u1) * cos(sigma) *
407  cos(azimuth)))));
408  lambda = atan2((sin(sigma) * sin(azimuth)), (cos(u1) * cos(sigma) -
409  sin(u1) * sin(sigma) * cos(azimuth)));
410  C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
411  omega = lambda - (1.0 - C) * spheroid->f * sin_alpha * (sigma + C * sin(sigma) *
412  (cos(two_sigma_m) + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)))));
413  lambda2 = r->lon + omega;
414  g->lat = lat2;
415  g->lon = lambda2;
416  return LW_SUCCESS;
417 }
char * r
Definition: cu_in_wkt.c:24
#define LW_SUCCESS
Definition: liblwgeom.h:80
#define POW2(x)
Definition: lwgeodetic.h:47
Datum distance(PG_FUNCTION_ARGS)
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 b
Definition: liblwgeom.h:317
double f
Definition: liblwgeom.h:318

References SPHEROID::b, distance(), SPHEROID::f, GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, LW_SUCCESS, POW2, r, spheroid_big_a(), spheroid_big_b(), and spheroid_mu2().

Referenced by lwgeom_project_spheroid(), and ptarray_area_spheroid().

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