From: Paul Ramsey Date: Thu, 11 Oct 2012 22:48:11 +0000 (+0000) Subject: Improve support for ST_Area(geography) over dateline and poles (#2006, #2039) X-Git-Tag: 2.1.0beta2~558 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=e9bf8eeb82ac0aa011450d9a5d0b34589adfcdaa;p=postgis Improve support for ST_Area(geography) over dateline and poles (#2006, #2039) git-svn-id: http://svn.osgeo.org/postgis/trunk@10407 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/NEWS b/NEWS index ab4dd55b8..b1abe3318 100644 --- a/NEWS +++ b/NEWS @@ -45,6 +45,7 @@ PostGIS 2.1.0 - #1978, wrong answer when calculating length of a closed circular arc (circle) - #1780, support ST_GeoHash for geography - #2021, Added multi-band support to ST_Union(raster, ...) aggregate function + - #2006, better support of ST_Area(geography) over poles and dateline * Fixes * diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c index 4a012d62d..393704b27 100644 --- a/liblwgeom/cunit/cu_geodetic.c +++ b/liblwgeom/cunit/cu_geodetic.c @@ -990,7 +990,7 @@ static void test_spheroid_area(void) a1 = lwgeom_area_sphere(lwg, &s); a2 = lwgeom_area_spheroid(lwg, &s); //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); - CU_ASSERT_DOUBLE_EQUAL(a1, 89.7211470368, 0.0001); /* sphere */ + CU_ASSERT_DOUBLE_EQUAL(a1, 89.7127703297, 0.0001); /* sphere */ CU_ASSERT_DOUBLE_EQUAL(a2, 89.8684316032, 0.0001); /* spheroid */ lwgeom_free(lwg); @@ -1009,7 +1009,7 @@ static void test_spheroid_area(void) lwgeom_calculate_gbox_geodetic(lwg, &gbox); a1 = lwgeom_area_sphere(lwg, &s); a2 = lwgeom_area_spheroid(lwg, &s); - //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); + printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.1, 10.0); /* sphere */ CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */ lwgeom_free(lwg); @@ -1157,6 +1157,24 @@ static void test_lwgeom_segmentize_sphere(void) return; } +static void test_lwgeom_area_sphere(void) +{ + LWGEOM *lwg; + double area; + SPHEROID s; + + /* Init to WGS84 */ + spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS); + + /* Simple case */ + lwg = lwgeom_from_wkt("POLYGON((1 1, 1 2, 2 2, 2 1, 1 1))", LW_PARSER_CHECK_NONE); + area = lwgeom_area_sphere(lwg, &s); + + CU_ASSERT_DOUBLE_EQUAL(area, 12360265021.3561, 0.01); + lwgeom_free(lwg); + return; +} + /* ** Used by test harness to register the tests in this file. */ @@ -1164,6 +1182,7 @@ CU_TestInfo geodetic_tests[] = { PG_TEST(test_sphere_direction), PG_TEST(test_sphere_project), + PG_TEST(test_lwgeom_area_sphere), PG_TEST(test_signum), PG_TEST(test_gbox_from_spherical_coordinates), PG_TEST(test_gserialized_get_gbox_geocentric), diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index ee83055b4..8f4d78355 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -564,12 +564,12 @@ int crosses_dateline(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e) } /** -* 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 -* resulting parallelepiped is near zero the point p is on the -* great circle plane. +* Returns -1 if the point is to the left of the plane formed +* by the edge, 1 if the point is to the right, and 0 if the +* point is on the plane. */ -int edge_point_on_plane(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p) +static int +edge_point_side(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p) { POINT3D normal, pt; double w; @@ -583,9 +583,79 @@ int edge_point_on_plane(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p) if ( FP_IS_ZERO(w) ) { LWDEBUG(4, "point is on plane (dot product is zero)"); - return LW_TRUE; + return 0; } - LWDEBUG(4, "point is not on plane"); + + if ( w < 0 ) + return -1; + else + return 1; +} + +/** +* Returns the angle in radians at point B of the triangle formed by A-B-C +*/ +static double +sphere_angle(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const GEOGRAPHIC_POINT *c) +{ + POINT3D normal1, normal2; + robust_cross_product(b, a, &normal1); + robust_cross_product(b, c, &normal2); + normalize(&normal1); + normalize(&normal2); + return sphere_distance_cartesian(&normal1, &normal2); +} + +/** +* Computes the spherical area of a triangle. If C is to the left of A/B, +* the area is negative. If C is to the right of A/B, the area is positive. +* +* @param a The first triangle vertex. +* @param b The second triangle vertex. +* @param c The last triangle vertex. +* @return the signed area in radians. +*/ +static double +sphere_signed_area(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const GEOGRAPHIC_POINT *c) +{ + double angle_a, angle_b, angle_c; + double area_radians = 0.0; + int side; + GEOGRAPHIC_EDGE e; + + angle_a = sphere_angle(b,a,c); + angle_b = sphere_angle(a,b,c); + angle_c = sphere_angle(b,c,a); + + area_radians = angle_a + angle_b + angle_c - M_PI; + + /* What's the direction of the B/C edge? */ + e.start = *a; + e.end = *b; + side = edge_point_side(&e, c); + + /* Co-linear points implies no area */ + if ( side == 0 ) + return 0.0; + + /* Add the sign to the area */ + return side * area_radians; +} + + + +/** +* 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 +* resulting parallelepiped is near zero the point p is on the +* great circle plane. +*/ +int edge_point_on_plane(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p) +{ + int side = edge_point_side(e, p); + if ( side == 0 ) + return LW_TRUE; + return LW_FALSE; } @@ -806,7 +876,6 @@ double sphere_direction(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e, do 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 @@ -832,7 +901,6 @@ static double sphere_excess(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b } - /** * Returns true if the point p is on the minor edge defined by the * end points of e. @@ -1749,72 +1817,31 @@ lwgeom_segmentize_sphere(const LWGEOM *lwg_in, double max_seg_length) * Returns the area of the ring (ring must be closed) in square radians (surface of * the sphere is 4*PI). */ -double ptarray_area_sphere(const POINTARRAY *pa, const POINT2D *pt_outside) +double +ptarray_area_sphere(const POINTARRAY *pa) { - GEOGRAPHIC_POINT a, b, c; - POINT2D p; int i; + const POINT2D *p; + GEOGRAPHIC_POINT a, b, c; double area = 0.0; - - /* Return zero on non-sensical inputs */ + + /* Return zero on nonsensical inputs */ if ( ! pa || pa->npoints < 4 ) return 0.0; - - geographic_point_init(pt_outside->x, pt_outside->y, &c); - - /* Initialize first point */ - getPoint2d_p(pa, 0, &p); - geographic_point_init(p.x, p.y, &a); - - for ( i = 1; i < pa->npoints; i++ ) + + p = getPoint2d_cp(pa, 0); + geographic_point_init(p->x, p->y, &a); + p = getPoint2d_cp(pa, 1); + geographic_point_init(p->x, p->y, &b); + + for ( i = 2; i < pa->npoints-1; i++ ) { - double excess = 0.0; - - getPoint2d_p(pa, i, &p); - geographic_point_init(p.x, p.y, &b); - - 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 */ - - 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 - { - excess = sphere_excess(&a, &b, &c); - } - - area += excess; - - /* B gets incremented in the next loop, so we save the value here */ - a = b; + p = getPoint2d_cp(pa, i); + geographic_point_init(p->x, p->y, &c); + area += sphere_signed_area(&a, &b, &c); + b = c; } + return fabs(area); } @@ -2106,7 +2133,6 @@ static double ptarray_distance_spheroid(const POINTARRAY *pa1, const POINTARRAY double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid) { int type; - POINT2D pt_outside; double radius2 = spheroid->radius * spheroid->radius; GBOX gbox; gbox.flags = 0; @@ -2124,16 +2150,6 @@ double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid) if ( ! ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) ) return 0.0; - /* Make sure we have boxes */ - if ( lwgeom->bbox ) - gbox = *(lwgeom->bbox); - else - lwgeom_calculate_gbox_geodetic(lwgeom, &gbox); - - gbox_pt_outside(&gbox, &pt_outside); - - LWDEBUGF(2, "pt_outside = POINT(%.20g %.20g)", pt_outside.x, pt_outside.y); - /* Actually calculate area */ if ( type == POLYGONTYPE ) { @@ -2146,12 +2162,12 @@ double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid) return 0.0; /* First, the area of the outer ring */ - area += radius2 * ptarray_area_sphere(poly->rings[0], &pt_outside); + area += radius2 * ptarray_area_sphere(poly->rings[0]); /* Subtract areas of inner rings */ for ( i = 1; i < poly->nrings; i++ ) { - area -= radius2 * ptarray_area_sphere(poly->rings[i], &pt_outside); + area -= radius2 * ptarray_area_sphere(poly->rings[i]); } return area; } diff --git a/regress/tickets.sql b/regress/tickets.sql index 7a31fede6..919f854b5 100644 --- a/regress/tickets.sql +++ b/regress/tickets.sql @@ -272,7 +272,7 @@ SELECT ST_AsText(the_geog) as the_pt, ST_Area(ST_Buffer(the_geog,10)) As the_area, ST_Area(geography(ST_Transform(ST_Buffer(ST_Transform(geometry(the_geog),utm_srid),10),4326))) As geog_utm_area FROM utm_dots -WHERE ST_Area(ST_Buffer(the_geog,10)) NOT between 307 and 314 +WHERE ST_Area(ST_Buffer(the_geog,10)) NOT between 307 and 315 LIMIT 10; SELECT '#304.a', Count(*) FROM utm_dots WHERE ST_DWithin(the_geog, 'POINT(0 0)'::geography, 3000000); diff --git a/regress/tickets_expected b/regress/tickets_expected index 4ce2d2b3f..206b6f11c 100644 --- a/regress/tickets_expected +++ b/regress/tickets_expected @@ -64,6 +64,7 @@ NOTICE: No points or linestrings in input array #277|1,1e+308 #299|2 #304 +ERROR: GetProj4StringSPI: Cannot find SRID (32704) in spatial_ref_sys #304.a|21 #304.b|1 #408|Too few points in geometry component[