Calculating path sinuosity in postgis

Question

Sinuosity is a fairly straightforward concept. How can this be calculated in postgis given a series of linestrings that represent the path?

Answer 1

Using the postgis documentation as a starting point, I have come up with the below solution, which is a function that takes a set of linestrings, and a segment length and returns a set of segments with associated sinuosity.

I'm relatively new to postgis, so any advice or criticism of the function or method welcomed.

create type SinuosityRecord as (id int, the_geom geometry, actualdistance double precision, lineardistance double precision, sinuositypercentage double precision, sinuosityindex double precision);

/*
    GetSinuosityIndex
      sql: expects an sql statement with a linestring geometry ordered sequentially
        linestringfieldname: the name of the field in the sql statement that represents the linestring geometry
        segmentlength: the length of segments you wish to use to calculate the sinuosity
        returns: 
            id : an id representing each segment
            the_geom: a linestring geometry for each segment
            actualdistance:  the actual length of the segment
            lineardistance: the linear distance between the segment start and end
            sinuositypercentage: the percentage of actual length of the difference in actual length and lineardistance
            sinuosityindex: the ratio of actual length to linear distance
*/
create or replace function GetSinuosityIndex ( sql text, linestringfieldname text, segmentlength numeric ) returns setof sinuosityrecord as
$BODY$
declare
    r sinuosityrecord;                          -- used to return the result set
    totallength double precision;       -- used to store the total length of the path
    serieslength int;                               -- used to store the number of segments that will be needed
begin

    -- calculates the total length of the path
    execute $$
    SELECT
        ST_LENGTH(the_geom) as totallength
    FROM
    (
        SELECT 
            ST_LINEMERGE(s.$$ || linestringfieldname || $$) AS the_geom
        FROM
            (
                $$ || sql || $$
            ) s
    ) x $$
    INTO
        totallength;

    -- determine the total number of segements that will be required
    serieslength := ceiling(totallength / segmentlength) + 1;

    for r in execute $$
    SELECT
        x.n as id,
        x.the_geom, 
        ST_LENGTH(x.the_geom) AS actualdistance, 
        ST_DISTANCE(ST_STARTPOINT(x.the_geom), ST_ENDPOINT(x.the_geom)) AS lineardistance, 
        (ST_LENGTH(x.the_geom) - ST_DISTANCE(ST_STARTPOINT(x.the_geom), ST_ENDPOINT(x.the_geom))) / ST_LENGTH(x.the_geom) * 100::double precision AS sinuositypercentage, 
        ST_LENGTH(x.the_geom) / ST_DISTANCE(ST_STARTPOINT(x.the_geom), ST_ENDPOINT(x.the_geom)) AS sinuosityindex
    FROM 
        (
            SELECT 
                n.n, 
                ST_LINE_SUBSTRING
                (
                    t.the_geom,
                    ($$ || segmentlength || $$ * n.n::numeric)::double precision / t.length, 
                    CASE
                        WHEN ($$ || segmentlength || $$ * (n.n + 1)::numeric)::double precision < t.length THEN ($$ || segmentlength || $$ * (n.n + 1)::numeric)::double precision / t.length
                        ELSE 1::double precision
                    END
                ) AS the_geom
            FROM 
            (
                SELECT 
                    ST_LINEMERGE(s.$$ || linestringfieldname || $$) AS the_geom, 
                    ST_LENGTH(s.$$ || linestringfieldname || $$) AS length
                FROM
                (
                    $$ || sql || $$
                ) s
        ) t
        CROSS JOIN 
        generate_series(0, $$ || serieslength || $$) n(n)
    WHERE
    ((n.n::numeric * $$ || segmentlength || $$)::double precision / t.length) < 1::double precision) x $$

    loop
        return next r;
    end loop;

    return;

end
$BODY$
language plpgsql;

An example of how this can then be called is:

SELECT
    *
FROM
    GetSinuosityIndex 
    ( 
        'SELECT 
            st_makeline(network_table.point_field) AS the_geom
        FROM
            network_table
            '
        ,'the_geom'
        ,600 
    )