Extrapoler une ligne dans PostGIS

J'essaie d'extrapoler à partir d'un segment de ligne pour trouver un point sur la ligne mais un tiers du chemin «retour», c'est-à-dire essayer de trouver un point new, des points donnés Aet Bci - dessous:

Étant donné une ligne, je peux l'interpoler pour obtenir une position à un pourcentage particulier le long de celle-ci:

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

J'ai essayé d'entrer un nombre négatif pour trouver un point le long de la ligne dans la direction opposée, mais cela échoue car l'argument d'interpolation doit être dans la plage [0, 1]

J'ai d'abord pensé à mettre la ligne à l'échelle, mais cela n'utilise pas le centre de la ligne comme origine, donc inutile pour mes besoins.

Une autre façon dont j'ai résolu un problème similaire précédemment est de le décomposer en étapes suivantes.

-- 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: correction de l'affectation de newlength de sorte qu'elle soit de 1 1/3 de la longueur, plutôt que de 1/3 et commutation A & B pour correspondre au diagramme.

Je l'ai résolu avec:

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))
          )

Explication:

(2-d) Mettez le point de départ à l'échelle par un facteur de 1,3333, en prenant le milieu du segment de ligne comme origine pour la mise à l'échelle.

Sortez le papier millimétré!

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

J'ai écrit une fonction pour cela:

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
-- /gis//a/104451/44921
-- /gis//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
);

Notez que cela donne la chaîne de caractères la plus longue mais pas le point de terminaison.

Code sur GitHub Gist (si vous votez ici, j'apprécierais une étoile là aussi)

Description des paramètres (au cas où vous les auriez manqués dans le commentaire de la fonction sql):

• La longueur du premier segment sera original_length * head_rate + head_constant.
• Si vous voulez qu'il soit doublé, le taux de tête est de 2, la constante est de 0.
• En Hongrie, nous utilisons normalement une projection EOV qui est basée sur un compteur. Donc, si je veux ajouter 2 mètres à la fin de la ligne, je mets le tail: rate à 1 et le tail_constant à 2.