9

Can anyone provide an example SQL for smoothing linestrings from postgis table using Bezier curves or Iterative Averaging (SIA) algorithm?

6

I created a small, naive script which converts input LineStrings to CompoundCurves based on some heuristics.

What it does:

  • Cuts down sharp corners to create a visually more appealing results than the original data.
  • Uses plpgsql. No additional extensions required.
  • Accepts an optional "smoothing factor" between 0 and 100 besides a geometry.

What it doesn't do:

  • Processes MultiLineStrings. For any other geometry type, it simply returns the input.
  • Uses Z and M values. It simply drops them. Use this only for 2D cartographic purposes.
  • Creates mathematically correct results. The results are far from correct, and might be even visually unaesthetic in some cases (e.g. sharp corners). I didn't test it thoroughly. Always review the results!
  • Runs fast. I'm sure it can be rewritten to a far more optimal form.
  • Does real smoothing. There are much better algorithms (e.g. Chaiken or the ones mentioned in the question) to use for real smoothing. This answer targets people like me searching for a pure PostGIS approach automatically creating some kind of curved lines from real data.

The script:

CREATE OR REPLACE FUNCTION CreateCurve(geom geometry, percent int DEFAULT 40)
    RETURNS geometry AS
$$
DECLARE
    result text;
    p0 geometry;
    p1 geometry;
    p2 geometry;
    intp geometry;
    tempp geometry;
    geomtype text := ST_GeometryType(geom);
    factor double precision := percent::double precision / 200;
    i integer;
BEGIN
    IF percent < 0 OR percent > 100 THEN
        RAISE EXCEPTION 'Smoothing factor must be between 0 and 100';
    END IF;
    IF geomtype != 'ST_LineString' OR factor = 0 THEN
        RETURN geom;
    END IF;
    result := 'COMPOUNDCURVE((';
    p0 := ST_PointN(geom, 1);
    IF ST_NPoints(geom) = 2 THEN
        p1:= ST_PointN(geom, 2);
        result := result || ST_X(p0) || ' ' || ST_Y(p0) || ',' || ST_X(p1) || ' ' || ST_Y(p1) || '))';
    ELSE
        FOR i IN 2..(ST_NPoints(geom) - 1) LOOP
            p1 := ST_PointN(geom, i);
            p2 := ST_PointN(geom, i + 1);
            result := result || ST_X(p0) || ' ' || ST_Y(p0) || ',';
            tempp := ST_Line_Interpolate_Point(ST_MakeLine(p1, p0), factor);
            p0 := ST_Line_Interpolate_Point(ST_MakeLine(p1, p2), factor);
            intp := ST_Line_Interpolate_Point(
                ST_MakeLine(
                    ST_Line_Interpolate_Point(ST_MakeLine(p0, p1), 0.5),
                    ST_Line_Interpolate_Point(ST_MakeLine(tempp, p1), 0.5)
                ), 0.5);
            result := result || ST_X(tempp) || ' ' || ST_Y(tempp) || '),CIRCULARSTRING(' || ST_X(tempp) || ' ' || ST_Y(tempp) || ',' || ST_X(intp) || ' ' ||
            ST_Y(intp) || ',' || ST_X(p0) || ' ' || ST_Y(p0) || '),(';
        END LOOP;
        result := result || ST_X(p0) || ' ' || ST_Y(p0) || ',' || ST_X(p2) || ' ' || ST_Y(p2) || '))';
    END IF;
    RETURN ST_SetSRID(result::geometry, ST_SRID(geom));
END;
$$
LANGUAGE 'plpgsql' IMMUTABLE;

As it returns curves in a geometry type, if you would like to use it in a GIS such as QGIS, you have to wrap it into PostGIS functions converting them around. The intended use syntax is:

SELECT ST_AsText(ST_CurveToLine(CreateCurve(geom))) AS geom FROM linestringtable;
  • This was a lifesaver! Thank you for the script. It seems like Chaikin smoothing will become available as a function from postgis 2.5 onwards, which I'm looking forward to. – she_weeds May 15 '18 at 8:19
1

This is still an open issue in PostGIS (and other GIS tools) as stated in the book "PostGIS in Action" in chapter 2.2.6 "Curved geometries".

Here are some references to algorithms and code:

0

You can try to convert your linestrings to curves with ST_LineToCurve and then back to linestrings with ST_CurveToLine.

You can set the number of segments per quarter circle you want in ST_CurveToLine.

  • LineToCurve is built to handing the outputs of CurveToLine, not to extract curves from arbitrary input. – Paul Ramsey Sep 3 '13 at 19:31
  • @PaulRamsey would bezier smoothing be added in next Postgis versions? I've been thinking on something like this for instance: webhelp.esri.com/arcgisdesktop/9.2/… – Gery Jan 2 '14 at 10:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.