lwt_GetFaceByPoint X
lwt_GetNodeByPoint X
lwt_GetEdgeByPoint X
-lwt_TopoGeo_AddPoint
+lwt_TopoGeo_AddPoint X
lwt_TopoGeo_AddLineString
lwt_TopoGeo_AddPolygon
*
*******************************************************************/
-
-/** Element of a TopoGeometry */
-typedef struct LWT_TOPOELEMENT_T {
- LWT_ELEMID id; /* primitive or topogeometry id */
- /** primitive type (0:node, 1:edge, 2:face) or layer id */
- int type;
-} LWT_TOPOELEMENT;
-
/**
* Topology errors type
*/
*/
LWT_ELEMID lwt_GetFaceByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol);
+
/*******************************************************************
*
* Topology population (non-ISO)
* @param point the point to add
* @param tol snap tolerance, the topology tolerance will be used if 0
*
- * @return identifier of added (or pre-existing) node
+ * @return identifier of added (or pre-existing) node or -1 on error
+ * (liblwgeom error handler will be invoked with error message)
*/
LWT_ELEMID lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol);
* @param line the line to add
* @param tol snap tolerance, the topology tolerance will be used if 0
* @param nedges output parameter, will be set to number of edges the
- * line was split into
+ * line was split into, or -1 on error
+ * (liblwgeom error handler will be invoked with error message)
*
* @return an array of <nedges> edge identifiers that sewed togheter
* will build up the input linestring (after snapping). Caller
- * will need to free the array using lwfree()
+ * will need to free the array using lwfree(), if not null.
*/
LWT_ELEMID* lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol,
int* nedges);
* @param poly the polygon to add
* @param tol snap tolerance, the topology tolerance will be used if 0
* @param nfaces output parameter, will be set to number of faces the
- * polygon was split into
+ * polygon was split into, or -1 on error
+ * (liblwgeom error handler will be invoked with error message)
*
* @return an array of <nfaces> face identifiers that sewed togheter
* will build up the input polygon (after snapping). Caller
- * will need to free the array using lwfree()
+ * will need to free the array using lwfree(), if not null.
*/
LWT_ELEMID* lwt_AddPolygon(LWT_TOPOLOGY* topo, LWPOLY* point, double tol,
int* nfaces);
-/**
- * Adds a geometry to the topology
- *
- * The given geometry will snap to existing nodes or
- * edges within given tolerance.
- * Existing edges or faces may be split by the operation.
- *
- * @param topo the topology to operate on
- * @param geom the geometry to add
- * @param tol snap tolerance, the topology tolerance will be used if 0
- * @param nelems output parameter, will be set to number of primitive
- * elements the geometry was split into.
- *
- * @return an array of <nelems> topoelements that taken togheter
- * will build up the input geometry (after snapping). Caller
- * will need to free the array using lwfree()
- *
- * @see lwt_CreateTopoGeom
- */
-LWT_TOPOELEMENT* lwt_AddGeometry(LWT_TOPOLOGY* topo, LWGEOM* geom,
- double tol, int* nelems);
-
/*******************************************************************
*
* ISO signatures here
* (coincident nodes, crossing edges,
* actual face containement)
*
- * @return ID of the newly added node
+ * @return ID of the newly added node, or -1 on error
+ * (liblwgeom error handler will be invoked with error message)
*
*/
LWT_ELEMID lwt_AddIsoNode(LWT_TOPOLOGY* topo, LWT_ELEMID face,
#include <stdio.h>
#include <inttypes.h> /* for PRId64 */
#include <errno.h>
+#include <math.h>
#ifdef WIN32
# define LWTFMT_ELEMID "lld"
return id;
}
+
+/* Return the smallest delta that can perturbate
+ * the maximum absolute value of a geometry ordinate
+ */
+static double
+_lwt_minTolerance( LWGEOM *g )
+{
+ const GBOX* gbox;
+ double max;
+
+ gbox = lwgeom_get_bbox(g);
+ if ( ! gbox ) return 0; /* empty */
+ max = FP_ABS(gbox->xmin);
+ if ( max < FP_ABS(gbox->xmax) ) max = FP_ABS(gbox->xmax);
+ if ( max < FP_ABS(gbox->ymin) ) max = FP_ABS(gbox->ymin);
+ if ( max < FP_ABS(gbox->ymax) ) max = FP_ABS(gbox->ymax);
+
+ return 3.6 * pow(10, - ( 15 - log(max?max:1.0) ) );
+}
+
+#define _LWT_MINTOLERANCE( topo, geom ) ( \
+ topo->precision ? topo->precision : _lwt_minTolerance(geom) )
+
+LWT_ELEMID
+lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol)
+{
+ int num, i;
+ double mindist;
+ LWT_ISO_NODE *nodes;
+ LWT_ISO_EDGE *edges;
+ LWGEOM *pt = lwpoint_as_lwgeom(point);
+ int flds;
+ LWT_ELEMID id = 0;
+
+ if ( ! tol ) tol = _LWT_MINTOLERANCE( topo, pt );
+
+ /*
+ -- 1. Check if any existing node is closer than the given precision
+ -- and if so pick the closest
+ TODO: use WithinBox2D
+ */
+ flds = LWT_COL_NODE_NODE_ID|LWT_COL_NODE_GEOM;
+ nodes = lwt_be_getNodeWithinDistance2D(topo, point, 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_NODE *n = &(nodes[i]);
+ LWGEOM *g = lwpoint_as_lwgeom(n->geom);
+ double dist = lwgeom_mindistance2d(g, pt);
+ if ( dist >= tol ) continue; /* must be closer than tolerated */
+ if ( ! id || dist < mindist )
+ {
+ id = n->node_id;
+ mindist = dist;
+ }
+ }
+ if ( id )
+ {
+ /* found an existing node */
+ if ( nodes ) _lwt_release_nodes(nodes, num);
+ return id;
+ }
+
+ initGEOS(lwnotice, lwgeom_geos_error);
+
+ /*
+ -- 2. Check if any existing edge falls within tolerance
+ -- and if so split it by a point projected on it
+ TODO: use WithinBox2D
+ */
+ flds = LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM;
+ edges = lwt_be_getEdgeWithinDistance2D(topo, point, tol, &num, flds, 1);
+ if ( num == -1 )
+ {
+ lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
+ return -1;
+ }
+ if ( num )
+ {{
+ /* The point is on or near an edge, split the edge */
+
+ LWT_ISO_EDGE *e = &(edges[0]);
+ LWGEOM *g = lwline_as_lwgeom(e->geom);
+ LWGEOM *prj;
+ int contains;
+ GEOSGeometry *prjg, *gg;
+
+ LWDEBUGF(1, "Splitting edge %" LWTFMT_ELEMID, e->edge_id);
+
+ /* project point to line, split edge by point */
+ prj = lwgeom_closest_point(g, pt);
+ if ( lwgeom_has_z(pt) )
+ {{
+ /*
+ -- This is a workaround for ClosestPoint lack of Z support:
+ -- http://trac.osgeo.org/postgis/ticket/2033
+ */
+ LWGEOM *tmp;
+ double z;
+ POINT4D p4d;
+ LWPOINT *prjpt;
+ /* add Z to "prj" */
+ tmp = lwgeom_force_3dz(prj);
+ prjpt = lwgeom_as_lwpoint(tmp);
+ getPoint4d_p(point->point, 0, &p4d);
+ z = p4d.z;
+ getPoint4d_p(prjpt->point, 0, &p4d);
+ p4d.z = z;
+ ptarray_set_point4d(prjpt->point, 0, &p4d);
+ lwgeom_free(prj);
+ prj = tmp;
+ }}
+ prjg = LWGEOM2GEOS(prj, 0);
+ if ( ! prjg ) {
+ lwgeom_free(prj);
+ _lwt_release_edges(edges, num);
+ lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
+ return -1;
+ }
+ gg = LWGEOM2GEOS(g, 0);
+ if ( ! gg ) {
+ lwgeom_free(prj);
+ _lwt_release_edges(edges, num);
+ GEOSGeom_destroy(prjg);
+ lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
+ return -1;
+ }
+ contains = GEOSContains(gg, prjg);
+ GEOSGeom_destroy(prjg);
+ GEOSGeom_destroy(gg);
+ if ( contains == 2 )
+ {
+ lwgeom_free(prj);
+ _lwt_release_edges(edges, num);
+ lwerror("GEOS exception on Contains: %s", lwgeom_geos_errmsg);
+ return -1;
+ }
+ if ( ! contains )
+ {{
+ double snaptol;
+ LWGEOM *snapedge;
+ LWLINE *snapline;
+ POINT4D p1, p2;
+
+ LWDEBUGF(1, "Edge %" LWTFMT_ELEMID
+ " does not contain projected point to it",
+ e->edge_id);
+
+ /*
+ -- The tolerance must be big enough for snapping to happen
+ -- and small enough to snap only to the projected point.
+ -- Unfortunately ST_Distance returns 0 because it also uses
+ -- a projected point internally, so we need another way.
+ */
+ snaptol = _lwt_minTolerance(prj);
+ snapedge = lwgeom_snap(g, prj, snaptol);
+ snapline = lwgeom_as_lwline(snapedge);
+
+ LWDEBUGF(1, "Edge snapped with tolerance %g", snaptol);
+
+ /* TODO: check if snapping did anything ? */
+#if POSTGIS_DEBUG_LEVEL > 0
+ {
+ size_t sz;
+ char *wkt1 = lwgeom_to_wkt(g, WKT_ISO, 15, &sz);
+ char *wkt2 = lwgeom_to_wkt(snapedge, WKT_ISO, 15, &sz);
+ LWDEBUGF(1, "Edge %s snapped became %s", wkt1, wkt2);
+ lwfree(wkt1);
+ lwfree(wkt2);
+ }
+#endif
+
+
+ /*
+ -- Snapping currently snaps the first point below tolerance
+ -- so may possibly move first point. See ticket #1631
+ */
+ getPoint4d_p(e->geom->points, 0, &p1);
+ getPoint4d_p(snapline->points, 0, &p2);
+ LWDEBUGF(1, "Edge first point is %g %g, "
+ "snapline first point is %g %g",
+ p1.x, p1.y, p2.x, p2.y);
+ if ( p1.x != p2.x || p1.y != p2.y )
+ {
+ LWDEBUG(1, "Snapping moved first point, re-adding it");
+ if ( LW_SUCCESS != ptarray_insert_point(snapline->points, &p1, 0) )
+ {
+ lwgeom_free(prj);
+ lwgeom_free(snapedge);
+ _lwt_release_edges(edges, num);
+ lwerror("GEOS exception on Contains: %s", lwgeom_geos_errmsg);
+ return -1;
+ }
+#if POSTGIS_DEBUG_LEVEL > 0
+ {
+ size_t sz;
+ char *wkt1 = lwgeom_to_wkt(g, WKT_ISO, 15, &sz);
+ LWDEBUGF(1, "Tweaked snapline became %s", wkt1);
+ lwfree(wkt1);
+ }
+#endif
+ }
+#if POSTGIS_DEBUG_LEVEL > 0
+ else {
+ LWDEBUG(1, "Snapping did not move first point");
+ }
+#endif
+
+ if ( -1 == lwt_ChangeEdgeGeom( topo, e->edge_id, snapline ) )
+ {
+ /* TODO: should have invoked lwerror already, leaking memory */
+ lwgeom_free(prj);
+ lwgeom_free(snapedge);
+ _lwt_release_edges(edges, num);
+ lwerror("lwt_ChangeEdgeGeom failed");
+ return -1;
+ }
+ lwgeom_free(snapedge);
+ }}
+#if POSTGIS_DEBUG_LEVEL > 0
+ else
+ {{
+ size_t sz;
+ char *wkt1 = lwgeom_to_wkt(g, WKT_ISO, 15, &sz);
+ char *wkt2 = lwgeom_to_wkt(prj, WKT_ISO, 15, &sz);
+ LWDEBUGF(1, "Edge %s contains projected point %s", wkt1, wkt2);
+ lwfree(wkt1);
+ lwfree(wkt2);
+ }}
+#endif
+
+ /* TODO: pass 1 as last argument (skipChecks) ? */
+ id = lwt_ModEdgeSplit( topo, e->edge_id, lwgeom_as_lwpoint(prj), 0 );
+ if ( -1 == id )
+ {
+ /* TODO: should have invoked lwerror already, leaking memory */
+ lwgeom_free(prj);
+ _lwt_release_edges(edges, num);
+ lwerror("lwt_ModEdgeSplit failed");
+ return -1;
+ }
+
+ lwgeom_free(prj);
+ _lwt_release_edges(edges, num);
+ }}
+ else
+ {
+ /* The point is isolated, add it as such */
+ /* TODO: pass 1 as last argument (skipChecks) ? */
+ id = lwt_AddIsoNode(topo, -1, point, 0);
+ if ( -1 == id )
+ {
+ /* should have invoked lwerror already, leaking memory */
+ lwerror("lwt_AddIsoNode failed");
+ return -1;
+ }
+ }
+
+ return id;
+}
+
+LWT_ELEMID*
+lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges)
+{
+ *nedges = -1;
+ lwerror("Not implemented yet");
+ return NULL;
+}
+
+LWT_ELEMID*
+lwt_AddPolygon(LWT_TOPOLOGY* topo, LWPOLY* point, double tol, int* nfaces)
+{
+ *nfaces = -1;
+ lwerror("Not implemented yet");
+ return NULL;
+}
* doring operations.
*/
bool data_changed;
+
+ int topoLoadFailMessageFlavor; /* 0:sql, 1:AddPoint */
};
LWT_BE_DATA be_data;
{
pfree(sqldata.data);
//cberror(be, "no topology named '%s' was found", name);
- cberror(be, "SQL/MM Spatial exception - invalid topology name");
+ if ( be->topoLoadFailMessageFlavor == 1 ) {
+ cberror(be, "No topology with name \"%s\" in topology.topology", name);
+ } else {
+ cberror(be, "SQL/MM Spatial exception - invalid topology name");
+ }
return NULL;
}
if ( SPI_processed > 1 )
* is valid for the whole backend lifetime */
old_context = MemoryContextSwitchTo( TopMemoryContext );
+ /* initialize backend data */
+ be_data.data_changed = false;
+ be_data.topoLoadFailMessageFlavor = 0;
+
/* register callbacks against liblwgeom-topo */
be_iface = lwt_CreateBackendIface(&be_data);
lwt_BackendIfaceRegisterCallbacks(be_iface, &be_callbacks);
SPI_finish();
PG_RETURN_INT32(node_id);
}
+
+/* TopoGeo_AddPoint(atopology, point, tolerance) */
+Datum TopoGeo_AddPoint(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(TopoGeo_AddPoint);
+Datum TopoGeo_AddPoint(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 ) {{
+ char buf[32];
+ char *ptr;
+ snprintf(buf, 32, "%s", lwtype_name(lwgeom_get_type(lwgeom)));
+ buf[31] = '\0';
+ ptr = buf;
+ while (*ptr) {
+ *ptr = toupper(*ptr);
+ ++ptr;
+ }
+ lwgeom_free(lwgeom);
+ PG_FREE_IF_COPY(geom, 1);
+ lwpgerror("Invalid geometry type (%s) passed to TopoGeo_AddPoint"
+ ", expected POINT", buf );
+ 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;
+
+ {
+ int pre = be_data.topoLoadFailMessageFlavor;
+ be_data.topoLoadFailMessageFlavor = 1;
+ topo = lwt_LoadTopology(be_iface, toponame);
+ be_data.topoLoadFailMessageFlavor = pre;
+ }
+ 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_AddPoint");
+ node_id = lwt_AddPoint(topo, pt, tol);
+ POSTGIS_DEBUG(1, "lwt_AddPoint 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);
+}
--
CREATE OR REPLACE FUNCTION topology.TopoGeo_AddPoint(atopology varchar, apoint geometry, tolerance float8 DEFAULT 0)
RETURNS int AS
-$$
-DECLARE
- id integer;
- rec RECORD;
- sql text;
- prj GEOMETRY;
- snapedge GEOMETRY;
- snaptol FLOAT8;
- tol FLOAT8;
- z FLOAT8;
-BEGIN
-
- -- 0. Check arguments
- IF geometrytype(apoint) != 'POINT' THEN
- RAISE EXCEPTION 'Invalid geometry type (%) passed to TopoGeo_AddPoint, expected POINT', geometrytype(apoint);
- END IF;
-
- -- Get tolerance, if 0 was given
- tol := COALESCE( NULLIF(tolerance, 0), topology._st_mintolerance(atopology, apoint) );
-
- -- 1. Check if any existing node is closer than the given precision
- -- and if so pick the closest
- sql := 'SELECT a.node_id FROM '
- || quote_ident(atopology)
- || '.node as a WHERE ST_DWithin(a.geom,$1,'
- || tol || ') AND ST_Distance($1, a.geom) < ' || tol
- || ' ORDER BY ST_Distance($1, a.geom) LIMIT 1';
-#ifdef POSTGIS_TOPOLOGY_DEBUG
- RAISE DEBUG '%', sql;
-#endif
- EXECUTE sql INTO id USING apoint;
- IF id IS NOT NULL THEN
- RETURN id;
- END IF;
-
-#ifdef POSTGIS_TOPOLOGY_DEBUG
- RAISE DEBUG 'No existing node within tolerance distance';
-#endif
-
- -- 2. Check if any existing edge falls within tolerance
- -- and if so split it by a point projected on it
- sql := 'SELECT a.edge_id, a.geom FROM '
- || quote_ident(atopology)
- || '.edge as a WHERE ST_DWithin(a.geom,$1,'
- || tol || ') ORDER BY ST_Distance($1, a.geom) LIMIT 1';
-#ifdef POSTGIS_TOPOLOGY_DEBUG
- RAISE DEBUG '%', sql;
-#endif
- EXECUTE sql INTO rec USING apoint;
- IF rec IS NOT NULL THEN
- -- project point to line, split edge by point
- prj := ST_ClosestPoint(rec.geom, apoint);
- -- This is a workaround for ClosestPoint lack of Z support:
- -- http://trac.osgeo.org/postgis/ticket/2033
- z := ST_Z(apoint);
- IF z IS NOT NULL THEN
- prj := ST_Translate(ST_Force_3DZ(prj), 0, 0, z); -- no ST_SetZ ...
- END IF;
-#ifdef POSTGIS_TOPOLOGY_DEBUG
- RAISE DEBUG 'Splitting edge % with closest point %', rec.edge_id, ST_AsText(prj);
-#endif
- IF NOT ST_Contains(rec.geom, prj) THEN
-#ifdef POSTGIS_TOPOLOGY_DEBUG
- RAISE DEBUG ' Snapping edge to contain closest point';
-#endif
- -- The tolerance must be big enough for snapping to happen
- -- and small enough to snap only to the projected point.
- -- Unfortunately ST_Distance returns 0 because it also uses
- -- a projected point internally, so we need another way.
- snaptol := topology._st_mintolerance(prj);
-#ifdef POSTGIS_TOPOLOGY_DEBUG
- RAISE DEBUG 'Tolerance for snapping to point % = %', ST_AsText(prj), snaptol;
-#endif
- snapedge := ST_Snap(rec.geom, prj, snaptol);
-
- -- Snapping currently snaps the first point below tolerance
- -- so may possibly move first point. See ticket #1631
- IF NOT ST_Equals(ST_StartPoint(rec.geom), ST_StartPoint(snapedge))
- THEN
-#ifdef POSTGIS_TOPOLOGY_DEBUG
- RAISE WARNING 'Snapping moved first edge vertex, fixing';
-#endif
- snapedge := ST_MakeLine(ST_StartPoint(rec.geom), snapedge);
- END IF;
-
-#ifdef POSTGIS_TOPOLOGY_DEBUG
- IF NOT ST_Contains(snapedge, prj) THEN -- or if equal ?
- RAISE WARNING 'Edge within % distance from node still does not contain the node after snapping to it with tolerance %', tol, snaptol;
- END IF;
-#endif
- PERFORM topology.ST_ChangeEdgeGeom(atopology, rec.edge_id, snapedge);
- END IF;
- id := topology.ST_ModEdgeSplit(atopology, rec.edge_id, prj);
- ELSE
-#ifdef POSTGIS_TOPOLOGY_DEBUG
- RAISE DEBUG 'No existing edge within tolerance distance';
-#endif
- id := topology.ST_AddIsoNode(atopology, NULL, apoint);
- END IF;
-
- RETURN id;
-END
-$$
-LANGUAGE 'plpgsql' VOLATILE;
+ 'MODULE_PATHNAME', 'TopoGeo_AddPoint'
+ LANGUAGE 'c' VOLATILE;
--} TopoGeo_AddPoint
--{
projects / postgis / commitdiff