Как мне найти широту / длину, которая находится в х км к северу от данной широты / длины?

Question

У меня есть код C #, который генерирует карты Google. Этот код просматривает все точки, которые мне нужно нанести на карту, и затем определяет границы прямоугольника, чтобы включить эти точки. Затем он передает эти границы в API Карт Google, чтобы соответствующим образом установить уровень масштабирования, чтобы показать все точки на карте.

Этот код работает нормально, но у меня есть новое требование.

С одной из точек может быть связана точность. Если это так, то я рисую круг вокруг точки с радиусом, установленным на значение точности. Опять же, это работает нормально, однако моя проверка границ теперь не делает то, что я хочу. Я хочу, чтобы в ограничительной рамке был полный круг.

Для этого требуется алгоритм, чтобы взять точку x и вычислить точку y, которая будет в z метрах к северу от x, а также в z метрах к югу от x.

У кого-нибудь есть этот алгоритм, желательно на C #. Я нашел общий алгоритм здесь , но, похоже, я не реализовал это правильно, так как ответы, которые я получаю, - это тысячи километров по течению.

Это общий пример

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 для расчета новой южной позиции.

Может ли кто-нибудь увидеть, что я сделал неправильно, или, возможно, предоставить известный алгоритм работы?

Answer 1

У меня очень похожий кусок кода. Это дало мне очень близкие результаты по сравнению с другой реализацией.

Я думаю, что проблема с вами в том, что вы используете «расстояние» как линейное расстояние в метрах вместо углового расстояния в радианах.

/// <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:)

Answer 2

Для ленивых (как и я;)) решение для копирования-вставки, версия Эриха Мирабала с очень незначительными изменениями:

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

Answer 3

Я не уверен, что я что-то здесь упускаю, но я думаю, что вопрос можно перефразировать следующим образом: «У меня есть точка широты и долготы, и я хочу найти точку х метров севернее и х метров южнее этот момент. "

Если это вопрос, то вам не нужно искать новую долготу (которая делает вещи проще), вам просто нужна новая широта. Степень широты составляет примерно 60 морских миль в длину в любом месте на Земле, а морская миля составляет 1852 метра. Итак, для новых широт х метров на север и юг:

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

Answer 4

Если у вас есть заданная широта и долгота, вы можете рассчитать правильную широту и долготу изменения широты на 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.

То же самое можно применить к долготе. Если у вас есть общее расстояние плюс изменение, вы можете рассчитать общее количество градусов аналогичным образом.

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 (широта) или 40074.cos (широта)

Answer 5

Есть проблемы с двумя уравнениями на Довольно удивительный сайт Эда Уильяма ... но я не анализировал их, чтобы понять почему.

Третье уравнение, которое я нашел здесь , похоже, дает правильные результаты.

Вот тестовый пример в 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";
?>

Примечание, которое я получил по электронной почте от автора (Эд Уильямс) первых двух уравнений:

Из моих «заметок о реализации»:

Обратите внимание на функцию мода. По-видимому, это реализовано по-другому в разных языках, с различными соглашениями о том, является ли знак результат следует за знаком делителя или дивидендом. (Мы хотим, чтобы знак следовать за делителем или быть евклидовым. Fmod и Java% C не работают.) В этом документе Mod (y, x) является остатком от деления y на x и всегда лежит в диапазоне 0 <= mod <x. Например: мод (2.3,2.) = 0,3 и мод (-2.3,2.) = 1,7 </p>

Если у вас есть функция floor (int в Excel), она возвращает floor (x) = "наибольшее целое число меньше или равно х", например этаж (-2,3) = - 3 а этаж (2.3) = 2

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

Следующее должно работать при отсутствии функции пола - независимо от "int" усекает или округляет вниз:

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

php похож на fmod в C и делает это «неправильно» для моих целей.

Answer 6

Точнее, если сначала спроецировать его на UTM , а затем проверить расстояние.

Надеюсь, это поможет