From 6dd39806d5b759d9ab8d8a11e7bbbf43efd19457 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Sat, 31 Oct 2009 00:01:45 +0000 Subject: [PATCH] Add ST_Length() implementation on spheroid and rationalize the sphere/spheroid implementations into a smaller shared set of functions. git-svn-id: http://svn.osgeo.org/postgis/trunk@4710 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/cunit/cu_geodetic.c | 48 ++-- liblwgeom/cunit/cu_geodetic.h | 2 +- liblwgeom/libgeom.h | 14 +- liblwgeom/liblwgeom.h | 1 + liblwgeom/lwgeodetic.c | 108 ++++++--- liblwgeom/lwgeodetic.h | 6 +- liblwgeom/lwspheroid.c | 386 ++------------------------------ liblwgeom/lwspheroid.h | 17 -- postgis/geography.sql.in.c | 27 ++- postgis/geography_measurement.c | 97 ++++---- 10 files changed, 222 insertions(+), 484 deletions(-) delete mode 100644 liblwgeom/lwspheroid.h diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c index 1c78ed177..5b7881450 100644 --- a/liblwgeom/cunit/cu_geodetic.c +++ b/liblwgeom/cunit/cu_geodetic.c @@ -480,6 +480,11 @@ void test_lwgeom_distance_sphere(void) LWGEOM *lwg1, *lwg2; GBOX gbox1, gbox2; double d; + SPHEROID s; + + /* Init and force spherical */ + spheroid_init(&s, 6378137.0, 6356752.314245179498); + s.a = s.b = s.radius; gbox1.flags = gflags(0, 0, 1); gbox2.flags = gflags(0, 0, 1); @@ -489,8 +494,8 @@ void test_lwgeom_distance_sphere(void) lwg2 = lwgeom_from_ewkt("LINESTRING(-10 -5, -5 0, 5 0, 10 -5)", PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg1, &gbox1); lwgeom_calculate_gbox_geodetic(lwg2, &gbox2); - d = lwgeom_distance_sphere(lwg1, lwg2, gbox1, gbox2, 0.0); - CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001); + d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0); + CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 180.0, 0.00001); lwgeom_free(lwg1); lwgeom_free(lwg2); @@ -499,7 +504,7 @@ void test_lwgeom_distance_sphere(void) lwg2 = lwgeom_from_ewkt("LINESTRING(-10 -5, -5 20, 5 0, 10 -5)", PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg1, &gbox1); lwgeom_calculate_gbox_geodetic(lwg2, &gbox2); - d = lwgeom_distance_sphere(lwg1, lwg2, gbox1, gbox2, 0.0); + d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0); CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001); lwgeom_free(lwg1); lwgeom_free(lwg2); @@ -509,8 +514,8 @@ void test_lwgeom_distance_sphere(void) lwg2 = lwgeom_from_ewkt("LINESTRING(-10 -5, -5 0, 5 0, 10 -5)", PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg1, &gbox1); lwgeom_calculate_gbox_geodetic(lwg2, &gbox2); - d = lwgeom_distance_sphere(lwg1, lwg2, gbox1, gbox2, 0.0); - CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001); + d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0); + CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 180.0, 0.00001); lwgeom_free(lwg1); lwgeom_free(lwg2); @@ -518,8 +523,8 @@ void test_lwgeom_distance_sphere(void) lwg2 = lwgeom_from_ewkt("POINT(-4 -1)", PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg1, &gbox1); lwgeom_calculate_gbox_geodetic(lwg2, &gbox2); - d = lwgeom_distance_sphere(lwg1, lwg2, gbox1, gbox2, 0.0); - CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 90.0, 0.00001); + d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0); + CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 90.0, 0.00001); lwgeom_free(lwg1); lwgeom_free(lwg2); @@ -528,7 +533,7 @@ void test_lwgeom_distance_sphere(void) lwg2 = lwgeom_from_ewkt("POINT(-1 -1)", PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg1, &gbox1); lwgeom_calculate_gbox_geodetic(lwg2, &gbox2); - d = lwgeom_distance_sphere(lwg1, lwg2, gbox1, gbox2, 0.0); + d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0); CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001); lwgeom_free(lwg1); lwgeom_free(lwg2); @@ -538,8 +543,8 @@ void test_lwgeom_distance_sphere(void) lwg2 = lwgeom_from_ewkt("POINT(-1 -1)", PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg1, &gbox1); lwgeom_calculate_gbox_geodetic(lwg2, &gbox2); - d = lwgeom_distance_sphere(lwg1, lwg2, gbox1, gbox2, 0.0); - CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001); + d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0); + CU_ASSERT_DOUBLE_EQUAL(d, 111178.142466, 0.1); lwgeom_free(lwg1); lwgeom_free(lwg2); @@ -548,7 +553,7 @@ void test_lwgeom_distance_sphere(void) lwg2 = lwgeom_from_ewkt("POINT(2 2)", PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg1, &gbox1); lwgeom_calculate_gbox_geodetic(lwg2, &gbox2); - d = lwgeom_distance_sphere(lwg1, lwg2, gbox1, gbox2, 0.0); + d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0); CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001); lwgeom_free(lwg1); lwgeom_free(lwg2); @@ -558,8 +563,8 @@ void test_lwgeom_distance_sphere(void) lwg2 = lwgeom_from_ewkt("0106000020E61000000100000001030000000100000007000000280EC3FB8CCA5EC0A5CDC747233C45402787C8F58CCA5EC0659EA2761E3C45400CED58DF8FCA5EC0C37FAE6E1E3C4540AE97B8E08FCA5EC00346F58B1F3C4540250359FD8ECA5EC05460628E1F3C45403738F4018FCA5EC05DC84042233C4540280EC3FB8CCA5EC0A5CDC747233C4540", PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg1, &gbox1); lwgeom_calculate_gbox_geodetic(lwg2, &gbox2); - d = lwgeom_distance_sphere(lwg1, lwg2, gbox1, gbox2, 0.0); - CU_ASSERT_DOUBLE_EQUAL(d * 6371009.0, 23630.8003, 0.1); + d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0); + CU_ASSERT_DOUBLE_EQUAL(d, 23630.8003, 0.1); lwgeom_free(lwg1); lwgeom_free(lwg2); @@ -567,40 +572,41 @@ void test_lwgeom_distance_sphere(void) void test_spheroid_distance(void) { - /* WGS 84 values */ - static double a = 6378137.0; - static double b = 6356752.314245179498; GEOGRAPHIC_POINT g1, g2; double d; + SPHEROID s; + + /* Init to WGS84 */ + spheroid_init(&s, 6378137.0, 6356752.314245179498); /* One vertical degree */ point_set(0.0, 0.0, &g1); point_set(0.0, 1.0, &g2); - d = spheroid_distance(g1, g2, a, b); + d = spheroid_distance(g1, g2, s); CU_ASSERT_DOUBLE_EQUAL(d, 110574.388615329, 0.001); /* Ten horizontal degrees */ point_set(-10.0, 0.0, &g1); point_set(0.0, 0.0, &g2); - d = spheroid_distance(g1, g2, a, b); + d = spheroid_distance(g1, g2, s); CU_ASSERT_DOUBLE_EQUAL(d, 1113194.90793274, 0.001); /* One horizonal degree */ point_set(-1.0, 0.0, &g1); point_set(0.0, 0.0, &g2); - d = spheroid_distance(g1, g2, a, b); + d = spheroid_distance(g1, g2, s); CU_ASSERT_DOUBLE_EQUAL(d, 111319.490779, 0.001); /* Around world w/ slight bend */ point_set(-180.0, 0.0, &g1); point_set(0.0, 1.0, &g2); - d = spheroid_distance(g1, g2, a, b); + d = spheroid_distance(g1, g2, s); CU_ASSERT_DOUBLE_EQUAL(d, 19893357.0704483, 0.001); /* Up to pole */ point_set(-180.0, 0.0, &g1); point_set(0.0, 90.0, &g2); - d = spheroid_distance(g1, g2, a, b); + d = spheroid_distance(g1, g2, s); CU_ASSERT_DOUBLE_EQUAL(d, 10001965.7295318, 0.001); } \ No newline at end of file diff --git a/liblwgeom/cunit/cu_geodetic.h b/liblwgeom/cunit/cu_geodetic.h index 7f055774c..102957a25 100644 --- a/liblwgeom/cunit/cu_geodetic.h +++ b/liblwgeom/cunit/cu_geodetic.h @@ -15,7 +15,7 @@ #include #include "CUnit/Basic.h" -#include "lwspheroid.h" +#include "lwgeodetic.h" #include "cu_tester.h" /*********************************************************************** diff --git a/liblwgeom/libgeom.h b/liblwgeom/libgeom.h index ef2e06020..43618d190 100644 --- a/liblwgeom/libgeom.h +++ b/liblwgeom/libgeom.h @@ -371,28 +371,28 @@ extern int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox); extern int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox); /** -* Calculate the geodetic distance from lwgeom1 to lwgeom2 on the unit sphere. -* Pass in a tolerance in radians. +* Initialize a spheroid object for use in geodetic functions. */ -extern double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, double tolerance); +extern void spheroid_init(SPHEROID *s, double a, double b); /** * Calculate the geodetic distance from lwgeom1 to lwgeom2 on the spheroid. -* Pass in a major axis, minor axis and tolerance using the same units for each (meters, generally). +* A spheroid with major axis == minor axis will be treated as a sphere. +* Pass in a tolerance in spheroid units. */ -extern double lwgeom_distance_spheroid(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, double a, double b, double tolerance); +extern double lwgeom_distance_spheroid(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, SPHEROID spheroid, double tolerance); /** * Calculate the geodetic area of a lwgeom on the unit sphere. The result * will have to by multiplied by the real radius^2 to get the real area. */ -extern double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox); +extern double lwgeom_area_spheroid(LWGEOM *lwgeom, GBOX gbox, SPHEROID s); /** * Calculate the geodetic length of a lwgeom on the unit sphere. The result * will have to by multiplied by the real radius to get the real length. */ -extern double lwgeom_length_sphere(LWGEOM *geom); +extern double lwgeom_length_spheroid(LWGEOM *geom, SPHEROID s); /** * Calculate covers predicate for two lwgeoms on the sphere. Currently diff --git a/liblwgeom/liblwgeom.h b/liblwgeom/liblwgeom.h index a69fc2bf6..bdec43cdd 100644 --- a/liblwgeom/liblwgeom.h +++ b/liblwgeom/liblwgeom.h @@ -251,6 +251,7 @@ typedef struct double f; /* flattening */ double e; /* eccentricity (first) */ double e_sq; /* eccentricity (first), squared */ + double radius; /* (2*a+b)/3 spherical average radius */ char name[20]; /* name of ellipse */ } SPHEROID; diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index 96d15ad3e..d80b6e358 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -1463,10 +1463,11 @@ int ptarray_point_in_ring(POINTARRAY *pa, POINT2D pt_outside, POINT2D pt_to_test } -static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double tolerance, int check_intersection) +static double ptarray_distance_spheroid(POINTARRAY *pa1, POINTARRAY *pa2, SPHEROID s, double tolerance, int check_intersection) { GEOGRAPHIC_EDGE e1, e2; GEOGRAPHIC_POINT g1, g2; + GEOGRAPHIC_POINT nearest1, nearest2; POINT2D p; double distance; int i, j; @@ -1485,7 +1486,12 @@ static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double t geographic_point_init(p.x, p.y, &g1); getPoint2d_p(pa2, 0, &p); geographic_point_init(p.x, p.y, &g2); - return sphere_distance(g1, g2); + /* Sphere special case, axes equal */ + if( s.a == s.b ) + distance = s.radius * sphere_distance(g1, g2); + else + distance = spheroid_distance(g1, g2, s); + return distance; } /* Handle point/line case here */ @@ -1520,14 +1526,39 @@ static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double t double d; getPoint2d_p(pa_many, i, &p); geographic_point_init(p.x, p.y, &(e1.end)); - d = edge_distance_to_point(e1, g1, 0); + /* Get the spherical distance between point and edge */ + d = s.radius * edge_distance_to_point(e1, g1, &g2); + /* New shortest distance! Record this distance / location */ if( d < distance ) + { distance = d; + nearest2 = g2; + } + /* We've gotten closer than the tolerance... */ if( d < tolerance ) - return distance; + { + /* Working on a sphere? The answer is correct, return */ + if( s.a == s.b ) + { + return distance; + } + /* On a spheroid? Confirm that we are *actually* closer than tolerance */ + else + { + d = spheroid_distance(g1, nearest2, s); + /* Yes, closer than tolerance, return! */ + if( d < tolerance ) + return d; + } + } e1.start = e1.end; } - return distance; + /* On sphere, return answer */ + if( s.a == s.b ) + return distance; + /* On spheroid, calculate final answer based on closest approach */ + else + return spheroid_distance(g1, nearest2, s); } /* Initialize start of line 1 */ @@ -1548,7 +1579,6 @@ static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double t for( j = 1; j < pa2->npoints; j++ ) { double d; - GEOGRAPHIC_POINT g; getPoint2d_p(pa2, j, &p); geographic_point_init(p.x, p.y, &(e2.end)); @@ -1558,18 +1588,33 @@ static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double t LWDEBUGF(4, "e2.start == GPOINT(%.6g %.6g) ", e2.start.lat, e2.start.lon); LWDEBUGF(4, "e2.end == GPOINT(%.6g %.6g) ", e2.end.lat, e2.end.lon); - if ( check_intersection && edge_intersection(e1, e2, &g) ) + if ( check_intersection && edge_intersection(e1, e2, &g1) ) { LWDEBUG(4,"edge intersection! returning 0.0"); return 0.0; } - d = edge_distance_to_edge(e1, e2, 0, 0); + d = s.radius * edge_distance_to_edge(e1, e2, &g1, &g2); LWDEBUGF(4,"got edge_distance_to_edge %.8g", d); if( d < distance ) + { distance = d; + nearest1 = g1; + nearest2 = g2; + } if( d < tolerance ) - return distance; + { + if( s.a == s.b ) + { + return d; + } + else + { + d = spheroid_distance(nearest1, nearest2, s); + if( d < tolerance ) + return d; + } + } /* Copy end to start to allow a new end value in next iteration */ e2.start = e2.end; @@ -1581,7 +1626,10 @@ static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double t } LWDEBUGF(4,"finished all loops, returning %.8g", distance); - return distance; + if( s.a == s.b ) + return distance; + else + return spheroid_distance(nearest1, nearest2, s); } @@ -1591,10 +1639,11 @@ static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double t * calculate external ring area and subtract internal ring area. A GBOX is * required to calculate an outside point. */ -double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox) +double lwgeom_area_spheroid(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid) { int type; POINT2D pt_outside; + double radius2 = spheroid.radius * spheroid.radius; assert(lwgeom); @@ -1625,12 +1674,12 @@ double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox) return 0.0; /* First, the area of the outer ring */ - area += ptarray_area_sphere(poly->rings[0], pt_outside); + area += radius2 * ptarray_area_sphere(poly->rings[0], pt_outside); /* Subtract areas of inner rings */ for( i = 1; i < poly->nrings; i++ ) { - area -= ptarray_area_sphere(poly->rings[i], pt_outside); + area -= radius2 * ptarray_area_sphere(poly->rings[i], pt_outside); } return area; } @@ -1644,7 +1693,7 @@ double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox) for( i = 0; i < col->ngeoms; i++ ) { - area += lwgeom_area_sphere(col->geoms[i], gbox); + area += lwgeom_area_spheroid(col->geoms[i], gbox, spheroid); } return area; } @@ -1659,7 +1708,7 @@ double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox) * longitude and latitude. Return immediately when the calulated distance drops * below the tolerance (useful for dwithin calculations). */ -double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, double tolerance) +double lwgeom_distance_spheroid(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, SPHEROID spheroid, double tolerance) { int type1, type2; int check_intersection = LW_FALSE; @@ -1699,7 +1748,7 @@ double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX else pa2 = ((LWLINE*)lwgeom2)->points; - return ptarray_distance_sphere(pa1, pa2, tolerance, check_intersection); + return ptarray_distance_spheroid(pa1, pa2, spheroid, tolerance, check_intersection); } /* Point/Polygon cases, if point-in-poly, return zero, else return distance. */ @@ -1734,7 +1783,7 @@ double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX /* Not inside, so what's the actual distance? */ for( i = 0; i < lwpoly->nrings; i++ ) { - double ring_distance = ptarray_distance_sphere(lwpoly->rings[i], lwpt->point, tolerance, check_intersection); + double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwpt->point, spheroid, tolerance, check_intersection); if( ring_distance < distance ) distance = ring_distance; if( distance < tolerance ) @@ -1779,7 +1828,7 @@ double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX /* Not contained, so what's the actual distance? */ for( i = 0; i < lwpoly->nrings; i++ ) { - double ring_distance = ptarray_distance_sphere(lwpoly->rings[i], lwline->points, tolerance, check_intersection); + double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwline->points, spheroid, tolerance, check_intersection); LWDEBUGF(4, "ring[%d] ring_distance = %.8g", i, ring_distance); if( ring_distance < distance ) distance = ring_distance; @@ -1816,7 +1865,7 @@ double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX { for( j = 0; j < lwpoly2->nrings; j++ ) { - double ring_distance = ptarray_distance_sphere(lwpoly1->rings[i], lwpoly2->rings[j], tolerance, check_intersection); + double ring_distance = ptarray_distance_spheroid(lwpoly1->rings[i], lwpoly2->rings[j], spheroid, tolerance, check_intersection); if( ring_distance < distance ) distance = ring_distance; if( distance < tolerance ) @@ -1835,7 +1884,7 @@ double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX for( i = 0; i < col->ngeoms; i++ ) { - double geom_distance = lwgeom_distance_sphere(col->geoms[i], lwgeom2, gbox1, gbox2, tolerance); + double geom_distance = lwgeom_distance_spheroid(col->geoms[i], lwgeom2, gbox1, gbox2, spheroid, tolerance); if( geom_distance < distance ) distance = geom_distance; if( distance < tolerance ) @@ -1853,7 +1902,7 @@ double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX for( i = 0; i < col->ngeoms; i++ ) { - double geom_distance = lwgeom_distance_sphere(lwgeom1, col->geoms[i], gbox1, gbox2, tolerance); + double geom_distance = lwgeom_distance_spheroid(lwgeom1, col->geoms[i], gbox1, gbox2, spheroid, tolerance); if( geom_distance < distance ) distance = geom_distance; if( distance < tolerance ) @@ -2270,7 +2319,7 @@ int lwgeom_check_geodetic(const LWGEOM *geom) return LW_FALSE; } -double ptarray_length_sphere(POINTARRAY *pa) +double ptarray_length_spheroid(POINTARRAY *pa, SPHEROID s) { GEOGRAPHIC_POINT a, b; POINT2D p; @@ -2290,7 +2339,12 @@ double ptarray_length_sphere(POINTARRAY *pa) getPoint2d_p(pa, i, &p); geographic_point_init(p.x, p.y, &b); - length += sphere_distance(a, b); + /* Special sphere case */ + if( s.a == s.b ) + length += s.radius * sphere_distance(a, b); + /* Spheroid case */ + else + length += spheroid_distance(a, b, s); /* B gets incremented in the next loop, so we save the value here */ a = b; @@ -2298,7 +2352,7 @@ double ptarray_length_sphere(POINTARRAY *pa) return length; } -double lwgeom_length_sphere(LWGEOM *geom) +double lwgeom_length_spheroid(LWGEOM *geom, SPHEROID s) { int type; int i = 0; @@ -2316,14 +2370,14 @@ double lwgeom_length_sphere(LWGEOM *geom) return 0.0; if ( type == LINETYPE ) - return ptarray_length_sphere(((LWLINE*)geom)->points); + return ptarray_length_spheroid(((LWLINE*)geom)->points, s); if ( type == POLYGONTYPE ) { LWPOLY *poly = (LWPOLY*)geom; for( i = 0; i < poly->nrings; i++ ) { - length += ptarray_length_sphere(poly->rings[i]); + length += ptarray_length_spheroid(poly->rings[i], s); } return length; } @@ -2334,7 +2388,7 @@ double lwgeom_length_sphere(LWGEOM *geom) for( i = 0; i < col->ngeoms; i++ ) { - length += lwgeom_length_sphere(col->geoms[i]); + length += lwgeom_length_spheroid(col->geoms[i], s); } return length; } diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h index 8facea5d0..5738642c1 100644 --- a/liblwgeom/lwgeodetic.h +++ b/liblwgeom/lwgeodetic.h @@ -77,5 +77,9 @@ int ptarray_point_in_ring(POINTARRAY *pa, POINT2D pt_outside, POINT2D pt_to_test double ptarray_area_sphere(POINTARRAY *pa, POINT2D pt_outside); double latitude_degrees_normalize(double lat); double longitude_degrees_normalize(double lon); -double ptarray_length_sphere(POINTARRAY *pa); +double ptarray_length_spheroid(POINTARRAY *pa, SPHEROID s); int geographic_point_equals(GEOGRAPHIC_POINT g1, GEOGRAPHIC_POINT g2); +/* +** Prototypes for spheroid functions. +*/ +double spheroid_distance(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, SPHEROID spheroid); diff --git a/liblwgeom/lwspheroid.c b/liblwgeom/lwspheroid.c index 49840724e..8c3eaaf21 100644 --- a/liblwgeom/lwspheroid.c +++ b/liblwgeom/lwspheroid.c @@ -1,11 +1,20 @@ -#include "lwspheroid.h" +#include "lwgeodetic.h" #define POW2(x) ((x)*(x)) -static double spheroid_mu2(double alpha, double a, double b) +void spheroid_init(SPHEROID *s, double a, double b) { - double b2 = POW2(b); - return POW2(cos(alpha)) * (POW2(a) - b2) / b2; + s->a = a; + s->b = b; + s->f = (a - b) / a; + s->e_sq = (a*a - b*b)/(a*a); + s->radius = (2.0 * a + b ) / 3.0; +} + +static double spheroid_mu2(double alpha, SPHEROID s) +{ + double b2 = POW2(s.b); + return POW2(cos(alpha)) * (POW2(s.a) - b2) / b2; } static double spheroid_big_a(double u2) @@ -18,11 +27,11 @@ static double spheroid_big_b(double u2) return (u2 / 1024.0) * (256.0 + u2 * (-128.0 + u2 * (74.0 - 47.0 * u2))); } -double spheroid_distance(GEOGRAPHIC_POINT r, GEOGRAPHIC_POINT s, double a, double b) +double spheroid_distance(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, SPHEROID spheroid) { - double lambda = (s.lon - r.lon); - double f = (a - b)/a; - double omf = 1 - f; + double lambda = (b.lon - a.lon); + double f = spheroid.f; + double omf = 1 - spheroid.f; double u1, u2; double cos_u1, cos_u2; double sin_u1, sin_u2; @@ -34,15 +43,15 @@ double spheroid_distance(GEOGRAPHIC_POINT r, GEOGRAPHIC_POINT s, double a, doubl int i = 0; /* Same point => zero distance */ - if( geographic_point_equals(r, s) ) + if( geographic_point_equals(a, b) ) { return 0.0; } - u1 = atan(omf * tan(r.lat)); + u1 = atan(omf * tan(a.lat)); cos_u1 = cos(u1); sin_u1 = sin(u1); - u2 = atan(omf * tan(s.lat)); + u2 = atan(omf * tan(b.lat)); cos_u2 = cos(u2); sin_u2 = sin(u2); @@ -83,368 +92,21 @@ double spheroid_distance(GEOGRAPHIC_POINT r, GEOGRAPHIC_POINT s, double a, doubl } while ( i < 999 && lambda != 0.0 && fabs((last_lambda - lambda)/lambda) > 1.0e-9 ); - u2 = spheroid_mu2(alpha, a, b); + u2 = spheroid_mu2(alpha, spheroid); big_a = spheroid_big_a(u2); big_b = spheroid_big_b(u2); delta_sigma = big_b * sin_sigma * (cos2_sigma_m + (big_b / 4.0) * (cos_sigma * (-1.0 + 2.0 * POW2(cos2_sigma_m)) - (big_b / 6.0) * cos2_sigma_m * (-3.0 + 4.0 * sqrsin_sigma) * (-3.0 + 4.0 * POW2(cos2_sigma_m)))); - distance = b * big_a * (sigma - delta_sigma); + distance = spheroid.b * big_a * (sigma - delta_sigma); /* Algorithm failure, distance == NaN, fallback to sphere */ if( distance != distance ) { - lwerror("spheroid_distance returned NaN: (%.20g %.20g) (%.20g %.20g) a = %.20g b = %.20g",r.lat, r.lon, s.lat, s.lon, a, b); - return (2.0 * a + b) * sphere_distance(r, s) / 3.0; + lwerror("spheroid_distance returned NaN: (%.20g %.20g) (%.20g %.20g) a = %.20g b = %.20g",a.lat, a.lon, b.lat, b.lon, spheroid.a, spheroid.b); + return spheroid.radius * sphere_distance(a, b); } return distance; } -double ptarray_distance_spheroid(POINTARRAY *pa1, POINTARRAY *pa2, double a, double b, double tolerance, int check_intersection) -{ - GEOGRAPHIC_EDGE e1, e2; - GEOGRAPHIC_POINT g1, g2; - GEOGRAPHIC_POINT nearest1, nearest2; - POINT2D p; - double radius = (2.0*a + b)/3.0; - double distance; - int i, j; - - /* Make result really big, so that everything will be smaller than it */ - distance = MAXFLOAT; - - /* Empty point arrays? Return negative */ - if ( pa1->npoints == 0 || pa1->npoints == 0 ) - return -1.0; - - /* Handle point/point case here */ - if ( pa1->npoints == 1 && pa2->npoints == 1 ) - { - getPoint2d_p(pa1, 0, &p); - geographic_point_init(p.x, p.y, &g1); - getPoint2d_p(pa2, 0, &p); - geographic_point_init(p.x, p.y, &g2); - return spheroid_distance(g1, g2, a, b); - } - - /* Handle point/line case here */ - if ( pa1->npoints == 1 || pa2->npoints == 1 ) - { - /* Handle one/many case here */ - int i; - POINTARRAY *pa_one, *pa_many; - - if( pa1->npoints == 1 ) - { - pa_one = pa1; - pa_many = pa2; - } - else - { - pa_one = pa2; - pa_many = pa1; - } - - /* Initialize our point */ - getPoint2d_p(pa_one, 0, &p); - geographic_point_init(p.x, p.y, &g1); - - /* Initialize start of line */ - getPoint2d_p(pa_many, 0, &p); - geographic_point_init(p.x, p.y, &(e1.start)); - - /* Iterate through the edges in our line */ - for( i = 1; i < pa_many->npoints; i++ ) - { - double d; - getPoint2d_p(pa_many, i, &p); - geographic_point_init(p.x, p.y, &(e1.end)); - d = radius * edge_distance_to_point(e1, g1, &g2); - if( d < distance ) - { - distance = d; - nearest2 = g2; - } - if( d < tolerance ) - { - double sd = spheroid_distance(g1, nearest2, a, b); - if( sd < tolerance ) - return sd; - } - e1.start = e1.end; - } - return spheroid_distance(g1, nearest2, a, b); - } - - /* Initialize start of line 1 */ - getPoint2d_p(pa1, 0, &p); - geographic_point_init(p.x, p.y, &(e1.start)); - - - /* Handle line/line case */ - for( i = 1; i < pa1->npoints; i++ ) - { - getPoint2d_p(pa1, i, &p); - geographic_point_init(p.x, p.y, &(e1.end)); - - /* Initialize start of line 2 */ - getPoint2d_p(pa2, 0, &p); - geographic_point_init(p.x, p.y, &(e2.start)); - - for( j = 1; j < pa2->npoints; j++ ) - { - double d; - GEOGRAPHIC_POINT g; - - getPoint2d_p(pa2, j, &p); - geographic_point_init(p.x, p.y, &(e2.end)); - - LWDEBUGF(4, "e1.start == GPOINT(%.6g %.6g) ", e1.start.lat, e1.start.lon); - LWDEBUGF(4, "e1.end == GPOINT(%.6g %.6g) ", e1.end.lat, e1.end.lon); - LWDEBUGF(4, "e2.start == GPOINT(%.6g %.6g) ", e2.start.lat, e2.start.lon); - LWDEBUGF(4, "e2.end == GPOINT(%.6g %.6g) ", e2.end.lat, e2.end.lon); - - if ( check_intersection && edge_intersection(e1, e2, &g) ) - { - LWDEBUG(4,"edge intersection! returning 0.0"); - return 0.0; - } - d = radius * edge_distance_to_edge(e1, e2, &g1, &g2); - LWDEBUGF(4,"got edge_distance_to_edge %.8g", d); - - if( d < distance ) - { - distance = d; - nearest1 = g1; - nearest2 = g2; - } - if( d < tolerance ) - { - double sd = spheroid_distance(nearest1, nearest2, a, b); - if( sd < tolerance ) - return sd; - } - - /* Copy end to start to allow a new end value in next iteration */ - e2.start = e2.end; - } - - /* Copy end to start to allow a new end value in next iteration */ - e1.start = e1.end; - - } - - return spheroid_distance(nearest1, nearest2, a, b); -} - - -double lwgeom_distance_spheroid(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, double a, double b, double tolerance) -{ - int type1, type2; - int check_intersection = LW_FALSE; - - assert(lwgeom1); - assert(lwgeom2); - - LWDEBUGF(4, "entered function, tolerance %.8g", tolerance); - - /* What's the distance to an empty geometry? We don't know. */ - if( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) ) - { - return 0.0; - } - - type1 = TYPE_GETTYPE(lwgeom1->type); - type2 = TYPE_GETTYPE(lwgeom2->type); - - - /* If the boxes aren't disjoint, we have to check for edge intersections */ - if( gbox_overlaps(gbox1, gbox2) ) - check_intersection = LW_TRUE; - - /* Point/line combinations can all be handled with simple point array iterations */ - if( ( type1 == POINTTYPE || type1 == LINETYPE ) && - ( type2 == POINTTYPE || type2 == LINETYPE ) ) - { - POINTARRAY *pa1, *pa2; - - if( type1 == POINTTYPE ) - pa1 = ((LWPOINT*)lwgeom1)->point; - else - pa1 = ((LWLINE*)lwgeom1)->points; - - if( type2 == POINTTYPE ) - pa2 = ((LWPOINT*)lwgeom2)->point; - else - pa2 = ((LWLINE*)lwgeom2)->points; - - return ptarray_distance_spheroid(pa1, pa2, a, b, tolerance, check_intersection); - } - - /* Point/Polygon cases, if point-in-poly, return zero, else return distance. */ - if( ( type1 == POLYGONTYPE && type2 == POINTTYPE ) || - ( type2 == POLYGONTYPE && type1 == POINTTYPE ) ) - { - POINT2D p; - LWPOLY *lwpoly; - LWPOINT *lwpt; - GBOX gbox; - double distance = MAXFLOAT; - int i; - - if( type1 == POINTTYPE ) - { - lwpt = (LWPOINT*)lwgeom1; - lwpoly = (LWPOLY*)lwgeom2; - gbox = gbox2; - } - else - { - lwpt = (LWPOINT*)lwgeom2; - lwpoly = (LWPOLY*)lwgeom1; - gbox = gbox1; - } - getPoint2d_p(lwpt->point, 0, &p); - - /* Point in polygon implies zero distance */ - if( lwpoly_covers_point2d(lwpoly, gbox, p) ) - return 0.0; - - /* Not inside, so what's the actual distance? */ - for( i = 0; i < lwpoly->nrings; i++ ) - { - double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwpt->point, a, b, tolerance, check_intersection); - if( ring_distance < distance ) - distance = ring_distance; - if( distance < tolerance ) - return distance; - } - return distance; - } - - /* Line/polygon case, if start point-in-poly, return zero, else return distance. */ - if( ( type1 == POLYGONTYPE && type2 == LINETYPE ) || - ( type2 == POLYGONTYPE && type1 == LINETYPE ) ) - { - POINT2D p; - LWPOLY *lwpoly; - LWLINE *lwline; - GBOX gbox; - double distance = MAXFLOAT; - int i; - - if( type1 == LINETYPE ) - { - lwline = (LWLINE*)lwgeom1; - lwpoly = (LWPOLY*)lwgeom2; - gbox = gbox2; - } - else - { - lwline = (LWLINE*)lwgeom2; - lwpoly = (LWPOLY*)lwgeom1; - gbox = gbox1; - } - getPoint2d_p(lwline->points, 0, &p); - - LWDEBUG(4, "checking if a point of line is in polygon"); - - /* Point in polygon implies zero distance */ - if( lwpoly_covers_point2d(lwpoly, gbox, p) ) - return 0.0; - - LWDEBUG(4, "checking ring distances"); - - /* Not contained, so what's the actual distance? */ - for( i = 0; i < lwpoly->nrings; i++ ) - { - double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwline->points, a, b, tolerance, check_intersection); - LWDEBUGF(4, "ring[%d] ring_distance = %.8g", i, ring_distance); - if( ring_distance < distance ) - distance = ring_distance; - if( distance < tolerance ) - return distance; - } - LWDEBUGF(4, "all rings checked, returning distance = %.8g", distance); - return distance; - - } - - /* Polygon/polygon case, if start point-in-poly, return zero, else return distance. */ - if( ( type1 == POLYGONTYPE && type2 == POLYGONTYPE ) || - ( type2 == POLYGONTYPE && type1 == POLYGONTYPE ) ) - { - POINT2D p; - LWPOLY *lwpoly1 = (LWPOLY*)lwgeom1; - LWPOLY *lwpoly2 = (LWPOLY*)lwgeom2; - double distance = MAXFLOAT; - int i, j; - - /* Point of 2 in polygon 1 implies zero distance */ - getPoint2d_p(lwpoly1->rings[0], 0, &p); - if( lwpoly_covers_point2d(lwpoly2, gbox2, p) ) - return 0.0; - - /* Point of 1 in polygon 2 implies zero distance */ - getPoint2d_p(lwpoly2->rings[0], 0, &p); - if( lwpoly_covers_point2d(lwpoly1, gbox1, p) ) - return 0.0; - - /* Not contained, so what's the actual distance? */ - for( i = 0; i < lwpoly1->nrings; i++ ) - { - for( j = 0; j < lwpoly2->nrings; j++ ) - { - double ring_distance = ptarray_distance_spheroid(lwpoly1->rings[i], lwpoly2->rings[j], a, b, tolerance, check_intersection); - if( ring_distance < distance ) - distance = ring_distance; - if( distance < tolerance ) - return distance; - } - } - return distance; - } - - /* Recurse into collections */ - if( lwgeom_is_collection(type1) ) - { - int i; - double distance = MAXFLOAT; - LWCOLLECTION *col = (LWCOLLECTION*)lwgeom1; - - for( i = 0; i < col->ngeoms; i++ ) - { - double geom_distance = lwgeom_distance_spheroid(col->geoms[i], lwgeom2, gbox1, gbox2, a, b, tolerance); - if( geom_distance < distance ) - distance = geom_distance; - if( distance < tolerance ) - return distance; - } - return distance; - } - - /* Recurse into collections */ - if( lwgeom_is_collection(type2) ) - { - int i; - double distance = MAXFLOAT; - LWCOLLECTION *col = (LWCOLLECTION*)lwgeom2; - - for( i = 0; i < col->ngeoms; i++ ) - { - double geom_distance = lwgeom_distance_spheroid(lwgeom1, col->geoms[i], gbox1, gbox2, a, b, tolerance); - if( geom_distance < distance ) - distance = geom_distance; - if( distance < tolerance ) - return distance; - } - return distance; - } - - - lwerror("arguments include unsupported geometry type (%s, %s)", lwgeom_typename(type1), lwgeom_typename(type1)); - return -1.0; - -} - diff --git a/liblwgeom/lwspheroid.h b/liblwgeom/lwspheroid.h deleted file mode 100644 index f59a57a11..000000000 --- a/liblwgeom/lwspheroid.h +++ /dev/null @@ -1,17 +0,0 @@ -#include "lwgeodetic.h" - -/* -** Useful information about spheroids: -** -** a = semi-major axis -** b = semi-minor axis = a - af -** f = flattening = (a - b)/a -** e = eccentricity (first) -** e_sq = eccentricity (first) squared = (a*a - b*b)/(a*a) -** omf = 1 - f -*/ - -double spheroid_distance(GEOGRAPHIC_POINT r, GEOGRAPHIC_POINT s, double a, double b); -double ptarray_distance_spheroid(POINTARRAY *pa1, POINTARRAY *pa2, double a, double b, double tolerance, int check_intersection); - - diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c index 0e53035be..233532158 100644 --- a/postgis/geography.sql.in.c +++ b/postgis/geography.sql.in.c @@ -437,7 +437,7 @@ CREATE OR REPLACE FUNCTION ST_AsGeoJson(int4, geography, int4, int4) -- Availability: 1.5.0 CREATE OR REPLACE FUNCTION _ST_Distance(geography, geography, float8, boolean) RETURNS float8 - AS 'MODULE_PATHNAME','geography_distance_sphere' + AS 'MODULE_PATHNAME','geography_distance' LANGUAGE 'C' IMMUTABLE STRICT; -- Availability: 1.5.0 @@ -460,6 +460,12 @@ CREATE OR REPLACE FUNCTION _ST_Expand(geography, float8) AS 'MODULE_PATHNAME','geography_expand' LANGUAGE 'C' IMMUTABLE STRICT; +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION ST_DWithin(geography, geography, float8, boolean) + RETURNS boolean + AS 'SELECT $1 && _ST_Expand($2,$3) AND $2 && _ST_Expand($1,$3) AND _ST_Distance($1, $2, $3, $4) < $3' + LANGUAGE 'SQL' IMMUTABLE; + -- Currently defaulting to spheroid calculations -- Availability: 1.5.0 CREATE OR REPLACE FUNCTION ST_DWithin(geography, geography, float8) @@ -467,17 +473,30 @@ CREATE OR REPLACE FUNCTION ST_DWithin(geography, geography, float8) AS 'SELECT $1 && _ST_Expand($2,$3) AND $2 && _ST_Expand($1,$3) AND _ST_Distance($1, $2, $3, true) < $3' LANGUAGE 'SQL' IMMUTABLE; +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION ST_Area(geography, boolean) + RETURNS float8 + AS 'MODULE_PATHNAME','geography_area' + LANGUAGE 'C' IMMUTABLE STRICT; + +-- Currently defaulting to spheroid calculations -- Availability: 1.5.0 CREATE OR REPLACE FUNCTION ST_Area(geography) RETURNS float8 - AS 'MODULE_PATHNAME','geography_area_sphere' + AS 'SELECT ST_Area($1, true)' + LANGUAGE 'SQL' IMMUTABLE STRICT; + +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION ST_Length(geography, boolean) + RETURNS float8 + AS 'MODULE_PATHNAME','geography_length' LANGUAGE 'C' IMMUTABLE STRICT; -- Availability: 1.5.0 CREATE OR REPLACE FUNCTION ST_Length(geography) RETURNS float8 - AS 'MODULE_PATHNAME','geography_length_sphere' - LANGUAGE 'C' IMMUTABLE STRICT; + AS 'SELECT ST_Length($1, true)' + LANGUAGE 'SQL' IMMUTABLE; -- Availability: 1.5.0 CREATE OR REPLACE FUNCTION ST_PointOutside(geography) diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c index 330978c09..5f8282036 100644 --- a/postgis/geography_measurement.c +++ b/postgis/geography_measurement.c @@ -23,9 +23,9 @@ #include "lwgeom_pg.h" /* For debugging macros. */ #include "geography.h" /* For utility functions. */ -Datum geography_distance_sphere(PG_FUNCTION_ARGS); -Datum geography_area_sphere(PG_FUNCTION_ARGS); -Datum geography_length_sphere(PG_FUNCTION_ARGS); +Datum geography_distance(PG_FUNCTION_ARGS); +Datum geography_area(PG_FUNCTION_ARGS); +Datum geography_length(PG_FUNCTION_ARGS); Datum geography_expand(PG_FUNCTION_ARGS); Datum geography_point_outside(PG_FUNCTION_ARGS); Datum geography_covers(PG_FUNCTION_ARGS); @@ -35,8 +35,8 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS); ** geography_distance_sphere(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance) ** returns double distance in meters */ -PG_FUNCTION_INFO_V1(geography_distance_sphere); -Datum geography_distance_sphere(PG_FUNCTION_ARGS) +PG_FUNCTION_INFO_V1(geography_distance); +Datum geography_distance(PG_FUNCTION_ARGS) { LWGEOM *lwgeom1 = NULL; LWGEOM *lwgeom2 = NULL; @@ -47,11 +47,25 @@ Datum geography_distance_sphere(PG_FUNCTION_ARGS) double tolerance; double distance; bool use_spheroid; + SPHEROID s; /* Get our geometry objects loaded into memory. */ g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); g2 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + /* Read our tolerance value. */ + tolerance = PG_GETARG_FLOAT8(2); + + /* Read our calculation type. */ + use_spheroid = PG_GETARG_BOOL(3); + + /* Initialize spheroid */ + spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS); + + /* Set to sphere if requested */ + if( ! use_spheroid ) + s.a = s.b = s.radius; + /* We need the bounding boxes in case of polygon calculations, which requires them to generate a stab-line to test point-in-polygon. */ if( ! gbox_from_gserialized(g1, &gbox1) || @@ -64,34 +78,14 @@ Datum geography_distance_sphere(PG_FUNCTION_ARGS) lwgeom1 = lwgeom_from_gserialized(g1); lwgeom2 = lwgeom_from_gserialized(g2); - /* Read our tolerance value. */ - tolerance = PG_GETARG_FLOAT8(2); - - /* Read our calculation type. */ - use_spheroid = PG_GETARG_BOOL(3); - - /* Sphere returns in radians, spheroid returns in meters. - Convert to radians for sphere. */ - if( ! use_spheroid ) - tolerance = tolerance / WGS84_RADIUS; - - /* Calculate the distance */ - if( use_spheroid ) - distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, gbox1, gbox2, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS, tolerance); - else - distance = lwgeom_distance_sphere(lwgeom1, lwgeom2, gbox1, gbox2, tolerance); + distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, gbox1, gbox2, s, tolerance); /* Something went wrong... should already be eloged */ if( distance < 0.0 ) { - elog(ERROR, "lwgeom_distance_sphere(oid) returned < 0.0"); + elog(ERROR, "lwgeom_distance_spheroid returned < 0.0"); PG_RETURN_NULL(); } - - /* Sphere returns in radians, spheroid returns in meters. - Convert from radians for sphere. */ - if( ! use_spheroid ) - distance = distance * WGS84_RADIUS; /* Clean up, but not all the way to the point arrays */ lwgeom_release(lwgeom1); @@ -147,19 +141,31 @@ Datum geography_expand(PG_FUNCTION_ARGS) } /* -** geography_area_sphere(GSERIALIZED *g) +** geography_area(GSERIALIZED *g) ** returns double area in meters square */ -PG_FUNCTION_INFO_V1(geography_area_sphere); -Datum geography_area_sphere(PG_FUNCTION_ARGS) +PG_FUNCTION_INFO_V1(geography_area); +Datum geography_area(PG_FUNCTION_ARGS) { LWGEOM *lwgeom = NULL; GBOX gbox; GSERIALIZED *g = NULL; double area; - + bool use_spheroid = LW_TRUE; + SPHEROID s; + /* Get our geometry object loaded into memory. */ g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + /* Read our calculation type */ + use_spheroid = PG_GETARG_BOOL(1); + + /* Initialize spheroid */ + spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS); + + /* User requests spherical calculation, turn our spheroid into a sphere */ + if( ! use_spheroid ) + s.a = s.b = s.radius; /* We need the bounding box to get an outside point for area algorithm */ if( ! gbox_from_gserialized(g, &gbox) ) @@ -171,7 +177,7 @@ Datum geography_area_sphere(PG_FUNCTION_ARGS) lwgeom = lwgeom_from_gserialized(g); /* Calculate the area */ - area = lwgeom_area_sphere(lwgeom, gbox); + area = lwgeom_area_spheroid(lwgeom, gbox, s); /* Something went wrong... */ if( area < 0.0 ) @@ -180,10 +186,6 @@ Datum geography_area_sphere(PG_FUNCTION_ARGS) PG_RETURN_NULL(); } - /* Currently normalizing with a fixed WGS84 radius, in future this - should be the average radius of the SRID in play */ - area = area * WGS84_RADIUS * WGS84_RADIUS; - /* Clean up, but not all the way to the point arrays */ lwgeom_release(lwgeom); @@ -195,20 +197,31 @@ Datum geography_area_sphere(PG_FUNCTION_ARGS) ** geography_length_sphere(GSERIALIZED *g) ** returns double length in meters */ -PG_FUNCTION_INFO_V1(geography_length_sphere); -Datum geography_length_sphere(PG_FUNCTION_ARGS) +PG_FUNCTION_INFO_V1(geography_length); +Datum geography_length(PG_FUNCTION_ARGS) { LWGEOM *lwgeom = NULL; GSERIALIZED *g = NULL; double length; + bool use_spheroid = LW_TRUE; + SPHEROID s; /* Get our geometry object loaded into memory. */ g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - lwgeom = lwgeom_from_gserialized(g); + + /* Read our calculation type */ + use_spheroid = PG_GETARG_BOOL(1); + + /* Initialize spheroid */ + spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS); + + /* User requests spherical calculation, turn our spheroid into a sphere */ + if( ! use_spheroid ) + s.a = s.b = s.radius; /* Calculate the length */ - length = lwgeom_length_sphere(lwgeom); + length = lwgeom_length_spheroid(lwgeom, s); /* Something went wrong... */ if( length < 0.0 ) @@ -217,10 +230,6 @@ Datum geography_length_sphere(PG_FUNCTION_ARGS) PG_RETURN_NULL(); } - /* Currently normalizing with a fixed WGS84 radius, in future this - should be the average radius of the SRID in play */ - length = length * WGS84_RADIUS; - /* Clean up, but not all the way to the point arrays */ lwgeom_release(lwgeom); -- 2.50.1