Python: Flächeninhalt berechnen

!i! · 2013-07-09 10:11:23

Hi,

ich möchte gerne von einer durch WGS84 Koordinaten (genauer: Boundingboxen) den Flächeninhalt automatisiert per Skript berechnen.
Ich habe recherchiert, aber mir ist nur das hier begegnet:
http://stackoverflow.com/questions/4681 … ing-python
Allerdings ist mir dabei das Zusammenspiel nicht klar, insbesondere wie die Koordinaten für PROJ zustande kommen

Hat jemand von euch sowas vielleicht schon mal gemacht?
Irgendwie habe ich gerade ein Brett vorm Kopf

wambacher · 2013-07-09 10:26:09

!i! wrote:

Irgendwie habe ich gerade ein Brett vorm Kopf

mag ich mich nicht zu äußern

Die Frage ist - zumindest für mich - ein wenig unklar formuliert.
Willst du die Flächen exakt an der BBOX abschneiden und nur die Größe der sichtbaren Bereiche berechnen? Oder aber die totalen Größen auch über den Rand hinaus?

Und wo kommen denn deine Daten her? Klar, aus Osm, aber wie?.

Gruss
walter

p.s. ich denke da an PostGis, aber es hängt halt von deinem Umfeld ab. Python sollte da keine Probleme machen.

Last edited by wambacher (2013-07-09 10:27:02)

!i! · 2013-07-09 10:41:21

Wie gesagt, ich will aus min/max lat/lon welche die Boundingbox beschreiben den Flächeninhalt bestimmen. Ich denke das ist, was du mit "sichtbaren" Bereich meinst.

Meinem Verständnis nach ist das ein Rechteck auf der projezierten ausgebreiteten Fläche. Da dort der Elipsoid mit eine Rolle spielt, müsste man das umprojezieren, damit man die Lage der BBOX mit einbezieht und eine längentreue Abbildung hat. Und dann sollte man eigentlich doch mit a*b weiterkommen?

Ich nutze ja bereits pyPROJ für die Distanz und würde gerne auch dabei bleiben:

def getDistance(p1lat,p1lon,p2lat,p2lon):
	g = Geod(ellps='WGS84', units="m")
	az12,az21,dist = g.inv(p1lon,p1lat,p2lon,p2lat) #fuer grosse und kleine Distanzen ok
	return (dist/1000) #in km

maxbe · 2013-07-09 10:43:23

!i! wrote:

Allerdings ist mir dabei das Zusammenspiel nicht klar, insbesondere wie die Koordinaten für PROJ zustande kommen

Das Polygon ist in Längen- und Breitengraden angegeben. Das ist recht unpraktisch, wenn man Flächen berechnen will, weil die Abstände der Längengrade nach Norden abnehmen. Also verwenden die pyproj, um das Polygon erstmal in eine Projektion umzurechnen, die echte Meter und Quadratmeter verwendet.

Diese Wahl der Projektion (+proj=aea +lat_1=37.0 +lat_2=41.0 +lat_0=39.0 +lon_0=-106.55) ist speziell für die eine Gegend in den USA gut ("aea" ist eine Albers-Kegelprojektion, die nehmen 106° West als Mittelmeridian zwischen 37° und 41° Nord).
Wenn Du hier in der Gegend was machst, solltest Du lieber irgendwas nehmen, was unsere Gegend korrekter abbildet (aea mit 10 Grad Ost und 47 bis 51 Nord z.B., oder UTM oder GK).

Grüße, Max

maxbe · 2013-07-09 10:58:12

!i! wrote:

Da dort der Elipsoid mit eine Rolle spielt, müsste man das umprojezieren, damit man die Lage der BBOX mit einbezieht und eine längentreue Abbildung hat. Und dann sollte man eigentlich doch mit a*b weiterkommen?

Wenn bei g.inv() richtige Meter rauskommen, reicht a*b, solange Du in der Grössenordnung eines Stadtplans bleibst und nicht auf 1/1000 genau rechnen willst. Bei grösseren Flächen wird es immer ungenauer. Deine BBox mit Koordinaten in Längen- und Breitengrad ist in der richtigen kugeligen Welt ja kein Rechteck, sondern eher sowas wie ein Trapez (wenn ein Pol dabei ist, sogar ein Dreieck).

Last edited by maxbe (2013-07-09 10:59:08)

seichter · 2013-07-09 11:42:18

Keine ebene Projektion kann eine gekrümmte Oberfläche fehlerfrei abbilden.
Bei sehr großen Flächen müsste man das Kugelsegment aus dem Raumwinkel und dem Radius berechnen. Beim Ellipsoid oder gar Geoid artet das aber in Arbeit aus.
In der Praxis (kleinere Flächen bis Bundesland, nicht in Polnähe) kommt man aber mit allen Projektionen zurecht, wenn man nur die Distanzen parallel zu den Breitengraden korrigiert. Damit wird das (eigentlich sphärische) Trapez in ein Rechteck verformt. Gauß-Krüger und UTM können die Koordinaten deshalb in Meter angeben.

!i! · 2013-07-09 12:32:06

Danke schon mal für eure Hinweise das bestätigt mich darin, dass ich da wohl doch nicht so viele Bretter vor dem Kopf hatte...

Da ich aber weltweite Changeset BBoxen analysiere (ja, die meisten sind sehr klein, was wenig verwunderlich ist), müsste ich es irgendwie hinkriegen, automatisch die passende Projektion zu wählen. Hat da jemand eine Idee wie man da rangehen müsste? Es kommt dabei sicherlich erstmal nicht auf den Meter genau an. Eine andere Alternative wäre sonst einen lokalen Ausschnitt zu nehmen, vielleich ganz Deutschland, aber da leidet ja auch immer die Aussagekraft drunter

maxbe · 2013-07-09 13:03:43

!i! wrote:

ja, die meisten sind sehr klein, was wenig verwunderlich ist), müsste ich es irgendwie hinkriegen, automatisch die passende Projektion zu wählen. Hat da jemand eine Idee wie man da rangehen müsste?

Wenn die Boxen ein paar Kilometer gross sind, nimm spherical mercator (epsg 3857 oder 900913) und rechne a*b*cos(breitengrad)*cos(breitengrad).

Alternativ zu Fuß: (Breite_oben - Breite_unten)/360*Erdumfang * (Länge_rechts - Länge_links)/360*Erdumfang * cos((Breite_oben+Breite_unten)/2)

Ist global, einfach, ganz ohne Umrechnerei mit Ellipsoiden und passt genau genug für eine Changeset-Bbox.

Bei riesigen Boxen hast Du immer ein Problem. Auch wenn Du "aea" nehmen würdest (was theoretisch global anwendbar wäre, weil das ist überall flächentreu, nur am dicken Ende des Kegels recht stark in die Breite verzerrt), müsstest Du erstmal sicherstellen, dass eine Box dort richtig auf ein Stück Kegelmantel projiziert wird, und nicht nur die 4 Ecken mit Strichen verbunden. Damit haben GISe oft Probleme, die hangeln sich oft nur von Punkt zu Punkt.

Grüße, Max

seichter · 2013-07-09 13:08:17

Oh je. Die meisten Changeset-Boxen dürften zwar klein sein, aber ein paar Riesenexemplare dürften schon dabei sein, wenn jemand z.B. zwei weit auseinanderliegende Gegenden bearbeitet und gemeinsam hochgeladen hat.

Jetzt kommt es darauf an: Wenn solche Edits nicht berücksichtigt werden sollen, weil sie das Ergebnis (z.B. bearbeitete Fläche) massiv verfälschen, hat man das Problem umschifft.
Wenn diese Changesets auch gezählt werden sollen, wird es unschön. Es gibt zwar flächentreue Projektionen, die verbiegen aber über große Distanzen gerade Linien fürchterlich. Da ist es dann nichts mit Breite mal Höhe.

Oranger Assistent · 2013-07-09 14:49:47

Tach.

!i! wrote:

ich möchte gerne von einer durch WGS84 Koordinaten (genauer: Boundingboxen) den Flächeninhalt automatisiert per Skript berechnen.

Gegeben: MinLon, MaxLon, MinLat, MaxLat (in Grad).

Die vier Werte definieren vier Punkte auf einer Kugel. Wenn wir diese vier Punkte über Meridianbögen und Parallelkreisbögen verbinden, erhalten wir ein Viereck. Gesucht ist die Fläche dieses Vierecks.

Lösung auf einer Kugel mit Radius R:

λ1 = minLon / 180° * π (Bogenmaß)
λ2 = maxLon / 180° * π (Bogenmaß)

φ1 = minLat / 180° * π (Bogenmaß)
φ2 = maxLat / 180° * π (Bogenmaß)

R = 6371000m (zum Beispiel)

F = (λ2-λ1) × (sinφ2-sinφ1) × R²

Diese Lösung ist exakt. Auf dem Ellipsoid ist die Formel etwas aufwendiger.
Das wird aber für Deine Anwendung sicher nicht benötigt.

Hat jemand von euch sowas vielleicht schon mal gemacht?

Ja. Punkt 11.

Der Assistent
(immer hilfreich, immer orange)

Edit: Typo

Last edited by Oranger Assistent (2013-07-09 14:56:07)

!i! · 2013-07-09 19:55:09

Danke für die vielen Antworten. Ich werde mal eine Stichprobe nehmen und schauen ob mir die Näherung genügt.
Prinzipiell zeigen die Daten nur relativ wenige extreme Ausreißer, von daher denke ich damit schon typische Mapping-Muster belegen zu können (ganz kleine Attribut-Fixe, größere Areale bingen, Botläufe, ...). Ist erstmal nur ein Zwischenprodukt, durch das ich einige Parameter sicher ableiten möchte.

fx99 · 2013-07-09 20:52:01

Ich weiß zwar nicht, warum Du die Flächen brauchst, aber rechen doch einfach in Quadrat-Grad. Das ist das einfachste und für eine ORdnungsmetrik sicherlich ausreichend.

seichter · 2013-07-09 21:34:51

fx99 wrote:

Ich weiß zwar nicht, warum Du die Flächen brauchst, aber rechen doch einfach in Quadrat-Grad. Das ist das einfachste und für eine ORdnungsmetrik sicherlich ausreichend.

Ist im Prinzip dasselbe wie die Formel von "Oranger Assistent", nur Winkelrechteck mit Radiusquadrat multipliziert und Breitenausdehnung korrigiert.
(Flächenelement in Kugelkoordinaten ist R*R*dλ*sinφ*dφ, für geogr. Koordinaten muss φ um 90° verschoben werden. Ohne Breitenkorrektur werden Flächen in höheren Breiten zu stark bewertet.)

Announcement

#1 2013-07-09 10:11:23

Python: Flächeninhalt berechnen

#2 2013-07-09 10:26:09

Re: Python: Flächeninhalt berechnen

#3 2013-07-09 10:41:21

Re: Python: Flächeninhalt berechnen

#4 2013-07-09 10:43:23

Re: Python: Flächeninhalt berechnen

#5 2013-07-09 10:58:12

Re: Python: Flächeninhalt berechnen

#6 2013-07-09 11:42:18

Re: Python: Flächeninhalt berechnen

#7 2013-07-09 12:32:06

Re: Python: Flächeninhalt berechnen

#8 2013-07-09 13:03:43

Re: Python: Flächeninhalt berechnen

#9 2013-07-09 13:08:17

Re: Python: Flächeninhalt berechnen

#10 2013-07-09 14:49:47

Re: Python: Flächeninhalt berechnen

#11 2013-07-09 19:55:09

Re: Python: Flächeninhalt berechnen

#12 2013-07-09 20:52:01

Re: Python: Flächeninhalt berechnen

#13 2013-07-09 21:34:51

Re: Python: Flächeninhalt berechnen

Board footer