Sinuosity is a fairly straightforward concept. How can this be calculated in postgis given a series of linestrings that represent the path?
Answer
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
)
:
