根据开始和结束地理坐标以指定间隔划分地理线

Posted 2023-02-23

【中文标题】根据开始和结束地理坐标以指定间隔划分地理线

【英文标题】：Divide Geo Line in specified interval based on start and end geo-coordinates

【发布时间】：2019-03-12 05:20:17

【问题描述】：

我有两个地理坐标（起点和终点），我可以在它们之间创建一条线。假设这两个坐标相距 30 米，那么我需要以 3 米的间隔找到每个点的地理位置。所以有10个这样的点是必需的。

我可以通过一些公式找到点，但这些点与起点和终点形成的线的方向不同。

到目前为止我所做的如下......

using System;

namespace Test

public class AzimuthCalculator

    public const double range = 0.00186411F; // in Miles
    public const double rangeInMeter = 3F;
    public const double factor = 0.003F;
    public static void Main(String[] args)
    
        double sLatitude = 28.6187763214111F;
        double sLongitude = 77.2093048095703F;

        double eLatitude = 28.6191763153134F;
        double eLongitude = 77.2097146511078F;

        Console.WriteLine($"Start Point : sLatitude, sLongitude");
        Console.WriteLine($"End Point : eLatitude,eLongitude"); 

        double distance = CalculateDistance(sLatitude, sLongitude, eLatitude, eLongitude);
        Console.WriteLine($"Distance : distance miles, MilesToMeter(distance) meter, (distance * 1.60934) kilometer");
        distance = distance * 1.60934;

        double numberOfIDS = distance/factor;

        double azimuthAngle = DegreeBearing(sLatitude, sLongitude, eLatitude, eLongitude);
        Console.WriteLine($"Azimuth : azimuthAngle");
        Console.WriteLine($"IDS : numberOfIDS");

        double constantAzimuth = (azimuthAngle/numberOfIDS);

        azimuthAngle = constantAzimuth;
        Console.WriteLine($"Original Azimuth : azimuthAngle");

        double[] getAnotherPoint = pointRadialDistance(sLatitude, sLongitude, azimuthAngle, distance);
        Console.WriteLine($"End Point : getAnotherPoint[0],getAnotherPoint[1]");

        distance = 0.003;   // 3 meter

        getAnotherPoint = pointRadialDistance(sLatitude, sLongitude, azimuthAngle, distance);


        int totalIDS = (Int32)(numberOfIDS);

        for(int i=0;i<totalIDS;i++)
        
            sLatitude = getAnotherPoint[0];
            sLongitude = getAnotherPoint[1];
            distance = 0.003;
            Console.WriteLine($"new PointLatLng(getAnotherPoint[0],getAnotherPoint[1]),");
            getAnotherPoint = pointRadialDistance(sLatitude, sLongitude, azimuthAngle, distance);
        
        Console.ReadLine();
           

    static double rEarth = 6371.01F; // Earth radius in km
    static double epsilon = 0.000001F;

    public static double[] pointRadialDistance(double lat1, double lon1, double bearing, double distance)
    
        double rlat1 = ToRad(lat1);
        double rlon1 = ToRad(lon1);
        double rbearing = ToRad(bearing);
        double rdistance = distance / rEarth; // normalize linear distance to radian angle

        double rlat = Math.Asin( Math.Sin(rlat1) * Math.Cos(rdistance) + Math.Cos(rlat1) * Math.Sin(rdistance) * Math.Cos(rbearing));
        double rlon = 0.0F;

        if ( Math.Cos(rlat) == 0 || Math.Abs(Math.Cos(rlat)) < epsilon) // Endpoint a pole
            rlon=rlon1;
        else
            rlon = ((rlon1 - Math.Asin( Math.Sin(rbearing) * Math.Sin(rdistance) / Math.Cos(rlat) ) + Math.PI ) % (2*Math.PI) ) - Math.PI;

        double lat = ToDegrees(rlat);
        double lon = ToDegrees(rlon);

        double[] listNew = new double[2];
        listNew[0] = lat;
        listNew[1] = lon;
        return (listNew);
    

    public static GeoLocation FindPointAtDistanceFrom(GeoLocation startPoint, double initialBearingRadians, double distanceKilometres)
    
        const double radiusEarthKilometres = 6371.01;
        var distRatio = distanceKilometres / radiusEarthKilometres;
        var distRatiosine = Math.Sin(distRatio);
        var distRatioCosine = Math.Cos(distRatio);

        var startLatRad = DegreesToRadians(startPoint.Latitude);
        var startLonRad = DegreesToRadians(startPoint.Longitude);

        var startLatCos = Math.Cos(startLatRad);
        var startLatSin = Math.Sin(startLatRad);

        var endLatRads = Math.Asin((startLatSin * distRatioCosine) + (startLatCos * distRatioSine * Math.Cos(initialBearingRadians)));

        var endLonRads = startLonRad
            + Math.Atan2(
                Math.Sin(initialBearingRadians) * distRatioSine * startLatCos,
                distRatioCosine - startLatSin * Math.Sin(endLatRads));

        return new GeoLocation
        
            Latitude = RadiansToDegrees(endLatRads),
            Longitude = RadiansToDegrees(endLonRads)
        ;
    

    public struct GeoLocation
    
        public double Latitude  get; set; 
        public double Longitude  get; set; 
    

    public static double DegreesToRadians(double degrees)
    
        const double degToRadFactor = Math.PI / 180;
        return degrees * degToRadFactor;
    

    public static double RadiansToDegrees(double radians)
    
        const double radToDegFactor = 180 / Math.PI;
        return radians * radToDegFactor;
    

    public static double CalculateDistance(double sLatitude, double sLongitude, double eLatitude, double eLongitude)
    
        var radiansOverDegrees = (Math.PI / 180.0);

        var sLatitudeRadians = sLatitude * radiansOverDegrees;
        var sLongitudeRadians = sLongitude * radiansOverDegrees;
        var eLatitudeRadians = eLatitude * radiansOverDegrees;
        var eLongitudeRadians = eLongitude * radiansOverDegrees;

        var dLongitude = eLongitudeRadians - sLongitudeRadians;
        var dLatitude = eLatitudeRadians - sLatitudeRadians;

        var result1 = Math.Pow(Math.Sin(dLatitude / 2.0), 2.0) + Math.Cos(sLatitudeRadians) * Math.Cos(eLatitudeRadians) * Math.Pow(Math.Sin(dLongitude / 2.0), 2.0);

        // Using 3956 as the number of miles around the earth
        var result2 = 3956.0 * 2.0 * Math.Atan2(Math.Sqrt(result1), Math.Sqrt(1.0 - result1));

        return result2;
       

    static double DegreeBearing(double lat1, double lon1,double lat2, double lon2)
    
        var dLon = ToRad(lon2 - lon1);
        var dPhi = Math.Log(Math.Tan(ToRad(lat2) / 2 + Math.PI / 4) / Math.Tan(ToRad(lat1) / 2 + Math.PI / 4));
        if (Math.Abs(dLon) > Math.PI)
        dLon = dLon > 0 ? - (2 * Math.PI - dLon) : (2 * Math.PI + dLon);
        return ToBearing(Math.Atan2(dLon, dPhi));            
    

    public static double ToRad(double degrees)
    
        return degrees * (Math.PI / 180);
    

    public static double ToDegrees(double radians)
    
        return radians * 180 / Math.PI;
    

    public static double ToBearing(double radians)
    
        // convert radians to degrees (as bearing: 0...360)
        return (ToDegrees(radians) + 360) % 360;
    

    public static double MeterToMiles(double meter)
    
        return (meter / 1609.344);
    

    public static double MilesToMeter(double miles)
    
        return (miles * 1609.344);
    

【参考方案1】：

为什么你从这里constantAzimuth = (azimuthAngle/numberOfIDS);计算错误的方位，然后再使用它？

您可以使用方法described here 计算大圆路径上的中间点（本质上是 SLERP - 球面线性插值）

Formula:    
a = sin((1−f)⋅δ) / sin δ
b = sin(f⋅δ) / sin δ
x = a ⋅ cos φ1 ⋅ cos λ1 + b ⋅ cos φ2 ⋅ cos λ2
y = a ⋅ cos φ1 ⋅ sin λ1 + b ⋅ cos φ2 ⋅ sin λ2
z = a ⋅ sin φ1 + b ⋅ sin φ2
φi = atan2(z, √x² + y²)
λi = atan2(y, x)
where   
    f is fraction along great circle route (f=0 is point 1, f=1 is point 2), 
    δ is the angular distance d/R between the two points.

【讨论】：

但是当我没有写constantAzimuth = (azimuthAngle/numberOfIDS)时，我还是有同样的问题

所以你可以计算 Destination point given distance and bearing from start point 的距离 3,6,9... 米并给出 sLatitude/sLongitude 而无需修改这些坐标（我怀疑你为什么一次又一次地改变一切）