From f6b715cbb801d91cfed1ce9d394ecac84a5e3e18 Mon Sep 17 00:00:00 2001 From: Sandro Santilli Date: Fri, 26 Sep 2014 08:51:07 +0000 Subject: [PATCH] Add ST_ClipByBox2D and lwgeom_clip_by_rect (#2939) Includes testcases and documentation Requires GEOS-3.5.0+ git-svn-id: http://svn.osgeo.org/postgis/trunk@12999 b70326c6-7e19-0410-871a-916f4a2858ee --- NEWS | 1 + doc/reference_processing.xml | 51 +++++++++++++++++++++++ liblwgeom/cunit/Makefile.in | 1 + liblwgeom/cunit/cu_clip_by_rect.c | 56 +++++++++++++++++++++++++ liblwgeom/cunit/cu_tester.c | 2 + liblwgeom/liblwgeom.h.in | 1 + liblwgeom/lwgeom_geos.c | 68 ++++++++++++++++++++++++++++++- postgis/lwgeom_geos.c | 65 ++++++++++++++++++++++++++++- postgis/postgis.sql.in | 8 ++++ regress/Makefile.in | 7 ++++ regress/clipbybox2d.sql | 13 ++++++ regress/clipbybox2d_expected | 6 +++ 12 files changed, 277 insertions(+), 2 deletions(-) create mode 100644 liblwgeom/cunit/cu_clip_by_rect.c create mode 100644 regress/clipbybox2d.sql create mode 100644 regress/clipbybox2d_expected diff --git a/NEWS b/NEWS index 1db1cea0d..4b35d1127 100644 --- a/NEWS +++ b/NEWS @@ -17,6 +17,7 @@ PostGIS 2.2.0 * New Features * + - #2939, ST_ClipByBox2D (Sandro Santilli / CartoDB) - #2247, ST_Retile and ST_CreateOverview: in-db raster overviews creation (Sandro Santilli / Vizzuality) - #899, -m shp2pgsql attribute names mapping -m switch diff --git a/doc/reference_processing.xml b/doc/reference_processing.xml index 52e01e2a5..99058c924 100644 --- a/doc/reference_processing.xml +++ b/doc/reference_processing.xml @@ -365,6 +365,57 @@ FROM (SELECT ST_Buffer( this function with standard OGC interface + + + + ST_ClipByBox2D + Returns the portion of a geometry falling within a rectangle. + + + + + + geometry ST_ClipByBox2D + geometry geom + box2d box + + + + + + Description + + +Clips a geometry by a 2D box in a fast but possibly dirty way. The output +geometry is not guaranteed to be valid (self-intersections for a polygon +may be introduced). Topologically invalid input geometries do not result +in exceptions being thrown. + + + Performed by the GEOS module. + Requires GEOS 3.5.0+ + + Availability: 2.2.0. + + + + + Examples + +-- Rely on implicit cast from geometry to box2d for the second parameter +SELECT ST_ClipByBox2D(the_geom, ST_MakeEnvelope(0,0,10,10)) FROM mytab; + + + + See Also + +, +, + + + + + ST_Collect diff --git a/liblwgeom/cunit/Makefile.in b/liblwgeom/cunit/Makefile.in index f89826227..fc3345934 100644 --- a/liblwgeom/cunit/Makefile.in +++ b/liblwgeom/cunit/Makefile.in @@ -34,6 +34,7 @@ OBJS= \ cu_tree.o \ cu_measures.o \ cu_node.o \ + cu_clip_by_rect.o \ cu_libgeom.o \ cu_split.o \ cu_stringbuffer.o \ diff --git a/liblwgeom/cunit/cu_clip_by_rect.c b/liblwgeom/cunit/cu_clip_by_rect.c new file mode 100644 index 000000000..f681eb413 --- /dev/null +++ b/liblwgeom/cunit/cu_clip_by_rect.c @@ -0,0 +1,56 @@ +/********************************************************************** + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.net + * + * Copyright (C) 2014 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. + * + **********************************************************************/ + +#include "CUnit/Basic.h" +#include "cu_tester.h" + +#include "liblwgeom.h" +#include "liblwgeom_internal.h" + +static void test_lwgeom_clip_by_rect(void) +{ +#if POSTGIS_GEOS_VERSION >= 35 + LWGEOM *in, *out; + const char *wkt; + char *tmp; + + /* Because i don't trust that much prior tests... ;) */ + cu_error_msg_reset(); + + wkt = "LINESTRING(0 0, 5 5, 10 0)"; + in = lwgeom_from_wkt(wkt, LW_PARSER_CHECK_NONE); + out = lwgeom_clip_by_rect(in, 5, 0, 10, 10); + tmp = lwgeom_to_ewkt(out); + /* printf("%s\n", tmp); */ + CU_ASSERT_STRING_EQUAL("LINESTRING(5 5,10 0)", tmp) + lwfree(tmp); lwgeom_free(out); lwgeom_free(in); + + wkt = "LINESTRING EMPTY"; + in = lwgeom_from_wkt(wkt, LW_PARSER_CHECK_NONE); + out = lwgeom_clip_by_rect(in, 5, 0, 10, 10); + tmp = lwgeom_to_ewkt(out); + /* printf("%s\n", tmp); */ + CU_ASSERT_STRING_EQUAL(wkt, tmp) + lwfree(tmp); lwgeom_free(out); lwgeom_free(in); + +#endif /* POSTGIS_GEOS_VERSION >= 35 */ +} + +/* +** Used by test harness to register the tests in this file. +*/ +CU_TestInfo clip_by_rect_tests[] = +{ + PG_TEST(test_lwgeom_clip_by_rect), + CU_TEST_INFO_NULL +}; +CU_SuiteInfo clip_by_rect_suite = {"clip_by_rect", NULL, NULL, clip_by_rect_tests}; diff --git a/liblwgeom/cunit/cu_tester.c b/liblwgeom/cunit/cu_tester.c index bed263426..3629957ed 100644 --- a/liblwgeom/cunit/cu_tester.c +++ b/liblwgeom/cunit/cu_tester.c @@ -26,6 +26,7 @@ extern CU_SuiteInfo print_suite; extern CU_SuiteInfo algorithms_suite; extern CU_SuiteInfo buildarea_suite; extern CU_SuiteInfo clean_suite; +extern CU_SuiteInfo clip_by_rect_suite; extern CU_SuiteInfo misc_suite; extern CU_SuiteInfo ptarray_suite; extern CU_SuiteInfo measures_suite; @@ -72,6 +73,7 @@ int main(int argc, char *argv[]) algorithms_suite, buildarea_suite, clean_suite, + clip_by_rect_suite, measures_suite, node_suite, wkt_out_suite, diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in index 976081d78..f5d535e6f 100644 --- a/liblwgeom/liblwgeom.h.in +++ b/liblwgeom/liblwgeom.h.in @@ -1904,6 +1904,7 @@ LWGEOM *lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2); LWGEOM *lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2); LWGEOM *lwgeom_symdifference(const LWGEOM* geom1, const LWGEOM* geom2); LWGEOM *lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2); +LWGEOM *lwgeom_clip_by_rect(const LWGEOM *geom1, double x0, double y0, double x1, double y1); /** * Snap vertices and segments of a geometry to another using a given tolerance. diff --git a/liblwgeom/lwgeom_geos.c b/liblwgeom/lwgeom_geos.c index 8ccc048fa..4eebe51ca 100644 --- a/liblwgeom/lwgeom_geos.c +++ b/liblwgeom/lwgeom_geos.c @@ -3,7 +3,7 @@ * PostGIS - Spatial Types for PostgreSQL * http://postgis.net * - * Copyright 2011-2012 Sandro Santilli + * Copyright 2011-2014 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. @@ -763,6 +763,72 @@ lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2) return result; } +LWGEOM * +lwgeom_clip_by_rect(const LWGEOM *geom1, double x0, double y0, double x1, double y1) +{ +#if POSTGIS_GEOS_VERSION < 35 + lwerror("The GEOS version this postgis binary " + "was compiled against (%d) doesn't support " + "'GEOSClipByRect' function (3.3.5+ required)", + POSTGIS_GEOS_VERSION); + return NULL; +#else /* POSTGIS_GEOS_VERSION >= 35 */ + LWGEOM *result ; + GEOSGeometry *g1, *g3 ; + int is3d ; + + /* A.Intersection(Empty) == Empty */ + if ( lwgeom_is_empty(geom1) ) + return lwgeom_clone(geom1); + + is3d = FLAGS_GET_Z(geom1->flags); + + initGEOS(lwnotice, lwgeom_geos_error); + + LWDEBUG(3, "clip_by_rect() START"); + + g1 = LWGEOM2GEOS(geom1); + if ( 0 == g1 ) /* exception thrown at construction */ + { + lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); + return NULL ; + } + + LWDEBUG(3, " constructed geometrys - calling geos"); + LWDEBUGF(3, " g1 = %s", GEOSGeomToWKT(g1)); + /*LWDEBUGF(3, "g1 is valid = %i",GEOSisvalid(g1)); */ + + g3 = GEOSClipByRect(g1,x0,y0,x1,y1); + GEOSGeom_destroy(g1); + + LWDEBUG(3, " clip_by_rect finished"); + + if (g3 == NULL) + { + lwerror("Error performing rectangular clipping: %s", + lwgeom_geos_errmsg); + return NULL; /* never get here */ + } + + LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ; + + result = GEOS2LWGEOM(g3, is3d); + GEOSGeom_destroy(g3); + + if (result == NULL) + { + lwerror("Error performing intersection: GEOS2LWGEOM: %s", + lwgeom_geos_errmsg); + return NULL ; /* never get here */ + } + + result->srid = geom1->srid; + + return result ; +#endif /* POSTGIS_GEOS_VERSION >= 35 */ +} + + /* ------------ BuildArea stuff ---------------------------------------------------------------------{ */ typedef struct Face_t { diff --git a/postgis/lwgeom_geos.c b/postgis/lwgeom_geos.c index b1b32e782..3084c92d7 100644 --- a/postgis/lwgeom_geos.c +++ b/postgis/lwgeom_geos.c @@ -3,7 +3,7 @@ * PostGIS - Spatial Types for PostgreSQL * http://postgis.net * - * Copyright 2009-2012 Sandro Santilli + * Copyright 2009-2014 Sandro Santilli * Copyright 2008 Paul Ramsey * Copyright 2001-2003 Refractions Research Inc. * @@ -1635,6 +1635,69 @@ Datum centroid(PG_FUNCTION_ARGS) PG_RETURN_POINTER(result); } +Datum ST_ClipByBox2d(PG_FUNCTION_ARGS); +PG_FUNCTION_INFO_V1(ST_ClipByBox2d); +Datum ST_ClipByBox2d(PG_FUNCTION_ARGS) +{ +#if POSTGIS_GEOS_VERSION < 35 + + lwerror("The GEOS version this PostGIS binary " + "was compiled against (%d) doesn't support " + "'GEOSClipByRect' function (3.5.0+ required)", + POSTGIS_GEOS_VERSION); + PG_RETURN_NULL(); + +#else /* POSTGIS_GEOS_VERSION >= 35 */ + + GSERIALIZED *geom1; + GSERIALIZED *result; + LWGEOM *lwgeom1, *lwresult ; + const GBOX *bbox1; + const GBOX *bbox2; + + geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + lwgeom1 = lwgeom_from_gserialized(geom1) ; + + bbox1 = lwgeom_get_bbox(lwgeom1); + if ( ! bbox1 ) { + /* empty clips to empty, no matter rect */ + lwgeom_free(lwgeom1); + PG_RETURN_POINTER(geom1); + } + + bbox2 = (GBOX *)PG_GETARG_POINTER(1); + + /* If bbox1 outside of bbox, return empty */ + if ( LW_FALSE == gbox_overlaps_2d(bbox1, bbox2) ) { + lwresult = lwgeom_construct_empty(lwgeom1->type, lwgeom1->srid, 0, 0); + lwgeom_free(lwgeom1); + PG_FREE_IF_COPY(geom1, 1); + result = geometry_serialize(lwresult) ; + lwgeom_free(lwresult) ; + PG_RETURN_POINTER(result); + } + + /* if bbox1 is covered by bbox, return lwgeom1 */ + if ( bbox1->xmax <= bbox2->xmax && bbox1->xmin >= bbox2->xmin && + bbox1->ymax <= bbox2->ymax && bbox1->ymin >= bbox2->ymin ) + { + lwgeom_free(lwgeom1); + PG_RETURN_POINTER(geom1); + } + + lwresult = lwgeom_clip_by_rect(lwgeom1, bbox2->xmin, bbox2->ymin, + bbox2->xmax, bbox2->ymax); + lwgeom_free(lwgeom1) ; + + result = geometry_serialize(lwresult) ; + lwgeom_free(lwresult) ; + + PG_FREE_IF_COPY(geom1, 0); + + PG_RETURN_POINTER(result); +#endif /* POSTGIS_GEOS_VERSION >= 35 */ +} + /*---------------------------------------------*/ diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in index 6e2976660..c882dd7f0 100644 --- a/postgis/postgis.sql.in +++ b/postgis/postgis.sql.in @@ -3110,6 +3110,14 @@ CREATE OR REPLACE FUNCTION ST_RemoveRepeatedPoints(geometry) LANGUAGE 'c' IMMUTABLE STRICT COST 100; +-- Requires GEOS >= 3.5.0 +-- Availability: 2.2.0 +CREATE OR REPLACE FUNCTION ST_ClipByBox2d(geom geometry, box box2d) + RETURNS geometry + AS 'MODULE_PATHNAME', 'ST_ClipByBox2d' + LANGUAGE 'c' IMMUTABLE STRICT + COST 50; + -------------------------------------------------------------------------------- -- ST_CleanGeometry / ST_MakeValid -------------------------------------------------------------------------------- diff --git a/regress/Makefile.in b/regress/Makefile.in index b8f27679e..e6c034c11 100644 --- a/regress/Makefile.in +++ b/regress/Makefile.in @@ -161,6 +161,13 @@ ifeq ($(shell expr $(POSTGIS_GEOS_VERSION) ">=" 34),1) delaunaytriangles endif +ifeq ($(shell expr $(POSTGIS_GEOS_VERSION) ">=" 34),1) + # GEOS-3.5 adds: + # ST_ClipByBox2d + TESTS += \ + clipbybox2d +endif + ifeq ($(HAVE_JSON),yes) # JSON-C adds: # ST_GeomFromGeoJSON() diff --git a/regress/clipbybox2d.sql b/regress/clipbybox2d.sql new file mode 100644 index 000000000..8c44e7a06 --- /dev/null +++ b/regress/clipbybox2d.sql @@ -0,0 +1,13 @@ +-- Overlap +SELECT ST_ClipByBox2d(ST_MakeEnvelope(0,0,10,10),ST_MakeEnvelope(5,5,20,20))::box2d; +-- Geom covers rect +SELECT ST_ClipByBox2d(ST_MakeEnvelope(0,0,10,10),ST_MakeEnvelope(5,5,8,8))::box2d; +-- Geom within rect +SELECT ST_ClipByBox2d(ST_MakeEnvelope(2,2,8,8),ST_MakeEnvelope(0,0,10,10))::box2d; +-- Multipoint with point inside, point outside and point on boundary +SELECT ST_AsText(ST_ClipByBox2d('MULTIPOINT(-1 -1, 0 0, 2 2)'::geometry,ST_MakeEnvelope(0,0,10,10))); +-- Invalid polygon (bow-tie) -- ST_Intersection throws an exception here +SELECT ST_AsText(ST_ClipByBox2d('POLYGON((0 0, 10 10, 0 10, 10 0, 0 0))', ST_MakeEnvelope(2,2,8,8))); +-- Invalid polygon (lineal self-intersection) -- ST_Intersection returns a collection +SELECT ST_AsText(ST_ClipByBox2d('POLYGON((0 0,5 4,5 6,0 10,10 10,5 6,5 4,10 0,0 0))', ST_MakeEnvelope(2,2,10,5))); + diff --git a/regress/clipbybox2d_expected b/regress/clipbybox2d_expected new file mode 100644 index 000000000..54b5dd273 --- /dev/null +++ b/regress/clipbybox2d_expected @@ -0,0 +1,6 @@ +BOX(5 5,10 10) +BOX(5 5,8 8) +BOX(2 2,8 8) +POINT(2 2) +POLYGON((2 2,8 2,2 8,8 8,2 2)) +POLYGON((2.5 2,5 4,5 5,5 4,7.5 2,2.5 2)) -- 2.50.1