User:Baiuvarius/Topo 50/Flächenberechnung
Flächenberechnung von geschlossenen Wegen und Multipolygonen mit Python in Maperitive
Grundlagen
Formel: Gaußsche Trapezformel (Wikipedia)
Die Form der Erde wird als Kugel mit einem Radius von 6371 km angenähert. Der berechnete Flächeninhalt in Quadratkilometern wird in ein neues Tag namens "area" geschrieben.
Vorbereitung
maperipy, math und OSM-Daten in Python laden
from maperipy import *
from maperipy.osm import *
import math
osm = None
osm_layer= None
layer_index=0
for layer in Map.layers:
layer_index += 1
if layer.layer_type == "OsmLayer":
osm = layer.osm
break
if osm == None:
raise AssertionError("There are no OSM map souces.")
Geschlossene Wege
Die Berechnung des Inhalts der Fläche, die von einem geschlossenen Weg umschlossen wird, ist einfach: Man arbeitet die Liste der Knoten der Reihe nach ab.
for way in osm.find_ways(lambda x : x.is_closed and (x.has_tag("landuse") or x.has_tag("natural") or x.has_tag("place")) and x.has_tag("name")):
a=0.0
n=way.nodes_count
for i in range(n):
kminus1=osm.node(way.nodes[(i-1)%n])
k=osm.node(way.nodes[i])
kplus1=osm.node(way.nodes[(i+1)%n])
a+=math.radians(k.location.x)*math.cos(math.radians(k.location.y))*math.radians(kplus1.location.y-kminus1.location.y)
area=abs(a)/2*6371*6371
way.set_tag("area",str(area))
print way.get_tag("name"),area,"km²"
Multipolygone
Die Berechnung des Flächeninhalts von Multipolygonen ist etwas komplizierter: Die einzelnen Wege müssen erst geschlossenen Ringen zugeordnet werden. Dazu beginnt man mit einem beliebigen Weg, fügt ihn zu dem aktuellen Ring hinzu und sucht dann einen Weg, dessen erster oder letzter Knoten mit dem letzten Knoten des Ringes übereinstimmt. Das wiederholt man, bis der Ring geschlossen ist, und berechnet dann den Flächeninhalt des Rings. Man bildet so lange neue Ringe, bis alle Wege einem Ring zugeordnet sind. Schließlich kann man die Flächeninhalte der einzelnen Ringe addieren oder subtrahieren, je nachdem, ob es sich um innere oder äußere Begrenzungen handelt.
for relation in osm.find_relations(lambda x: x.has_tag("type","multipolygon") and (x.has_tag("landuse") or x.has_tag("natural") or x.has_tag("place")) and x.has_tag("name")):
incomplete=False
a=0.0
#assemble all member ways in "unassigned" lists
outer_ways=[]
inner_ways=[]
#list of nodes in current ring
ring=[]
ring_area=0.0
for element in relation.members:
if element.ref_type!=OsmReferenceType.WAY:
continue
#set dummy tag because way.nodes seems not to work on untagged ways
try:
osm.way(element.ref_id).set_tag("key","value")
except:
incomplete=True
break
if element.role=="inner":
inner_ways.append(element.ref_id)
else:
outer_ways.append(element.ref_id)
if incomplete:
print relation.get_tag("name"),"(id",relation.id, ") incomplete"
continue
while len(outer_ways)>0:
#assign one unassigned way to the current ring
way=osm.way(outer_ways[0])
for i in range(way.nodes_count):
ring.append(way.nodes[i])
del outer_ways[0]
#check if ring is closed
while ring[-1]!=ring[0]:
#look for way that starts or ends with end node of current ring
for way_id in outer_ways:
if osm.way(way_id).nodes[0]==ring[-1]:
way=osm.way(way_id)
for i in range(way.nodes_count):
ring.append(way.nodes[i])
outer_ways.remove(way_id)
break
elif osm.way(way_id).nodes[-1]==ring[-1]:
way=osm.way(way_id)
for i in range(way.nodes_count-1,-1,-1):
ring.append(way.nodes[i])
outer_ways.remove(way_id)
break
break
n=len(ring)
for i in range(n):
kminus1=osm.node(ring[(i-1)%n])
k=osm.node(ring[i])
kplus1=osm.node(ring[(i+1)%n])
ring_area+=math.radians(k.location.x)*math.cos(math.radians(k.location.y))*math.radians(kplus1.location.y-kminus1.location.y)
a+=abs(ring_area)
#reinitialise ring
ring=[]
ring_area=0.0
while len(inner_ways)>0:
#assign one unassigned way to the current ring
way=osm.way(inner_ways[0])
for i in range(way.nodes_count):
ring.append(way.nodes[i])
del inner_ways[0]
#check if ring is closed
while ring[-1]!=ring[0]:
#look for way that starts or ends with end node of current ring
for way_id in inner_ways:
if osm.way(way_id).nodes[0]==ring[-1]:
way=osm.way(way_id)
for i in range(way.nodes_count):
ring.append(way.nodes[i])
inner_ways.remove(way_id)
break
elif osm.way(way_id).nodes[-1]==ring[-1]:
way=osm.way(way_id)
for i in range(way.nodes_count-1,-1,-1):
ring.append(way.nodes[i])
inner_ways.remove(way_id)
break
break
n=len(ring)
for i in range(n):
kminus1=osm.node(ring[(i-1)%n])
k=osm.node(ring[i])
kplus1=osm.node(ring[(i+1)%n])
ring_area+=math.radians(k.location.x)*math.cos(math.radians(k.location.y))*math.radians(kplus1.location.y-kminus1.location.y)
a-=abs(ring_area)
#reinitialise ring
ring=[]
ring_area=0.0
area=a/2*6371*6371
relation.set_tag("area",str(area))
print relation.get_tag("name"),area,"km²"