PostGIS  2.2.8dev-r@@SVN_REVISION@@

◆ ptarray_area_spheroid()

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

Definition at line 480 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().

481 {
482  GEOGRAPHIC_POINT a, b;
483  POINT2D p;
484  int i;
485  double area = 0.0;
486  GBOX gbox2d;
487  int in_south = LW_FALSE;
488  double delta_lon_tolerance;
489  double latitude_min;
490 
491  gbox2d.flags = gflags(0, 0, 0);
492 
493  /* Return zero on non-sensical inputs */
494  if ( ! pa || pa->npoints < 4 )
495  return 0.0;
496 
497  /* Get the raw min/max values for the latitudes */
499 
500  if ( SIGNUM(gbox2d.ymin) != SIGNUM(gbox2d.ymax) )
501  lwerror("ptarray_area_spheroid: cannot handle ptarray that crosses equator");
502 
503  /* Geodetic bbox < 0.0 implies geometry is entirely in southern hemisphere */
504  if ( gbox2d.ymax < 0.0 )
505  in_south = LW_TRUE;
506 
507  LWDEBUGF(4, "gbox2d.ymax %.12g", gbox2d.ymax);
508 
509  /* Tolerance for strip area calculation */
510  if ( in_south )
511  {
512  delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymin) / 8.0) - 2.0) / 10000.0;
513  latitude_min = deg2rad(fabs(gbox2d.ymax));
514  }
515  else
516  {
517  delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymax) / 8.0) - 2.0) / 10000.0;
518  latitude_min = deg2rad(gbox2d.ymin);
519  }
520 
521  /* Initialize first point */
522  getPoint2d_p(pa, 0, &p);
523  geographic_point_init(p.x, p.y, &a);
524 
525  for ( i = 1; i < pa->npoints; i++ )
526  {
527  GEOGRAPHIC_POINT a1, b1;
528  double strip_area = 0.0;
529  double delta_lon = 0.0;
530  LWDEBUGF(4, "edge #%d", i);
531 
532  getPoint2d_p(pa, i, &p);
533  geographic_point_init(p.x, p.y, &b);
534 
535  a1 = a;
536  b1 = b;
537 
538  /* Flip into north if in south */
539  if ( in_south )
540  {
541  a1.lat = -1.0 * a1.lat;
542  b1.lat = -1.0 * b1.lat;
543  }
544 
545  LWDEBUGF(4, "in_south %d", in_south);
546 
547  LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) );
548 
549  if ( crosses_dateline(&a, &b) )
550  {
551  double shift;
552 
553  if ( a1.lon > 0.0 )
554  shift = (M_PI - a1.lon) + 0.088; /* About 5deg more */
555  else
556  shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */
557 
558  LWDEBUGF(4, "shift: %.8g", shift);
559  LWDEBUGF(4, "before shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
560  point_shift(&a1, shift);
561  point_shift(&b1, shift);
562  LWDEBUGF(4, "after shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
563 
564  }
565 
566 
567  delta_lon = fabs(b1.lon - a1.lon);
568 
569  LWDEBUGF(4, "a1(%.18g %.18g) b1(%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon);
570  LWDEBUGF(4, "delta_lon %.18g", delta_lon);
571  LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance);
572 
573  if ( delta_lon > 0.0 )
574  {
575  if ( delta_lon < delta_lon_tolerance )
576  {
577  strip_area = spheroid_striparea(&a1, &b1, latitude_min, spheroid);
578  LWDEBUGF(4, "strip_area %.12g", strip_area);
579  area += strip_area;
580  }
581  else
582  {
583  GEOGRAPHIC_POINT p, q;
584  double step = floor(delta_lon / delta_lon_tolerance);
585  double distance = spheroid_distance(&a1, &b1, spheroid);
586  double pDistance = 0.0;
587  int j = 0;
588  LWDEBUGF(4, "step %.18g", step);
589  LWDEBUGF(4, "distance %.18g", distance);
590  step = distance / step;
591  LWDEBUGF(4, "step %.18g", step);
592  p = a1;
593  while (pDistance < (distance - step * 1.01))
594  {
595  double azimuth = spheroid_direction(&p, &b1, spheroid);
596  j++;
597  LWDEBUGF(4, " iteration %d", j);
598  LWDEBUGF(4, " azimuth %.12g", azimuth);
599  pDistance = pDistance + step;
600  LWDEBUGF(4, " pDistance %.12g", pDistance);
601  spheroid_project(&p, spheroid, step, azimuth, &q);
602  strip_area = spheroid_striparea(&p, &q, latitude_min, spheroid);
603  LWDEBUGF(4, " strip_area %.12g", strip_area);
604  area += strip_area;
605  LWDEBUGF(4, " area %.12g", area);
606  p.lat = q.lat;
607  p.lon = q.lon;
608  }
609  strip_area = spheroid_striparea(&p, &b1, latitude_min, spheroid);
610  area += strip_area;
611  }
612  }
613 
614  /* B gets incremented in the next loop, so we save the value here */
615  a = b;
616  }
617  return fabs(area);
618 }
int ptarray_calculate_gbox_cartesian(const POINTARRAY *pa, GBOX *gbox)
Calculate box (x/y) and add values to gbox.
Definition: g_box.c:512
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:172
int npoints
Definition: liblwgeom.h:355
int crosses_dateline(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Definition: lwgeodetic.c:616
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:345
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:441
void point_shift(GEOGRAPHIC_POINT *p, double shift)
Shift a point around by a number of radians.
Definition: lwgeodetic.c:136
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:32
double x
Definition: liblwgeom.h:312
double ymin
Definition: liblwgeom.h:278
#define LW_FALSE
Definition: liblwgeom.h:62
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:61
double ymax
Definition: liblwgeom.h:279
double y
Definition: liblwgeom.h:312
int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
Definition: lwgeom_api.c:448
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:268
Datum distance(PG_FUNCTION_ARGS)
uint8_t flags
Definition: liblwgeom.h:275
uint8_t gflags(int hasz, int hasm, int geodetic)
Construct a new flags char.
Definition: g_util.c:130
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:156
#define deg2rad(d)
Conversion functions.
Definition: lwgeodetic.h:59
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:74
#define SIGNUM(n)
Macro that returns: -1 if n < 0, 1 if n > 0, 0 if n == 0.
Here is the call graph for this function:
Here is the caller graph for this function: