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²"