Creating Perpendicular Lines

[Elham Peiravian]

Hello everyone,

I have linestrings representing road segments. I need to find the midpoints of these segments and draw lines perpendicular to the segments going through the midpoint. I have this query but the output in incorrect. Any suggestions?

DROP TABLE if exists p_line;
create table if not exists p_line(
id integer NOT NULL
    PRIMARY KEY AUTOINCREMENT);
SELECT ADDGEOMETRYCOLUMN('p_line', 'line1', 4326, 'Linestring', 'XY' , 1);
SELECT ADDGEOMETRYCOLUMN('p_line', 'line2', 4326, 'Linestring', 'XY' , 1);

WITH
  midpoint AS (
    SELECT
      ST_Line_Interpolate_Point(geojson, 0.5) AS point
    FROM road_segmentation
  )
INSERT INTO p_line(line1, line2)
SELECT
  MakeLine(point, MakePoint(st_X(point) + 10, st_y(point))) AS line1,
  MakeLine(point, MakePoint(st_X(point), st_y(point) + 10)) AS line2
FROM midpoint;

[Jukka Rahkonen]

Hi,

There may be more clever ways, but I think that you could utilize OffsetCurve function and then connect the middle points of the original segments to the middle points of the segments with offset.

-Jukka Rahkonen-

[Uwe Schmitz]

I've done this several times and what you definitely need ist the direction of the line segment

in which your midpoint lies in. Convert this direction to a unit vector, swap x<->y and mutiply the

resulting vector with your required distance. You will get the two perpendicular points by

adding the resulting vector and its negative (multiply by -1) to the calculated midpoint.

That said... i don't know how to achieve this through SQL/Spatialite.

May be someone else can help here.

Regards

uwe

[Stefano Verugi]

Hi Uwe

it would be interesting to see the code in the language you are familiar with, I didn't fully understand the procedure in the way it was described

thanks for your help

regards

[Uwe Schmitz]

Hi Stefano,

here a quick example in the language I am most familiar with... Tcl https://www.tcl.tk/

It's a simple language and I think you can easily translate the idea to other languages.

I put some comments in to make it more readable.

Because Google groups tends to mess up the program code,

I put a copy here: http://paste.tclers.tk/5783

#=====================================================

# get direction of line segment given
# by start-/endpoint as unit vector [xu:yu]
proc pp2uv {x1 y1 x2 y2} {
   set dx [expr {$x2 - $x1}]
   set dy [expr {$y2 - $y1}]
   set ab [expr {sqrt ($dx * $dx + $dy * $dy)}]

   set xu [expr {$dx / $ab}]
   set yu [expr {$dy / $ab}]

   return [list $xu $yu]
}

# get perpedicular line on linesegment ls at
# position t with lenth l
# ls       line segment {x1 y1 x2 y2}
# t        positon on ls [0..1]
# l        length of result segment
# returns:
#   {px1 py1 px2 py2} start-/endpoint perpendicular linesegment
proc getPpline {ls t l} {
   lassign $ls x1 y1 x2 y2
   # calculate point at t
   set it [expr {1-$t}]
   set xp [expr {$x1 * $t + $x2 * $it}]
   set yp [expr {$y1 * $t + $y2 * $it}]
   # get unit vector of ls
   lassign [pp2uv $x1 $y1 $x2 $y2] xu yu
   
   set l2 [expr {$l / 2.}]

   # generate start point of perpedicular line
   set px1 [expr {$xp - $yu * $l2}]
   set py1 [expr {$yp + $xu * $l2}]

   # generate end point of perpedicular line
   set px2 [expr {$xp + $yu * $l2}]
   set py2 [expr {$yp - $xu * $l2}]

   return [list $px1 $py1 $px2 $py2]
}

set linesegment {1. 1. 9. 9.}
set ppl [getPpline $linesegment 0.5 2.0]

puts $ppl

#=====================================================

Regards

uwe

[a.fu...@lqt.it]

it seems a classic problem quite trivial to solve using any
programming language (C, Python, Java etc) but which is
practically impossible to deal with in pure SQL

bye Sandro

[Jukka Rahkonen]

Hi,

I am not sure if it is practically impossible. Tricky it may be, but I think that the followind approuch might work:

1. Create an offset line with desired distanse with OffsetCurve.

2. Dissolve segments with DissolveSegments.

3. Find centroids of original and dissolved segments with ST_Centroid.

4. Join the centroids of the original segments and corresponding offset segments.

However, I do not know if the lines that connect the centroids are guaranteed to be perpendicular to the original segments.

-Jukka Rahkonen-

[Totò Fiandaca]

Have you considered using the ST_ClosestPoint function?
I think it's for you.

Greetings

https://postgis.net/docs/ST_ClosestPoint.html

--
You received this message because you are subscribed to the Google Groups "SpatiaLite Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to spatialite-use...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/spatialite-users/dbd8a51f-fa48-4d22-8e6e-9a051d1d2fc0n%40googlegroups.com.

--

Ing. Salvatore Fiandaca
mobile.:+39 327.493.8955 
m: pigrecoin...@gmail.com

C.F.: FNDSVT71E29Z103G

P.IVA: 06597870820

membro QGIS Italia - http://qgis.it/
socio GFOSS.it - http://gfoss.it/

member OSGeo - https://www.osgeo.org/member/fiandaca/

blog: https://pigrecoinfinito.com/ 
FB: Co-admin - https://www.facebook.com/qgis.it/  

TW: https://twitter.com/totofiandaca

43°51'0.54"N  10°34'27.62"E - EPSG:4326

“Se la conoscenza deve essere aperta a tutti,

perchè mai limitarne l’accesso?” 

R. Stallman

Questo documento, allegati inclusi, contiene informazioni di proprietà di FIANDACA SALVATORE e deve essere utilizzato esclusivamente dal destinatario in relazione alle finalità per le quali è stato ricevuto. E' vietata qualsiasi forma di riproduzione o divulgazione senza l'esplicito consenso di FIANDACA SALVATORE. Qualora fosse stato ricevuto per errore si prega di informare tempestivamente il mittente e distruggere la copia in proprio possesso.

[Jukka Rahkonen]

I did bad thinking. If lines are segmented before offsetting then the result should be good.

-Jukka-

[Stefano Verugi]

Hi,

please consider the following:

select
makeline(

        st_project(c.geometry,offsetvalue,l."azm"+radians(90)),

st_project(c.geometry,offsetvalue,l."azm"-radians(90))) as geometry
 from center_point c, myline l;

/* 

center_point is from a table containing the midpoints of the line using the st_line_interpolate_point(linestring, fraction)

myline is a table with the line containing a field "azm" calculated with the function St_azimuth(startpoint, endpoint) of  the line itself

st_project(center_point, offsetvalue, azm), well, returns a projected point of center_point based on the other 2 parameters)works with latlong coordinates

*/

I can prepare something more complete next weekend but I thing the core is in the above lines

hope it helps

PS if you are using linestrings defined by more the 2 points the azimuth of the center_point needs to use the two adjacent vertices

[Uwe Schmitz]

Jukka,

Hi,

I am not sure if it is practically impossible. Tricky it may be, but I think that the followind approuch might work:

1. Create an offset line with desired distanse with OffsetCurve.

2. Dissolve segments with DissolveSegments.

3. Find centroids of original and dissolved segments with ST_Centroid.

4. Join the centroids of the original segments and corresponding offset segments.

However, I do not know if the lines that connect the centroids are guaranteed to be perpendicular to the original segments.

-Jukka Rahkonen-

I'm pretty sure your solution leads to the same result as my procedural algorithm above,

but with a little bit more overhead.

If I understand the OP right, he is looking for a general solution 

applicable to arbitrary LineStrings. My solution assumes that the LineSegment containing 

the "center" is already identified.

For a general solution on LineStrings some more problems come to my mind.

What if the center point lies exactly (or nearly) on a vertex? May be it's

best to calculate the bisector of the adjacent LineSegments?

Last not least I think Sandro is right here: Something very easy to implement in

a "classic" programming language, but very hard to do it in plain SQL.

Regards

uwe

[Bwertz M]

I recently stumbled across this thread that has a bit matured by now, but I cannot agree with Sandro's statement for this problem described above. 

The following SQL script calculates the center of arbitrary input linestrings and checks the northern orientation of the line segments that touch this center. Based on the azimuth, a left and right vertical to the center of the line is calculated, although these may vary in length.

Two different perpendicular line geometries are generated. One approach is based on simple trigonometric functions (geom_trig), and the other one on the ST_Project() function (geom_proj). Interesting though, both geometries vary from each other. "geom_proj" contains segments which are not axactly 90° on a plain map projected in a metric CRS, while the once from "geom_trig" do. It seems that ST_Project() computes the output on the Ellipsoid and not on the Great Circle.

with

--PARAMETERS

para ("key", "value") as (

   values(

     'offset_left', 500

   ),(

     'offset_right', 2500

   ),(

     'epsg', 25832

   )

)

--01> Linienmittelpunkt

--- ST_Line_Interpolate_Point( line Curve , fraction Double precision ) : Point

, cte1 as (

   select pk, ST_Line_Interpolate_Point(geom, 0.5) as g_ln_midpt

   from st_line

)

--02> Linienteil: Start - Mitte

--- ST_Line_Substring( line Curve , start_fraction Double precision , end_fraction Double precision ) : Curve

, cte2 as (

   select pk, ST_Line_Substring(geom, 0.0, 0.5) as g_ln_sub1

   from st_line

)

--- Segmente

, cte2b as (

   select pk, st_dissolveSegments(g_ln_sub1) as g_segs_sub1

   from cte2

)

--03> Linienteil: Mitte - Ende

--- ST_Line_Substring( line Curve , start_fraction Double precision , end_fraction Double precision ) : Curve

, cte3 as (

   select pk, ST_Line_Substring(geom, 0.5, 1.0) as g_ln_sub2

   from st_line

)

--- Segmente

, cte3b as (

   select pk, st_dissolveSegments(g_ln_sub2) as g_segs_sub2

   from cte3

)

--04> Nachbar-Segmente Linienmittelpunkt

, cte4 as (

   select a.pk

   -- letztes Segment der ersten Linienhälfte

   , st_GeometryN(a.g_segs_sub1, st_numGeometries(a.g_segs_sub1)) as g_seg1

   -- erstes Segment der zweiten Linienhälfte

   , st_GeometryN(b.g_segs_sub2, 1) as g_seg2

   from cte2b as a

   join cte3b as b

   using(pk)

)

--05> Azimuth: Segmente

--- ST_Azimuth( pt1 Geometry , pt2 Geometry ) : Double precision

, cte5 as (

   select pk, ST_Azimuth(st_startpoint(g_seg1), st_endpoint(g_seg2)) as azi

   from cte4

)

--06> Mittelsenkrechte

--- ST_Project( start_point Geometry, distance Double precision, azimuth Double precision ) : Point(EPSG:4326)

--- st_Project() verlangt nach WGS84 (EPSG:4326)

, cte6 as (

  -- rechtes Segment

  select a.pk, 'right' as off_dir

    , st_x(a.g_ln_midpt) as x0, st_y(a.g_ln_midpt) as y0

    , sin(azi + pi()/2) * c."value" as x1, cos(azi + pi()/2) * c."value" as y1

    , ST_Project(

    st_transform(g_ln_midpt,4326)

    , c."value"

    , b.azi + pi()/2 -- = +radians(90)

  ) as geom

  from cte1 as a

  join cte5 as b

  using (pk)

  join para as c

  where c."key" = 'offset_right'

  union

  -- linkes Segment

  select a.pk, 'left'

    , st_x(a.g_ln_midpt) as x0, st_y(a.g_ln_midpt) as y0

    , sin(azi - pi()/2) * c."value" as x1, cos(azi - pi()/2) * c."value" as y1

    , ST_Project(

    st_transform(g_ln_midpt,4326)

    , c."value"

    , b.azi - pi()/2

  )

  from cte1 as a

  join cte5 as b

  using (pk)

  join para as c

  where c."key" = 'offset_left'

)

select a.pk, off_dir

  , makeLine(makePoint(x0, y0, c."value"), makePoint(x0 + x1, y0 + y1, c."value")) as g_trig

  , st_transform(makeLine(st_transform(a.g_ln_midpt,4326), b.geom), c."value") as g_proj

from cte1 as a

join cte6 as b

using (pk)

join para as c

where c."key" = 'epsg'

;

Greatings

Bernd

[a.fu...@lqt.it]

On Mon, 11 Dec 2023 05:28:57 -0800 (PST), Bwertz M wrote:
> I recently stumbled across this thread that has a bit matured by now,
> but I cannot agree with Sandro's statement for this problem described
> above. 
>

Sandro simply said:
"it seems a classic problem quite trivial to solve using any
programming language (C, Python, Java etc) but which is
practically impossible to deal with in pure SQL"

> The following SQL script calculates the center of arbitrary input
> linestrings and checks the northern orientation of the line segments
> that touch this center. Based on the azimuth, a left and right
> vertical to the center of the line is calculated, although these may
> vary in length.
>
> Two different perpendicular line geometries are generated. One
> approach is based on simple trigonometric functions (geom_trig), and
> the other one on the ST_Project() function (geom_proj). Interesting
> though, both geometries vary from each other. "geom_proj" contains
> segments which are not axactly 90° on a plain map projected in a
> metric CRS, while the once from "geom_trig" do. It seems that
> ST_Project() computes the output on the Ellipsoid and not on the
> Great
> Circle.
>

I offer you my sincere congratulations.
A brilliant tour de force that solves the problem in terms of pure
Spatial SQL. BRAVO

that being said: it's one of the most complicated SQL scripts I've
ever come across.
It's so long and complex that I gave up before I finished reading it;
I can easily imagine the reactions of novice beginners.
Final consideration: it has every appearance of being deadly slow at
run time.

I confirm my previous assessment: it's a problem that can be solved
in an infinitely simpler, more elegant and quicker way using any
programming language.

Tackling it in pure SQL still remains an excellent intellectual
challenge to test one's personal skills on the most extreme and
imaginative use of Spatial SQL

bye Sandro