C#用lat/long数组计算多个曲率

C# Calculate multiple curvatures with lat/long array

我有一组纬度和经度构成了从起点到终点的整条路线,但我真的很难计算整条铁路路线的所有曲率,包括左转为负值和右转积极(或反之亦然,以更容易者为准)。

我也有代码可以给出方位角以及两个方位角之间的区别,但就我所知(下面的代码)。

从下图中的绿色标记处出发,向西(左)行驶,有:

本质上就是在计算上面每个转弯的转弯半径list/picture。

所需的输出对象

public EdgeCurve[] Curves {get; set;}

public class EdgeCurve
{
    public double StartLatitude {get;s et;}
    public double StartLongitude {get; set;}

    public double EndLatitude {get; set;}
    public double EndLongitude {get; set;}

    public double Curve {get; set;}

    //Optional
    public double Distance {get;set;}
}

轴承代码

public static class BearingHelper
{

    /// <summary>
    /// Angle between two points, where North is 0 and South is 180.
    /// </summary>
    /// <param name="lat0">Latitude for point A.</param>
    /// <param name="lon0">Longitude for point A.</param>
    /// <param name="lat1">Latitude for point B.</param>
    /// <param name="lon1">Longitude for point B.</param>
    /// <returns>Angle between points A and B, in degrees.</returns>
    public static double Degrees(double lat0, double lon0, double lat1, double lon1)
    {
        var dlon = (lon1 - lon0).ToRadians();
        var dphi = Math.Log(Math.Tan(lat1.ToRadians() / 2 + PI_d4) / Math.Tan(lat0.ToRadians() / 2 + PI_d4));

        if (Math.Abs(dlon) > Math.PI)
        {
            dlon = dlon > 0 ? -(PI_m2 - dlon) : (PI_m2 + dlon);
        }

        return (Math.Atan2(dlon, dphi).ToDegrees() + 360) % 360;
    }

    public static double AbsDifference(double bearing0, double bearing1)
    {
        double diff = Math.Abs(bearing0 - bearing1);

        if (diff > 180)
            diff -= 360;

        return Math.Abs(diff);
    }
}

public static class CoordinateMath
{
    /// <summary>
    /// PI multiplied by 2.
    /// </summary>
    public const double PI_m2 = 6.28318530718;

    /// <summary>
    /// PI divided by 4.
    /// </summary>
    public const double PI_d4 = 0.78539816339;

    /// <summary>
    /// Convert degrees to radians.
    /// </summary>
    public static double ToRadians(this double deg) => deg * 0.01745329252;

    /// <summary>
    /// Convert radians to degrees.
    /// </summary>
    public static double ToDegrees(this double rad) => rad * 57.295779513;
}

数组 Lat/Longs

[
  {
    "Latitude": 53.743647,
    "Longitude": -0.347531
  },
  {
    "Latitude": 53.743635,
    "Longitude": -0.348089
  },
  {
    "Latitude": 53.74365,
    "Longitude": -0.348648
  },
  {
    "Latitude": 53.74369,
    "Longitude": -0.349114
  },
  {
    "Latitude": 53.743754,
    "Longitude": -0.3496
  },
  {
    "Latitude": 53.743867,
    "Longitude": -0.35029
  },
  {
    "Latitude": 53.743941,
    "Longitude": -0.350679
  },
  {
    "Latitude": 53.744063,
    "Longitude": -0.351096
  },
  {
    "Latitude": 53.744518,
    "Longitude": -0.352272
  },
  {
    "Latitude": 53.744557,
    "Longitude": -0.352369
  },
  {
    "Latitude": 53.744888,
    "Longitude": -0.353199
  },
  {
    "Latitude": 53.744925,
    "Longitude": -0.353296
  },
  {
    "Latitude": 53.745078,
    "Longitude": -0.353724
  },
  {
    "Latitude": 53.745176,
    "Longitude": -0.354035
  },
  {
    "Latitude": 53.745339,
    "Longitude": -0.354672
  },
  {
    "Latitude": 53.745494,
    "Longitude": -0.355461
  },
  {
    "Latitude": 53.745591,
    "Longitude": -0.356196
  },
  {
    "Latitude": 53.745644,
    "Longitude": -0.356898
  },
  {
    "Latitude": 53.74566,
    "Longitude": -0.357678
  },
  {
    "Latitude": 53.745633,
    "Longitude": -0.358616
  },
  {
    "Latitude": 53.745577,
    "Longitude": -0.359361
  },
  {
    "Latitude": 53.745521,
    "Longitude": -0.359924
  },
  {
    "Latitude": 53.745439,
    "Longitude": -0.360563
  },
  {
    "Latitude": 53.745344,
    "Longitude": -0.361198
  },
  {
    "Latitude": 53.745192,
    "Longitude": -0.362056
  },
  {
    "Latitude": 53.745051,
    "Longitude": -0.362748
  },
  {
    "Latitude": 53.744811,
    "Longitude": -0.363821
  },
  {
    "Latitude": 53.744539,
    "Longitude": -0.364885
  },
  {
    "Latitude": 53.74428,
    "Longitude": -0.365809
  },
  {
    "Latitude": 53.744023,
    "Longitude": -0.366656
  },
  {
    "Latitude": 53.743862,
    "Longitude": -0.367141
  },
  {
    "Latitude": 53.743723,
    "Longitude": -0.367545
  },
  {
    "Latitude": 53.7434,
    "Longitude": -0.368412
  },
  {
    "Latitude": 53.743096,
    "Longitude": -0.369159
  },
  {
    "Latitude": 53.742528,
    "Longitude": -0.370427
  },
  {
    "Latitude": 53.742197,
    "Longitude": -0.371099
  },
  {
    "Latitude": 53.741909,
    "Longitude": -0.371663
  },
  {
    "Latitude": 53.741679,
    "Longitude": -0.372077
  },
  {
    "Latitude": 53.741462,
    "Longitude": -0.372498
  },
  {
    "Latitude": 53.74131,
    "Longitude": -0.372765
  },
  {
    "Latitude": 53.737842,
    "Longitude": -0.379223
  },
  {
    "Latitude": 53.737396,
    "Longitude": -0.380011
  },
  {
    "Latitude": 53.737152,
    "Longitude": -0.380388
  },
  {
    "Latitude": 53.737006,
    "Longitude": -0.380587
  },
  {
    "Latitude": 53.736736,
    "Longitude": -0.380912
  },
  {
    "Latitude": 53.736466,
    "Longitude": -0.381183
  },
  {
    "Latitude": 53.736294,
    "Longitude": -0.381331
  },
  {
    "Latitude": 53.736187,
    "Longitude": -0.381417
  },
  {
    "Latitude": 53.735827,
    "Longitude": -0.381659
  },
  {
    "Latitude": 53.735664,
    "Longitude": -0.381751
  },
  {
    "Latitude": 53.735378,
    "Longitude": -0.381882
  },
  {
    "Latitude": 53.735185,
    "Longitude": -0.381952
  },
  {
    "Latitude": 53.734976,
    "Longitude": -0.382014
  },
  {
    "Latitude": 53.734575,
    "Longitude": -0.382077
  },
  {
    "Latitude": 53.734372,
    "Longitude": -0.382083
  },
  {
    "Latitude": 53.734044,
    "Longitude": -0.382055
  },
  {
    "Latitude": 53.733148,
    "Longitude": -0.381884
  },
  {
    "Latitude": 53.732913,
    "Longitude": -0.381838
  },
  {
    "Latitude": 53.732526,
    "Longitude": -0.381762
  },
  {
    "Latitude": 53.732156,
    "Longitude": -0.381722
  },
  {
    "Latitude": 53.731976,
    "Longitude": -0.381718
  },
  {
    "Latitude": 53.731854,
    "Longitude": -0.381728
  }
]

假设您对将路径建模为一系列线段感到满意,您可以通过计算相邻线段之间的带符号角度来简单地计算曲率。

伪代码

for(var i = 1; i < points.Length - 1; i++){
    var a = (points[i] - points[i-1]).Normalized;
    var b = (points[i+1] - points[i]).Normalized;
    // calculate signed angle between vectors
    angle = atan2( a.x*b.y - a.y*b.x, a.x*b.x + a.y*b.y );
}

根据finding signed angle between vectors

计算符号角的代码

如果你想要平滑的插值,比如Bézier spline,你仍然可以创建它,把它分割成更短的线段,然后使用上面的方法。但是,如果您想要平滑段的单个曲率值,则需要定义该值应该是多少。因为平滑函数将具有连续的曲率范围。

请注意,这假设平面坐标,如果您有 lat/long,您可能想要重新投影到平面坐标系,可能每个线段对一个。