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

Posted

技术标签:

【中文标题】根据开始和结束地理坐标以指定间隔划分地理线【英文标题】: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 而无需修改这些坐标(我怀疑你为什么一次又一次地改变一切)

以上是关于根据开始和结束地理坐标以指定间隔划分地理线的主要内容,如果未能解决你的问题,请参考以下文章

地理空间坐标和以千米为单位的距离

如何将以网格间隔收集的点数据转换为 r 中的地理参考数据集?

地理空间坐标和距离(以公里为单位)

地理坐标系和投影坐标系

原始线要素类为地理坐标系,如何获取以米为单位的距离?

Arcgis 地理坐标系和投影坐标系