19

I'm trying to extrapolate from a line segment to find a point on the line but a 3rd of the way 'back', i.e. trying to find point new, given points A and B below:

enter image description here

Given a line, I can interpolate it to get a position at any particular percentage along it:

=# select st_line_interpolate_point(
   st_makeline('0101000020E6100000300DC347C49418C03EE8D9ACFAA44A40', 
               '0101000020E6100000FB743C66A03218C0CDCCCCCCCC7C4A40'), 
   0.333);
0101000020E6100000ED45B41D537718C069C6A2E9EC984A40

I tried entering a negative number to find a point along the line in the opposite direction, but that fails as the interpolation argument has to be in the range [0, 1]

I thought about first scaling the line, but that doesn't use the centre of the line as the origin, so is useless for my purposes.

21

Another way I have solved a similar problem previously is to break it down into the following steps.

-- get the points A and B given a line L
A := ST_STARTPOINT(L);
B := ST_ENDPOINT(L);

-- get the bearing from point B --> A
azimuth := ST_AZIMUTH(B,A);

-- get the length of the line A --> B
length := ST_DISTANCE(A,B);
newlength := length + (length * (1/3));   -- increase the line length by 1/3

-- create a new point 1/3 as far away from A as B is from A
newpoint := ST_TRANSLATE(A, sin(azimuth) * newlength, cos(azimuth) * newlength);

EDIT: fixed the assignment of newlength so that it is 1 1/3 the length, rather than 1/3 and switched A & B around to match diagram.

  • We have a winner! Much more comprehensible. – EoghanM Sep 10 '12 at 15:34
  • This is pretty cool – Nathan W Oct 17 '12 at 5:29
  • Thanks. I originally had this snippet from some work I was doing manually interpolating between contour lines - turned out it was a waste of time what I was trying to do with the contours, but glad this this snippet helped someone else out! :) – Jayden Oct 17 '12 at 6:29
6

Have solved it with:

F = 1.3333
st_affine(A, F, 0, 
             0, F, 
            (F-1)*-st_x(st_line_interpolate_point(st_makeline(A, B), 0.5)), 
            (F-1)*-st_y(st_line_interpolate_point(st_makeline(A, B), 0.5))
          )

Explanation:

(2-d) Scale the start point by a factor of 1.3333, taking the midpoint of the line segment as the origin for scaling.

Get out the graph paper!

http://en.wikipedia.org/wiki/Affine_transformation

2

I wrote a function for this:

CREATE OR REPLACE FUNCTION st_extend (
    geom geometry,
    head_rate double precision,
    head_constant double precision,
    tail_rate double precision,
    tail_constant double precision)
  RETURNS geometry AS
$BODY$
-- Extends a linestring.
-- First segment get extended by length * head_rate + head_constant.
-- Last segment get extended by length * tail_rate + tail_constant.
--
-- References:
-- http://blog.cleverelephant.ca/2015/02/breaking-linestring-into-segments.html
-- https://gis.stackexchange.com/a/104451/44921
-- https://gis.stackexchange.com/a/16701/44921
WITH segment_parts AS (
SELECT
(pt).path[1]-1 as segment_num
,
CASE
WHEN
  (nth_value((pt).path, 2) OVER ()) = (pt).path
AND
  (last_value((pt).path) OVER ()) = (pt).path
THEN
  3
WHEN
  (nth_value((pt).path, 2) OVER ()) = (pt).path
THEN
  1
WHEN
  (last_value((pt).path) OVER ()) = (pt).path
THEN
  2
ELSE
  0
END AS segment_flag
,
(pt).geom AS a
,
lag((pt).geom, 1, NULL) OVER () AS b
FROM ST_DumpPoints($1) pt
)
,
extended_segment_parts
AS
(
SELECT
  *
  ,
  ST_Azimuth(a,b) AS az1
  ,
  ST_Azimuth(b,a) AS az2
  ,
  ST_Distance(a,b) AS len
FROM
segment_parts
where b IS NOT NULL
)
,
expanded_segment_parts
AS
(
SELECT
  segment_num
  ,
  CASE
  WHEN
    bool(segment_flag & 2)
  THEN
    ST_Translate(b, sin(az2) * (len*tail_rate+tail_constant), cos(az2) * (len*tail_rate+tail_constant))
  ELSE
    a
  END
  AS a
  ,
  CASE
  WHEN
    bool(segment_flag & 1)
  THEN
    ST_Translate(a, sin(az1) * (len*head_rate+head_constant), cos(az1) * (len*head_rate+head_constant))
  ELSE
    b
  END
  AS b
FROM extended_segment_parts
)
,
expanded_segment_lines
AS
(
SELECT
  segment_num
  ,
  ST_MakeLine(a, b) as geom
FROM
expanded_segment_parts
)
SELECT
  ST_LineMerge(ST_Collect(geom ORDER BY segment_num)) AS geom
FROM expanded_segment_lines
;
$BODY$
LANGUAGE sql;

Usage:

SELECT st_extend(
st_makeline(
  '0101000020E6100000300DC347C49418C03EE8D9ACFAA44A40', 
  '0101000020E6100000FB743C66A03218C0CDCCCCCCCC7C4A40'
),
1.333::double precision,
0::double precision,
1::double precision,
0::double precision
);

Note that this yields the longer linestring but not the endpoint.

Code on GitHub Gist (if you upvote here I would appreciate a star there too)

Description of parameters (in case you missed them in the comment of the sql function):

  • First segment length will be original_length * head_rate + head_constant.
  • If you want it to get doubled then head rate is 2, constant is 0.
  • We in hungary normally use EOV projection which is meter based. So if I want to add 2 meters to the end of the line, i set the tail:rate to 1 and the tail_constant to 2.
  • This works very well. Can you add some info about head_rate, head_constant, tail_rate, and tail_constant? They aren't explained here or on your GitHub. I am assuming that head rate = scale factor for line extension after the endpoint and tail rate = scale factor for line extension before the startpoint. How do the constants work? What effect do they have? – jbalk Nov 16 at 0:09
  • Its in the comment. First segment length will be original_length * head_rate + head_constant. If you want it to get doubled then head rate is 2, constant is 0. We in hungary normally use EOV projection which is meter based. So if I want to add 2 meters to the end of the line, i set the tail:rate to 1 and the tail_constant to 2. – SzieberthAdam Nov 17 at 9:31
  • Thank you! And thank you for sharing this function. It works perfectly and it runs quickly. – jbalk Nov 18 at 18:15

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.