calcul de la superficie Polygon en utilisant Latitude et Longitude produite à partir de l'espace cartésien et un fichier mondial

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

Question

donné une série de paires de coordonnées GPS, il faut que pour calculer l'aire d'un polygone (n-gon). Ceci est relativement faible (pas plus de 50 000 pieds carrés). Les géocodes sont créés en appliquant une transformation de affines avec les données d'un fichier mondial.

J'ai essayé d'utiliser une approche en deux étapes en faisant convertir les géocodage en coordonnées cartésiennes: 

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

j'utilise un calcul de produit croisé pour déterminer la zone.

Le problème est que les résultats sont un peu hors de la précision (environ 1%). Est-ce que je peux examiner pour améliorer?

Merci.

Autres conseils

Je suis en train de modifier une carte Google de sorte qu'un utilisateur peut calculer la zone d'un polygone en cliquant sur les sommets. Il ne donnait pas correcte zones jusqu'à ce que je fait que le Math.cos (latAnchor) était en radians premier 

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

est devenu: 

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

où lon, lonAnchor et latAnchor sont en degrés. Fonctionne comme un charme maintenant.

J'ai vérifié sur Internet pour différentes formules de la zone polygone (ou code) mais n'a pas trouvé une bonne ou à mettre en œuvre facile.

Maintenant, je l'ai écrit l'extrait de code à la zone de calcul d'un polygone dessiné sur la surface de la terre. Le polygone peut avoir n sommets avec chaque sommet a ayant sa propre latitude longitude.

Quelques points importants

  1. L'entrée de tableau pour cette fonction sera « n + 1 » éléments. Le dernier élément aura les mêmes valeurs que celle du premier.
  2. J'ai écrit très basique code C #, de sorte que les gars peuvent également adapter dans d'autres langues.
  3. 6378137 est la valeur du rayon de la terre, en mètres.

  4. La zone de sortie aura unité de mètres carrés 

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

Sur la base de la solution par Risky Pathak est la solution ici pour SQL (Redshift) aux zones calculate pour GeoJSON multipolygones (avec l'hypothèse que linestring 0 est le polygone le plus extérieur) 

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
;

L'extrait de RiskyPathak Adapté à 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;
}

Merci Risky Pathak!

Dans l'esprit de partage, voici mon adaptation à Delphes: 

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;
Licencié sous: CC-BY-SA avec attribution
Non affilié à StackOverflow
scroll top