Path radius from lat/long data

37 次查看(过去 30 天)
John F
John F 2012-9-6
I am working on computing the approximate radius of curvature of a road. I have lat/long coordinates for points along the path, and I'd like to estimate the circular radius between two points that are maybe 25m apart. My data is taken at about 1m intervals.
I've been using a lot of homemade functions to do this, and now I have the Mapping Toolbox, which appears to have a lot more options for the calculations.
Any help would be greatly appreciated.

回答(4 个)

Sean de Wolski
Sean de Wolski 2012-9-6
A place to start:
We might be able to be of more assistance if you gave us sample values and expected results etc...
  2 个评论
John F
John F 2012-9-7
*distance * will only give the GC distance and the azimuth. If you know the start and end azimuth and the GC distance (the chord length across the circle), you can infer a turn radius as:
r = GC Distance / (2*sin(Az diff/2)).
I was hoping to have some sort of built-in way to do it.
John F
John F 2012-9-7
Below I've attached a series of lat/long pairs. This is a snipit of the data I'm using.
42.7835436729962 -86.2020406675281
42.7835518998064 -86.2020559300967
42.7835599541611 -86.2020712802969
42.783567857253 -86.202086588022
42.7835756962754 -86.2021016751859
42.7835836211166 -86.2021164228966
42.7835916232908 -86.2021307936418
42.7835997240812 -86.2021448481797
42.7836079311272 -86.2021587029057
42.7836162014008 -86.2021720493005
42.7836245485008 -86.2021851558504
42.7836329434407 -86.2021980596979
42.7836413419187 -86.202210804503
42.7836498275778 -86.2022234149829
42.7836583625339 -86.2022359816687
42.7836672670222 -86.2022493198003
42.7836765644761 -86.2022633572753
42.7836877230373 -86.2022789853265
42.7836986757104 -86.2022945067202
42.7837095207999 -86.2023097926742
42.7837214161949 -86.2023260410005
42.7837331903952 -86.2023419421105
42.7837455493554 -86.2023578947349
42.7837578268076 -86.2023734772032
42.7837699346272 -86.202388628398
42.7837820196009 -86.2024034746581
42.7837939990458 -86.2024181024752
42.7838058530985 -86.2024323509415
42.7838175083986 -86.2024463479894
42.7838289706887 -86.2024602067699
42.7838403166793 -86.2024737875164
42.7838515605653 -86.2024870203106
42.783862711647 -86.20250001056
42.7838738458047 -86.202512539887
42.7838849010244 -86.2025245684777
42.7838958519132 -86.2025362512888
42.7839067036926 -86.2025476768544
42.7839174906479 -86.2025588035517
42.7839281531605 -86.2025697049991
42.7839386387775 -86.2025803631241
42.7839489772986 -86.2025908361929
42.7839592766637 -86.2026013220483
42.7839695849305 -86.2026117110859
42.783980010619 -86.2026218405116
42.783990447273 -86.202631647871
42.7840008105244 -86.2026412330825
42.7840112381582 -86.2026506030633
42.7840222969203 -86.2026602966352
42.7840333852842 -86.2026698785234
42.78404440119 -86.2026793082072
42.7840553692155 -86.2026884508119
42.7840664132898 -86.2026972452513
42.7840773494053 -86.2027057215037
42.784087973405 -86.2027139793035
42.7840983260447 -86.2027222223533
42.7841086304375 -86.2027305310157
42.7841189644356 -86.2027388391499
42.7841293852338 -86.2027470545824
42.7841405758021 -86.2027554669283
42.784151768432 -86.2027637352618
42.7841629881604 -86.2027719170737
42.7841742126824 -86.2027799538251
42.7841860452533 -86.2027882467188
42.7841978855006 -86.2027963889979
42.7842097372356 -86.202804459438
42.7842215579688 -86.2028123149511
42.7842334079254 -86.2028200584785
42.7842452437391 -86.2028277962096
42.7842576318921 -86.2028359517141
42.7842698993575 -86.2028440850134
42.7842821653297 -86.2028520714827
42.7842984490032 -86.2028613056736
42.784314580663 -86.2028700991099
42.7843305427649 -86.2028786079394
42.7843509454516 -86.2028873954075
42.784370944074 -86.2028957057821
42.7843920927253 -86.202903994588
42.7844160647641 -86.2029123634886
42.7844393243509 -86.2029203055591
42.7844617478533 -86.2029279343462
42.7844834399407 -86.2029352613239
42.7845044371786 -86.2029423099627
42.7845247701247 -86.2029492953445
42.7845444971385 -86.2029560068923
42.7845636820415 -86.2029623828627
42.7845824483295 -86.2029684987096

请先登录,再进行评论。


Babak
Babak 2012-9-6
编辑:Babak 2012-9-6
It depends your data are in 2D (x-y plane) or 3D (xyz plane).
There are equations to find the radius of curvature for each case. For example for 2D data if you have a function y=f(x), or sampled data points from this function, the radius of curvature at any point x0, is mathematically proven to be:
rho = abs(1+f'(x0)^2))^(1.5)/abs(f"(x0))
where f'(x0) and f"(x0) are the first and second derivative of the function at the point of interest.
Now, if I were in your shoes, I would try to find a "finite difference approximation" for these derivatives with the data I have and use the above equation to find the radius of curvature.
Offcourse, if your data are in 3D you need to use the proper equation. See
for more details.
  1 个评论
John F
John F 2012-9-6
My data is 2D. So, that's not an issue.
The problem I see with this method is that my points are lat/long in degrees, and I need the radius of curvature in meters. If my grid were all meters, it would be simple.

请先登录,再进行评论。


Babak
Babak 2012-9-6
编辑:Babak 2012-9-6
You must be using the geographical coordinate system. It is same as the spherical coordinates. Here is the transformation to Cartesian coordinates:
x = r*sin(theta)*cos(phi)
y = r*sin(theta)*sin(phi)
z = r*cos(phi)
where theta is the angle of inclination from the plane of reference. So you have a bunch of theta and phi angles. You can use the above transformation and find a whole bunch of (x/r,y/r,z/r) = (a,b,c) data. You have the cord distance between two points of (x1,y1,z1) and (x2,y2,z2) and need to compute the arc distance btween these two points. The cord distance (25m) between the points you have is
d_cord_12 = sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2)
By factorizing r you can write
(x1,y1,z1) = r1*(a1,b1,c1)
and
(x1,y1,z1) = r2*(a2,b2,c2)
If the points are close to each other you can estimate r1=r2=r and reduce the d_cord_12 equation to
d_cord_12 = r * sqrt((a1-a2)^2+(b1-b2)^2+(c1-c2)^2)
by having (a1,b1,c1) and (a2,b2,c2) and d_cord_12 you can compute r as
r = d_cord_12 / sqrt((a1-a2)^2+(b1-b2)^2+(c1-c2)^2)

Star Strider
Star Strider 2012-9-7
编辑:Star Strider 2012-9-8
I'm not sure there is a completely built-in way to get the distances along the road you are interested in, but this seems to produce [x y z] in meters, since the results of [xchk ychk zchk] — delivering a y value of 60 nautical miles per degree longitude at the equator — appear to be accurate and appropriate. I used sph2cart in these calculations. I named your data array AzEl here:
d2r = @(a) pi*a/180;
ka = 6378.1370; % Equatorial radius (km)
kb = 6356.7523; % Polar radius (km)
a = ka*1000;
b = kb*1000;
phi = d2r(AzEl(:,1));
R = sqrt(((a^2*cos(phi)).^2 + (b^2*sin(phi)).^2) ./ ((a*cos(phi)).^2 + (b*sin(phi)).^2)); % Radius at geodedic latitude (phi)
Rm = median(R);
[xchk ychk zchk] = sph2cart(d2r(1), d2r(0), a);
ynm = ychk/1852; % Convert from meters to nautical miles
[x y z] = sph2cart(d2r(AzEl(:,2)), d2r(AzEl(:,1)), Rm);
A further check indicates the median of DistVct — the vector of distances between your road data points — to be about 1.193 meters:
DifMtx = [diff(x) diff(y) diff(z)];
DistVct = hypot(DifMtx(:,2), DifMtx(:,3));
The geodedic data and related calculations are from the Wikipedia article on Earth Radius.
If this does what you want it to do, you can code it as a function, providing an almost-built-in capability.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by