PostGIS  2.5.0beta1dev-r@@SVN_REVISION@@

◆ ptarray_area_spheroid()

static double ptarray_area_spheroid ( const POINTARRAY pa,
const SPHEROID spheroid 
)
static

Definition at line 494 of file lwspheroid.c.

References area(), crosses_dateline(), deg2rad, distance(), GBOX::flags, geographic_point_init(), getPoint2d_p(), gflags(), GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, LW_FALSE, LW_TRUE, LWDEBUGF, lwerror(), POINTARRAY::npoints, point_shift(), ptarray_calculate_gbox_cartesian(), SIGNUM, spheroid_direction(), spheroid_distance(), spheroid_project(), spheroid_striparea(), POINT2D::x, POINT2D::y, GBOX::ymax, and GBOX::ymin.

Referenced by lwgeom_area_spheroid(), and spheroid_big_b().

495 {
496  GEOGRAPHIC_POINT a, b;
497  POINT2D p;
498  int i;
499  double area = 0.0;
500  GBOX gbox2d;
501  int in_south = LW_FALSE;
502  double delta_lon_tolerance;
503  double latitude_min;
504 
505  gbox2d.flags = gflags(0, 0, 0);
506 
507  /* Return zero on non-sensical inputs */
508  if ( ! pa || pa->npoints < 4 )
509  return 0.0;
510 
511  /* Get the raw min/max values for the latitudes */
513 
514  if ( SIGNUM(gbox2d.ymin) != SIGNUM(gbox2d.ymax) )
515  lwerror("ptarray_area_spheroid: cannot handle ptarray that crosses equator");
516 
517  /* Geodetic bbox < 0.0 implies geometry is entirely in southern hemisphere */
518  if ( gbox2d.ymax < 0.0 )
519  in_south = LW_TRUE;
520 
521  LWDEBUGF(4, "gbox2d.ymax %.12g", gbox2d.ymax);
522 
523  /* Tolerance for strip area calculation */
524  if ( in_south )
525  {
526  delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymin) / 8.0) - 2.0) / 10000.0;
527  latitude_min = deg2rad(fabs(gbox2d.ymax));
528  }
529  else
530  {
531  delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymax) / 8.0) - 2.0) / 10000.0;
532  latitude_min = deg2rad(gbox2d.ymin);
533  }
534 
535  /* Initialize first point */
536  getPoint2d_p(pa, 0, &p);
537  geographic_point_init(p.x, p.y, &a);
538 
539  for ( i = 1; i < pa->npoints; i++ )
540  {
541  GEOGRAPHIC_POINT a1, b1;
542  double strip_area = 0.0;
543  double delta_lon = 0.0;
544  LWDEBUGF(4, "edge #%d", i);
545 
546  getPoint2d_p(pa, i, &p);
547  geographic_point_init(p.x, p.y, &b);
548 
549  a1 = a;
550  b1 = b;
551 
552  /* Flip into north if in south */
553  if ( in_south )
554  {
555  a1.lat = -1.0 * a1.lat;
556  b1.lat = -1.0 * b1.lat;
557  }
558 
559  LWDEBUGF(4, "in_south %d", in_south);
560 
561  LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) );
562 
563  if ( crosses_dateline(&a, &b) )
564  {
565  double shift;
566 
567  if ( a1.lon > 0.0 )
568  shift = (M_PI - a1.lon) + 0.088; /* About 5deg more */
569  else
570  shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */
571 
572  LWDEBUGF(4, "shift: %.8g", shift);
573  LWDEBUGF(4, "before shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
574  point_shift(&a1, shift);
575  point_shift(&b1, shift);
576  LWDEBUGF(4, "after shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
577 
578  }
579 
580 
581  delta_lon = fabs(b1.lon - a1.lon);
582 
583  LWDEBUGF(4, "a1(%.18g %.18g) b1(%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon);
584  LWDEBUGF(4, "delta_lon %.18g", delta_lon);
585  LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance);
586 
587  if ( delta_lon > 0.0 )
588  {
589  if ( delta_lon < delta_lon_tolerance )
590  {
591  strip_area = spheroid_striparea(&a1, &b1, latitude_min, spheroid);
592  LWDEBUGF(4, "strip_area %.12g", strip_area);
593  area += strip_area;
594  }
595  else
596  {
597  GEOGRAPHIC_POINT p, q;
598  double step = floor(delta_lon / delta_lon_tolerance);
599  double distance = spheroid_distance(&a1, &b1, spheroid);
600  double pDistance = 0.0;
601  int j = 0;
602  LWDEBUGF(4, "step %.18g", step);
603  LWDEBUGF(4, "distance %.18g", distance);
604  step = distance / step;
605  LWDEBUGF(4, "step %.18g", step);
606  p = a1;
607  while (pDistance < (distance - step * 1.01))
608  {
609  double azimuth = spheroid_direction(&p, &b1, spheroid);
610  j++;
611  LWDEBUGF(4, " iteration %d", j);
612  LWDEBUGF(4, " azimuth %.12g", azimuth);
613  pDistance = pDistance + step;
614  LWDEBUGF(4, " pDistance %.12g", pDistance);
615  spheroid_project(&p, spheroid, step, azimuth, &q);
616  strip_area = spheroid_striparea(&p, &q, latitude_min, spheroid);
617  LWDEBUGF(4, " strip_area %.12g", strip_area);
618  area += strip_area;
619  LWDEBUGF(4, " area %.12g", area);
620  p.lat = q.lat;
621  p.lon = q.lon;
622  }
623  strip_area = spheroid_striparea(&p, &b1, latitude_min, spheroid);
624  area += strip_area;
625  }
626  }
627 
628  /* B gets incremented in the next loop, so we save the value here */
629  a = b;
630  }
631  return fabs(area);
632 }
int ptarray_calculate_gbox_cartesian(const POINTARRAY *pa, GBOX *gbox)
Calculate box (x/y) and add values to gbox.
Definition: g_box.c:542
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.
Definition: lwspheroid.c:186
int crosses_dateline(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Definition: lwgeodetic.c:635
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...
Definition: lwspheroid.c:359
Datum area(PG_FUNCTION_ARGS)
static double spheroid_striparea(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, double latitude_min, const SPHEROID *spheroid)
This function doesn&#39;t work for edges crossing the dateline or in the southern hemisphere.
Definition: lwspheroid.c:455
void point_shift(GEOGRAPHIC_POINT *p, double shift)
Shift a point around by a number of radians.
Definition: lwgeodetic.c:151
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:47
double x
Definition: liblwgeom.h:330
double ymin
Definition: liblwgeom.h:296
#define LW_FALSE
Definition: liblwgeom.h:76
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
double ymax
Definition: liblwgeom.h:297
double y
Definition: liblwgeom.h:330
double spheroid_direction(const GEOGRAPHIC_POINT *r, const GEOGRAPHIC_POINT *s, const SPHEROID *spheroid)
Computes the direction of the geodesic joining two points on the spheroid.
Definition: lwspheroid.c:282
Datum distance(PG_FUNCTION_ARGS)
uint8_t flags
Definition: liblwgeom.h:293
uint8_t gflags(int hasz, int hasm, int geodetic)
Construct a new flags char.
Definition: g_util.c:145
int getPoint2d_p(const POINTARRAY *pa, uint32_t n, POINT2D *point)
Definition: lwgeom_api.c:338
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:171
#define deg2rad(d)
Conversion functions.
Definition: lwgeodetic.h:74
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
#define SIGNUM(n)
Macro that returns: -1 if n < 0, 1 if n > 0, 0 if n == 0.
uint32_t npoints
Definition: liblwgeom.h:373
Here is the call graph for this function:
Here is the caller graph for this function: