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

时间:2022-09-15 11:49:36

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 here but I appear to have not implemented this correctly as the answers I am getting are 1000s of km adrift.

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

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?

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

7 个解决方案

#1


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-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.

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)在任意给定的纬度调整到地球的周长。

#2


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.DegreesToRadians is :)

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

#3


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

#4


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

#5


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

在埃德威廉的相当令人敬畏的网站上有两个方程式存在问题...但我没有分析它们以了解原因。

A third equation that I found here seems 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

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

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

如果你有一个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:

以下应该在没有floor函数的情况下工作 - 无论“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一样,为我的目的而言它“错误”。

#6


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

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

Hope this helps

希望这可以帮助

#7


For what it's worth, I have an example in PHP which can do what the OP is requesting. In my example it is drawing a box around a starting lat/long coordinate, but the code can easily be used to get a single point, X number of km or miles away.

对于它的价值,我在PHP中有一个例子可以做OP请求的事情。在我的示例中,它是在起始纬度/长坐标周围绘制一个框,但代码可以很容易地用于获得单个点,X数千公里或数英里之外。

http://www.richardpeacock.com/blog/2011/11/draw-box-around-coordinate-google-maps-based-miles-or-kilometers

#1


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-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.

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)在任意给定的纬度调整到地球的周长。

#2


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.DegreesToRadians is :)

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

#3


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

#4


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

#5


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

在埃德威廉的相当令人敬畏的网站上有两个方程式存在问题...但我没有分析它们以了解原因。

A third equation that I found here seems 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

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

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

如果你有一个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:

以下应该在没有floor函数的情况下工作 - 无论“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一样,为我的目的而言它“错误”。

#6


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

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

Hope this helps

希望这可以帮助

#7


For what it's worth, I have an example in PHP which can do what the OP is requesting. In my example it is drawing a box around a starting lat/long coordinate, but the code can easily be used to get a single point, X number of km or miles away.

对于它的价值,我在PHP中有一个例子可以做OP请求的事情。在我的示例中,它是在起始纬度/长坐标周围绘制一个框,但代码可以很容易地用于获得单个点,X数千公里或数英里之外。

http://www.richardpeacock.com/blog/2011/11/draw-box-around-coordinate-google-maps-based-miles-or-kilometers