From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Sat, 10 Oct 2009 00:08:14 +0000 (+0000)
Subject: Fix ST_Area(geography) calculation to be more... correct.
X-Git-Tag: 1.5.0b1~379
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=40f0d6756c59b6e72be4088c49c4d87f35f34aaf;p=postgis

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 f00c1e67c..0cb1be269 100644
--- 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 64123c90f..808d8679e 100644
--- 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);