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);
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);
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);
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);
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);
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);
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);
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);
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);
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
#include <string.h>
#include "CUnit/Basic.h"
-#include "lwspheroid.h"
+#include "lwgeodetic.h"
#include "cu_tester.h"
/***********************************************************************
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
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;
}
-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;
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 */
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 */
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, "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;
}
LWDEBUGF(4,"finished all loops, returning %.8g", distance);
- return distance;
+ if( s.a == s.b )
+ return distance;
+ else
+ return spheroid_distance(nearest1, nearest2, s);
}
* 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);
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;
}
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;
}
* 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;
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. */
/* 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 )
/* 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;
{
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 )
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 )
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 )
return LW_FALSE;
}
-double ptarray_length_sphere(POINTARRAY *pa)
+double ptarray_length_spheroid(POINTARRAY *pa, SPHEROID s)
{
GEOGRAPHIC_POINT a, b;
POINT2D p;
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;
return length;
}
-double lwgeom_length_sphere(LWGEOM *geom)
+double lwgeom_length_spheroid(LWGEOM *geom, SPHEROID s)
{
int type;
int i = 0;
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;
}
for( i = 0; i < col->ngeoms; i++ )
{
- length += lwgeom_length_sphere(col->geoms[i]);
+ length += lwgeom_length_spheroid(col->geoms[i], s);
}
return length;
}
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);
-#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)
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;
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);
}
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;
-
-}
-
+++ /dev/null
-#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);
-
-
-- 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
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)
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)
#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);
** 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;
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) ||
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);
}
/*
-** 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) )
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 )
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);
** 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 )
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);
projects / postgis / commitdiff