Segment distance along path (imported from kml) using mapping toolbox

As an example - suppose I import a kml of a highway route and want to generate the lat/longs of where to place the highway marker signs. I need to find the distance traveled along the road at 1 mile incriments. Roads are not straight or polynomial curves so I can not simply interpolate (? I assume). I think the mapping toolbox must have something for this, I just don't know the magic phrase to help find it. Or a similar problem: I'm mapping a trip from NY to FL, assume I can import a KML of the roads, if I need to stop every 150 miles for gas, where am I stopping along that route? How can I solve the displacement along a path that is denoted by a set of lat/longs? It should be simple so I feel I am missing something obvious.

Answer 1

[Andres]

Maybe this somewhat naive approach gives you a start once your road coordinate resolution is high enough and you are okay with cartesian coordinates? Sorry I don't know the capabilities of the mapping toolbox.

% signpost distance
signpost_dist = 1;

% generate some road points
resolution = 21;
x1 = linspace(0, pi/2, resolution);
x2 = linspace(pi/2, pi/2, ceil(resolution/2));
x3 = linspace(pi/2, pi, resolution);

y1 = sin(x1);
y2 = linspace(1, 0, ceil(resolution/2));
y3 = cos(x3)/2;

x = [x1, x2, x3];
y = [y1, y2, y3];

% calculate curve length
t = 0:numel(y)-1;
s = cart2alc(x, y, t);

x_signpost = interp1(s, x, signpost_dist:signpost_dist:s(end));
y_signpost = interp1(s, y, signpost_dist:signpost_dist:s(end));

figure
plot(x,y,'.-',DisplayName='road')
hold on
plot(x_signpost, y_signpost, 'rs', DisplayName='signpost')
grid on
axis equal
xlabel('x (mi)')
ylabel('y (mi)')
title([num2str(s(end)) ' mile road'])

function [s,K] = cart2alc(x,y,t)

% CART2ALC cartesian to arc length and curvature coordinates
% 
% [s,K] = cart2alc(x,y)
% [s,K] = cart2alc(x,y,t)
%

if nargin < 3
    isParametric = false;
else
    isParametric = true;
end

if isParametric
    dxpdt = gradient(x)./gradient(t);
    dypdt = gradient(y)./gradient(t);

    d2xpdt2 = gradient(dxpdt)./gradient(t);
    d2ypdt2 = gradient(dypdt)./gradient(t);
    
    s = cumtrapz(t,hypot(dxpdt,dypdt));
else
    dxpdt = 1;
    dypdt = gradient(y)./gradient(x);
    
    d2xpdt2 = 0;
    d2ypdt2 = gradient(dypdt)./gradient(x);
    
    s = cumtrapz(x,hypot(dxpdt,dypdt));
end

K = (dxpdt.*d2ypdt2 - dypdt.*d2xpdt2)./hypot(dxpdt,dypdt).^3;

end

[Andres]

... and, as always, if there's a match on the file exchange from John D'Errico, try it first https://de.mathworks.com/matlabcentral/fileexchange/34871-arclength

Answer 2

[KSSV]

This should be useful for you: https://in.mathworks.com/matlabcentral/fileexchange/34874-interparc