Fix potential corner case in sphere area calculation (#451)

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

diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c
index 4b8fa95163cbd6ba04d0f8280852bf124d9cb648..c1a05d7baa3f8ede96a707299dccb206bd76497f 100644 (file)
--- a/liblwgeom/cunit/cu_geodetic.c
+++ b/liblwgeom/cunit/cu_geodetic.c
@@ -932,30 +932,29 @@ static void test_spheroid_area(void)
        /* One-degree square */
        lwg = lwgeom_from_ewkt("POLYGON((8.5 2,8.5 1,9.5 1,9.5 2,8.5 2))", PARSER_CHECK_NONE);
        lwgeom_calculate_gbox_geodetic(lwg, &gbox);
-       //a1 = lwgeom_area_sphere(lwg, &gbox, &s);
+       a1 = lwgeom_area_sphere(lwg, &gbox, &s);
        a2 = lwgeom_area_spheroid(lwg, &gbox, &s);
        //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
-       CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.1, 10.0); /* sphere */
+       CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.1, 10.0); /* sphere */
        CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
 
        /* One-degree square *near* dateline */
-       lwg = lwgeom_from_ewkt("POLYGON((178.5 2,178.5 1,179.5 1,179.5 2,178.5 2))", PARSER_CHECK_NONE);
+       lwg = lwgeom_from_ewkt("POLYGON((179.5 2,179.5 1,178.5 1,178.5 2,179.5 2))", PARSER_CHECK_NONE);
        lwgeom_calculate_gbox_geodetic(lwg, &gbox);
-       //a1 = lwgeom_area_sphere(lwg, &gbox, &s);
+       a1 = lwgeom_area_sphere(lwg, &gbox, &s);
        a2 = lwgeom_area_spheroid(lwg, &gbox, &s);
        //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
-       CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.1, 10.0); /* sphere */
+       CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.1, 10.0); /* sphere */
        CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
 
        /* One-degree square *across* dateline */
-       lwg = lwgeom_from_ewkt(" POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", PARSER_CHECK_NONE);
+       lwg = lwgeom_from_ewkt("POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", PARSER_CHECK_NONE);
        lwgeom_calculate_gbox_geodetic(lwg, &gbox);
        a1 = lwgeom_area_sphere(lwg, &gbox, &s);
        a2 = lwgeom_area_spheroid(lwg, &gbox, &s);
        //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
-       //CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.1, 10.0); /* sphere */
+       CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.3679, 10.0); /* sphere */
        CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
-
 }
 
 

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index a09ef0af6d2a0ad80fcfe964d7b43044ae7c76b2..7851d041bb480afeabfecdf8713905457ea994ed 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -1384,12 +1384,29 @@ double ptarray_area_sphere(const POINTARRAY *pa, const POINT2D *pt_outside)
                        else
                                shift = (M_PI - b.lon) + 0.088; /* About 5deg more */
 
+                       LWDEBUGF(4, "before shift a1(%.8g %.8g) b1(%.8g %.8g) c1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon, c1.lat, c1.lon);
                        point_shift(&a1, shift);
                        point_shift(&b1, shift);
                        point_shift(&c1, shift);
+                       LWDEBUGF(4, "after shift a1(%.8g %.8g) b1(%.8g %.8g) c1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon, c1.lat, c1.lon);
                        excess = sphere_excess(&a1, &b1, &c1);
                }
-               else
+               else if( crosses_dateline(&a, &c) )
+               {
+                       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 - c.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);
                }