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.
采纳的回答
Andres
2023-9-27
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
1 个评论
Andres
2023-9-27
... 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
更多回答（1 个）
KSSV
2023-9-27
1 个投票
类别
在 帮助中心 和 File Exchange 中查找有关 Geographic Plots 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
