PostGIS: Creating a line inside a polygon

I love working with the PostGIS API.

It is a deep and wide river to swim in.

But sometime having read, researched, experimented, and not found the right solution one can be forced to roll one’s own so to speak to get the job done.
For the problem below, if someone can point out a way to create what is needed using standard PostGIS API calls, then let me know.

The problem is this.

The customer wants to create transects along 20m segments chopped from a set of river centrelines.

For each 20m segment, a centroid must be created (at 10m), and an average bearing for the segment calculated from the first and last point.

A transect (straight line) must then be formed at the centroid with a direction at right angles to the segment bearing.
The transect must then be fitted inside a floodplain polygon so that it starts and ends by touching the boundary of the local part of the floodplain polygon.
The length of the created transect is compared to the intersection of the transect with a river polygon to compute a ratio of river width to floodplain width.

The floodplain polygon is very convoluted and extensive.

Floodplain

To be clear, what is not required is a set of homogeneous transects create from an envelope that covers the whole floodplain (whole to the part solution).

Envelope Transects

The required solution creates transects that are specific to river linestring segments that vary from segment to segment.
And for each of these, the transect must fit perfectly within that part of the floodplain polygon within which it resides, so that the transect line has its starting point touch one boundary, and its ending point to touch another.
So, the length of the transect varies for each and ever transect and its local area.

Techniques that guess a transect length using a width greater than the polygon are problematic as they can cross different parts of the polygon than the local area.

Transect crossing more than one part

A solution therefore is needed that can compute/create a transect perfectly inside the local area.

Final (Localised) Transect

While I tried many approaches (hoping not to write code), in the end I wrote a function called ST_InsideLine whose signature looks like this:

 CREATE OR REPLACE FUNCTION spdba.ST_InsideLine(
     p_point           IN geometry,
     p_direction_start IN NUMERIC,
     p_direction_end   IN NUMERIC,
     p_polygon         IN geometry,
     p_dIncrement      IN NUMERIC DEFAULT 5.0
 )
 RETURNS GEOMETRY
 ...

The function’s parameters are:

  • p_point (geometry) — Starting point which must fall inside p_polygon
  • p_direction_start (numeric) — Whole circle bearing in degrees defining direction of line from p_point to the start of the line
  • p_direction_end (numeric) — Whole circle bearing in degrees defining direction of line from p_point to the end of the line
  • p_polygon (geometry) — Polygon for which the inside line must be fitted.
  • p_dIncrement (numeric) — Line increment (distance) for incrementally extending a line, and testing line to find its boundary intersection point.

The function works as follows:

This function creates a line that lies inside p_polygon but whose start and end points touch p_polygon’s boundary.

Line is generated from a starting point and two bearings.
The algorithm first generates a line from the supplied point to the line’s start point by extending and testing the line by p_dIncrement until it finds a p_polygon boundary.
After finding point for first half of line, the second line is generated using p_direction_end.
If p_direction_end is null or the same as p_direction_start, a default direction of p_direction_start – 180.0 is used.
The algorithm uses a stepping approach: it first creates a line at p_direction_start for p_dIncrement distance. If the line does not touch a p_polygon boundary point it increases the line length by p_dIncrement and tests again.
The stepping process continues until the line touches or crosses the boundary.
Once the two halves are created, they are unioned together and the resulting line returned.

Here is an picture of the solution for a test polygon.

Line Inside

The following example shows how, by varying the bearings so they are different from each other, the function can create “bent” linestrings. Here is an example:

 SELECT ST_AsText(
          ST_SnapToGrid(
            spdba.ST_InsideLine(
              ST_SetSrid(ST_Point(82,60),28355),
              0.0::NUMERIC,
              generate_series(10,350,10)::NUMERIC,
              ST_GeomFromText('POLYGON((8.003 66.926, 11.164 70.086, 13.692 70.929, 19.171 70.929, 23.385 70.508, 26.546 70.297, 33.078 71.983, 36.871 74.301, 43.824 75.776, 51.199 75.986, 59.206 74.511, 62.788 71.772, 64.685 70.719, 73.535 71.351, 78.592 69.244, 83.649 64.187, 84.913 62.501, 86.178 57.022, 85.756 53.019, 85.124 49.226, 86.81 45.433, 87.863 40.376, 89.338 37.215, 89.338 32.58, 87.653 27.522, 83.438 18.462, 81.12 15.933, 74.799 17.619, 77.538 25.205, 80.067 30.472, 80.488 37.215, 78.381 41.219, 75.22 53.229, 72.06 60.394, 62.999 63.133, 52.463 65.451, 46.353 66.926, 37.714 63.344, 29.496 62.501, 20.646 61.447, 14.114 62.922, 9.899 61.447, 3.157 63.765, 3.367 64.187, 8.003 66.926))',28355),
              5.0
        ),
            0.001
        )
        ) AS geom ;
 -- results
 --
 LINESTRING(82.875 64.961,82 60,82 65.836)
 LINESTRING(83.557 64.279,82 60,82 65.836)
 LINESTRING(84.083 63.608,82 60,82 65.836)
 LINESTRING(84.529 63.014,82 60,82 65.836)
 LINESTRING(84.924 62.453,82 60,82 65.836)
 LINESTRING(85.08 61.778,82 60,82 65.836)
 LINESTRING(85.22 61.172,82 60,82 65.836)
 LINESTRING(85.354 60.591,82 60,82 65.836)
 LINESTRING(85.49 60,82 60,82 65.836)
 LINESTRING(85.639 59.358,82 60,82 65.836)
 LINESTRING(85.811 58.613,82 60,82 65.836)
 LINESTRING(86.027 57.675,82 60,82 65.836)
 LINESTRING(86.127 56.537,82 60,82 65.836)
 LINESTRING(85.991 55.244,82 60,82 65.836)
 LINESTRING(85.798 53.421,82 60,82 65.836)
 LINESTRING(85.374 50.729,82 60,82 65.836)
 LINESTRING(87.701 27.667,82 60,82 65.836)
 LINESTRING(82 16.893,82 60,82 65.836)
 LINESTRING(78.611 40.781,82 60,82 65.836)
 LINESTRING(77.031 46.348,82 60,82 65.836)
 LINESTRING(76.119 49.814,82 60,82 65.836)
 LINESTRING(75.482 52.232,82 60,82 65.836)
 LINESTRING(74.872 54.019,82 60,82 65.836)
 LINESTRING(74.216 55.506,82 60,82 65.836)
 LINESTRING(73.585 56.937,82 60,82 65.836)
 LINESTRING(72.938 58.402,82 60,82 65.836)
 LINESTRING(72.234 60,82 60,82 65.836)
 LINESTRING(58.022 64.228,82 60,82 65.836)
 LINESTRING(40.559 75.083,82 60,82 65.836)
 LINESTRING(55.77 75.144,82 60,82 65.836)
 LINESTRING(68.869 71.018,82 60,82 65.836)
 LINESTRING(72.535 71.28,82 60,82 65.836)
 LINESTRING(76.052 70.302,82 60,82 65.836)
 LINESTRING(78.66 69.176,82 60,82 65.836)
 LINESTRING(80.751 67.085,82 60,82 65.836)

Visually, the result looks like this:

Line Inside Polygon

The code is available as part of my package of PostGIS functions on my shop