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

StackOverflow https://stackoverflow.com/questions/1125144

  •  13-09-2019
  •  | 
  •  

Вопрос

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

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

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

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

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

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

Это было полезно?

Решение

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

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

Другие советы

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

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

/// <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 есть :)

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

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

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

Если это вопрос, то вам не нужно находить новую долготу (что упрощает задачу), вам просто нужна новая широта.Градус широты имеет длину примерно 60 морских миль в любой точке Земли, а морская миля равна 1852 метрам.Итак, для новых широт 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

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

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

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

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

Из моих "заметок по внедрению":

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

Если у вас есть функция floor (int в Excel), которая возвращает floor(x)= "наибольшее целое число, меньшее или равное x", напримерэтаж (-2,3)=-3 и этаж (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 похож на fmod в C и делает это "неправильно" для моих целей.

Это будет более точно, если вы сначала перепроектируете его на UTM а затем проверьте расстояние.

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

Лицензировано под: CC-BY-SA с атрибуция
Не связан с StackOverflow
scroll top