C# 如何找到给定纬度/经度以北 x 公里的纬度/经度？

How do I find the lat/long that is x km north of a given lat/long?

提问by Steve Weet

I have some C# code that generates google maps. This codes looks at all the Points I need to plot on the map and then works out the Bounds of a rectangle to include those points. It then passes this bounds to the Google Maps API to set the zoom level appropriately to show all of the points on the map.

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

This code is working fine however I have a new requirement.

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

One of the points may have a precision associated with it. If this is the case then I draw a circle around the point with the radius set to the precision value. Again this works fine however my bounds checking is now not doing what I want it to do. I want to have the bounding box include the complete circle.

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

This requires an algorithm to take a point x and calculate the point y that would be z metres north of x and also z metres south of x.

这需要一个算法来获取一个点 x 并计算在 x 以北 z 米处以及 x 以南 z 米处的点 y。

Does anyone have this algorithm, preferably in C#. I did find a generic algorithm herebut I appear to have not implemented this correctly as the answers I am getting are 1000s of km adrift.

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

This is the Generic example

这是通用示例

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

And this is my C# translation.

这是我的 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;
    }

I am calling this function with a bearing of 0 to calculate the new northerly position and a value of 180 to calculate the new southerly position.

我用 0 的方位角调用这个函数来计算新的北向位置，用 180 的值来计算新的南向位置。

Can anyone either see what I have done wrong or perhaps provide a known working algorithm?

谁能看到我做错了什么，或者提供一个已知的工作算法？

采纳答案by Sam152

If you have a given latitude and longitude you can calculate the correct latitude and longitude of an x-km change in latitude like so:

如果您有给定的纬度和经度，则可以像这样计算 x 公里纬度变化的正确纬度和经度：

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.

The same can apply to longitude. If you have the total distance plus the change you can calculate the total degrees in a similar fashion.

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

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.

Again, these calculations should work, but I'm running off pure intuition here, but the logic does seem to hold true.

同样，这些计算应该有效，但我在这里完全脱离了直觉，但逻辑似乎确实成立。

Edit: As pointed out by Skizz 40,075 needs to be adjusted to the circumference of the earth at any given latitude using 2.pi.r.cos(lat) or 40074.cos(lat)

编辑：正如 Skizz 所指出的，40,075 需要使用 2.pi.r.cos(lat) 或 40074.cos(lat) 调整到任何给定纬度的地球周长

回答by Pablo Cabrera

It is more accurate if you first reproject it to UTMand then check the distance.

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

Hope this helps

希望这可以帮助

回答by Erich Mirabal

I have a very similar piece of code. It got me very close results when compared to another implementation.

我有一段非常相似的代码。与另一个实现相比，它给了我非常接近的结果。

I think the problem with yours is that you are using "distance" as linear distance in meters instead of angular distance in radians.

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

/// <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);
}

Where

在哪里

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

and LatLonAlt is in degrees/meters (conversion takes place internally). Adjust as needed.

和 LatLonAlt 以度/米为单位（转换发生在内部）。根据需要进行调整。

I assume you can figure out what the value for UnitConstants.DegreesToRadiansis :)

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

回答by ryan_s

I'm not sure if I'm missing something here, but I think the question could be rephrased as, "I have a lat/lon point, and I want to find the point x meters north and x meters south of that point."

我不确定我是否在这里遗漏了什么，但我认为这个问题可以改写为，“我有一个纬度/经度点，我想找到该点以北 x 米和以南 x 米的点。 ”

If that's the question then you don't need to find a new longitude (which makes things simpler), you just need a new latitude. A degree of latitude is roughly 60 nautical miles long anywhere on Earth, and a nautical mile is 1,852 meters. So, for new latitudes x meters north and south:

如果这是问题，那么您不需要找到新的经度（这使事情更简单），您只需要一个新的纬度。在地球上的任何地方，纬度大约为 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)

This is not completely accurate because the Earth is not a perfect sphere with exactly 60 nautical miles between each degree of latitude. However, the other answers assume that lines of latitude are equidistant, so I'm assuming you don't care about that. If you're interested in how much error that might introduce, there is a nice table on Wikipedia that shows "Surface distance per 1° change in latitude" for different latitudes at this link:

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

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

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

回答by Sy Moen

There are problems with the two equations on Ed William's rather awesome site... but I didn't analyze them to see why.

Ed William 相当棒的网站上的两个方程存在问题......但我没有分析它们以了解原因。

A third equation that I found hereseems to give proper results.

我在这里找到的第三个方程似乎给出了正确的结果。

Here is the test case in php... the third equation is correct, the first two give wildly incorrect values for longitude.

这是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";
?>

Note I recieved by email from the author (Ed Williams) of the first two equations:

请注意，我收到了前两个方程的作者 (Ed Williams) 发来的电子邮件：

From my "implementation notes":

Note on the mod function. This appears to be implemented differently in different languages, with differing conventions on whether the sign of the result follows the sign of the divisor or the dividend. (We want the sign to follow the divisor or be Euclidean. C's fmod and Java's % do not work.) In this document, Mod(y,x) is the remainder on dividing y by x and always lies in the range 0 <= mod < x. For instance: mod(2.3,2.)=0.3 and mod(-2.3,2.)=1.7

If you have a floor function (int in Excel), that returns floor(x)= "largest integer less than or equal to x" e.g. floor(-2.3)=-3 and floor(2.3) =2

从我的“实施笔记”：

注意 mod 函数。这在不同的语言中似乎以不同的方式实现，对于结果的符号是跟随除数还是被除数的符号有不同的约定。（我们希望符号跟随除数或欧几里得。C 的 fmod 和 Java 的 % 不起作用。）在本文档中，Mod(y,x) 是 y 除以 x 的余数，并且始终位于 0 <= 范围内模 < 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)

The following should work in the absence of a floor function- regardless of whether "int" truncates or rounds downward:

以下应该在没有地板功能的情况下工作 - 无论“int”是截断还是向下舍入：

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

php is like fmod in C and does it "wrong" for my purposes.

php 就像 C 中的 fmod 并且它对我的目的来说是“错误的”。

回答by Zar Shardan

For lazy people, (like me ;) ) a copy-paste solution, Erich Mirabal's version with very minor changes:

对于懒惰的人，（像我一样；））复制粘贴解决方案，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;
}

Usage:

用法：

[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}");
    }
}