Расчет площади многоугольника с использованием широты и долготы, полученных из декартового пространства и файла мира

Question

Учитывая серию координатных пар GPS, мне нужно вычислить площадь многоугольника (n-gon). Это относительно мало (не более 50000 кв. Футов). Геокоды создаются путем применения аффинного преобразования с данными из файла мира.

Я попытался использовать двухэтапный подход путем преобразования геокодов в декартовы координаты:

double xPos = (lon-lonAnchor)*( Math.toRadians( 6378137 ) )*Math.cos( latAnchor );
double yPos = (lat-latAnchor)*( Math.toRadians( 6378137 ) );

затем я использую кросс-произведение для определения площади.

Проблема в том, что результаты немного погрешны (около 1%). Есть ли что-нибудь, что я могу изучить, чтобы улучшить это?

Спасибо.

Answer 1

Я проверил в интернете различные формулы (или код) многоугольников, но не нашел ни одного хорошего или простого в реализации.

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

Несколько важных моментов

1. Массив, введенный в эту функцию, будет иметь элементы "n + 1". Последний элемент будет иметь те же значения, что и первый.
2. Я написал очень простой код на C #, чтобы ребята также могли адаптировать его на другом языке.
3. 6378137 - значение радиуса Земли в метрах.

4. Выходная площадь будет иметь единицу квадратного метра

   private static double CalculatePolygonArea(IList<MapPoint> coordinates)
   {
       double area = 0;
   
       if (coordinates.Count > 2)
       {
           for (var i = 0; i < coordinates.Count - 1; i++)
           {
               MapPoint p1 = coordinates[i];
               MapPoint p2 = coordinates[i + 1];
               area += ConvertToRadian(p2.Longitude - p1.Longitude) * (2 + Math.Sin(ConvertToRadian(p1.Latitude)) + Math.Sin(ConvertToRadian(p2.Latitude)));
           }
   
           area = area * 6378137 * 6378137 / 2;
       }
   
       return Math.Abs(area);
   }
   
   private static double ConvertToRadian(double input)
   {
       return input * Math.PI / 180;
   }
   

Answer 2

Я изменяю карту Google, чтобы пользователь мог вычислить площадь многоугольника, щелкая вершины.Он не давал правильных областей, пока я не убедился, что Math.cos (latAnchor) был в радианах первым

Итак:

double xPos = (lon-lonAnchor)*( Math.toRadians( 6378137 ) )*Math.cos( latAnchor );

стало:

double xPos = (lon-lonAnchor)*( 6378137*PI/180 ) )*Math.cos( latAnchor*PI/180 );

где lon, lonAnchor и latAnchor указаны в градусах.Теперь работает как шарм.

Answer 3

1% ошибка кажется немного высокой только из-за вашего приближения. Вы сравниваете с фактическими измерениями или идеальным расчетом? Помните, что в GPS также есть ошибка, которая может внести свой вклад.

Если вам нужен более точный метод для этого, есть хороший ответ на этот вопрос. Если вы идете по более быстрому пути, вы можете использовать геоид WGS84 вместо своей базовой сферы для преобразования в декартовы координаты (ECEF). Вот вики-ссылка для этого преобразования.

Answer 4

Спасибо, Рискованный Патхак!

В духе обмена, вот моя адаптация в Delphi:

interface

uses 
  System.Math; 

TMapGeoPoint = record
  Latitude: Double;
  Longitude: Double;
end;


function AreaInAcres(AGeoPoints: TList<TMapGeoPoint>): Double;

implementation

function AreaInAcres(AGeoPoints: TList<TMapGeoPoint>): Double;
var
  Area: Double;
  i: Integer;
  P1, P2: TMapGeoPoint;
begin
 Area := 0;

 // We need at least 2 points
 if (AGeoPoints.Count > 2) then
 begin
   for I := 0 to AGeoPoints.Count - 1 do
   begin
     P1 := AGeoPoints[i];
     if i < AGeoPoints.Count - 1  then
       P2 := AGeoPoints[i + 1]
     else
       P2 := AGeoPoints[0];
     Area := Area + DegToRad(P2.Longitude - P1.Longitude) * (2 + 
        Sin(DegToRad(P1.Latitude)) + Sin(DegToRad(P2.Latitude)));
    end;

    Area := Area * 6378137 * 6378137 / 2;

  end;

  Area := Abs(Area); //Area (in sq meters)

  // 1 Square Meter = 0.000247105 Acres
  result := Area * 0.000247105;
end;

Answer 5

Адаптированный фрагмент RiskyPathak для PHP

function CalculatePolygonArea($coordinates) {
    $area = 0;
    $coordinatesCount = sizeof($coordinates);
    if ($coordinatesCount > 2) {
      for ($i = 0; $i < $coordinatesCount - 1; $i++) {
        $p1 = $coordinates[$i];
        $p2 = $coordinates[$i + 1];
        $p1Longitude = $p1[0];
        $p2Longitude = $p2[0];
        $p1Latitude = $p1[1];
        $p2Latitude = $p2[1];
        $area += ConvertToRadian($p2Longitude - $p1Longitude) * (2 + sin(ConvertToRadian($p1Latitude)) + sin(ConvertToRadian($p2Latitude)));
      }
    $area = $area * 6378137 * 6378137 / 2;
    }
    return abs(round(($area));
}

function ConvertToRadian($input) {
    $output = $input * pi() / 180;
    return $output;
}

Answer 6

На основе решения Риски Патхака здесь представлено решение для SQL (Redshift) для вычисления областей для мультиполигонов GeoJSON (с предположением, что линейная строка 0 является самым внешним многоугольником)

create or replace view geo_area_area as 
with points as (
    select ga.id as key_geo_area
    , ga.name, gag.linestring
    , gag.position
    , radians(gag.longitude) as x
    , radians(gag.latitude) as y
    from geo_area ga
    join geo_area_geometry gag on (gag.key_geo_area = ga.id)
)
, polygons as (
    select key_geo_area, name, linestring, position 
    , x
    , lag(x) over (partition by key_geo_area, linestring order by position) as prev_x
    , y
    , lag(y) over (partition by key_geo_area, linestring order by position) as prev_y
    from points
)
, area_linestrings as (
    select key_geo_area, name, linestring
    , abs( sum( (x - prev_x) * (2 + sin(y) + sin(prev_y)) ) ) * 6378137 * 6378137 / 2 / 10^6 as area_km_squared
    from polygons
    where position != 0
    group by 1, 2, 3
)
select key_geo_area, name
, sum(case when linestring = 0 then area_km_squared else -area_km_squared end) as area_km_squared
from area_linestrings
group by 1, 2
;