Fix ST_Area(geography) calculation to be more... correct.

git-svn-id: http://svn.osgeo.org/postgis/trunk@4636 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index f00c1e67c4aa699b45b87bbea404a180af4fc987..0cb1be26935947fa69d54dce94556068b3b7fde7 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -162,6 +162,18 @@ void point_rad2deg(GEOGRAPHIC_POINT *p)
        p->lon = rad2deg(p->lon);
 }
 
+/** 
+* Shift a point around by a number of radians
+*/
+static void point_shift(GEOGRAPHIC_POINT *p, double shift)
+{
+       double lon = p->lon + shift;
+       if( lon > M_PI )
+               p->lon = -1.0 * M_PI + (lon - M_PI);
+       return;
+}
+
+
 static int geographic_point_equals(GEOGRAPHIC_POINT g1, GEOGRAPHIC_POINT g2)
 {
        return FP_EQUALS(g1.lat, g2.lat) && FP_EQUALS(g1.lon, g2.lon);
@@ -353,6 +365,28 @@ void y_to_z(POINT3D *p)
 }
 
 
+static inline int crosses_dateline(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e)
+{
+       double sign_s = signum(s.lon);
+       double sign_e = signum(e.lon);
+       double ss = fabs(s.lon);
+       double ee = fabs(e.lon);
+       if( sign_s == sign_e )
+       {
+               return LW_FALSE;
+       }
+       else
+       {
+               double dl = ss + ee;
+               if( dl < M_PI )
+                       return LW_FALSE;
+               else if ( FP_EQUALS(dl, M_PI) )
+                       return LW_FALSE;
+               else
+                       return LW_TRUE;
+       }
+}
+
 /**
 * Returns true if the point p is on the great circle plane.
 * Forms the scalar triple product of A,B,p and if the volume of the
@@ -552,6 +586,27 @@ double sphere_distance_cartesian(POINT3D s, POINT3D e)
        return acos(dot_product(s, e));
 }
 
+/**
+* Given two points on a unit sphere, calculate the direction from s to e.
+*/
+static double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e, double d)
+{
+       double heading = 0.0;
+       
+       if( FP_IS_ZERO(cos(s.lat)) )
+               return (s.lat > 0.0) ? M_PI : 0.0;
+
+       if( sin(e.lon - s.lon) < 0.0 )
+       {
+               heading = acos((sin(e.lat) - sin(s.lat) * cos(d)) / (sin(d) * cos(s.lat)));
+       }
+       else
+       {
+               heading = -1.0 * acos((sin(e.lat) - sin(s.lat) * cos(d)) / (sin(d) * cos(s.lat)));
+       }
+       return heading;
+}
+
 /**
 * Computes the spherical excess of a spherical triangle defined by
 * the three vectices A, B, C. Computes on the unit sphere (i.e., divides
@@ -568,25 +623,15 @@ static double sphere_excess(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, GEOGRAPHIC_P
        double a_dist = sphere_distance(b, c);
        double b_dist = sphere_distance(c, a);
        double c_dist = sphere_distance(a, b);
-       double sign = (a.lon - b.lon) / fabs(a.lon - b.lon);
+       double hca = sphere_direction(c, a, b_dist);
+       double hcb = sphere_direction(c, b, a_dist);
+       double sign = signum(hcb-hca);
        double ss = (a_dist + b_dist + c_dist) / 2.0;
-       double E = tan(ss/2.0) * tan((ss-a_dist)/2.0) * tan((ss-b_dist)/2.0) * tan((ss-c_dist)/2.0);
+       double E = tan(ss/2.0)*tan((ss-a_dist)/2.0)*tan((ss-b_dist)/2.0)*tan((ss-c_dist)/2.0);
        return 4.0 * atan(sqrt(E)) * sign;
 }
 
-/**
-* Given two points on a unit sphere, calculate the direction from s to e.
-*/
-double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e)
-{
-       double latS = s.lat;
-       double latE = e.lon;
-       double dlng = e.lat - s.lon;
-       double heading = atan2(sin(dlng) * cos(latE),
-                              cos(latS) * sin(latE) -
-                              sin(latS) * cos(latE) * cos(dlng)) / M_PI;
-       return heading;
-}
+
 
 /**
 * Returns true if the point p is on the minor edge defined by the
@@ -1324,9 +1369,33 @@ double ptarray_area_sphere(POINTARRAY *pa, POINT2D pt_outside)
 
        for( i = 1; i < pa->npoints; i++ )
        {
+               double excess = 0.0;
+               
                getPoint2d_p(pa, i, &p);
                geographic_point_init(p.x, p.y, &b);
-               area += sphere_excess(a, b, c);
+               
+               if( crosses_dateline(a, b) )
+               {
+                       GEOGRAPHIC_POINT a1 = a, b1 = b, c1 = c;
+                       double shift;
+
+                       if( a.lon > 0.0 )
+                               shift = (M_PI - a.lon) + 0.088; /* About 5deg more */
+                       else
+                               shift = (M_PI - b.lon) + 0.088; /* About 5deg more */
+
+                       point_shift(&a1, shift);
+                       point_shift(&b1, shift);
+                       point_shift(&c1, shift);
+                       excess = sphere_excess(a1, b1, c1);
+               }
+               else
+               {
+                       excess = sphere_excess(a, b, c);
+               }
+               
+               area += excess;
+               
                /* B gets incremented in the next loop, so we save the value here */
                a = b;
        }

diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h
index 64123c90f8bd6fa61f5e736705c3b8e6e1f06fd1..808d8679e2f61f6c9caae020ad684e8377994daf 100644 (file)
--- a/liblwgeom/lwgeodetic.h
+++ b/liblwgeom/lwgeodetic.h
@@ -60,7 +60,6 @@ int clairaut_cartesian(POINT3D start, POINT3D end, GEOGRAPHIC_POINT *g_top, GEOG
 int clairaut_geographic(GEOGRAPHIC_POINT start, GEOGRAPHIC_POINT end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom);
 double sphere_distance(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e);
 double sphere_distance_cartesian(POINT3D s, POINT3D e);
-double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e);
 int sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPHIC_POINT *n);
 int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox);
 int edge_calculate_gbox_slow(GEOGRAPHIC_EDGE e, GBOX *gbox);