Polygonfläche berechnen

[Michi (Gast)]

Hallo Zusammen,

Ich möchte gerne die Fläche mehrerer Polygone berechnen. Konkret handelt 
es sich um eine Weltkarte, alle physischen Kontinente und Inseln liegen 
mir als einzelne Pfade (Liste aus Punkten) vor. Das sind ca. 3000 Pfade 
die alle geschlossen sind. Nun würde ich gerne jeweils deren Fläche 
ausrechnen. (Karten-Projektion bzw. Tatsächlich Fläche sind mir egal).

Ich habe also mal die Gaußsche Trapezformel angesetzt. Bei einfachen 
Polygonen funktioniert das auch wunderbar, bei anderen/komplexeren kommt 
da nur blödsinn raus. Meine Vermutung geht dahin, dass die Trapezformel 
eine "saubere" Mittellinie braucht ... das kann ich aber nicht 
gewährleisten bei den komplexen Formen der Kontinente. Leider bin ich 
kein Mathematiker. Kann mir das jemand erklären?

Jetzt suche ich also eine besser Variante für die Berechnung der 
Polygonfläche. Hat jemand Links für mich oder sogar pseudocode / Python 
?


Danke

[Peter Lustig (Gast)]

Moin,

dann mach es so wie alle... Nimm Dreiecke dafür.
FEM macht das so, weil es schnell und einfach ist.

Viel Erfolg.

[H-G S. (haenschen)]

Hach waren das Zeiten ... als mein Objekt durch die selbst erstellte 
Höhenkarte lief!

Doch plötzlich versank sie an genau einer Stelle ein Drittel im Boden 
... seither zweifele ich ein wenig an den Berechnungen mit Dreiecken 
:-)

[Yalu X. (yalu) (Moderator)]

Ich nehme an, die Koordinaten der Eckpunkte sind gegeben. Dann ist es
ganz einfach und in zwei oder drei Codezeilen programmiert:

  https://de.wikipedia.org/wiki/Polygon#Fl.C3.A4che

Edit:

Ah, die hast du ja schon ohne Erfolg programmiert. Dann hast du entweder
etwas falsch gemacht, oder es sind so viele Punkte, dass sich
Rundungsfehler bemerkbar machen.

Hast du ein Beispiel für so ein Polygon und dem von dir berechneten
Ergebnis? Dann kann ich schauen, ob ich etwas anderes herausbekomme.

Edit 2:

Michi schrieb:
> Meine Vermutung geht dahin, dass die Trapezformel
> eine "saubere" Mittellinie braucht ... das kann ich aber nicht
> gewährleisten bei den komplexen Formen der Kontinente.

Nein, die Formel funktioniert für beliebige Polygone, wenn man sie
richtig anwendet, sogar für Polygone mit Löchern.

[Vektoranfänger (Gast)]

> Konkret handelt
> es sich um eine Weltkarte, alle physischen Kontinente und Inseln liegen
> mir als einzelne Pfade (Liste aus Punkten) vor.

Ich hätte auch gern solche Daten, am liebsten in verschiedenen 
Auflösungen: z.b. grob mal 500..1000km pro Segment (o. kl. Inseln)
Dann vllt. 10x/100x/1000x besser aufgelöst (bloss Richtangaben, bin 
flexibel. Koordinatensystem übrigens auch egal).

Wo gibt es solche Perimeterdaten frei?

[Michi (Gast)]

Yalu X. schrieb:
> Ah, die hast du ja schon ohne Erfolg programmiert. Dann hast du entweder
> etwas falsch gemacht, oder es sind so viele Punkte, dass sich
> Rundungsfehler bemerkbar machen.

Meine Implementierung wäre

areas_list = list() #neue leere Liste
for path in path_list:
 if len(path) > 0:
  a = 0.0
  for pntID, pnt in enumerate(path):
   a += (pnt.y + path[(pntID + 1)%len(path)].y) * (pnt.x - path[(pntID + 1)%len(path)].x)
   a = abs(a/2)
  areas_list.append(a)
 else:
  areas_list.append(0.0)

Wegen Beispielen komme ich gleich nochmal zurück - wie hättest du die 
Daten denn gerne? CSV? Irgendwie speziell formatiert?

[Yalu X. (yalu) (Moderator)]

Michi schrieb:
> wie hättest du die
> Daten denn gerne? CSV? Irgendwie speziell formatiert?

Egal. Irgendwie so, dass man sie auch mit den Augen lesen kann.

[123 (Gast)]

Yalu X. schrieb:
> Michi schrieb:
>> wie hättest du die
>> Daten denn gerne? CSV? Irgendwie speziell formatiert?
>
> Egal. Irgendwie so, dass man sie auch mit den Augen lesen kann.

Habs dir angehängt -- vielen Dank für deine Hilfe

[Yalu X. (yalu) (Moderator)]

Ich bekomme mit deinem Code (der korrekt aussieht)

  16987844.326731376


heraus, was durchaus plausibel ist. Bist du sicher, dass du die 
Zahlenwerte richtig einliest und auch das richtige Ergebnis ausgibst (es 
werden ja mehrere Flächen berechnet).

[T.roll (Gast)]

Vektoranfänger schrieb:
> Ich hätte auch gern solche Daten

Rohdaten gibts bei Openstreetmap (aufgeteilt downloadbar z.B. bei 
http://download.geofabrik.de/).

Damit kannst du dich beliebig austoben. Gibt genug Programme die 
Filtern, anpassen, usw.:
http://wiki.openstreetmap.org/wiki/Category:OSM_processing

[Wolfgang (Gast)]

Vektoranfänger schrieb:
> Wo gibt es solche Perimeterdaten frei?

Suche mal nach world+vector+shorelines

[Michi (Gast)]

Also ich habe meine von http://indiemapper.io/

[Michi (Gast)]

Ich korrigiere mal meinen Python Code von weiter oben. War natürlich ein 
Fehler in der Programmierung meinerseits. Es ist eine Einrückung, die 
Halbierung der Fläche darf nur am Ende der Summe gemacht werden und 
nicht bei jedem neuen Punkt.


Richtig also:

areas_list = list() #neue leere Liste
for path in path_list:
 if len(path) > 0:
  a = 0.0
  for pntID, pnt in enumerate(path):
   a += (pnt.y + path[(pntID + 1)%len(path)].y) * (pnt.x - path[(pntID + 1)%len(path)].x)
  a = abs(a/2)
  areas_list.append(a)
 else:
  areas_list.append(0.0)

[Can-Tak (Gast)]

Michi schrieb:
> Es ist eine Einrückung

Python ... man hasst es oder verabscheut es.

[Yalu X. (yalu) (Moderator)]

Yalu X. schrieb:
> Ich bekomme mit deinem Code (der korrekt aussieht)
>
>   16987844.326731376

Michi schrieb:
> Ich korrigiere mal meinen Python Code von weiter oben. War natürlich ein
> Fehler in der Programmierung meinerseits. Es ist eine Einrückung, ...

Upps, das hätte ich eigentlich beim Austesten bemerken müssen.

Dass ich den Fehler übersehen habe, liegt daran, dass ich beim Kopieren
deines Codes unbewusst die Einrückungstiefe an meine Gewohnheiten (2
Leerzeichen pro Ebene) angepasst und dabei ich intuitiv die falsche
Einrückung der betreffenden Zeile korrigiert habe.