[SCM] PostGIS branch stable-3.3 updated. 3.3.6-12-g96145b118
git at osgeo.org
git at osgeo.org
Thu Mar 14 18:21:19 PDT 2024
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".
The branch, stable-3.3 has been updated
via 96145b118098e1e4723438c4bafef11826c7ec87 (commit)
from 237be25f4e3d2f51816633b61f43552f4bdb3677 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 96145b118098e1e4723438c4bafef11826c7ec87
Author: Regina Obe <lr at pcorp.us>
Date: Thu Mar 14 18:48:18 2024 -0400
Use geographiclib also for sphere calcs.
Closes #5671 for PostGIS 3.3.7
diff --git a/NEWS b/NEWS
index ead91bc60..2ec6557bf 100644
--- a/NEWS
+++ b/NEWS
@@ -11,6 +11,8 @@ xxxx/xx/xx
set to off (Sandro Santilli)
- #5589, ST_3DDistance error for shared first point (Paul Ramsey)
- #5686, ST_NumInteriorRings and Triangle crash (Paul Ramsey)
+ - #5671, Bug in ST_Area function with use_spheroid=false
+ (Paul Ramsey, Regina Obe)
PostGIS 3.3.6
diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c
index 12506476d..5e89aef3b 100644
--- a/liblwgeom/cunit/cu_geodetic.c
+++ b/liblwgeom/cunit/cu_geodetic.c
@@ -1564,8 +1564,9 @@ static void test_lwgeom_area_sphere(void)
double area;
SPHEROID s;
- /* Init to WGS84 */
- spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
+ /* Init to Sphere */
+ spheroid_init(&s, WGS84_RADIUS, WGS84_RADIUS);
+ s.a = s.b = s.radius;
/* Simple case */
lwg = lwgeom_from_wkt("POLYGON((1 1, 1 2, 2 2, 2 1, 1 1))", LW_PARSER_CHECK_NONE);
@@ -1577,17 +1578,17 @@ static void test_lwgeom_area_sphere(void)
/* Robustness tests, from ticket #3393 */
lwg = lwgeom_from_wkt("POLYGON((0 78.703946026663,0 0,179.999997913235 0,179.999997913235 -33.0888306884702,0 78.703946026663))", LW_PARSER_CHECK_NONE);
area = lwgeom_area_sphere(lwg, &s);
- CU_ASSERT_DOUBLE_EQUAL(area, 127516467322130, 1.0);
+ CU_ASSERT_DOUBLE_EQUAL(area, 127516466960671, 1.0);
lwgeom_free(lwg);
lwg = lwgeom_from_wkt("POLYGON((0 78.703946026662,0 0,179.999997913235 0,179.999997913235 -33.0888306884702,0 78.703946026662))", LW_PARSER_CHECK_NONE);
area = lwgeom_area_sphere(lwg, &s);
- CU_ASSERT_DOUBLE_EQUAL(area, 127516467322130, 1.0);
+ CU_ASSERT_DOUBLE_EQUAL(area, 127516466960671, 1.0);
lwgeom_free(lwg);
lwg = lwgeom_from_wkt("POLYGON((0 78.703946026664,0 0,179.999997913235 0,179.999997913235 -33.0888306884702,0 78.703946026664))", LW_PARSER_CHECK_NONE);
area = lwgeom_area_sphere(lwg, &s);
- CU_ASSERT_DOUBLE_EQUAL(area, 127516467322130, 1.0);
+ CU_ASSERT_DOUBLE_EQUAL(area, 127516466960671, 1.0);
lwgeom_free(lwg);
/* end #3393 */
}
diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index c2b2f7833..d739a9b44 100644
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -714,58 +714,6 @@ edge_point_side(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
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
@@ -1805,39 +1753,6 @@ 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)
-{
- uint32_t i;
- const POINT2D *p;
- GEOGRAPHIC_POINT a, b, c;
- double area = 0.0;
-
- /* Return zero on nonsensical inputs */
- if ( ! pa || pa->npoints < 4 )
- return 0.0;
-
- 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++ )
- {
- 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);
-}
-
-
static double ptarray_distance_spheroid(const POINTARRAY *pa1, const POINTARRAY *pa2, const SPHEROID *s, double tolerance, int check_intersection)
{
GEOGRAPHIC_EDGE e1, e2;
@@ -2029,67 +1944,14 @@ static double ptarray_distance_spheroid(const POINTARRAY *pa1, const POINTARRAY
/**
-* Calculate the area of an LWGEOM. Anything except POLYGON, MULTIPOLYGON
-* and GEOMETRYCOLLECTION return zero immediately. Multi's recurse, polygons
-* calculate external ring area and subtract internal ring area. A GBOX is
-* required to calculate an outside point.
+* Delegate to the spheroid function with a spherically
+* parameterized spheroid.
*/
double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
{
- int type;
- double radius2 = spheroid->radius * spheroid->radius;
-
- assert(lwgeom);
-
- /* No area in nothing */
- if ( lwgeom_is_empty(lwgeom) )
- return 0.0;
-
- /* Read the geometry type number */
- type = lwgeom->type;
-
- /* Anything but polygons and collections returns zero */
- if ( ! ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) )
- return 0.0;
-
- /* Actually calculate area */
- if ( type == POLYGONTYPE )
- {
- LWPOLY *poly = (LWPOLY*)lwgeom;
- uint32_t i;
- double area = 0.0;
-
- /* Just in case there's no rings */
- if ( poly->nrings < 1 )
- return 0.0;
-
- /* First, the area of the outer ring */
- 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]);
- }
- return area;
- }
-
- /* Recurse into sub-geometries to get area */
- if ( type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE )
- {
- LWCOLLECTION *col = (LWCOLLECTION*)lwgeom;
- uint32_t i;
- double area = 0.0;
-
- for ( i = 0; i < col->ngeoms; i++ )
- {
- area += lwgeom_area_sphere(col->geoms[i], spheroid);
- }
- return area;
- }
-
- /* Shouldn't get here. */
- return 0.0;
+ SPHEROID s;
+ spheroid_init(&s, WGS84_RADIUS, WGS84_RADIUS);
+ return lwgeom_area_spheroid(lwgeom, &s);
}
diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h
index 6364db210..39053ef5f 100644
--- a/liblwgeom/lwgeodetic.h
+++ b/liblwgeom/lwgeodetic.h
@@ -128,7 +128,6 @@ int lwline_covers_lwpoint(const LWLINE* lwline, const LWPOINT* lwpoint);
int lwpoly_intersects_line(const LWPOLY* lwpoly, const POINTARRAY* line);
int lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside);
int ptarray_point_in_ring(const POINTARRAY *pa, const POINT2D *pt_outside, const POINT2D *pt_to_test);
-double ptarray_area_sphere(const POINTARRAY *pa);
double latitude_degrees_normalize(double lat);
double longitude_degrees_normalize(double lon);
double ptarray_length_spheroid(const POINTARRAY *pa, const SPHEROID *s);
diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c
index 1911472c3..9a941888c 100644
--- a/postgis/geography_measurement.c
+++ b/postgis/geography_measurement.c
@@ -540,10 +540,7 @@ Datum geography_area(PG_FUNCTION_ARGS)
s.a = s.b = s.radius;
/* Calculate the area */
- if ( use_spheroid )
- area = lwgeom_area_spheroid(lwgeom, &s);
- else
- area = lwgeom_area_sphere(lwgeom, &s);
+ area = lwgeom_area_spheroid(lwgeom, &s);
/* Clean up */
lwgeom_free(lwgeom);
-----------------------------------------------------------------------
Summary of changes:
NEWS | 2 +
liblwgeom/cunit/cu_geodetic.c | 11 +--
liblwgeom/lwgeodetic.c | 148 ++--------------------------------------
liblwgeom/lwgeodetic.h | 1 -
postgis/geography_measurement.c | 5 +-
5 files changed, 14 insertions(+), 153 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list