Fix for st_area(geography) over the dateline (#450)

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

diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c
index 688bd5bc2158a2982a97c9f0692cd58f0610a765..4b8fa95163cbd6ba04d0f8280852bf124d9cb648 100644 (file)
--- a/liblwgeom/cunit/cu_geodetic.c
+++ b/liblwgeom/cunit/cu_geodetic.c
@@ -920,7 +920,6 @@ static void test_spheroid_area(void)
        CU_ASSERT_DOUBLE_EQUAL(a1, 89.7211470368, 0.0001); /* sphere */
        CU_ASSERT_DOUBLE_EQUAL(a2, 89.8684316032, 0.0001); /* spheroid */
 
-
        /* Big-ass polygon */
        lwg = lwgeom_from_ewkt("POLYGON((-2 3, -2 4, -1 4, -1 3, -2 3))", PARSER_CHECK_NONE);
        lwgeom_calculate_gbox_geodetic(lwg, &gbox);
@@ -930,6 +929,33 @@ static void test_spheroid_area(void)
        CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.1, 10.0); /* sphere */
        CU_ASSERT_DOUBLE_EQUAL(a2, 12286574431.9, 10.0); /* spheroid */
 
+       /* 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);
+       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(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);
+       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(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);
+       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(a2, 12304814950.073, 100.0); /* spheroid */
+
 }
 
 

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index c62abef36fae38fcf59bf28863192bc51c443e4f..a09ef0af6d2a0ad80fcfe964d7b43044ae7c76b2 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -131,6 +131,8 @@ 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);
+       else
+               p->lon = lon;
        return;
 }
 

diff --git a/liblwgeom/lwspheroid.c b/liblwgeom/lwspheroid.c
index 176a1e0bec22b275014a17ac4159b022ced37247..039f9f7ff1007c3a973ef08829ffe571ccb3fe33 100644 (file)
--- a/liblwgeom/lwspheroid.c
+++ b/liblwgeom/lwspheroid.c
@@ -430,6 +430,8 @@ static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *sphero
 
                LWDEBUGF(4, "in_south %d", in_south);
 
+               LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) );
+
                if ( crosses_dateline(&a, &b) )
                {
                        double shift;
@@ -439,15 +441,18 @@ static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *sphero
                        else
                                shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */
 
+                       LWDEBUGF(4, "shift: %.8g", shift);
+                       LWDEBUGF(4, "before shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
                        point_shift(&a1, shift);
                        point_shift(&b1, shift);
+                       LWDEBUGF(4, "after shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
+                       
                }
 
-               LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) );
 
                delta_lon = fabs(b1.lon - a1.lon);
 
-               LWDEBUGF(4, "(%.18g %.18g) (%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon);
+               LWDEBUGF(4, "a1(%.18g %.18g) b1(%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon);
                LWDEBUGF(4, "delta_lon %.18g", delta_lon);
                LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance);