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);
+
+ /* Point in polygon implies zero distance */
+ if( lwpoly_covers_point2d(lwpoly, gbox, p) )
+ return 0.0;
+
+ /* 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);
+ if( ring_distance < distance )
+ distance = ring_distance;
+ if( distance < tolerance )
+ return distance;
+ }
+ return distance;
+
}
/* Polygon/polygon case, if start point-in-poly, return zero, else return distance. */
if( ( type1 == POLYGONTYPE && type2 == LINETYPE ) ||
( type2 == POLYGONTYPE && type1 == LINETYPE ) )
{
- }
+ 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_sphere(lwpoly1->rings[i], lwpoly2->rings[j], tolerance);
+ if( ring_distance < distance )
+ distance = ring_distance;
+ if( distance < tolerance )
+ return distance;
+ }
+ }
+ return distance;
+ }
+
+ /* Recurse into collections */
+ if( lwgeom_contains_subgeoms(type1) )
+ {
+ int i;
+ double distance = MAXFLOAT;
+ LWCOLLECTION *col = (LWCOLLECTION*)lwgeom1;
+
+ for( i = 0; i < col->ngeoms; i++ )
+ {
+ double geom_distance = lwgeom_distance_sphere(col->geoms[i], lwgeom2, gbox1, gbox2, tolerance);
+ if( geom_distance < distance )
+ distance = geom_distance;
+ if( distance < tolerance )
+ return distance;
+ }
+ return distance;
+ }
+
+ /* Recurse into collections */
+ if( lwgeom_contains_subgeoms(type2) )
+ {
+ int i;
+ double distance = MAXFLOAT;
+ LWCOLLECTION *col = (LWCOLLECTION*)lwgeom2;
+
+ for( i = 0; i < col->ngeoms; i++ )
+ {
+ double geom_distance = lwgeom_distance_sphere(lwgeom1, col->geoms[i], gbox1, gbox2, 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;
}
lwgeom_rtree.o \
geography_inout.o \
geography_gist.o \
- geography_estimate.o
+ geography_estimate.o \
+ geography_distance.o
# Objects to build using PGXS
OBJS=$(PG_OBJS)
*
**********************************************************************/
+
+/**********************************************************************
+** GIDX structure.
+**
+** This is an n-dimensional (practically, the
+** implementation is 2-4 dimensions) box used for index keys. The
+** coordinates are anonymous, so we can have any number of dimensions.
+** The sizeof a GIDX is 1 + 2 * ndims * sizeof(float).
+** The order of values in a GIDX is
+** xmin,xmax,ymin,ymax,zmin,zmax...
+*/
+typedef struct
+{
+ int32 varsize;
+ float c[1];
+} GIDX;
+
+/*
+** For some GiST support functions, it is more efficient to allocate
+** memory for our GIDX off the stack and cast that memory into a GIDX.
+** But, GIDX is variable length, what to do? We'll bake in the assumption
+** that no GIDX is more than 4-dimensional for now, and ensure that much
+** space is available.
+** 4 bytes varsize + 4 dimensions * 2 ordinates * 4 bytes float size = 36 bytes
+*/
+#define GIDX_MAX_SIZE 36
+
+
/* Useful functions for all GEOGRAPHY handlers. */
GSERIALIZED* geography_serialize(LWGEOM *lwgeom);
void geography_valid_typmod(LWGEOM *lwgeom, int32 typmod);
void geography_valid_type(uchar type);
+int geography_datum_gidx(Datum geography_datum, GIDX *gidx);
+GIDX* gidx_new(int ndims);
+void gbox_from_gidx(GIDX *gidx, GBOX *gbox);
--
CREATE OR REPLACE FUNCTION _ST_AsGeoJson(int4, geography, int4, int4)
- RETURNS text
+ RETURNS text
AS 'MODULE_PATHNAME','geography_as_geojson'
LANGUAGE 'C' IMMUTABLE STRICT;
-- ST_AsGeoJson(geography, precision) / version=1 options=0
CREATE OR REPLACE FUNCTION ST_AsGeoJson(geography, int4)
- RETURNS TEXT
- AS 'SELECT _ST_AsGeoJson(1, $1, $2, 0)'
- LANGUAGE 'SQL' IMMUTABLE STRICT;
+ RETURNS TEXT
+ AS 'SELECT _ST_AsGeoJson(1, $1, $2, 0)'
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
-- ST_AsGeoJson(geography) / precision=15 version=1 options=0
CREATE OR REPLACE FUNCTION ST_AsGeoJson(geography)
- RETURNS TEXT
- AS 'SELECT _ST_AsGeoJson(1, $1, 15, 0)'
- LANGUAGE 'SQL' IMMUTABLE STRICT;
+ RETURNS TEXT
+ AS 'SELECT _ST_AsGeoJson(1, $1, 15, 0)'
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
-- ST_AsGeoJson(version, geography) / precision=15 options=0
CREATE OR REPLACE FUNCTION ST_AsGeoJson(int4, geography)
- RETURNS TEXT
- AS 'SELECT _ST_AsGeoJson($1, $2, 15, 0)'
- LANGUAGE 'SQL' IMMUTABLE STRICT;
+ RETURNS TEXT
+ AS 'SELECT _ST_AsGeoJson($1, $2, 15, 0)'
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
-- ST_AsGeoJson(version, geography, precision) / options=0
CREATE OR REPLACE FUNCTION ST_AsGeoJson(int4, geography, int4)
- RETURNS TEXT
- AS 'SELECT _ST_AsGeoJson($1, $2, $3, 0)'
- LANGUAGE 'SQL' IMMUTABLE STRICT;
+ RETURNS TEXT
+ AS 'SELECT _ST_AsGeoJson($1, $2, $3, 0)'
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
-- ST_AsGeoJson(geography, precision, options) / version=1
CREATE OR REPLACE FUNCTION ST_AsGeoJson(geography, int4, int4)
- RETURNS TEXT
- AS 'SELECT _ST_AsGeoJson(1, $1, $2, $3)'
- LANGUAGE 'SQL' IMMUTABLE STRICT;
+ RETURNS TEXT
+ AS 'SELECT _ST_AsGeoJson(1, $1, $2, $3)'
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
-- ST_AsGeoJson(version, geography, precision,options)
CREATE OR REPLACE FUNCTION ST_AsGeoJson(int4, geography, int4, int4)
- RETURNS TEXT
- AS 'SELECT _ST_AsGeoJson($1, $2, $3, $4)'
- LANGUAGE 'SQL' IMMUTABLE STRICT;
+ RETURNS TEXT
+ AS 'SELECT _ST_AsGeoJson($1, $2, $3, $4)'
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+-- ---------- ---------- ---------- ---------- ---------- ---------- ----------
+-- Measurement Functions
+-- Availability: 1.5.0
+-- ---------- ---------- ---------- ---------- ---------- ---------- ----------
+-- Stop calculation and return once distance is less than tolerance
+CREATE OR REPLACE FUNCTION _ST_Distance(geography, geography, float8)
+ RETURNS float8
+ AS 'MODULE_PATHNAME','geography_distance_sphere'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION ST_Distance(geography, geography)
+ RETURNS float8
+ AS 'SELECT _ST_Distance($1, $2, 0.0)'
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+
+
+-- ---------- ---------- ---------- ---------- ---------- ---------- ----------
COMMIT;
--- /dev/null
+/**********************************************************************
+ * $Id: geography_inout.c 4535 2009-09-28 18:16:21Z colivier $
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "postgres.h"
+
+#include "../postgis_config.h"
+
+#include <math.h>
+#include <float.h>
+#include <string.h>
+#include <stdio.h>
+#include <errno.h>
+
+#include "libgeom.h" /* For standard geometry types. */
+#include "lwgeom_pg.h" /* For debugging macros. */
+#include "geography.h" /* For utility functions. */
+
+Datum geography_distance_sphere(PG_FUNCTION_ARGS);
+
+
+/*
+** geography_in(cstring) returns *GSERIALIZED
+*/
+PG_FUNCTION_INFO_V1(geography_distance_sphere);
+Datum geography_distance_sphere(PG_FUNCTION_ARGS)
+{
+ LWGEOM *lwgeom1 = NULL;
+ LWGEOM *lwgeom2 = NULL;
+ GBOX gbox1;
+ GBOX gbox2;
+ GIDX *gidx1 = gidx_new(3);
+ GIDX *gidx2 = gidx_new(3);
+ GSERIALIZED *g1 = NULL;
+ GSERIALIZED *g2 = NULL;
+ double tolerance;
+ double distance;
+
+ /* We need the bounding boxes in case of polygon calculations,
+ which requires them to generate a stab-line to test point-in-polygon. */
+ geography_datum_gidx(PG_GETARG_DATUM(0), gidx1);
+ geography_datum_gidx(PG_GETARG_DATUM(0), gidx2);
+ gbox_from_gidx(gidx1, &gbox1);
+ gbox_from_gidx(gidx2, &gbox2);
+ pfree(gidx1);
+ pfree(gidx2);
+
+ /* 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));
+ lwgeom1 = lwgeom_from_gserialized(g1);
+ lwgeom2 = lwgeom_from_gserialized(g2);
+
+ /* Read our tolerance value. */
+ tolerance = PG_GETARG_FLOAT8(2);
+
+ /* Calculate the distance */
+ distance = lwgeom_distance_sphere(lwgeom1, lwgeom2, &gbox1, &gbox2, tolerance);
+
+ /* Something went wrong... should already be eloged */
+ if( distance < 0.0 )
+ {
+ elog(ERROR, "Error in geography_distance_sphere calculation.");
+ PG_RETURN_NULL();
+ }
+
+ /* Clean up, but not all the way to the point arrays */
+ lwgeom_release(lwgeom1);
+ lwgeom_release(lwgeom2);
+
+ PG_RETURN_FLOAT8(distance);
+
+}
+
#include "lwgeom_pg.h" /* For debugging macros. */
#include "geography.h" /* For utility functions. */
-/*********************************************************************************
-** GIDX structure.
-**
-** This is an n-dimensional (practically, the
-** implementation is 2-4 dimensions) box used for index keys. The
-** coordinates are anonymous, so we can have any number of dimensions.
-** The sizeof a GIDX is 1 + 2 * ndims * sizeof(float).
-** The order of values in a GIDX is
-** xmin,xmax,ymin,ymax,zmin,zmax...
-*/
-typedef struct
-{
- int32 varsize;
- float c[1];
-} GIDX;
-/*
-** For some GiST support functions, it is more efficient to allocate
-** memory for our GIDX off the stack and cast that memory into a GIDX.
-** But, GIDX is variable length, what to do? We'll bake in the assumption
-** that no GIDX is more than 4-dimensional for now, and ensure that much
-** space is available.
-** 4 bytes varsize + 4 dimensions * 2 ordinates * 4 bytes float size = 36 bytes
-*/
-#define GIDX_MAX_SIZE 36
-
/*
** When is a node split not so good? If more than 90% of the entries
** end up in one of the children.
/* Allocates a new GIDX on the heap of the requested dimensionality */
-static GIDX* gidx_new(int ndims)
+GIDX* gidx_new(int ndims)
{
size_t size = GIDX_SIZE(ndims);
GIDX *g = (GIDX*)palloc(size);
}
+void gbox_from_gidx(GIDX *a, GBOX *gbox)
+{
+ gbox->xmin = (double)GIDX_GET_MIN(a,0);
+ gbox->ymin = (double)GIDX_GET_MIN(a,1);
+ gbox->zmin = (double)GIDX_GET_MIN(a,2);
+ gbox->xmax = (double)GIDX_GET_MAX(a,0);
+ gbox->ymax = (double)GIDX_GET_MAX(a,1);
+ gbox->zmax = (double)GIDX_GET_MAX(a,2);
+}
+
+
/*
** Overlapping GIDX box test.
**
** full geography and return the box based on that. If no box is available,
** return G_FAILURE, otherwise G_SUCCESS.
*/
-static int geography_datum_gidx(Datum geography_datum, GIDX *gidx)
+int geography_datum_gidx(Datum geography_datum, GIDX *gidx)
{
GSERIALIZED *gpart;
int result = G_SUCCESS;
projects / postgis / commitdiff