From: Paul Ramsey Date: Fri, 30 Oct 2009 20:45:50 +0000 (+0000) Subject: Add in spheroid calculations for ST_Distance and ST_DWithin. X-Git-Tag: 1.5.0b1~321 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=260b853a5879566f91aaa1f22335e12a1fe6f583;p=postgis Add in spheroid calculations for ST_Distance and ST_DWithin. git-svn-id: http://svn.osgeo.org/postgis/trunk@4709 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/libgeom.h b/liblwgeom/libgeom.h index 5714ab92c..ef2e06020 100644 --- a/liblwgeom/libgeom.h +++ b/liblwgeom/libgeom.h @@ -376,6 +376,12 @@ extern int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox); */ extern double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, double tolerance); +/** +* 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). +*/ +extern double lwgeom_distance_spheroid(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, double a, double b, 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. diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index 0e5b2ec29..96d15ad3e 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -1350,7 +1350,7 @@ void gbox_pt_outside(GBOX gbox, POINT2D *pt_outside) * Returns the area of the ring (ring must be closed) in square radians (surface of * the sphere is 4*PI). */ -double ptarray_area_sphere(POINTARRAY *pa, POINT2D pt_outside) +double ptarray_area_sphere(POINTARRAY *pa, POINT2D pt_outside) { GEOGRAPHIC_POINT a, b, c; POINT2D p; diff --git a/liblwgeom/lwspheroid.c b/liblwgeom/lwspheroid.c index 1af59c497..49840724e 100644 --- a/liblwgeom/lwspheroid.c +++ b/liblwgeom/lwspheroid.c @@ -99,4 +99,352 @@ double spheroid_distance(GEOGRAPHIC_POINT r, GEOGRAPHIC_POINT s, double a, doubl } return distance; -} \ No newline at end of file +} + +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 index 023fabd66..f59a57a11 100644 --- a/liblwgeom/lwspheroid.h +++ b/liblwgeom/lwspheroid.h @@ -12,5 +12,6 @@ */ 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.h b/postgis/geography.h index ec601c5d8..dc8d6d406 100644 --- a/postgis/geography.h +++ b/postgis/geography.h @@ -21,6 +21,7 @@ #define WGS84_MINOR_AXIS 6356752.314245179498 #define WGS84_RADIUS ((2.0 * WGS84_MAJOR_AXIS + WGS84_MINOR_AXIS ) / 3.0) + /********************************************************************** ** Useful functions for all GEOGRAPHY handlers. */ diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c index 28bfe3c2c..0e53035be 100644 --- a/postgis/geography.sql.in.c +++ b/postgis/geography.sql.in.c @@ -435,15 +435,22 @@ CREATE OR REPLACE FUNCTION ST_AsGeoJson(int4, geography, int4, int4) -- Stop calculation and return immediately once distance is less than tolerance -- Availability: 1.5.0 -CREATE OR REPLACE FUNCTION _ST_Distance(geography, geography, float8) +CREATE OR REPLACE FUNCTION _ST_Distance(geography, geography, float8, boolean) RETURNS float8 AS 'MODULE_PATHNAME','geography_distance_sphere' LANGUAGE 'C' IMMUTABLE STRICT; +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION ST_Distance(geography, geography, boolean) + RETURNS float8 + AS 'SELECT _ST_Distance($1, $2, 0.0, $3)' + LANGUAGE 'SQL' IMMUTABLE STRICT; + +-- Currently defaulting to spheroid calculations -- Availability: 1.5.0 CREATE OR REPLACE FUNCTION ST_Distance(geography, geography) RETURNS float8 - AS 'SELECT _ST_Distance($1, $2, 0.0)' + AS 'SELECT _ST_Distance($1, $2, 0.0, true)' LANGUAGE 'SQL' IMMUTABLE STRICT; -- Only expands the bounding box, the actual geometry will remain unchanged, use with care. @@ -453,10 +460,11 @@ CREATE OR REPLACE FUNCTION _ST_Expand(geography, float8) AS 'MODULE_PATHNAME','geography_expand' LANGUAGE 'C' IMMUTABLE STRICT; +-- Currently defaulting to spheroid calculations -- Availability: 1.5.0 CREATE OR REPLACE FUNCTION ST_DWithin(geography, geography, float8) RETURNS boolean - AS 'SELECT $1 && _ST_Expand($2,$3) AND $2 && _ST_Expand($1,$3) AND _ST_Distance($1, $2, $3) < $3' + 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 diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c index ad73e5176..330978c09 100644 --- a/postgis/geography_measurement.c +++ b/postgis/geography_measurement.c @@ -46,6 +46,7 @@ Datum geography_distance_sphere(PG_FUNCTION_ARGS) GSERIALIZED *g2 = NULL; double tolerance; double distance; + bool use_spheroid; /* Get our geometry objects loaded into memory. */ g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); @@ -65,21 +66,32 @@ Datum geography_distance_sphere(PG_FUNCTION_ARGS) /* Read our tolerance value. */ tolerance = PG_GETARG_FLOAT8(2); - tolerance = tolerance / WGS84_RADIUS; + + /* 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 */ - distance = lwgeom_distance_sphere(lwgeom1, lwgeom2, gbox1, gbox2, tolerance); + 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); /* Something went wrong... should already be eloged */ if( distance < 0.0 ) { - elog(ERROR, "lwgeom_distance_sphere returned < 0.0"); + elog(ERROR, "lwgeom_distance_sphere(oid) returned < 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 */ - distance = distance * WGS84_RADIUS; + /* 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);