From 3febad808557de504f5b68fc09cf92c15c1f95d8 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Mon, 6 Feb 2012 23:25:00 +0000 Subject: [PATCH] Have ST_Union aggregate use UnaryUnion from GEOS-3.0.0 (#922) git-svn-id: http://svn.osgeo.org/postgis/trunk@9056 b70326c6-7e19-0410-871a-916f4a2858ee --- postgis/lwgeom_geos.c | 187 ++++++++++++++++++++ topology/test/regress/legacy_query_expected | 2 +- 2 files changed, 188 insertions(+), 1 deletion(-) diff --git a/postgis/lwgeom_geos.c b/postgis/lwgeom_geos.c index fd424c777..80c04360f 100644 --- a/postgis/lwgeom_geos.c +++ b/postgis/lwgeom_geos.c @@ -250,7 +250,193 @@ Datum hausdorffdistancedensify(PG_FUNCTION_ARGS) */ PG_FUNCTION_INFO_V1(pgis_union_geometry_array); Datum pgis_union_geometry_array(PG_FUNCTION_ARGS) +#if POSTGIS_GEOS_VERSION >= 33 { + /* + ** For GEOS >= 3.3, use the new UnaryUnion functionality to merge the + ** terminal collection from the ST_Union aggregate + */ + Datum datum; + ArrayType *array; + + int is3d = LW_FALSE, gotsrid = LW_FALSE; + int nelems = 0, i = 0, geoms_size = 0, curgeom = 0; + + GSERIALIZED *gser_out = NULL; + + GEOSGeometry *g = NULL; + GEOSGeometry *g_union = NULL; + GEOSGeometry **geoms = NULL; + + int srid = SRID_UNKNOWN; + + size_t offset = 0; + bits8 *bitmap; + int bitmask; + int empty_type = 0; + + datum = PG_GETARG_DATUM(0); + + /* Null array, null geometry (should be empty?) */ + if ( (Pointer *)datum == NULL ) PG_RETURN_NULL(); + + array = DatumGetArrayTypeP(datum); + + /* How many things in our array? */ + nelems = ArrayGetNItems(ARR_NDIM(array), ARR_DIMS(array)); + + /* PgSQL supplies a bitmap of which array entries are null */ + bitmap = ARR_NULLBITMAP(array); + + /* Empty array? Null return */ + if ( nelems == 0 ) PG_RETURN_NULL(); + + /* One-element union is the element itself! */ + if ( nelems == 1 ) + { + /* If the element is a NULL then we need to handle it separately */ + if (bitmap && (*bitmap & 1) == 0) + PG_RETURN_NULL(); + else + PG_RETURN_POINTER((GSERIALIZED *)(ARR_DATA_PTR(array))); + } + + /* Ok, we really need GEOS now ;) */ + initGEOS(lwnotice, lwgeom_geos_error); + + /* + ** Collect the non-empty inputs and stuff them into a GEOS collection + */ + geoms_size = nelems; + geoms = palloc( sizeof(GEOSGeometry*) * geoms_size ); + + /* + ** We need to convert the array of GSERIALIZED into a GEOS collection. + ** First make an array of GEOS geometries. + */ + offset = 0; + bitmap = ARR_NULLBITMAP(array); + bitmask = 1; + for ( i = 0; i < nelems; i++ ) + { + /* Only work on non-NULL entries in the array */ + if ((bitmap && (*bitmap & bitmask) != 0) || !bitmap) + { + GSERIALIZED *gser_in = (GSERIALIZED *)(ARR_DATA_PTR(array)+offset); + + /* Check for SRID mismatch in array elements */ + if ( gotsrid ) + { + error_if_srid_mismatch(srid, gserialized_get_srid(gser_in)); + } + else + { + /* Initialize SRID/dimensions info */ + srid = gserialized_get_srid(gser_in); + is3d = gserialized_has_z(gser_in); + gotsrid = 1; + } + + /* Don't include empties in the union */ + if ( gserialized_is_empty(gser_in) ) + { + int gser_type = gserialized_get_type(gser_in); + if (gser_type > empty_type) + { + empty_type = gser_type; + } + } + else + { + g = (GEOSGeometry *)POSTGIS2GEOS(gser_in); + + /* Uh oh! Exception thrown at construction... */ + if ( ! g ) + { + lwerror("One of the geometries in the set " + "could not be converted to GEOS: %s", lwgeom_geos_errmsg); + PG_RETURN_NULL(); + } + + /* Ensure we have enough space in our storage array */ + if ( curgeom == geoms_size ) + { + geoms_size *= 2; + geoms = repalloc( geoms, sizeof(GEOSGeometry*) * geoms_size ); + } + + offset += INTALIGN(VARSIZE(gser_in)); + geoms[curgeom] = g; + curgeom++; + } + } + + /* Advance NULL bitmap */ + if (bitmap) + { + bitmask <<= 1; + if (bitmask == 0x100) + { + bitmap++; + bitmask = 1; + } + } + } + + /* + ** Take our GEOS geometries and turn them into a GEOS collection, + ** then pass that into cascaded union. + */ + if (curgeom > 0) + { + g = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, geoms, curgeom); + if ( ! g ) + { + lwerror("Could not create GEOS COLLECTION from geometry array: %s", lwgeom_geos_errmsg); + PG_RETURN_NULL(); + } + + g_union = GEOSUnaryUnion(g); + GEOSGeom_destroy(g); + if ( ! g_union ) + { + lwerror("GEOSUnaryUnion: %s", + lwgeom_geos_errmsg); + PG_RETURN_NULL(); + } + + GEOSSetSRID(g_union, srid); + gser_out = GEOS2POSTGIS(g_union, is3d); + GEOSGeom_destroy(g_union); + } + /* No real geometries in our array, any empties? */ + else + { + /* If it was only empties, we'll return the largest type number */ + if ( empty_type > 0 ) + { + PG_RETURN_POINTER(geometry_serialize(lwgeom_construct_empty(empty_type, srid, is3d, 0))); + } + /* Nothing but NULL, returns NULL */ + else + { + PG_RETURN_NULL(); + } + } + + if ( ! gser_out ) + { + /* Union returned a NULL geometry */ + PG_RETURN_NULL(); + } + + PG_RETURN_POINTER(gser_out); +} + +#else +{ +/* For GEOS < 3.3, use the old CascadedUnion function for polygons and + brute force two-by-two for other types. */ Datum datum; ArrayType *array; int is3d = 0; @@ -564,6 +750,7 @@ Datum pgis_union_geometry_array(PG_FUNCTION_ARGS) PG_RETURN_POINTER(result); } +#endif /* POSTGIS_GEOS_VERSION >= 33 */ /** * @example ST_UnaryUnion {@link #geomunion} SELECT ST_UnaryUnion( diff --git a/topology/test/regress/legacy_query_expected b/topology/test/regress/legacy_query_expected index e92a2402f..1f045a1ae 100644 --- a/topology/test/regress/legacy_query_expected +++ b/topology/test/regress/legacy_query_expected @@ -41,7 +41,7 @@ E20E19|6|MULTILINESTRING((21 6,21 14,21 22)) E25|7|MULTILINESTRING((9 35,13 35)) R1a|8|MULTILINESTRING((9 14,21 14,35 14)) S1S2|MULTIPOINT(21 14,35 14) -R1R2|MULTILINESTRING((36 38,38 35,41 34,42 33,45 32,47 28,50 28,52 32,57 33,57 36,59 39,61 38,62 41,47 42,45 40,41 40),(9 14,21 14,35 14)) +R1R2|MULTILINESTRING((9 14,21 14,35 14),(36 38,38 35,41 34,42 33,45 32,47 28,50 28,52 32,57 33,57 36,59 39,61 38,62 41,47 42,45 40,41 40)) R4|MULTILINESTRING((25 30,25 35)) P1P2|MULTIPOLYGON(((21 6,9 6,9 14,9 22,21 22,35 22,35 14,35 6,21 6))) P3P4|MULTIPOLYGON(((47 14,47 6,35 6,35 14,35 22,47 22,47 14)),((25 30,17 30,17 40,31 40,31 30,25 30))) -- 2.40.0