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


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
    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;
    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) || '))';
        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_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));

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

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:


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.

  • 1
    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

My take on smoothing - Smoothing Iterative Averaging (SIA) + ST_SimplifyPreserveTopology using PostGIS:

Keep in mind for good results: ss < si, sl = 0 if ss > 1

sl is my addition to algorithm to not change longer segments as I use it for smoothing out GPS accuracy issues in line tracks. After SIA I also use St_SimplifyPreserveTopology to remove tiny nuances and downsize number of vertices. I don't care about 3D and have only Line Strings. Shouldn't be hard to adapt if you use multiline strings and 3D.


update table_name set geom = ST_SimplifyPreserveTopology(ST_SIA(geom, 1, 3, 20), 1);


-- Smoothing Iterative Averaging (SIA)
-- Parameters:
-- Smoothing Sensitivity (ss) - positive integer
-- Smoothing Iteration (si) - positive integer
-- Smoothing Limit (sl) - 0+ numeric - maximum distance between neighbouring vertices, 0 = do not apply limit
create or replace function ST_SIA(geom geometry, ss int default 1, si int default 5, sl numeric default 0)
    returns geometry as $$
    p_prev geometry;
    pi geometry;
    p_next geometry;
    new_line geometry;
    s integer;
    i integer;
    i_prev integer;
    i_next integer;
    p_count integer;
    if ss < 1 or si < 1 or sl < 0 then
        raise exception 'Smoothing sensitivity and smoothing iteration must be positive integer! Smoothing Limit 0+!';
    end if;
    if ST_GeometryType(geom) != 'ST_LineString' or ST_NPoints(geom) <= ss*2+1 then
        return geom;
    end if;
    new_line := geom;
    p_count := st_npoints(geom);
    for s in 1..si loop
        for i in 2..p_count-1 loop
            i_prev := i-ss;
            i_next := i+ss;
            if i <= ss then
                i_prev := 1;
            end if;
            if i+ss > p_count then
                i_next := p_count;
            end if;
            p_prev := ST_PointN(geom, i_prev);
            pi := ST_PointN(geom, i);
            p_next := ST_PointN(geom, i_next);
            if sl = 0 or ST_Length(ST_MakeLine(p_prev, p_next)) < sl then
                new_line := ST_SetPoint(new_line, i-1, ST_Centroid(ST_MakeLine(array[p_prev,pi,p_next])));
            end if;
        end loop;
        geom := new_line;
    end loop;
    return geom;
$$ language plpgsql;

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.