31

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

这段代码运行良好,但是我有一个新要求。

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

这需要一种算法来获取一个点 x 并计算点 y,该点位于 x 以北 z 米处和 x 以南 z 米处。

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

这是通用示例

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 的值来计算新的南向位置。

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

4

7 回答 7

26

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.

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

于 2009-07-14T13:34:32.600 回答
14

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

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}");
    }
}
于 2016-07-06T21:12:04.227 回答
7

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

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:

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:

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

于 2009-07-18T09:23:53.733 回答
6

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

new-lat = ((old-km-north + x-km-change)/40,075) * 360)
           ^ is the ratio of the                  ^ times the ratio of the circle
           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
           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)

于 2009-07-14T12:55:56.723 回答
4

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
            $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:

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(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:

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.

于 2011-04-11T20:20:36.323 回答
0

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

Hope this helps

于 2009-07-14T13:21:00.760 回答
0

For people who want a java version Eirch's code

/**
 * move latlng point by rang and bearing
 *
 * @param latLng  point
 * @param range   range in meters
 * @param bearing bearing in degrees
 * @return new LatLng
 */
public static LatLng moveLatLng(LatLng latLng, double range, double bearing) {
    double EarthRadius = 6378137.0;
    double DegreesToRadians = Math.PI / 180.0;
    double RadiansToDegrees = 180.0 / Math.PI;

    final double latA = latLng.latitude * DegreesToRadians;
    final double lonA = latLng.longitude * DegreesToRadians;
    final double angularDistance = range / EarthRadius;
    final double trueCourse = bearing * DegreesToRadians;

    final double lat = Math.asin(
            Math.sin(latA) * Math.cos(angularDistance) +
                    Math.cos(latA) * Math.sin(angularDistance) * Math.cos(trueCourse));

    final double dlon = Math.atan2(
            Math.sin(trueCourse) * Math.sin(angularDistance) * Math.cos(latA),
            Math.cos(angularDistance) - Math.sin(latA) * Math.sin(lat));

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

    return new LatLng(lat * RadiansToDegrees, lon * RadiansToDegrees);
}
于 2020-09-16T06:51:37.183 回答