lwt_ModEdgeHeal X
lwt_NewEdgeHeal X
+lwt_GetFaceByPoint X
+lwt_GetNodeByPoint X
+lwt_GetEdgeByPoint X
+lwt_TopoGeo_AddPoint
+lwt_TopoGeo_AddLineString
+lwt_TopoGeo_AddPolygon
*/
void lwt_FreeTopology(LWT_TOPOLOGY* topo);
+/**
+ * Retrieve the id of a node at a point location
+ *
+ * @param topo the topology to operate on
+ * @param point the point to use for query
+ * @param tol max distance around the given point to look for a node
+ * @return a node identifier if one is found, 0 if none is found, -1
+ * on error (multiple nodes within distance).
+ * The liblwgeom error handler will be invoked in case of error.
+ */
+LWT_ELEMID lwt_GetNodeByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol);
+
+/**
+ * Find the edge-id of an edge that intersects a given point
+ *
+ * @param topo the topology to operate on
+ * @param point the point to use for query
+ * @param tol max distance around the given point to look for an
+ * intersecting edge
+ * @return an edge identifier if one is found, 0 if none is found, -1
+ * on error (multiple edges within distance).
+ * The liblwgeom error handler will be invoked in case of error.
+ */
+LWT_ELEMID lwt_GetEdgeByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol);
+
+/**
+ * Find the face-id of a face containing a given point
+ *
+ * @param topo the topology to operate on
+ * @param point the point to use for query
+ * @param tol max distance around the given point to look for a
+ * containing face
+ * @return a face identifier if one is found (0 if universe), -1
+ * on error (multiple faces within distance or point on node
+ * or edge).
+ * The liblwgeom error handler will be invoked in case of error.
+ */
+LWT_ELEMID lwt_GetFaceByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol);
+
/*******************************************************************
*
* Topology population (non-ISO)
#include "../postgis_config.h"
-#define POSTGIS_DEBUG_LEVEL 1
+/*#define POSTGIS_DEBUG_LEVEL 1*/
#include "lwgeom_log.h"
#include "liblwgeom_internal.h"
{
return _lwt_HealEdges( topo, e1, e2, 0 );
}
+
+LWT_ELEMID
+lwt_GetNodeByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol)
+{
+ LWT_ISO_NODE *elem;
+ int num;
+ int flds = LWT_COL_NODE_NODE_ID|LWT_COL_NODE_GEOM; /* geom not needed */
+ LWT_ELEMID id = 0;
+ POINT2D qp; /* query point */
+
+ if ( ! getPoint2d_p(pt->point, 0, &qp) )
+ {
+ lwerror("Empty query point");
+ return -1;
+ }
+ elem = lwt_be_getNodeWithinDistance2D(topo, pt, tol, &num, flds, 0);
+ if ( num == -1 )
+ {
+ lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
+ return -1;
+ }
+ else if ( num )
+ {
+ if ( num > 1 )
+ {
+ _lwt_release_nodes(elem, num);
+ lwerror("Two or more nodes found");
+ return -1;
+ }
+ id = elem[0].node_id;
+ }
+
+ return id;
+}
+
+LWT_ELEMID
+lwt_GetEdgeByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol)
+{
+ LWT_ISO_EDGE *elem;
+ int num, i;
+ int flds = LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM; /* GEOM is not needed */
+ LWT_ELEMID id = 0;
+ LWGEOM *qp = lwpoint_as_lwgeom(pt); /* query point */
+
+ if ( lwgeom_is_empty(qp) )
+ {
+ lwerror("Empty query point");
+ return -1;
+ }
+ elem = lwt_be_getEdgeWithinDistance2D(topo, pt, tol, &num, flds, 0);
+ if ( num == -1 )
+ {
+ lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
+ return -1;
+ }
+ for (i=0; i<num;++i)
+ {
+ LWT_ISO_EDGE *e = &(elem[i]);
+#if 0
+ LWGEOM* geom;
+ double dist;
+
+ if ( ! e->geom )
+ {
+ _lwt_release_edges(elem, num);
+ lwnotice("Corrupted topology: edge %" LWTFMT_ELEMID
+ " has null geometry", e->edge_id);
+ continue;
+ }
+
+ /* Should we check for intersection not being on an endpoint
+ * as documented ? */
+ geom = lwline_as_lwgeom(e->geom);
+ dist = lwgeom_mindistance2d_tolerance(geom, qp, tol);
+ if ( dist > tol ) continue;
+#endif
+
+ if ( id )
+ {
+ _lwt_release_edges(elem, num);
+ lwerror("Two or more edges found");
+ return -1;
+ }
+ else id = e->edge_id;
+ }
+
+ return id;
+}
+
+LWT_ELEMID
+lwt_GetFaceByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol)
+{
+ LWT_ELEMID id = 0;
+ LWT_ISO_EDGE *elem;
+ int num, i;
+ int flds = LWT_COL_EDGE_EDGE_ID |
+ LWT_COL_EDGE_GEOM |
+ LWT_COL_EDGE_FACE_LEFT |
+ LWT_COL_EDGE_FACE_RIGHT;
+ LWGEOM *qp = lwpoint_as_lwgeom(pt);
+
+ id = lwt_be_getFaceContainingPoint(topo, pt);
+ if ( id == -2 ) {
+ lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
+ return -1;
+ }
+
+ if ( id > 0 )
+ {
+ return id;
+ }
+ id = 0; /* or it'll be -1 for not found */
+
+ LWDEBUG(1, "No face properly contains query point,"
+ " looking for edges");
+
+ /* Not in a face, may be in universe or on edge, let's check
+ * for distance */
+ /* NOTE: we never pass a tolerance of 0 to avoid ever using
+ * ST_Within, which doesn't include endpoints matches */
+ elem = lwt_be_getEdgeWithinDistance2D(topo, pt, tol?tol:1e-5, &num, flds, 0);
+ if ( num == -1 )
+ {
+ lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
+ return -1;
+ }
+ for (i=0; i<num; ++i)
+ {
+ LWT_ISO_EDGE *e = &(elem[i]);
+ LWT_ELEMID eface = 0;
+ LWGEOM* geom;
+ double dist;
+
+ if ( ! e->geom )
+ {
+ _lwt_release_edges(elem, num);
+ lwnotice("Corrupted topology: edge %" LWTFMT_ELEMID
+ " has null geometry", e->edge_id);
+ continue;
+ }
+
+ /* don't consider dangling edges */
+ if ( e->face_left == e->face_right )
+ {
+ LWDEBUGF(1, "Edge %" LWTFMT_ELEMID
+ " is dangling, won't consider it", e->edge_id);
+ continue;
+ }
+
+ geom = lwline_as_lwgeom(e->geom);
+ dist = lwgeom_mindistance2d_tolerance(geom, qp, tol);
+
+ LWDEBUGF(1, "Distance from edge %" LWTFMT_ELEMID
+ " is %g (tol=%g)", e->edge_id, dist, tol);
+
+ /* we won't consider edges too far */
+ if ( dist > tol ) continue;
+ if ( e->face_left == 0 ) {
+ eface = e->face_right;
+ }
+ else if ( e->face_right == 0 ) {
+ eface = e->face_left;
+ }
+ else {
+ _lwt_release_edges(elem, num);
+ lwerror("Two or more faces found");
+ return -1;
+ }
+
+ if ( id && id != eface )
+ {
+ _lwt_release_edges(elem, num);
+ lwerror("Two or more faces found"
+#if 0 /* debugging */
+ " (%" LWTFMT_ELEMID
+ " and %" LWTFMT_ELEMID ")", id, eface
+#endif
+ );
+ return -1;
+ }
+ else id = eface;
+ }
+ if ( num ) _lwt_release_edges(elem, num);
+
+ return id;
+}
#include "liblwgeom_internal.h" /* for gbox_clone */
#include "liblwgeom_topo.h"
-#define POSTGIS_DEBUG_LEVEL 1
+
+/*#define POSTGIS_DEBUG_LEVEL 1*/
#include "lwgeom_log.h"
#include "lwgeom_pg.h"
PG_RETURN_INT32(ret);
}
+
+/* GetNodeByPoint(atopology, point, tolerance) */
+Datum GetNodeByPoint(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(GetNodeByPoint);
+Datum GetNodeByPoint(PG_FUNCTION_ARGS)
+{
+ text* toponame_text;
+ char* toponame;
+ double tol;
+ LWT_ELEMID node_id;
+ GSERIALIZED *geom;
+ LWGEOM *lwgeom;
+ LWPOINT *pt;
+ LWT_TOPOLOGY *topo;
+
+ toponame_text = PG_GETARG_TEXT_P(0);
+ toponame = text2cstring(toponame_text);
+ PG_FREE_IF_COPY(toponame_text, 0);
+
+ geom = PG_GETARG_GSERIALIZED_P(1);
+ lwgeom = lwgeom_from_gserialized(geom);
+ pt = lwgeom_as_lwpoint(lwgeom);
+ if ( ! pt ) {
+ lwgeom_free(lwgeom);
+ PG_FREE_IF_COPY(geom, 1);
+ lwpgerror("Node geometry must be a point");
+ PG_RETURN_NULL();
+ }
+
+ tol = PG_GETARG_FLOAT8(2);
+ if ( tol < 0 )
+ {
+ PG_FREE_IF_COPY(geom, 1);
+ lwpgerror("Tolerance must be >=0");
+ PG_RETURN_NULL();
+ }
+
+ if ( SPI_OK_CONNECT != SPI_connect() ) {
+ lwpgerror("Could not connect to SPI");
+ PG_RETURN_NULL();
+ }
+ be_data.data_changed = false;
+
+ topo = lwt_LoadTopology(be_iface, toponame);
+ pfree(toponame);
+ if ( ! topo ) {
+ /* should never reach this point, as lwerror would raise an exception */
+ SPI_finish();
+ PG_RETURN_NULL();
+ }
+
+ POSTGIS_DEBUG(1, "Calling lwt_GetNodeByPoint");
+ node_id = lwt_GetNodeByPoint(topo, pt, tol);
+ POSTGIS_DEBUG(1, "lwt_GetNodeByPoint returned");
+ lwgeom_free(lwgeom);
+ PG_FREE_IF_COPY(geom, 1);
+ lwt_FreeTopology(topo);
+
+ if ( node_id == -1 ) {
+ /* should never reach this point, as lwerror would raise an exception */
+ SPI_finish();
+ PG_RETURN_NULL();
+ }
+
+ SPI_finish();
+ PG_RETURN_INT32(node_id);
+}
+
+/* GetEdgeByPoint(atopology, point, tolerance) */
+Datum GetEdgeByPoint(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(GetEdgeByPoint);
+Datum GetEdgeByPoint(PG_FUNCTION_ARGS)
+{
+ text* toponame_text;
+ char* toponame;
+ double tol;
+ LWT_ELEMID node_id;
+ GSERIALIZED *geom;
+ LWGEOM *lwgeom;
+ LWPOINT *pt;
+ LWT_TOPOLOGY *topo;
+
+ toponame_text = PG_GETARG_TEXT_P(0);
+ toponame = text2cstring(toponame_text);
+ PG_FREE_IF_COPY(toponame_text, 0);
+
+ geom = PG_GETARG_GSERIALIZED_P(1);
+ lwgeom = lwgeom_from_gserialized(geom);
+ pt = lwgeom_as_lwpoint(lwgeom);
+ if ( ! pt ) {
+ lwgeom_free(lwgeom);
+ PG_FREE_IF_COPY(geom, 1);
+ lwpgerror("Node geometry must be a point");
+ PG_RETURN_NULL();
+ }
+
+ tol = PG_GETARG_FLOAT8(2);
+ if ( tol < 0 )
+ {
+ PG_FREE_IF_COPY(geom, 1);
+ lwpgerror("Tolerance must be >=0");
+ PG_RETURN_NULL();
+ }
+
+ if ( SPI_OK_CONNECT != SPI_connect() ) {
+ lwpgerror("Could not connect to SPI");
+ PG_RETURN_NULL();
+ }
+ be_data.data_changed = false;
+
+ topo = lwt_LoadTopology(be_iface, toponame);
+ pfree(toponame);
+ if ( ! topo ) {
+ /* should never reach this point, as lwerror would raise an exception */
+ SPI_finish();
+ PG_RETURN_NULL();
+ }
+
+ POSTGIS_DEBUG(1, "Calling lwt_GetEdgeByPoint");
+ node_id = lwt_GetEdgeByPoint(topo, pt, tol);
+ POSTGIS_DEBUG(1, "lwt_GetEdgeByPoint returned");
+ lwgeom_free(lwgeom);
+ PG_FREE_IF_COPY(geom, 1);
+ lwt_FreeTopology(topo);
+
+ if ( node_id == -1 ) {
+ /* should never reach this point, as lwerror would raise an exception */
+ SPI_finish();
+ PG_RETURN_NULL();
+ }
+
+ SPI_finish();
+ PG_RETURN_INT32(node_id);
+}
+
+/* GetFaceByPoint(atopology, point, tolerance) */
+Datum GetFaceByPoint(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(GetFaceByPoint);
+Datum GetFaceByPoint(PG_FUNCTION_ARGS)
+{
+ text* toponame_text;
+ char* toponame;
+ double tol;
+ LWT_ELEMID node_id;
+ GSERIALIZED *geom;
+ LWGEOM *lwgeom;
+ LWPOINT *pt;
+ LWT_TOPOLOGY *topo;
+
+ toponame_text = PG_GETARG_TEXT_P(0);
+ toponame = text2cstring(toponame_text);
+ PG_FREE_IF_COPY(toponame_text, 0);
+
+ geom = PG_GETARG_GSERIALIZED_P(1);
+ lwgeom = lwgeom_from_gserialized(geom);
+ pt = lwgeom_as_lwpoint(lwgeom);
+ if ( ! pt ) {
+ lwgeom_free(lwgeom);
+ PG_FREE_IF_COPY(geom, 1);
+ lwpgerror("Node geometry must be a point");
+ PG_RETURN_NULL();
+ }
+
+ tol = PG_GETARG_FLOAT8(2);
+ if ( tol < 0 )
+ {
+ PG_FREE_IF_COPY(geom, 1);
+ lwpgerror("Tolerance must be >=0");
+ PG_RETURN_NULL();
+ }
+
+ if ( SPI_OK_CONNECT != SPI_connect() ) {
+ lwpgerror("Could not connect to SPI");
+ PG_RETURN_NULL();
+ }
+ be_data.data_changed = false;
+
+ topo = lwt_LoadTopology(be_iface, toponame);
+ pfree(toponame);
+ if ( ! topo ) {
+ /* should never reach this point, as lwerror would raise an exception */
+ SPI_finish();
+ PG_RETURN_NULL();
+ }
+
+ POSTGIS_DEBUG(1, "Calling lwt_GetFaceByPoint");
+ node_id = lwt_GetFaceByPoint(topo, pt, tol);
+ POSTGIS_DEBUG(1, "lwt_GetFaceByPoint returned");
+ lwgeom_free(lwgeom);
+ PG_FREE_IF_COPY(geom, 1);
+ lwt_FreeTopology(topo);
+
+ if ( node_id == -1 ) {
+ /* should never reach this point, as lwerror would raise an exception */
+ SPI_finish();
+ PG_RETURN_NULL();
+ }
+
+ SPI_finish();
+ PG_RETURN_INT32(node_id);
+}
-- if near the point there are two or more edges it throw an exception.
--
CREATE OR REPLACE FUNCTION topology.GetEdgeByPoint(atopology varchar, apoint geometry, tol1 float8)
- RETURNS int
-AS
-$$
-DECLARE
- sql text;
- idedge int;
-BEGIN
- --
- -- Atopology and apoint are required
- --
- IF atopology IS NULL OR apoint IS NULL THEN
- RAISE EXCEPTION 'Invalid null argument';
- END IF;
-
- --
- -- Apoint must be a point
- --
- IF substring(geometrytype(apoint), 1, 5) != 'POINT'
- THEN
- RAISE EXCEPTION 'Node geometry must be a point';
- END IF;
-
- --
- -- Tolerance must be >= 0
- --
- IF tol1 < 0
- THEN
- RAISE EXCEPTION 'Tolerance must be >=0';
- END IF;
-
-
- if tol1 = 0 then
- sql := 'SELECT edge_id FROM '
- || quote_ident(atopology)
- || '.edge_data WHERE ST_Intersects(geom, $1)';
- else
- sql := 'SELECT edge_id FROM '
- || quote_ident(atopology)
- || '.edge_data WHERE ST_DWithin(geom, $1, $2)';
- end if;
-
- BEGIN
- EXECUTE sql INTO STRICT idedge USING apoint, tol1;
- EXCEPTION
- WHEN NO_DATA_FOUND THEN
- idedge = 0;
- WHEN TOO_MANY_ROWS THEN
- RAISE EXCEPTION 'Two or more edges found';
- END;
-
- RETURN idedge;
-
-END
-$$
-LANGUAGE 'plpgsql' STABLE STRICT;
+ RETURNS int AS
+ 'MODULE_PATHNAME', 'GetEdgeByPoint'
+ LANGUAGE 'c' STABLE STRICT;
--} GetEdgeByPoint
-- if near the point there are two or more faces it throw an exception.
--
CREATE OR REPLACE FUNCTION topology.GetFaceByPoint(atopology varchar, apoint geometry, tol1 float8)
- RETURNS int
-AS
-$$
-DECLARE
- sql text;
- idface int;
-BEGIN
-
- idface := -1;
-
- --
- -- Atopology and apoint are required
- --
- IF atopology IS NULL OR apoint IS NULL THEN
- RAISE EXCEPTION 'Invalid null argument';
- END IF;
-
- --
- -- Apoint must be a point
- --
- IF substring(geometrytype(apoint), 1, 5) != 'POINT'
- THEN
- RAISE EXCEPTION 'Node geometry must be a point';
- END IF;
-
- --
- -- Tolerance must be >= 0
- --
- IF tol1 < 0
- THEN
- RAISE EXCEPTION 'Tolerance must be >=0';
- END IF;
- --
- -- first test is to check if there is inside an mbr
- --
- if tol1 = 0 then
- sql := 'SELECT face_id FROM '
- || quote_ident(atopology)
- || '.face WHERE mbr && $1 LIMIT 1';
- else
- sql := 'SELECT face_id FROM '
- || quote_ident(atopology)
- || '.face as a WHERE ST_DWithin(mbr, $1, $2) LIMIT 1';
- end if;
-
- BEGIN
- EXECUTE sql INTO STRICT idface USING apoint, tol1;
- EXCEPTION
- WHEN NO_DATA_FOUND THEN
- idface = 0;
- END;
-
- if idface > 0 then
- --
- -- probably there is something so now check the exact test
- --
-
- if tol1 = 0 then
- sql := 'SELECT e.face_id FROM ('
- || 'SELECT d.face_id,ST_BuildArea(ST_Union(geom)) as geom FROM ('
- || 'SELECT b.edge_id as edge_id,b.left_face as face_id,b.geom as geom FROM '
- || quote_ident(atopology) || '.edge_data as b,'
- || '(SELECT a.face_id FROM '
- || quote_ident(atopology) || '.face as a '
- || 'WHERE ST_Intersects(a.mbr,$1)=true'
- || ') as c '
- || 'WHERE (b.left_face = c.face_id) '
- || ' UNION ALL '
- || 'SELECT b.edge_id as edge_id, b.right_face as face_id, b.geom as geom FROM '
- || quote_ident(atopology) || '.edge_data as b,'
- || '(SELECT a.face_id FROM '
- || quote_ident(atopology) || '.face as a '
- || 'WHERE ST_Intersects(a.mbr,$1)=true'
- || ') as c '
- || 'WHERE (b.right_face = c.face_id) '
- || ') as d '
- || 'GROUP BY face_id '
- || ') as e '
- || 'WHERE ST_Intersects(e.geom, $1)';
- else
- sql := 'SELECT e.face_id FROM ('
- || 'SELECT d.face_id,ST_BuildArea(ST_Union(geom)) as geom FROM ('
- || 'SELECT b.edge_id as edge_id,b.left_face as face_id,b.geom as geom FROM '
- || quote_ident(atopology) || '.edge_data as b,'
- || '(SELECT a.face_id FROM '
- || quote_ident(atopology) || '.face as a '
- || 'WHERE ST_DWithin(a.mbr,$1,$2)'
- || ') as c '
- || 'WHERE (b.left_face = c.face_id) '
- || ' UNION ALL '
- || 'SELECT b.edge_id as edge_id, b.right_face as face_id, b.geom as geom FROM '
- || quote_ident(atopology) || '.edge_data as b,'
- || '(SELECT a.face_id FROM '
- || quote_ident(atopology) || '.face as a '
- || 'WHERE ST_DWithin(a.mbr,$1,$2)'
- || ') as c '
- || 'WHERE (b.right_face = c.face_id) '
- || ') as d '
- || 'GROUP BY face_id '
- || ') as e '
- || 'WHERE ST_DWithin(e.geom, $1, $2)';
- end if;
-
- RAISE DEBUG ' ==> %',sql;
-
- BEGIN
- EXECUTE sql INTO STRICT idface USING apoint, tol1;
- EXCEPTION
- WHEN NO_DATA_FOUND THEN
- idface = 0;
- WHEN TOO_MANY_ROWS THEN
- RAISE EXCEPTION 'Two or more faces found';
- END;
-
- end if;
-
- RETURN idface;
-
-END
-$$
-LANGUAGE 'plpgsql' STABLE STRICT;
+ RETURNS int AS
+ 'MODULE_PATHNAME', 'GetFaceByPoint'
+ LANGUAGE 'c' STABLE STRICT;
--} GetFaceByPoint
-- if near the point there are two or more nodes it throw an exception.
--
CREATE OR REPLACE FUNCTION topology.GetNodeByPoint(atopology varchar, apoint geometry, tol1 float8)
- RETURNS int
-AS
-$$
-DECLARE
- sql text;
- idnode int;
-BEGIN
- --
- -- Atopology and apoint are required
- --
- IF atopology IS NULL OR apoint IS NULL THEN
- RAISE EXCEPTION 'Invalid null argument';
- END IF;
-
- --
- -- Apoint must be a point
- --
- IF substring(geometrytype(apoint), 1, 5) != 'POINT'
- THEN
- RAISE EXCEPTION 'Node geometry must be a point';
- END IF;
-
- --
- -- Tolerance must be >= 0
- --
- IF tol1 < 0
- THEN
- RAISE EXCEPTION 'Tolerance must be >=0';
- END IF;
-
-
- if tol1 = 0 then
- sql := 'SELECT node_id FROM '
- || quote_ident(atopology)
- || '.node WHERE ST_Intersects(geom, $1)';
- else
- sql := 'SELECT node_id FROM '
- || quote_ident(atopology)
- || '.node WHERE ST_DWithin(geom, $1, $2)';
- end if;
-
- BEGIN
- EXECUTE sql INTO STRICT idnode USING apoint, tol1;
- EXCEPTION
- WHEN NO_DATA_FOUND THEN
- idnode = 0;
- WHEN TOO_MANY_ROWS THEN
- RAISE EXCEPTION 'Two or more nodes found';
- END;
-
- RETURN idnode;
-
-END
-$$
-LANGUAGE 'plpgsql' STABLE STRICT;
+ RETURNS int AS
+ 'MODULE_PATHNAME', 'GetNodeByPoint'
+ LANGUAGE 'c' STABLE STRICT;
--} GetNodeByPoint
projects / postgis / commitdiff