[SCM] PostGIS branch master updated. 3.4.0rc1-998-gb7990ecf8
git at osgeo.org
git at osgeo.org
Fri Mar 8 13:14:09 PST 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, master has been updated
via b7990ecf8b6981f4a20e7121d0e75222075df6e9 (commit)
from 7bd482a696fd4ac419bf76f4483143112b0396d5 (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 b7990ecf8b6981f4a20e7121d0e75222075df6e9
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Fri Mar 8 13:13:38 2024 -0800
Calculate sphere area using the spheroid code,
but with a spherical initialization. Remove
no longer used code. References #5671
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 982e4e4af..689d4c070 100644
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -1805,39 +1805,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 +1996,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 dc8bde060..2181d3c3c 100644
--- a/postgis/geography_measurement.c
+++ b/postgis/geography_measurement.c
@@ -539,10 +539,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:
liblwgeom/cunit/cu_geodetic.c | 11 ++---
liblwgeom/lwgeodetic.c | 96 +++--------------------------------------
liblwgeom/lwgeodetic.h | 1 -
postgis/geography_measurement.c | 5 +--
4 files changed, 12 insertions(+), 101 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list