pggeodb.nancy.inra.fr/db_robot - db_robot on pggeodb.nancy.inra.fr

Function: st_concavehull

Schema

public

Owner

albenard

Descriptions

args: geomA, target_percent, allow_holes=false - The concave hull of a geometry represents a possibly concave geometry that encloses all geometries within the set. You can think of it as shrink wrapping.

Options

Definition

CREATE OR REPLACE FUNCTION public.st_concavehull (
 param_geom public.geometry,
 param_pctconvex double precision,
 param_allow_holes boolean = false
)
RETURNS public.geometry AS
$span$
DECLARE
var_convhull geometry := public.ST_ConvexHull(param_geom);
var_param_geom geometry := param_geom;
var_initarea float := public.ST_Area(var_convhull);
var_newarea float := var_initarea;
var_div integer := 6;
var_tempgeom geometry;
var_tempgeom2 geometry;
var_cent geometry;
var_geoms geometry[4];
var_enline geometry;
var_resultgeom geometry;
var_atempgeoms geometry[];
var_buf float := 1;
BEGIN
-- We start with convex hull as our base
var_resultgeom := var_convhull;

IF param_pctconvex = 1 THEN
return var_resultgeom;
ELSIF public.ST_GeometryType(var_param_geom) = 'ST_Polygon' THEN -- it is as concave as it is going to get
IF param_allow_holes THEN -- leave the holes
RETURN var_param_geom;
ELSE -- remove the holes
var_resultgeom := public.ST_MakePolygon(public.ST_ExteriorRing(var_param_geom));
RETURN var_resultgeom;
END IF;
END IF;
IF public.ST_Dimension(var_resultgeom) > 1 AND param_pctconvex BETWEEN 0 and 0.98 THEN
-- get linestring that forms envelope of geometry
var_enline := public.ST_Boundary(public.ST_Envelope(var_param_geom));
var_buf := public.ST_Length(var_enline)/1000.0;
IF public.ST_GeometryType(var_param_geom) = 'ST_MultiPoint' AND public.ST_NumGeometries(var_param_geom) BETWEEN 4 and 200 THEN
-- we make polygons out of points since they are easier to cave in.
-- Note we limit to between 4 and 200 points because this process is slow and gets quadratically slow
var_buf := sqrt(public.ST_Area(var_convhull)*0.8/(public.ST_NumGeometries(var_param_geom)*public.ST_NumGeometries(var_param_geom)));
var_atempgeoms := ARRAY(SELECT geom FROM public.ST_DumpPoints(var_param_geom));
-- 5 and 10 and just fudge factors
var_tempgeom := public.ST_Union(ARRAY(SELECT geom
FROM (
-- fuse near neighbors together
SELECT DISTINCT ON (i) i,  public.ST_Distance(var_atempgeoms[i],var_atempgeoms[j]), public.ST_Buffer(public.ST_MakeLine(var_atempgeoms[i], var_atempgeoms[j]) , var_buf*5, 'quad_segs=3') As geom
FROM generate_series(1,array_upper(var_atempgeoms, 1)) As i
INNER JOIN generate_series(1,array_upper(var_atempgeoms, 1)) As j
ON (
NOT public.ST_Intersects(var_atempgeoms[i],var_atempgeoms[j])
AND public.ST_DWithin(var_atempgeoms[i],var_atempgeoms[j], var_buf*10)
)
UNION ALL
-- catch the ones with no near neighbors
SELECT i, 0, public.ST_Buffer(var_atempgeoms[i] , var_buf*10, 'quad_segs=3') As geom
FROM generate_series(1,array_upper(var_atempgeoms, 1)) As i
LEFT JOIN generate_series(ceiling(array_upper(var_atempgeoms,1)/2)::integer,array_upper(var_atempgeoms, 1)) As j
ON (
NOT public.ST_Intersects(var_atempgeoms[i],var_atempgeoms[j])
AND public.ST_DWithin(var_atempgeoms[i],var_atempgeoms[j], var_buf*10)
)
WHERE j IS NULL
ORDER BY 1, 2
) As foo ) );
IF public.ST_IsValid(var_tempgeom) AND public.ST_GeometryType(var_tempgeom) = 'ST_Polygon' THEN
var_tempgeom := public.ST_ForceSFS(public.ST_Intersection(var_tempgeom, var_convhull));
IF param_allow_holes THEN
var_param_geom := var_tempgeom;
ELSE
var_param_geom := public.ST_MakePolygon(public.ST_ExteriorRing(var_tempgeom));
END IF;
return var_param_geom;
ELSIF public.ST_IsValid(var_tempgeom) THEN
var_param_geom := public.ST_ForceSFS(public.ST_Intersection(var_tempgeom, var_convhull));
END IF;
END IF;

IF public.ST_GeometryType(var_param_geom) = 'ST_Polygon' THEN
IF NOT param_allow_holes THEN
var_param_geom := public.ST_MakePolygon(public.ST_ExteriorRing(var_param_geom));
END IF;
return var_param_geom;
END IF;
           var_cent := public.ST_Centroid(var_param_geom);
           IF (public.ST_XMax(var_enline) - public.ST_XMin(var_enline) ) > var_buf AND (public.ST_YMax(var_enline) - public.ST_YMin(var_enline) ) > var_buf THEN
                   IF public.ST_Dwithin(public.ST_Centroid(var_convhull) , public.ST_Centroid(public.ST_Envelope(var_param_geom)), var_buf/2) THEN
               -- If the geometric dimension is > 1 and the object is symettric (cutting at centroid will not work -- offset a bit)
                       var_cent := public.ST_Translate(var_cent, (public.ST_XMax(var_enline) - public.ST_XMin(var_enline))/1000,  (public.ST_YMAX(var_enline) - public.ST_YMin(var_enline))/1000);
                   ELSE
                       -- uses closest point on geometry to centroid. I can't explain why we are doing this
                       var_cent := public.ST_ClosestPoint(var_param_geom,var_cent);
                   END IF;
                   IF public.ST_DWithin(var_cent, var_enline,var_buf) THEN
                       var_cent := public.ST_centroid(public.ST_Envelope(var_param_geom));
                   END IF;
                   -- break envelope into 4 triangles about the centroid of the geometry and returned the clipped geometry in each quadrant
                   FOR i in 1 .. 4 LOOP
                      var_geoms[i] := public.ST_MakePolygon(public.ST_MakeLine(ARRAY[public.ST_PointN(var_enline,i), public.ST_PointN(var_enline,i+1), var_cent, public.ST_PointN(var_enline,i)]));
                      var_geoms[i] := public.ST_ForceSFS(public.ST_Intersection(var_param_geom, public.ST_Buffer(var_geoms[i],var_buf)));
                      IF public.ST_IsValid(var_geoms[i]) THEN

                      ELSE
                           var_geoms[i] := public.ST_BuildArea(public.ST_MakeLine(ARRAY[public.ST_PointN(var_enline,i), public.ST_PointN(var_enline,i+1), var_cent, public.ST_PointN(var_enline,i)]));
                      END IF;
                   END LOOP;
                   var_tempgeom := public.ST_Union(ARRAY[public.ST_ConvexHull(var_geoms[1]), public.ST_ConvexHull(var_geoms[2]) , public.ST_ConvexHull(var_geoms[3]), public.ST_ConvexHull(var_geoms[4])]);
                   --RAISE NOTICE 'Curr vex % ', public.ST_AsText(var_tempgeom);
                   IF public.ST_Area(var_tempgeom) <= var_newarea AND public.ST_IsValid(var_tempgeom)  THEN --AND public.ST_GeometryType(var_tempgeom) ILIKE '%Polygon'

                       var_tempgeom := public.ST_Buffer(public.ST_ConcaveHull(var_geoms[1],least(param_pctconvex + param_pctconvex/var_div),true),var_buf, 'quad_segs=2');
                       FOR i IN 1 .. 4 LOOP
                           var_geoms[i] := public.ST_Buffer(public.ST_ConcaveHull(var_geoms[i],least(param_pctconvex + param_pctconvex/var_div),true), var_buf, 'quad_segs=2');
                           IF public.ST_IsValid(var_geoms[i]) Then
                               var_tempgeom := public.ST_Union(var_tempgeom, var_geoms[i]);
                           ELSE
                               RAISE NOTICE 'Not valid % %', i, public.ST_AsText(var_tempgeom);
                               var_tempgeom := public.ST_Union(var_tempgeom, public.ST_ConvexHull(var_geoms[i]));
                           END IF;
                       END LOOP;

                       --RAISE NOTICE 'Curr concave % ', public.ST_AsText(var_tempgeom);
                       IF public.ST_IsValid(var_tempgeom) THEN
                           var_resultgeom := var_tempgeom;
                       END IF;
                       var_newarea := public.ST_Area(var_resultgeom);
                   ELSIF public.ST_IsValid(var_tempgeom) THEN
                       var_resultgeom := var_tempgeom;
                   END IF;

                   IF public.ST_NumGeometries(var_resultgeom) > 1  THEN
                       var_tempgeom := public._ST_ConcaveHull(var_resultgeom);
                       IF public.ST_IsValid(var_tempgeom) AND public.ST_GeometryType(var_tempgeom) ILIKE 'ST_Polygon' THEN
                           var_resultgeom := var_tempgeom;
                       ELSE
                           var_resultgeom := public.ST_Buffer(var_tempgeom,var_buf, 'quad_segs=2');
                       END IF;
                   END IF;
                   IF param_allow_holes = false THEN
                   -- only keep exterior ring since we do not want holes
                       var_resultgeom := public.ST_MakePolygon(public.ST_ExteriorRing(var_resultgeom));
                   END IF;
               ELSE
                   var_resultgeom := public.ST_Buffer(var_resultgeom,var_buf);
               END IF;
               var_resultgeom := public.ST_ForceSFS(public.ST_Intersection(var_resultgeom, public.ST_ConvexHull(var_param_geom)));
           ELSE
               -- dimensions are too small to cut
               var_resultgeom := public._ST_ConcaveHull(var_param_geom);
           END IF;
           RETURN var_resultgeom;
END;
$span$
LANGUAGE 'plpgsql'
IMMUTABLE
RETURNS NULL ON NULL INPUT
SECURITY INVOKER
COST 100;

COMMENT ON FUNCTION public.st_concavehull(param_geom public.geometry, param_pctconvex double precision, param_allow_holes boolean)
IS 'args: geomA, target_percent, allow_holes=false - The concave hull of a geometry represents a possibly concave geometry that encloses all geometries within the set. You can think of it as shrink wrapping.';

This file was generated with SQL Manager for PostgreSQL (www.pgsqlmanager.com) at 26/02/2014 11:51