Coordinate offset and heading

[giancan]

Dear friends,
I need our help for a new task.
I have a gps track (sequence of lat long and altitude) and I would like to know if there is a way to make an offset (or a parallel line if you prefer) of this line to its right and one on the left (which means also up and down when the line is not north-south and so on) few meters apart.
Is there any library which allows me to do it easily?
Thanks in advance,
G.

PS: coordinate system is not important now, meaning that I can change it to the one that better suits to the specific task.

[rrashkin]

Let's forget the altitude component for now since I don't think you're really mapping in 3-d here. So you have a series of lat/long. Let's say you know how to convert that to a more convenient x/y system. Now, if your track is a straight line, then you can easily turn it into a parametric representation: y=mx+b where m is the slope (δy/δx) and b is the y-intercept. Then a parallel line has the same m and an offset (+ or -) of b. Is your track straight? If not, is it piece-wise straight?

[giancan]

Dear rrashkin,
unfortunately the track is not straight. Quite on the contrary, most of the times is even orbiting... but I don't need the new line to be as "dense" as the input. I mean, if the input has point every meter, the offset line may have every t2 or 3, I don't care. I just need to know with one look which is the left and which is the right of my line (or the driver and the passenger if you prefer).
Yes, as you say, let's ignore the altitude or let's assume it is 0 or whatever is more convenient.
Thanks,
G.

[rrashkin]

Well, I'm sure with polynomial recursion you can make it as complicated as you need to but if we take every 2 points on the reference track to define a line, then you can define a parallel line that is offset from that as I described above. Then you just do that for every point pair. It's tedious but that's why we have computers.

[Gribouillis]

I have formulas, if (x0, y0) and (x1, y1) are two points on the track. The line joining them has the equation

(y1-y0) * x - (x1-x0) * y = x0 * y1 - x1 * y0

The parallel line on the right at a distance r has the equation

(y1-y0) * x - (x1-x0) * y = x0 * y1 - x1 * y0 - r * sqrt((x1-x0)**2 + (y1-y0)**2)

The parallel line on the left has the same equation with +r instead of -r.

[giancan]

Wow, it looks good... but unfortunately me maths skills are quite low...it will take me a couple of hours to convert it for my purposes. How will I get to the new line coordinate?

[Gribouillis]

Try to built points shown in this picture

[Gribouillis]

Here is a function to compute the point M1 prime in the picture above, given the coordinates of M0, M1, M2. Change r into -r to compute the point on the left side

#!/usr/bin/env python
# -*-coding: utf8-*-
from __future__ import unicode_literals, print_function, division
from math import sqrt

def right_M(r, M0, M1, M2):
    x, y = zip(*(M0, M1, M2))
    dx, dy = ((t[1]-t[0], t[2]-t[1]) for t in (x, y))
    norm = tuple(sqrt(dx[i] ** 2 + dy[i] ** 2) for i in (0, 1))
    u, v = ((-dy[i]/norm[i], dx[i]/norm[i]) for i in (0,1))
    p = 1 + u[0]*v[0]+u[1]*v[1]
    w = ((u[0]+v[0])/p, (u[1]+v[1])/p)
    return (x[1] - r * w[0], y[1] - r * w[1])

if __name__ == "__main__":
    M = [(0, 0), (4, 6), (8, 9), (14, 10.5)]
    print(right_M(1.5, *M[0:3]))
    print(right_M(1.5, *M[1:4]))

""" my output -->
    (5.10555127546399, 4.954163456597992)
    (8.657670780786754, 7.618253085590066)    
"""

[giancan]

OK, as I said, my maths skills are bad.
May I ask you to make an example?
Say that we have the following:
x0 , y0 = 10.15 , 24.32
x1 , y1 = 12.31 , 23.73

What px0 (parallel x0), py0 and px1,py1 will be with an offset of 1 meter on the right?
Thanks and sorry ;)

[Gribouillis]

This should be correct

>>> x0 , y0 = 10.15 , 24.32
>>> x1 , y1 = 12.31 , 23.73
>>> dx = x1 - x0
>>> dy = y1 - y0
>>> from math import sqrt
>>> qq = sqrt(dx**2 + dy**2)
>>> ux, uy = -dy/qq, dx/qq # the vector (ux, uy) is orthogonal to the line and normalized
>>> r = 1 # the 1 meter
>>> px0, py0 = x0 - r * ux, y0 - r * uy
>>> px1, py1 = x1 - r * ux, y1 - r * uy
>>> print(px0, py0)
(9.886504720249329, 23.3553393148111)
>>> print(px1, py1)
(12.04650472024933, 22.7653393148111)

[giancan]

I see. Thanks a lot!
But now should I take x1,y1 and x2,y2 or x2,y2 and x3,y3?
Because if I take x1,y1 and x2,y2 I get a new px1,py1... should I make an average between the two couple of values or should I take 1st and 2nd; 3rd and 4th; 5th and 6th; ...; ...? I don't need a great accuracy.
I copy here a list of points from "real life" :)
Thanks
//////////////////////////////////////////////////////
2761066.458,4209083.553
2760885.230,4209203.495
2760702.221,4209317.843
2760506.633,4209436.833
2760316.499,4209562.510
2760122.803,4209691.190
2759937.790,4209831.062
2756662.771,4210809.258
2756436.792,4210791.790
2756209.923,4210754.669
2755984.390,4210688.616
2755775.888,4210602.640
2755571.728,4210508.202
2755359.331,4210408.171
2755148.937,4210293.538
2754594.566,4209838.703
2754493.599,4209652.845
2754385.285,4209454.435
2754286.211,4209273.905
2754294.560,4209039.479
//////////////////////////////////////////////////////

[Gribouillis]

Well, here is a right path computed with the previous function. The first and last points are missing, you can add them using the px0, py0 method

#!/usr/bin/env python
# -*-coding: utf8-*-
from __future__ import unicode_literals, print_function, division
from math import sqrt

data = """
2761066.458,4209083.553
2760885.230,4209203.495
2760702.221,4209317.843
2760506.633,4209436.833
2760316.499,4209562.510
2760122.803,4209691.190
2759937.790,4209831.062
2756662.771,4210809.258
2756436.792,4210791.790
2756209.923,4210754.669
2755984.390,4210688.616
2755775.888,4210602.640
2755571.728,4210508.202
2755359.331,4210408.171
2755148.937,4210293.538
2754594.566,4209838.703
2754493.599,4209652.845
2754385.285,4209454.435
2754286.211,4209273.905
2754294.560,4209039.479
""".strip().split("\n")
M = list(tuple(float(v) for v in line.split(',')) for line in data)

def right_M(r, M0, M1, M2):
    x, y = zip(*(M0, M1, M2))
    dx, dy = ((t[1]-t[0], t[2]-t[1]) for t in (x, y))
    norm = tuple(sqrt(dx[i] ** 2 + dy[i] ** 2) for i in (0, 1))
    u, v = ((-dy[i]/norm[i], dx[i]/norm[i]) for i in (0,1))
    p = 1 + u[0]*v[0]+u[1]*v[1]
    w = ((u[0]+v[0])/p, (u[1]+v[1])/p)
    return (x[1] - r * w[0], y[1] - r * w[1])

if __name__ == "__main__":
    r = 1 # distance (same unit as the points)
    for i in xrange(len(M) - 2):
        print(right_M(r, *M[i:i+3]))

""" my output -->
   (2760885.7709899466, 4209204.336130894)
   (2760702.7458359217, 4209318.69422428)
   (2760507.1687698085, 4209437.677572424)
   (2760317.0513888276, 4209563.343587389)
   (2760123.3817485413, 4209692.006075578)
   (2759938.2491084263, 4209831.96852491)
   (2756662.8791700928, 4210810.269344599)
   (2756436.672511555, 4210792.7837467715)
   (2756209.7008925797, 4210755.645955941)
   (2755984.0579236434, 4210689.560748778)
   (2755775.487301093, 4210603.556452203)
   (2755571.305043517, 4210509.108156547)
   (2755358.8783515077, 4210409.063172728)
   (2755148.3756782985, 4210294.370962217)
   (2754593.7805975624, 4209839.352114648)
   (2754492.72078066, 4209653.32325918)
   (2754384.4078044216, 4209454.915134511)
   (2754285.201830749, 4209274.144655606)   
"""

All this should be made faster by using numpy.

(image added for r=40)

[giancan]

Great! Thanks really a lot