Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question

1 Answer 1

up vote 5 down vote accepted

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 
    )
share|improve this answer
    
Welcome to our site! Thanks for sharing this. –  whuber Sep 7 '12 at 22:10

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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