PostGIS  2.1.10dev-r@@SVN_REVISION@@
static double ptarray_area_spheroid ( const POINTARRAY pa,
const SPHEROID spheroid 
)
static

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

368 {
369  GEOGRAPHIC_POINT a, b;
370  POINT2D p;
371  int i;
372  double area = 0.0;
373  GBOX gbox2d;
374  int in_south = LW_FALSE;
375  double delta_lon_tolerance;
376  double latitude_min;
377 
378  gbox2d.flags = gflags(0, 0, 0);
379 
380  /* Return zero on non-sensical inputs */
381  if ( ! pa || pa->npoints < 4 )
382  return 0.0;
383 
384  /* Get the raw min/max values for the latitudes */
386 
387  if ( signum(gbox2d.ymin) != signum(gbox2d.ymax) )
388  lwerror("ptarray_area_spheroid: cannot handle ptarray that crosses equator");
389 
390  /* Geodetic bbox < 0.0 implies geometry is entirely in southern hemisphere */
391  if ( gbox2d.ymax < 0.0 )
392  in_south = LW_TRUE;
393 
394  LWDEBUGF(4, "gbox2d.ymax %.12g", gbox2d.ymax);
395 
396  /* Tolerance for strip area calculation */
397  if ( in_south )
398  {
399  delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymin) / 8.0) - 2.0) / 10000.0;
400  latitude_min = deg2rad(fabs(gbox2d.ymax));
401  }
402  else
403  {
404  delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymax) / 8.0) - 2.0) / 10000.0;
405  latitude_min = deg2rad(gbox2d.ymin);
406  }
407 
408  /* Initialize first point */
409  getPoint2d_p(pa, 0, &p);
410  geographic_point_init(p.x, p.y, &a);
411 
412  for ( i = 1; i < pa->npoints; i++ )
413  {
414  GEOGRAPHIC_POINT a1, b1;
415  double strip_area = 0.0;
416  double delta_lon = 0.0;
417  LWDEBUGF(4, "edge #%d", i);
418 
419  getPoint2d_p(pa, i, &p);
420  geographic_point_init(p.x, p.y, &b);
421 
422  a1 = a;
423  b1 = b;
424 
425  /* Flip into north if in south */
426  if ( in_south )
427  {
428  a1.lat = -1.0 * a1.lat;
429  b1.lat = -1.0 * b1.lat;
430  }
431 
432  LWDEBUGF(4, "in_south %d", in_south);
433 
434  LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) );
435 
436  if ( crosses_dateline(&a, &b) )
437  {
438  double shift;
439 
440  if ( a1.lon > 0.0 )
441  shift = (M_PI - a1.lon) + 0.088; /* About 5deg more */
442  else
443  shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */
444 
445  LWDEBUGF(4, "shift: %.8g", shift);
446  LWDEBUGF(4, "before shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
447  point_shift(&a1, shift);
448  point_shift(&b1, shift);
449  LWDEBUGF(4, "after shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
450 
451  }
452 
453 
454  delta_lon = fabs(b1.lon - a1.lon);
455 
456  LWDEBUGF(4, "a1(%.18g %.18g) b1(%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon);
457  LWDEBUGF(4, "delta_lon %.18g", delta_lon);
458  LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance);
459 
460  if ( delta_lon > 0.0 )
461  {
462  if ( delta_lon < delta_lon_tolerance )
463  {
464  strip_area = spheroid_striparea(&a1, &b1, latitude_min, spheroid);
465  LWDEBUGF(4, "strip_area %.12g", strip_area);
466  area += strip_area;
467  }
468  else
469  {
470  GEOGRAPHIC_POINT p, q;
471  double step = floor(delta_lon / delta_lon_tolerance);
472  double distance = spheroid_distance(&a1, &b1, spheroid);
473  double pDistance = 0.0;
474  int j = 0;
475  LWDEBUGF(4, "step %.18g", step);
476  LWDEBUGF(4, "distance %.18g", distance);
477  step = distance / step;
478  LWDEBUGF(4, "step %.18g", step);
479  p = a1;
480  while (pDistance < (distance - step * 1.01))
481  {
482  double azimuth = spheroid_direction(&p, &b1, spheroid);
483  j++;
484  LWDEBUGF(4, " iteration %d", j);
485  LWDEBUGF(4, " azimuth %.12g", azimuth);
486  pDistance = pDistance + step;
487  LWDEBUGF(4, " pDistance %.12g", pDistance);
488  spheroid_project(&p, spheroid, step, azimuth, &q);
489  strip_area = spheroid_striparea(&p, &q, latitude_min, spheroid);
490  LWDEBUGF(4, " strip_area %.12g", strip_area);
491  area += strip_area;
492  LWDEBUGF(4, " area %.12g", area);
493  p.lat = q.lat;
494  p.lon = q.lon;
495  }
496  strip_area = spheroid_striparea(&p, &b1, latitude_min, spheroid);
497  area += strip_area;
498  }
499  }
500 
501  /* B gets incremented in the next loop, so we save the value here */
502  a = b;
503  }
504  return fabs(area);
505 }
int ptarray_calculate_gbox_cartesian(const POINTARRAY *pa, GBOX *gbox)
Calculate box (x/y) and add values to gbox.
Definition: g_box.c:471
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:59
int npoints
Definition: liblwgeom.h:327
int crosses_dateline(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Definition: lwgeodetic.c:615
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:232
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't work for edges crossing the dateline or in the southern hemisphere.
Definition: lwspheroid.c:328
#define signum(a)
Ape a java function.
Definition: lwgeodetic.h:66
void point_shift(GEOGRAPHIC_POINT *p, double shift)
Shift a point around by a number of radians.
Definition: lwgeodetic.c:137
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:33
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
double x
Definition: liblwgeom.h:284
double ymin
Definition: liblwgeom.h:250
#define LW_FALSE
Definition: liblwgeom.h:52
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:51
double ymax
Definition: liblwgeom.h:251
double y
Definition: liblwgeom.h:284
int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
Definition: lwgeom_api.c:434
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:155
Datum distance(PG_FUNCTION_ARGS)
uint8_t flags
Definition: liblwgeom.h:247
uint8_t gflags(int hasz, int hasm, int geodetic)
Construct a new flags char.
Definition: g_util.c:131
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:157
#define deg2rad(d)
Conversion functions.
Definition: lwgeodetic.h:60
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55

Here is the call graph for this function:

Here is the caller graph for this function: