diff --git a/postgis/Makefile.in b/postgis/Makefile.in index e9fa5c1c0..06a41ec48 100644 --- a/postgis/Makefile.in +++ b/postgis/Makefile.in @@ -40,6 +40,7 @@ PG_OBJS=lwgeom_pg.o \ lwgeom_geos.o \ lwgeom_geos_prepared.o \ lwgeom_geos_clean.o \ + lwgeom_geos_split.o \ lwgeom_export.o \ lwgeom_in_gml.o \ lwgeom_in_kml.o \ diff --git a/postgis/lwgeom_geos_split.c b/postgis/lwgeom_geos_split.c new file mode 100644 index 000000000..b8ed1f99c --- /dev/null +++ b/postgis/lwgeom_geos_split.c @@ -0,0 +1,186 @@ +/********************************************************************** + * $Id: lwgeom_geos.c 5258 2010-02-17 21:02:49Z strk $ + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.refractions.net + * + * Copyright 2009-2010 Sandro Santilli + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU General Public Licence. See the COPYING file. + * + ********************************************************************** + * + * ST_SplitGeometry + * + * Split polygon by line, line by line, line by point. + * Returns components as a collection + * + * Author: Sandro Santilli + * + * Work done for Faunalia (http://www.faunalia.it) with fundings + * from Regione Toscana - Sistema Informativo per il Governo + * del Territorio e dell'Ambiente (RT-SIGTA). + * + * Thanks to the PostGIS community for sharing poly/line ideas [1] + * + * [1] http://trac.osgeo.org/postgis/wiki/UsersWikiSplitPolygonWithLineString + * + * + **********************************************************************/ + +#include "lwgeom_geos.h" +#include "funcapi.h" + +#include +#include + +/* #define POSTGIS_DEBUG_LEVEL 4 */ + +/* Initializes GEOS internally */ +static LWGEOM* lwline_split_by_point(LWLINE* lwgeom_in, LWPOINT* blade_in); +static LWGEOM* +lwline_split_by_point(LWLINE* lwline_in, LWPOINT* blade_in) +{ + GEOSGeometry* line; + GEOSGeometry* point; + double loc, dist; + int intersects; + POINT2D pt; + POINTARRAY* pa1; + POINTARRAY* pa2; + LWGEOM** components; + LWCOLLECTION* out; + + /* Possible outcomes: + * + * 1. The point is not on the line + * -> Return original geometry (cloned) + * 2. The point is on the line + * -> Return a multiline with all elements + */ + + getPoint2d_p(blade_in->point, 0, &pt); + loc = ptarray_locate_point(lwline_in->points, &pt, &dist); + + /* lwnotice("Location: %g -- Distance: %g", loc, dist); */ + + if ( dist > 0 ) { /* TODO: accept a tolerance ? */ + /* No intersection */ + return (LWGEOM*)lwline_clone(lwline_in); + } + + /* There is an intersection, let's get two substrings */ + + if ( loc == 0 || loc == 1 ) + { + /* Intersection is on the boundary or outside */ + return (LWGEOM*)lwline_clone(lwline_in); + } + + pa1 = ptarray_substring(lwline_in->points, 0, loc); + pa2 = ptarray_substring(lwline_in->points, loc, 1); + + /* TODO: check if either pa1 or pa2 are empty ? */ + + components = lwalloc(sizeof(LWGEOM*)*2); + components[0] = (LWGEOM*)lwline_construct(lwline_in->SRID, NULL, pa1); + components[1] = (LWGEOM*)lwline_construct(lwline_in->SRID, NULL, pa2); + + out = lwcollection_construct(MULTILINETYPE, lwline_in->SRID, + NULL, 2, components); + + /* That's all folks */ + + return (LWGEOM*)out; +} + +static LWGEOM* lwline_split(LWLINE* lwgeom_in, LWGEOM* blade_in); +static LWGEOM* +lwline_split(LWLINE* lwline_in, LWGEOM* blade_in) +{ + switch (TYPE_GETTYPE(blade_in->type)) + { + case POINTTYPE: + return lwline_split_by_point(lwline_in, (LWPOINT*)blade_in); + + case LINETYPE: + default: + lwerror("Splitting a Line by a %s is unsupported", + lwgeom_typename(TYPE_GETTYPE(blade_in->type))); + return NULL; + } + return NULL; +} + +static LWGEOM* lwpoly_split(LWPOLY* lwgeom_in, LWGEOM* blade_in); +static LWGEOM* +lwpoly_split(LWPOLY* lwgeom_in, LWGEOM* blade_in) +{ + switch (TYPE_GETTYPE(blade_in->type)) + { + case LINETYPE: + default: + lwerror("Splitting a Polygon by a %s is unsupported", + lwgeom_typename(TYPE_GETTYPE(blade_in->type))); + return NULL; + } + return NULL; +} + +static LWGEOM* lwgeom_split(LWGEOM* lwgeom_in, LWGEOM* blade_in); +static LWGEOM* +lwgeom_split(LWGEOM* lwgeom_in, LWGEOM* blade_in) +{ + switch (TYPE_GETTYPE(lwgeom_in->type)) + { + case LINETYPE: + return lwline_split((LWLINE*)lwgeom_in, blade_in); + + case POLYGONTYPE: + return lwpoly_split((LWPOLY*)lwgeom_in, blade_in); + + default: + lwerror("Splitting of %d geometries is unsupported", + lwgeom_typename(TYPE_GETTYPE(lwgeom_in->type))); + return NULL; + } + +} + + +Datum ST_SplitGeometry(PG_FUNCTION_ARGS); +PG_FUNCTION_INFO_V1(ST_SplitGeometry); +Datum ST_SplitGeometry(PG_FUNCTION_ARGS) +{ +#if 0 && POSTGIS_GEOS_VERSION < 33 + elog(ERROR, "You need GEOS-3.3.0 or up for ST_CleanGeometry"); + PG_RETURN_NULL(); +#else /* POSTGIS_GEOS_VERSION >= 33 */ + + PG_LWGEOM *in, *blade_in, *out; + LWGEOM *lwgeom_in, *lwblade_in, *lwgeom_out; + + in = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + lwgeom_in = lwgeom_deserialize(SERIALIZED_FORM(in)); + + blade_in = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + lwblade_in = lwgeom_deserialize(SERIALIZED_FORM(blade_in)); + + lwgeom_out = lwgeom_split(lwgeom_in, lwblade_in); + if ( ! lwgeom_out ) { + PG_FREE_IF_COPY(in, 0); + PG_FREE_IF_COPY(blade_in, 1); + PG_RETURN_NULL(); + } + + out = pglwgeom_serialize(lwgeom_out); + + PG_FREE_IF_COPY(in, 0); + PG_FREE_IF_COPY(blade_in, 1); + + PG_RETURN_POINTER(out); + +#endif /* POSTGIS_GEOS_VERSION >= 33 */ +} + diff --git a/postgis/postgis.sql.in.c b/postgis/postgis.sql.in.c index 598e1c42d..7d143bf81 100644 --- a/postgis/postgis.sql.in.c +++ b/postgis/postgis.sql.in.c @@ -4082,6 +4082,26 @@ CREATE OR REPLACE FUNCTION ST_CleanGeometry(geometry) LANGUAGE 'C' IMMUTABLE STRICT COST 100; +-------------------------------------------------------------------------------- +-- ST_SplitGeometry +-------------------------------------------------------------------------------- + +-- ST_SplitGeometry(in geometry, blade geometry) +-- +-- Split a geometry in parts after cutting it with given blade. +-- Returns a collection containing all parts. +-- +-- Note that multi-part geometries will be returned exploded, +-- no matter relation to blade. +-- +-- Availability: 2.0.0 +-- +CREATE OR REPLACE FUNCTION ST_SplitGeometry(geometry, geometry) + RETURNS geometry + AS 'MODULE_PATHNAME', 'ST_SplitGeometry' + LANGUAGE 'C' IMMUTABLE STRICT + COST 100; + -------------------------------------------------------------------------------- -- Aggregates and their supporting functions -------------------------------------------------------------------------------- diff --git a/regress/Makefile.in b/regress/Makefile.in index 48ec0791d..94aee7fca 100644 --- a/regress/Makefile.in +++ b/regress/Makefile.in @@ -68,7 +68,8 @@ TESTS = \ dumppoints \ wmsservers \ tickets \ - remove_repeated_points + remove_repeated_points \ + split # Styled buffer only if GEOS >= 3.2 ifeq ($(shell expr $(POSTGIS_GEOS_VERSION) ">=" 32),1) diff --git a/regress/split.sql b/regress/split.sql new file mode 100644 index 000000000..7c9005d84 --- /dev/null +++ b/regress/split.sql @@ -0,0 +1,9 @@ +-- Point on line +select '1',st_astext(st_splitgeometry('LINESTRING(0 0, 10 0)', 'POINT(5 0)')); + +-- Point on line boundary +select '2',st_astext(st_splitgeometry('LINESTRING(0 0, 10 0)', 'POINT(10 0)')); + +-- Point off line +select '3',st_astext(st_splitgeometry('LINESTRING(0 0, 10 0)', 'POINT(5 1)')); + diff --git a/regress/split_expected b/regress/split_expected new file mode 100644 index 000000000..f6cc7bec9 --- /dev/null +++ b/regress/split_expected @@ -0,0 +1,3 @@ +1|MULTILINESTRING((0 0,5 0),(5 0,10 0)) +2|LINESTRING(0 0,10 0) +3|LINESTRING(0 0,10 0)