Pobieranie danych z pliku OSM w Dynamo/Python

Serwis mapowy OpenStreetMap umożliwia eksportowanie danych wektorowych pożądanego obszaru do pliku XML. W bibliotece Dynamo istnieją gotowe zestawy węzłów umożliwiające ekstrakcję tych danych i zamianę na geometrię Dynamo, jak Elk lub DynaMaps. Czasem jednak zachodzi potrzeba przetworzenia tych danych w sposób, którego nie przewidzieli autorzy gotowych rozwiązań. Poniżej przedstawię prostą metodę na pobranie danych z pliku OSM przy pomocy Pythona.

Do obsługi pliku XML użyję standardowej biblioteki xml.etree.ElementTree. Przykładowo pokażę jak wydobyć obrysy budynków wraz z liczbą kondygnacji. Co prawda OSM podaje niekiedy dokładną wysokość budynku, jednak częściej zdarzało mi się znajdować właśnie informację o piętrach.

Typowa sekcja opisująca budynek w pliku OSM wygląda tak:

 <way id="28095147" visible="true" version="4" changeset="49823984" timestamp="2017-06-25T22:16:38Z" user="kocio" uid="52087">
  <nd ref="308595334"/>
  <nd ref="308595359"/>
  <nd ref="308595338"/>
  <nd ref="308595345"/>
  <nd ref="4936195428"/>
  <nd ref="308595331"/>
  <nd ref="308595334"/>
  <tag k="addr:city" v="Warszawa"/>
  <tag k="addr:housenumber" v="8"/>
  <tag k="addr:postcode" v="00-020"/>
  <tag k="addr:street" v="Chmielna"/>
  <tag k="building" v="yes"/>
  <tag k="building:levels" v="3"/>
 </way>

Sekcje obiektów, takich jak budynki, opisane są tagiem <way>. Wewnątrz takiej sekcji znajdują się m.in. elementy <tag>, zawierające dwa atrybuty „k” i „v” (key i value). Jeśli dany obiekt jest budynkiem to znajdziemy w nim<tag> zawierający atrybut „k” o wartości „building”. W ten sposób możemy przeszukać plik XML pod kątem budynków.

Innym istotnym elementem sekcji są pola <nd>. Zawierają one jeden tylko atrybut, będący referencją do innych obiektów w pliku. Pola te, węzły, opisują narożniki obrysu budynku.

I tak, w tymże samym pliku znajdziemy sekcję <node> o odpowiedniej wartości „id” dla każdej z powyższych referencji. Przykładowo:

<node id="308595334" visible="true" version="2" changeset="51689947" timestamp="2017-09-03T10:59:36Z" user="tkedt" uid="2284681" lat="52.2331940" lon="21.0179163"/>

Najważniejszymi dla nas atrybutami węzła są „lat” i „lon” opisujące współrzędne geograficzne punktu. Aby były one użyteczne w Dynamo musimy jest przekonwertować na układ współrzędnych x,y – czego nie będę szczegółowo tu opisywał. Kompletny kod na końcu posta zawiera funkcję konwertującą.

Ostatnim elementem budynku, który będzie nas w naszym przykładzie interesował jest pole <tag> o wartości atrybutu „k” równej „building:levels”. Określa ona oczywiście liczbę kondygnacji budynku – w atrybucie „v”.

Jako, że pobierając dane dla budynków będziemy potrzebowali informacji o węzłach, pierwszym krokiem w skrypcie będzie utworzenie listy (a w zasadzie słownika) wszystkich węzłów z pliku. W ten sposób będziemy mogli łatwo odszukać współrzędne interesującego nas wierzchołka.

root = xml.parse(file_path).getroot()
children = root.getchildren()

nodes = {}

for c in children:
	if c.tag == "node":
		nodes[c.attrib["id"]] = [c.attrib["lat"], c.attrib["lon"]]

Powyższy kod tworzy słownik „nodes”, w którym kluczem jest „id” węzła, a wartością lista zawierająca współrzędne geograficzne.

W kodzie całej funkcji znajdziemy też fragment przetwarzający tag o nazwie <bounds>. Dotyczy on konwersji współrzędnych i jak wyżej znów nie będę się w niego wgłębiał.

Następnie możemy już przejść do właściwej części kodu, czyli ekstrakcji budynków.

output = []

for c in children:
	if c.tag == "way":
		outline_nodes = []
		is_building = False
		levels = 0
		for t in c:
			if t.tag == "tag":
				if t.attrib["k"] == "building":
					is_building = True
				if t.attrib["k"] == "building:levels":
					levels = int(t.attrib["v"])
			if t.tag == "nd":
				outline_nodes.append(t.attrib["ref"])
		if is_building:
			outline = []
			for on in outline_nodes:
				lat = float(nodes[on][0])
				lon = float(nodes[on][1])
				xy = ll_to_xy(lat,lon)
				p = Point.ByCoordinates(xy[0],xy[1])
				outline.append(p)
			output.append([outline,levels])

Powyższy kod dla każdego znalezionego taga <way> sprawdza czy zawiera on element tag, mówiący o tym, że dana sekcja opisuje budynek. Wyszukuje również elementy opisujące liczbę kondygnacji lub odniesienia do węzłów.

Jeśli dana sekcja okaże się budynkiem, to w drugiej części pętli tworzona jest struktura danych zawierająca obrys budynku (utworzony z węzłów) oraz liczba kondygnacji. Następnie jest ona przekazywana do listy wyników naszej funkcji.

Poniżej prezentuję cały kod, który na wejściu przyjmuje ściężkę do pliku OSM. Zwracana do Dynamo jest lista budynków, z uzyskanymi przez nas danymi o nich.

import sys
sys.path.append(r'C:\Program Files (x86)\IronPython 2.7\DLLs')
sys.path.append(r'C:\Program Files (x86)\IronPython 2.7\Lib')
import xml.etree.ElementTree as xml
import math

import clr
clr.AddReference('ProtoGeometry')
from Autodesk.DesignScript.Geometry import *

r_in_meters = 6371000
lat_center = 0
cos_lat = 0
base_lat = 0
base_lon = 0

def ll_to_xy(lat,lon):
	return [r_in_meters * math.radians(lon-base_lon) * cos_lat,r_in_meters * math.radians(lat-base_lat)]
	

#The inputs to this node will be stored as a list in the IN variables.
file_path = IN[0]

root = xml.parse(file_path).getroot()
children = root.getchildren()

output = []
nodes = {}

for c in children:
	if c.tag == "node":
		nodes[c.attrib["id"]] = [c.attrib["lat"], c.attrib["lon"]]
	if c.tag == "bounds":
		minLat = float(c.attrib["minlat"])
		maxLat = float(c.attrib["maxlat"])
		minLon = float(c.attrib["minlon"])
		maxLon = float(c.attrib["maxlon"])
		
		base_lat = minLat + (maxLat-minLat)/2
		base_lon = minLon + (maxLon-minLon)/2
		cos_lat = math.cos(math.radians(base_lat))

for c in children:
	if c.tag == "way":
		outline_nodes = []
		is_building = False
		levels = 0
		for t in c:
			if t.tag == "tag":
				if t.attrib["k"] == "building":
					is_building = True
				if t.attrib["k"] == "building:levels":
					levels = int(t.attrib["v"])
			if t.tag == "nd":
				outline_nodes.append(t.attrib["ref"])
		if is_building:
			outline = []
			for on in outline_nodes:
				lat = float(nodes[on][0])
				lon = float(nodes[on][1])
				xy = ll_to_xy(lat,lon)
				p = Point.ByCoordinates(xy[0],xy[1])
				outline.append(p)
			output.append([outline,levels])
		

#Assign your output to the OUT variable.
OUT = output
Posts created 32

Dodaj komentarz

Twój adres email nie zostanie opublikowany.

Related Posts

Begin typing your search term above and press enter to search. Press ESC to cancel.

Back To Top