Path radius from lat/long data

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.

Answer 1

[Sean de Wolski]

A place to start:

distance

We might be able to be of more assistance if you gave us sample values and expected results etc...

[John F]

*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]

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

Answer 2

[Babak]

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

http://en.wikipedia.org/wiki/Radius_of_curvature_(applications)

for more details.

[John F]

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.

Answer 3

[Babak]

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)

http://en.wikipedia.org/wiki/Spherical_coordinate_system

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)

Answer 4

[Star Strider]

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.