如何找到给定纬度/长度以北x km的纬度/经度?

问题描述 投票:26回答:6

我有一些生成谷歌地图的C#代码。这些代码查看我需要在地图上绘制的所有点,然后计算出矩形的边界以包含这些点。然后,它会将此边界传递给Google Maps API,以适当地设置缩放级别,以显示地图上的所有点。

这段代码工作正常,但我有一个新的要求。

其中一个点可能具有与之相关的精度。如果是这种情况,那么我在半径设置为精度值的点周围绘制一个圆。再次,这工作正常,但我的边界检查现在没有做我想要它做的事情。我希望边界框包含完整的圆圈。

这需要算法采用点x并计算将在x以北z米以及x以南z米的点y。

有没有人有这个算法,最好是在C#中。我确实找到了一个通用的算法here,但我似乎没有正确实现这个,因为我得到的答案是千米的漂移。

这是通用示例

Lat/lon given radial and distance

A point {lat,lon} is a distance d out on the tc radial from point 1 if:

     lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
     IF (cos(lat)=0)
        lon=lon1      // endpoint a pole
     ELSE
        lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
     ENDIF

这是我的C#翻译。

  // Extend a Point North/South by the specified distance
    public static Point ExtendPoint(Point _pt, int _distance, int _bearing )
    {
        Decimal lat = 0.0;
        Decimal lng = 0.0;

        lat = Math.Asin(Math.Sin(_pt.Lat) * Math.Cos(_distance) + Math.Cos(_pt.Lat) * 
            Math.Sin(_distance) * Math.Cos(_bearing));

         if (Math.Cos(lat) == 0)
         {
            lng = _pt.Lng;      // endpoint a pole
         }
         else 
         {
             lng = (
                 (_pt.Lng - Math.Asin(Math.Sin(_bearing) * Math.Sin(_distance) / Math.Cos(lat)) 
                 + Math.PI) % (2 * Math.PI)) - Math.PI;
         }

         ret = new Point(lat,lng);
         return ret;
    }

我调用此函数的轴承为0来计算新的北偏位置,并将值180设置为计算新的南风位置。

任何人都可以看到我做错了或者可能提供一个已知的工作算法吗?

c# geometry gis
6个回答
6
投票

如果您有一个给定的纬度和经度,您可以计算纬度x-km变化的正确纬度和经度,如下所示:

new-lat = ((old-km-north + x-km-change)/40,075) * 360)
           ^ is the ratio of the                  ^ times the ratio of the circle
           of the earth the change                by 360 to get the total ratio 
           covers.                                covered in degrees.

这同样适用于经度。如果您有总距离加上更改,您可以以类似的方式计算总度数。

new-long = ((old-km-east + x-km-change)/40,075) * 360)
           ^ is the ratio of the                  ^ times the ratio of the circle
           of the earth the change                by 360 to get the total ratio 
           covers.                                covered in degrees.

同样,这些计算应该有效,但我在这里运行纯粹的直觉,但逻辑似乎确实如此。

编辑:正如Skizz 40,075指出的那样,需要使用2.pi.r.cos(lat)或40074.cos(lat)在任意给定的纬度调整到地球的周长。


21
投票

我有一段非常相似的代码。与其他实现相比,它让我得到了非常接近的结果。

我认为与你的问题在于你使用“距离”作为以米为单位的线性距离而不是以弧度表示的角距离。

/// <summary>
/// Calculates the end-point from a given source at a given range (meters) and bearing (degrees).
/// This methods uses simple geometry equations to calculate the end-point.
/// </summary>
/// <param name="source">Point of origin</param>
/// <param name="range">Range in meters</param>
/// <param name="bearing">Bearing in degrees</param>
/// <returns>End-point from the source given the desired range and bearing.</returns>
public static LatLonAlt CalculateDerivedPosition(LatLonAlt source, double range, double bearing)
{
    double latA = source.Latitude * UnitConstants.DegreesToRadians;
    double lonA = source.Longitude * UnitConstants.DegreesToRadians;
    double angularDistance = range / GeospatialConstants.EarthRadius;
    double trueCourse = bearing * UnitConstants.DegreesToRadians;

    double lat = Math.Asin(
        Math.Sin(latA) * Math.Cos(angularDistance) + 
        Math.Cos(latA) * Math.Sin(angularDistance) * Math.Cos(trueCourse));

    double dlon = Math.Atan2(
        Math.Sin(trueCourse) * Math.Sin(angularDistance) * Math.Cos(latA), 
        Math.Cos(angularDistance) - Math.Sin(latA) * Math.Sin(lat));

    double lon = ((lonA + dlon + Math.PI) % UnitConstants.TwoPi) - Math.PI;

    return new LatLonAlt(
        lat * UnitConstants.RadiansToDegrees, 
        lon * UnitConstants.RadiansToDegrees, 
        source.Altitude);
}

哪里

public const double EarthRadius = 6378137.0;   //  WGS-84 ellipsoid parameters

和LatLonAlt以度/米为单位(转换在内部进行)。根据需要调整。

我假设你可以弄清楚UnitConstants.DegreesToRadians的价值是什么:)


9
投票

对于懒惰的人(像我一样))复制粘贴解决方案,Erich Mirabal的版本有很小的变化:

using System.Device.Location; // add reference to System.Device.dll
public static class GeoUtils
{
    /// <summary>
    /// Calculates the end-point from a given source at a given range (meters) and bearing (degrees).
    /// This methods uses simple geometry equations to calculate the end-point.
    /// </summary>
    /// <param name="source">Point of origin</param>
    /// <param name="range">Range in meters</param>
    /// <param name="bearing">Bearing in degrees</param>
    /// <returns>End-point from the source given the desired range and bearing.</returns>
    public static GeoCoordinate CalculateDerivedPosition(this GeoCoordinate source, double range, double bearing)
    {
        var latA = source.Latitude * DegreesToRadians;
        var lonA = source.Longitude * DegreesToRadians;
        var angularDistance = range / EarthRadius;
        var trueCourse = bearing * DegreesToRadians;

        var lat = Math.Asin(
            Math.Sin(latA) * Math.Cos(angularDistance) +
            Math.Cos(latA) * Math.Sin(angularDistance) * Math.Cos(trueCourse));

        var dlon = Math.Atan2(
            Math.Sin(trueCourse) * Math.Sin(angularDistance) * Math.Cos(latA),
            Math.Cos(angularDistance) - Math.Sin(latA) * Math.Sin(lat));

        var lon = ((lonA + dlon + Math.PI) % (Math.PI*2)) - Math.PI;

        return new GeoCoordinate(
            lat * RadiansToDegrees,
            lon * RadiansToDegrees,
            source.Altitude);
    }

    private const double DegreesToRadians = Math.PI/180.0;
    private const double RadiansToDegrees = 180.0/ Math.PI;
    private const double EarthRadius = 6378137.0;
}

用法:

[TestClass]
public class CalculateDerivedPositionUnitTest
{
    [TestMethod]
    public void OneDegreeSquareAtEquator()
    {
        var center = new GeoCoordinate(0, 0);
        var radius = 111320;
        var southBound = center.CalculateDerivedPosition(radius, -180);
        var westBound = center.CalculateDerivedPosition(radius, -90);
        var eastBound = center.CalculateDerivedPosition(radius, 90);
        var northBound = center.CalculateDerivedPosition(radius, 0);

        Console.Write($"leftBottom: {southBound.Latitude} , {westBound.Longitude} rightTop: {northBound.Latitude} , {eastBound.Longitude}");
    }
}

7
投票

我不确定我是否在这里遗漏了一些东西,但我认为这个问题可以改为:“我有一个纬度/经度点,我想要找到该点以北x米和该点以南x米的点。 “

如果这是问题,那么你不需要找到新的经度(这会使事情变得更简单),你只需要一个新的纬度。地球上任何地方的纬度大约为60海里,海里的长度为1,852米。因此,对于新的纬度x南北:

north_lat = lat + x / (1852 * 60)
north_lat = min(north_lat, 90)

south_lat = lat - x / (1852 * 60)
south_lat = max(south_lat, -90)

这并不完全准确,因为地球不是一个完美的球体,每个纬度之间恰好有60海里。然而,其他答案假设纬度线是等距的,所以我假设你不关心它。如果您对可能引入的错误感兴趣,维基百科上有一个很好的表格,显示了此链接中不同纬度的“纬度每1°变化的表面距离”:

http://en.wikipedia.org/wiki/Latitude#Degree_length


4
投票

Ed William's rather awesome site上的两个方程存在问题......但我没有分析它们以了解原因。

I found here的第三个等式似乎给出了正确的结果。

这是php中的测试用例......第三个等式是正确的,前两个给经度不正确的值。

<?php
            $lon1 = -108.553412; $lat1 = 35.467155; $linDistance = .5; $bearing = 170;
            $lon1 = deg2rad($lon1); $lat1 = deg2rad($lat1);
            $distance = $linDistance/6371;  // convert dist to angular distance in radians
            $bearing = deg2rad($bearing);

            echo "lon1: " . rad2deg($lon1) . " lat1: " . rad2deg($lat1) . "<br>\n";

// doesn't work
            $lat2 = asin(sin($lat1) * cos($distance) + cos($lat1) * sin($distance) * cos($bearing) );
            $dlon = atan2(sin($bearing) * sin($distance) * cos($lat1), cos($distance) - sin($lat1) * sin($lat2));
            $lon2 = (($lon1 - $dlon + M_PI) % (2 * M_PI)) - M_PI;  // normalise to -180...+180

            echo "lon2: " . rad2deg($lon2) . " lat2: " . rad2deg($lat2) . "<br>\n";

// same results as above
            $lat3 = asin( (sin($lat1) * cos($distance)) + (cos($lat1) * sin($distance) * cos($bearing)));
            $lon3 = (($lon1 - (asin(sin($bearing) * sin($distance) / cos($lat3))) + M_PI) % (2 * M_PI)) - M_PI;

            echo "lon3: " . rad2deg($lon3) . " lat3: " . rad2deg($lat3) . "<br>\n";

// gives correct answer... go figure
            $lat4 = asin(sin($lat1) * cos($linDistance/6371) + cos($lat1) * sin($linDistance/6371) * cos($bearing) );
            $lon4 = $lon1 + atan2( (sin($bearing) * sin($linDistance/6371) * cos($lat1) ), (cos($linDistance/6371) - sin($lat1) * sin($lat2)));

            echo "lon4: " . rad2deg($lon4) . " lat4: " . rad2deg($lat4) . "<br>\n";
?>

注意我收到了作者(Ed Williams)发送的前两个方程式的电子邮件:

从我的“实施说明”:

关于mod功能的注意事项。这似乎在不同的语言中以不同的方式实现,对结果的符号是否遵循除数或被除数的符号有不同的约定。 (我们希望符号遵循除数或者是欧几里德.C的fmod和Java的%不起作用。)在本文档中,Mod(y,x)是将y除以x并且始终位于0 <=的范围内的余数。 mod <x。例如:mod(2.3,2。)= 0.3和mod(-2.3,2。)= 1.7

如果你有一个floor函数(Excel中的int),则返回floor(x)=“小于或等于x的最大整数”,例如floor(-2.3)= - 3和floor(2.3)= 2

mod(y,x) = y - x*floor(y/x)

以下应该在没有floor函数的情况下工作 - 无论“int”是否向下截断或向下舍入:

mod=y - x * int(y/x)
if ( mod < 0) mod = mod + x

php就像C中的fmod一样,为我的目的而言它“错误”。


0
投票

如果您首先将其重新投影到UTM然后检查距离,则更准确。

希望这可以帮助

© www.soinside.com 2019 - 2024. All rights reserved.