Определение, находится ли координата GPS в пределах многоугольника на карте

Question

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

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

Есть пара вопросов, подобных этойздесь, на StackOverflow, но они описывают решение только в общих терминах и для 2D-плоскости, в то время как я ищу существующую реализацию, которая поддерживает полигоны, определенные парами широта / долгота в WGS 84 .

Какие API или сервисы существуют для таких тестов на столкновение?

Answer 1

Вот Java-программа, которая использует функцию, которая будет возвращать true, если широта / долгота найдена внутри многоугольника, определенного списком широт / долгот, с демонстрацией состояния Флориды.

Я не уверен, имеет ли это дело с тем, что широта / долгота системы GPS не является координатной плоскостью x / y. Для моего использования я продемонстрировал, что это работает (я думаю, что если вы укажете достаточное количество точек в ограничительной рамке, это размывает эффект того, что земля является сферой, а прямые линии между двумя точками на земле не являются прямой линией стрелки). .

Сначала укажите точки, которые составляют угловые точки многоугольника, он может иметь вогнутые и выпуклые углы. Координаты, которые я использую ниже, отслеживают периметр штата Флорида.

метод coordinate_is_inside_polygon использует алгоритм, который я не совсем понимаю. Вот официальное объяснение от источника, откуда я его взял:

"... решение, предложенное Филиппом Реверди, состоит в том, чтобы вычислить сумму углов, образованных между контрольной точкой и каждой парой точек, составляющих многоугольник. Если эта сумма равна 2pi, то точка является внутренней точкой, если 0 тогда точка является внешней точкой. Это также работает для полигонов с отверстиями, если полигон определен путем, состоящим из совпадающих ребер в и из отверстия, как это принято во многих пакетах САПР. "

Мои модульные тесты показывают, что он работает надежно, даже когда ограничивающий прямоугольник имеет форму 'C' или даже имеет форму Torus . (Мои модульные тесты проверяют множество точек во Флориде и проверяют, возвращает ли функция значение true. И я выбираю несколько координат в любой точке мира и проверяю, что она возвращает false. Я выбираю места по всему миру, которые могут его запутать.

Я не уверен, что это сработает, если ограничивающий прямоугольник пересекает экватор, меридиан или любую область, где координаты меняются от -180 -> 180, -90 -> 90. Или ваш многоугольник вращается вокруг Земли вокруг северного / южного полюса. Для меня это нужно только для работы по периметру Флориды. Если вам нужно определить многоугольник, который охватывает землю или пересекает эти линии, вы можете обойти его, создав два многоугольника, один из которых представляет область на одной стороне меридиана, а другой - область на другой стороне и проверит, является ли ваша точка находится в любой из этих точек.

Вот где я нашел этот алгоритм: Определение, лежит ли точка внутри многоугольника - Решение 2

Запустите его для себя, чтобы дважды проверить его.

Поместите это в файл с именем Runner.java

import java.util.ArrayList;
public class Runner
{
    public static double PI = 3.14159265;
    public static double TWOPI = 2*PI;
    public static void main(String[] args) {
    ArrayList<Double> lat_array = new ArrayList<Double>();
    ArrayList<Double> long_array = new ArrayList<Double>();

    //This is the polygon bounding box, if you plot it, 
    //you'll notice it is a rough tracing of the parameter of 
    //the state of Florida starting at the upper left, moving 
    //clockwise, and finishing at the upper left corner of florida.

    ArrayList<String> polygon_lat_long_pairs = new ArrayList<String>();
    polygon_lat_long_pairs.add("31.000213,-87.584839");  
    //lat/long of upper left tip of florida.
    polygon_lat_long_pairs.add("31.009629,-85.003052");
    polygon_lat_long_pairs.add("30.726726,-84.838257");
    polygon_lat_long_pairs.add("30.584962,-82.168579");
    polygon_lat_long_pairs.add("30.73617,-81.476441");  
    //lat/long of upper right tip of florida.
    polygon_lat_long_pairs.add("29.002375,-80.795288");
    polygon_lat_long_pairs.add("26.896598,-79.938355");
    polygon_lat_long_pairs.add("25.813738,-80.059204");
    polygon_lat_long_pairs.add("24.93028,-80.454712");
    polygon_lat_long_pairs.add("24.401135,-81.817017");
    polygon_lat_long_pairs.add("24.700927,-81.959839");
    polygon_lat_long_pairs.add("24.950203,-81.124878");
    polygon_lat_long_pairs.add("26.0015,-82.014771");
    polygon_lat_long_pairs.add("27.833247,-83.014527");
    polygon_lat_long_pairs.add("28.8389,-82.871704");
    polygon_lat_long_pairs.add("29.987293,-84.091187");
    polygon_lat_long_pairs.add("29.539053,-85.134888");
    polygon_lat_long_pairs.add("30.272352,-86.47522");
    polygon_lat_long_pairs.add("30.281839,-87.628784");

    //Convert the strings to doubles.       
    for(String s : polygon_lat_long_pairs){
        lat_array.add(Double.parseDouble(s.split(",")[0]));
        long_array.add(Double.parseDouble(s.split(",")[1]));
    }

   //prints TRUE true because the lat/long passed in is
    //inside the bounding box.
    System.out.println(coordinate_is_inside_polygon(
            25.7814014D,-80.186969D,
            lat_array, long_array));

    //prints FALSE because the lat/long passed in 
    //is Not inside the bounding box.
    System.out.println(coordinate_is_inside_polygon(
            25.831538D,-1.069338D, 
            lat_array, long_array));

}
public static boolean coordinate_is_inside_polygon(
    double latitude, double longitude, 
    ArrayList<Double> lat_array, ArrayList<Double> long_array)
{       
       int i;
       double angle=0;
       double point1_lat;
       double point1_long;
       double point2_lat;
       double point2_long;
       int n = lat_array.size();

       for (i=0;i<n;i++) {
          point1_lat = lat_array.get(i) - latitude;
          point1_long = long_array.get(i) - longitude;
          point2_lat = lat_array.get((i+1)%n) - latitude; 
          //you should have paid more attention in high school geometry.
          point2_long = long_array.get((i+1)%n) - longitude;
          angle += Angle2D(point1_lat,point1_long,point2_lat,point2_long);
       }

       if (Math.abs(angle) < PI)
          return false;
       else
          return true;
}

public static double Angle2D(double y1, double x1, double y2, double x2)
{
   double dtheta,theta1,theta2;

   theta1 = Math.atan2(y1,x1);
   theta2 = Math.atan2(y2,x2);
   dtheta = theta2 - theta1;
   while (dtheta > PI)
      dtheta -= TWOPI;
   while (dtheta < -PI)
      dtheta += TWOPI;

   return(dtheta);
}

public static boolean is_valid_gps_coordinate(double latitude, 
    double longitude)
{
    //This is a bonus function, it's unused, to reject invalid lat/longs.
    if (latitude > -90 && latitude < 90 && 
            longitude > -180 && longitude < 180)
    {
        return true;
    }
    return false;
}
}

Магия демона должна быть проверена юнитом. Поместите это в файл с именем MainTest.java, чтобы убедиться, что он работает для вас

import java.util.ArrayList;
import org.junit.Test;
import static org.junit.Assert.*;

public class MainTest {
@Test
public void test_lat_long_in_bounds(){
    Runner r = new Runner();
    //These make sure the lat/long passed in is a valid gps 
    //lat/long coordinate.  These should be valid. 
    assertTrue(r.is_valid_gps_coordinate(25, -82));
    assertTrue(r.is_valid_gps_coordinate(-25, -82));
    assertTrue(r.is_valid_gps_coordinate(25, 82));
    assertTrue(r.is_valid_gps_coordinate(-25, 82));
    assertTrue(r.is_valid_gps_coordinate(0, 0));
    assertTrue(r.is_valid_gps_coordinate(89, 179));
    assertTrue(r.is_valid_gps_coordinate(-89, -179));
    assertTrue(r.is_valid_gps_coordinate(89.999, 179));
    //If your bounding box crosses the equator or prime meridian, 
    then you have to test for those situations still work.
}
@Test
public void realTest_for_points_inside()
{
    ArrayList<Double> lat_array = new ArrayList<Double>();
    ArrayList<Double> long_array = new ArrayList<Double>();
    ArrayList<String> polygon_lat_long_pairs = new ArrayList<String>();
    //upper left tip of florida.
    polygon_lat_long_pairs.add("31.000213,-87.584839");
    polygon_lat_long_pairs.add("31.009629,-85.003052");
    polygon_lat_long_pairs.add("30.726726,-84.838257");
    polygon_lat_long_pairs.add("30.584962,-82.168579");
    polygon_lat_long_pairs.add("30.73617,-81.476441");  
    //upper right tip of florida.
    polygon_lat_long_pairs.add("29.002375,-80.795288");
    polygon_lat_long_pairs.add("26.896598,-79.938355");
    polygon_lat_long_pairs.add("25.813738,-80.059204");
    polygon_lat_long_pairs.add("24.93028,-80.454712");
    polygon_lat_long_pairs.add("24.401135,-81.817017");
    polygon_lat_long_pairs.add("24.700927,-81.959839");
    polygon_lat_long_pairs.add("24.950203,-81.124878");
    polygon_lat_long_pairs.add("26.0015,-82.014771");
    polygon_lat_long_pairs.add("27.833247,-83.014527");
    polygon_lat_long_pairs.add("28.8389,-82.871704");
    polygon_lat_long_pairs.add("29.987293,-84.091187");
    polygon_lat_long_pairs.add("29.539053,-85.134888");
    polygon_lat_long_pairs.add("30.272352,-86.47522");
    polygon_lat_long_pairs.add("30.281839,-87.628784");

    for(String s : polygon_lat_long_pairs){
        lat_array.add(Double.parseDouble(s.split(",")[0]));
        long_array.add(Double.parseDouble(s.split(",")[1]));
    }

    Runner r = new Runner();
    ArrayList<String> pointsInside = new ArrayList<String>();
    pointsInside.add("30.82112,-87.255249");
    pointsInside.add("30.499804,-86.8927");
    pointsInside.add("29.96826,-85.036011");
    pointsInside.add("30.490338,-83.981323");
    pointsInside.add("29.825395,-83.344116");
    pointsInside.add("30.215406,-81.828003");
    pointsInside.add("29.299813,-82.728882");
    pointsInside.add("28.540135,-81.212769");
    pointsInside.add("27.92065,-82.619019");
    pointsInside.add("28.143691,-81.740113");
    pointsInside.add("27.473186,-80.718384");
    pointsInside.add("26.769154,-81.729126");
    pointsInside.add("25.853292,-80.223999");
    pointsInside.add("25.278477,-80.707398");
    pointsInside.add("24.571105,-81.762085");   //bottom tip of keywest
    pointsInside.add("24.900388,-80.663452");
    pointsInside.add("24.680963,-81.366577");

    for(String s : pointsInside)
    {
        assertTrue(r.coordinate_is_inside_polygon(
            Double.parseDouble(s.split(",")[0]), 
            Double.parseDouble(s.split(",")[1]), 
            lat_array, long_array));
    }
}

@Test
public void realTest_for_points_outside()
{
    ArrayList<Double> lat_array = new ArrayList<Double>();
    ArrayList<Double> long_array = new ArrayList<Double>();

    ArrayList<String> polygon_lat_long_pairs = new ArrayList<String>();
    //upper left tip, florida.
    polygon_lat_long_pairs.add("31.000213,-87.584839");
    polygon_lat_long_pairs.add("31.009629,-85.003052");
    polygon_lat_long_pairs.add("30.726726,-84.838257");
    polygon_lat_long_pairs.add("30.584962,-82.168579");
    polygon_lat_long_pairs.add("30.73617,-81.476441");
    //upper right tip, florida.
    polygon_lat_long_pairs.add("29.002375,-80.795288");
    polygon_lat_long_pairs.add("26.896598,-79.938355");
    polygon_lat_long_pairs.add("25.813738,-80.059204");
    polygon_lat_long_pairs.add("24.93028,-80.454712");
    polygon_lat_long_pairs.add("24.401135,-81.817017");
    polygon_lat_long_pairs.add("24.700927,-81.959839");
    polygon_lat_long_pairs.add("24.950203,-81.124878");
    polygon_lat_long_pairs.add("26.0015,-82.014771");
    polygon_lat_long_pairs.add("27.833247,-83.014527");
    polygon_lat_long_pairs.add("28.8389,-82.871704");
    polygon_lat_long_pairs.add("29.987293,-84.091187");
    polygon_lat_long_pairs.add("29.539053,-85.134888");
    polygon_lat_long_pairs.add("30.272352,-86.47522");
    polygon_lat_long_pairs.add("30.281839,-87.628784");

    for(String s : polygon_lat_long_pairs)
    {
        lat_array.add(Double.parseDouble(s.split(",")[0]));
        long_array.add(Double.parseDouble(s.split(",")[1]));
    }

    Runner r = new Runner();

    ArrayList<String> pointsOutside = new ArrayList<String>();
    pointsOutside.add("31.451159,-87.958374");
    pointsOutside.add("31.319856,-84.607544");
    pointsOutside.add("30.868282,-84.717407");
    pointsOutside.add("31.338624,-81.685181");
    pointsOutside.add("29.452991,-80.498657");
    pointsOutside.add("26.935783,-79.487915");
    pointsOutside.add("25.159207,-79.916382");
    pointsOutside.add("24.311058,-81.17981");
    pointsOutside.add("25.149263,-81.838989");
    pointsOutside.add("27.726326,-83.695679");
    pointsOutside.add("29.787263,-87.024536");
    pointsOutside.add("29.205877,-62.102052");
    pointsOutside.add("14.025751,-80.690919");
    pointsOutside.add("29.029276,-90.805666");
    pointsOutside.add("-12.606032,-70.151369");
    pointsOutside.add("-56.520716,-172.822269");
    pointsOutside.add("-75.89666,9.082024");
    pointsOutside.add("-24.078567,142.675774");
    pointsOutside.add("84.940737,177.480462");
    pointsOutside.add("47.374545,9.082024");
    pointsOutside.add("25.831538,-1.069338");
    pointsOutside.add("0,0");

    for(String s : pointsOutside){
        assertFalse(r.coordinate_is_inside_polygon(
            Double.parseDouble(s.split(",")[0]),
            Double.parseDouble(s.split(",")[1]), lat_array, long_array));
    }
}
}
//The list of lat/long inside florida bounding box all return true.
//The list of lat/long outside florida bounding box all return false.

Я использовал eclipse IDE для запуска java с использованием java 1.6.0. Для меня все юнит-тесты пройдены. Вам нужно включить jar-файл junit 4 в ваш путь к классам или импортировать его в Eclipse.

Answer 2

Сначала я думал так же, как Шэб (его предложение называется Алгоритм каста лучей ), но у меня были мысли, похожие на Spacedman:

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

Я реализовал и протестировал математически правильный способ сделать это, например, пересекая большие круги и определяя, находится ли одна из двух пересекающихся точек на обеих дугах.(Примечание: я выполнил шаги, описанные здесь , но обнаружил несколько ошибок: функция sign отсутствует в конце шага 6 (непосредственно перед arcsin), а окончательный тест представляет собой числовой мусор(поскольку вычитание плохо обусловлено), используйте скорее L_1T >= max(L_1a, L_1b), чтобы проверить, находится ли S1 на первой дуге и т. д.)

Это также чрезвычайно медленный и числовой кошмар (оценивает ~ 100тригонометрические функции, между прочим);он оказался непригодным для использования в наших встроенных системах.

Однако есть хитрость : если рассматриваемая область достаточно мала, просто выполните стандартную картографическую проекцию, например, сферическая проекция Меркатора каждой точки:

// latitude, longitude in radians
x = longitude;
y = log(tan(pi/4 + latitude/2));

Затем можно применить приведение лучей, где пересечение дуг проверяется этой функцией:

public bool ArcsIntersecting(double x1, double y1, double x2, double y2, 
  double x3, double y3, double x4, double y4)
    {

    double vx1 = x2 - x1;
    double vy1 = y2 - y1;

    double vx2 = x4 - x3;
    double vy2 = y4 - y3;

    double denom = vx1 * vy2 - vx2 * vy1;

    if (denom == 0) { return false; } // edges are parallel

    double t1 = (vx2 * (y1 - y3) - vy2 * (x1 - x3)) / denom;

    double t2;

    if (vx2 != 0) { t2 = (x1 - x3 + t1 * vx1) / vx2; }
    else if (vy2 != 0) { t2 = (y1 - y3 + t1 * vy1) / vy2; }
    else { return false; } // edges are matching

    return min(t1, t2) >= 0 && max(t1, t2) <= 1;
}

Answer 3

Если у вас есть координаты WGS84 на сфере, то ваш многоугольник делит сферу на две области - как мы узнаем, какая область находится «внутри», а какая «снаружи» многоугольника?Вопрос по сути бессмысленный!

Например, предположим, что многоугольник образовал линию экватора - северное полушарие «внутри» или «снаружи»?

Answer 4

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

Answer 5

Версия JavaScript -

{
const PI = 3.14159265;
const TWOPI = 2*PI;
function isCoordinateInsidePitch(latitude, longitude, latArray, longArray)
{       
       let angle=0;
       let p1Lat;
       let p1Long;
       let p2Lat;
       let p2Long;
       let n = latArray.length;

       for (let i = 0; i < n; i++) {
          p1Lat = latArray[i] - latitude;
          p1Long = longArray[i] - longitude;
          p2Lat = latArray[(i+1)%n] - latitude;
          p2Long = longArray[(i+1)%n] - longitude;
          angle += angle2D(p1Lat,p1Long,p2Lat,p2Long);
       }

       return !(Math.abs(angle) < PI);
}

function angle2D(y1, x1, y2, x2)
{
   let dtheta,theta1,theta2;

   theta1 = Math.atan2(y1,x1);
   theta2 = Math.atan2(y2,x2);
   dtheta = theta2 - theta1;
   while (dtheta > PI)
      dtheta -= TWOPI;
   while (dtheta < -PI)
      dtheta += TWOPI;

   return dtheta;
}

function isValidCoordinate(latitude,longitude)
{
    return (
    latitude !== '' && longitude !== '' && !isNaN(latitude) 
     && !isNaN(longitude) && latitude > -90 &&
     latitude < 90 && longitude > -180 && longitude < 180
     )
}
let latArray = [32.10458, 32.10479, 32.1038, 32.10361];
let longArray = [34.86448, 34.86529, 34.86563, 34.86486];
// true
console.log(isCoordinateInsidePitch(32.104447, 34.865108,latArray, longArray));
// false
// isCoordinateInsidePitch(32.104974, 34.864576,latArray, longArray);
// true
// isValidCoordinate(0, 0)
// true
// isValidCoordinate(32.104974, 34.864576)
}

Answer 6

код Runner.Java в VB.NET

В интересах пользователей .NET тот же код помещен в VB.NET. Попробовал и довольно быстро. Пробовал с 350000 записей, он заканчивается всего за несколько минут. Но, как сказал автор, мне еще предстоит протестировать сценарии, пересекающие экватор, мультизоны и т. Д.

Использование

If coordinate_is_inside_polygon(CurLat, CurLong, Lat_Array, Long_Array) Then
    MsgBox("Location " & CurLat & "," & CurLong & " is within polygon boundary")
Else
    MsgBox("Location " & CurLat & "," & CurLong & " is NOT within polygon boundary")
End If

'Функции

Public Function coordinate_is_inside_polygon(ByVal latitude As Double, ByVal longitude As Double, ByVal lat_array() As Double, ByVal long_array() As Double) As Boolean
    Dim i As Integer
    Dim angle As Double = 0
    Dim point1_lat As Double
    Dim point1_long As Double
    Dim point2_lat As Double
    Dim point2_long As Double
    Dim n As Integer = lat_array.Length()
    For i = 0 To n - 1
        point1_lat = lat_array(i) - latitude
        point1_long = long_array(i) - longitude
        point2_lat = lat_array((i + 1) Mod n) - latitude
        point2_long = long_array((i + 1) Mod n) - longitude
        angle += Angle2D(point1_lat, point1_long, point2_lat, point2_long)
    Next

    If Math.Abs(angle) < PI Then Return False Else Return True
End Function

Public Function Angle2D(ByVal y1 As Double, ByVal x1 As Double, ByVal y2 As Double, ByVal x2 As Double) As Double
    Dim dtheta, theta1, theta2 As Double
    theta1 = Math.Atan2(y1, x1)
    theta2 = Math.Atan2(y2, x2)
    dtheta = theta2 - theta1
    While dtheta > PI
         dtheta -= TWOPI
    End While

    While dtheta < -PI
        dtheta += TWOPI
    End While
    Return (dtheta)
End Function



 Public Function is_valid_gps_coordinate(ByVal latitude As Double, ByVal longitude As Double) As Boolean
        If latitude > -90 AndAlso latitude < 90 AndAlso longitude > -180 AndAlso longitude < 180 Then
            Return True
        End If

        Return False
End Function

Answer 7

Вот алгоритм, написанный на Go: Он принимает координаты точки в формате [lat, long] и многоугольник в формате [[lat, long], [lat, long] ...]. Алгоритм соединит первую и последнюю точку в срезе полигона

import "math"

// ContainsLocation determines whether the point is inside the polygon
func ContainsLocation(point []float64, polygon [][]float64, geodesic 
   bool) bool {
    size := len(polygon)
    if size == 0 {
        return false
    }

    var (
        lat2, lng2, dLng3 float64
    )

    lat3 := toRadians(point[0])
    lng3 := toRadians(point[1])
    prev := polygon[size-1]
    lat1 := toRadians(prev[0])
    lng1 := toRadians(prev[1])
    nIntersect := 0

    for _, v := range polygon {
        dLng3 = wrap(lng3-lng1, -math.Pi, math.Pi)
        // Special case: point equal to vertex is inside.
        if lat3 == lat1 && dLng3 == 0 {
            return true
        }
        lat2 = toRadians(v[0])
        lng2 = toRadians(v[1])
        // Offset longitudes by -lng1.
        if intersects(lat1, lat2, wrap(lng2-lng1, -math.Pi, math.Pi), lat3, dLng3, geodesic) {
            nIntersect++
        }
        lat1 = lat2
        lng1 = lng2
    }

    return (nIntersect & 1) != 0
}

func toRadians(p float64) float64 {
    return p * (math.Pi / 180.0)
}

func wrap(n, min, max float64) float64 {
    if n >= min && n < max {
        return n
    }
    return mod(n-min, max-min) + min
}

func mod(x, m float64) float64 {
    return math.Remainder(math.Remainder(x, m)+m, m)
}

func intersects(lat1, lat2, lng2, lat3, lng3 float64, geodesic bool) bool {
    // Both ends on the same side of lng3.
    if (lng3 >= 0 && lng3 >= lng2) || (lng3 < 0 && lng3 < lng2) {
        return false
    }
    // Point is South Pole.
    if lat3 <= -math.Pi/2 {
        return false
    }
    // Any segment end is a pole.
    if lat1 <= -math.Pi/2 || lat2 <= -math.Pi/2 || lat1 >= math.Pi/2 || lat2 >= math.Pi/2 {
        return false
    }
    if lng2 <= -math.Pi {
        return false
    }
    linearLat := (lat1*(lng2-lng3) + lat2*lng3) / lng2
    // Northern hemisphere and point under lat-lng line.
    if lat1 >= 0 && lat2 >= 0 && lat3 < linearLat {
        return false
    }
    // Southern hemisphere and point above lat-lng line.
    if lat1 <= 0 && lat2 <= 0 && lat3 >= linearLat {
        return true
    }
    // North Pole.
    if lat3 >= math.Pi/2 {
        return true
    }

    // Compare lat3 with latitude on the GC/Rhumb segment corresponding to lng3.
    // Compare through a strictly-increasing function (tan() or mercator()) as convenient.
    if geodesic {
        return math.Tan(lat3) >= tanLatGC(lat1, lat2, lng2, lng3)
    }
    return mercator(lat3) >= mercatorLatRhumb(lat1, lat2, lng2, lng3)
}

func tanLatGC(lat1, lat2, lng2, lng3 float64) float64 {
    return (math.Tan(lat1)*math.Sin(lng2-lng3) + math.Tan(lat2)*math.Sin(lng3)) / math.Sin(lng2)
}

func mercator(lat float64) float64 {
    return math.Log(math.Tan(lat*0.5 + math.Pi/4))
}

func mercatorLatRhumb(lat1, lat2, lng2, lng3 float64) float64 {
    return (mercator(lat1)*(lng2-lng3) + mercator(lat2)*lng3) / lng2
}

Answer 8

Предполагая, что вы справляетесь со случаем обтекания меридиана и пересечения экватора (путем добавления смещений) - вы не можете просто рассматривать это как простую двухмерную точку в многоугольнике?