self-hosted/postgis
self-hosted/postgis
1
0
Fork 0
mirror of https://git.osgeo.org/gitea/postgis/postgis synced 2024-10-24 09:02:37 +00:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse source 
git-svn-id: http://svn.osgeo.org/postgis/trunk@5345 b70326c6-7e19-0410-871a-916f4a2858ee
This commit is contained in:
Paul Ramsey 2010-02-25 15:13:48 +00:00
parent 735f89558a
commit d8095c6034
2 changed files with 30 additions and 14 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

25
liblwgeom/cunit/cu_geodetic.c
View file

@ -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);
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);
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);
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);
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, 12360265021.3679, 10.0); /* sphere */
CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
}

19
liblwgeom/lwgeodetic.c
View file

@ -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 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
else
{
excess = sphere_excess(&a, &b, &c);
}