Bauforum-Logo

Offenes Forum Bauingenieurwesen

log in | registrieren

zurück zum Forum
  Mix-Ansicht

3D-Modell jedes beliebigen Grundstücks in NRW (Software)

verfasst von Martin Vogel Homepage E-Mail, Dortmund / Bochum, 16.03.2017, 15:14 Uhr

Die Qualität der seit Januar 2017 frei verfügbaren Geländehöhendaten für NRW kann einem den Atem stocken lassen. Hier werfen wir einen Blick auf das Gelände der in einen ehemaligen Steinbruch an der Volme hineingebauten Stadthalle Hagen. Visualisiert wurden die aufbereiteten Geländedaten mit dem kostenlosen CAD-Programm FreeCAD. Die horizontale Auflösung beträgt einen Meter, die vertikale Auflösung der Höhendaten beträgt sogar nur einen Zentimeter:

[image]
Datenquelle: https://www.opengeodata.nrw.de/produkte/geobasis/dgm/dgm1/ ; Land NRW (2017) ; Datenlizenz Deutschland - Namensnennung - Version 2.0 (http://www.govdata.de/dl-de/by-2-0)

Auf den Servern des Opengeodata-Projektes gibt es sogar noch höher aufgelöste Daten, diese sind aber etwas unhandlicher zu verarbeiten, da sie nicht in einem festen Raster, sondern als ungeordnete Punktewolke mit etwa 4 Punkten pro Quadratmeter vorliegen.

Es ist nun für jedermann ganz einfach, sich ein 3D-Modell seines Grundstücks oder eines beliebigen anderen Geländeausschnittes selber zu erstellen. Sämtliche benötigte Software dazu ist kostenlos verfügbar. Bis zur 3D-druckfähigen Geländemodelldatei sind es nur wenige Schritte:

1. Von der Webseite https://www.opengeodata.nrw.de/produkte/geobasis/dgm/dgm1/ die Dateien für die gewünschten Ortschaften herunterladen.
Die ZIP-Dateien sind im Schnitt ein paar hundert Megabyte groß.

2. Diese ZIP-Dateien in einen beliebigen Ordner entpacken.
Es werden jeweils dutzende XYZ-Dateien angelegt, jede ist 124 MB groß und enthält 4 Millionen Höhenwerte einer 2000x2000 m² großen Geländekachel.

Wer in die XYZ-Dateien hineinschaut, stellt fest, dass die Koordinaten in einem für Nichtvermesser recht ungewöhnlichen Format vorliegen. Dieses Datenformat mit dem etwas sperrigen Namen „EPSG:25832 ETRS89/UTM 32.Zone“ lässt sich aber in gewöhnliche GPS-Koordinaten mit Längen- und Breitengrad umrechnen. Ebenso können bekannte Geokoordinaten in das EPSG-Format konvertiert werden.

3. Zur Umrechnung verwenden wir ein selbstgeschriebenes Python-3-Programm, welches das Python-Modul "utm" von Tobias Bieniek verwendet. Laden Sie die folgende ZIP-Datei herunter, entpacken Sie sie und starten Sie das Programm "Gelaendemodell.py".
[image]ZIP-Archiv_8QQNV90P1.zip

Falls Sie Python 3 noch nicht installiert haben, finden Sie es auf https://www.python.org/downloads/ – im Moment ist Version 3.6.0 aktuell.

4. Das Programm führt Sie in einem Dialog schrittweise zum gewünschten Ergebnis. Sie benötigen lediglich noch die Koordinaten von zwei Punkten, mit denen Sie ein Rechteck definieren, für welches unser Python-Programm das Geländemodell erzeugen soll. Diese Koordinaten können Sie direkt aus Google Maps oder einem anderen Kartendienst kopieren oder abschreiben. Für die Stadthalle Hagen oben waren das zum Beispiel die Punkte „51.350383, 7.482606“ und „51.356366, 7.491436“. Welche Eckpunkte des Rechtecks Sie aus der Karte abgreifen, ist egal, solange diese nur diagonal gegenüberliegen und sich innerhalb von Nordrheinwestfalen befinden.

5. Nach Angabe der gewünschten Auflösung legt das Programm drei Dateien neu an: eine XYZ-Datei enthält einfach nur eine Liste aller gefundenen Höhendaten des ausgewählten Bereichs, eine SCR-Datei enthält Befehle für die CAD-Programme AutoCAD und BricsCAD zur Erstellung eines Volumenkörpers, über den Sie sehr simpel Massenberechnungen zu Aushub und Auftrag anstellen können und schließlich eine Stereolithographiedatei im STL-Format, aus der Sie mit einem 3D-Drucker ein greifbares Modell erschaffen können. Die Schrittweite für die Auflösung sollte anfangs auf einen recht hohen Wert eingestellt werden, um die erzeugten Dateien klein und damit schnell bearbeitbar zu halten. Ein 10-Meter-Raster oder sogar ein 100-Meter-Raster zu verwenden, ist meistens gar keine schlechte Idee. Mehr als vielleicht 1000 Punkte pro Grundstück werden oft gar nicht benötigt.

Zum Ansehen und Weiterbearbeiten der STL-Datei empfiehlt sich das kostenlose Programm FreeCAD, in dessen erstaunlich mächtigen Funktionsumfang man sich allerdings erst einarbeiten muss. Kurzanleitung: STL-Datei öffnen, Inhalt mit dem Lupen-Icon bildfüllend darstellen, Verschieben der Ansicht mit der mittleren Maustaste, Drehen mit zusätzlichem Drücken der linken Maustaste (Bedienungsanleitung).

Viel Spaß damit![image]

#! /usr/bin/env python3

# Python-Skript zur Erstellung von Höhenmodellen

# Ein durch zwei Koordinaten im WGS84-Dezimalgradformat (Google Maps)
# aufgespanntes Rechteck gibt die Eckpunkte einer Geländeoberfläche an,
# für die jeweils ein Höhenmodell im STL-Format für 3D-Drucker, eine
# auf die Auswahl reduzierte XYZ-Datei und ein AutoCAD-Befehlsskript (SCR)
# erzeugt werden. Das CAD-Skript enthält die Befehle zur Erzeugung eines
# Prismenfeldes aus 3D-Körpern, sodass Volumenbefehle angewendet werden
# können.

# Das Programm wurde für die öffentlichen Höhendaten für digitale
# Geländemodelle in NRW geschrieben, lässt sich aber mit geringen Anpassungen
# auch für andere DGM-Daten nutzen.

# Datenquelle:
# http://www.bezreg-koeln.nrw.de/
#        brk_internet/geobasis/hoehenmodelle/gelaendemodelle/
# bzw. https://www.opengeodata.nrw.de/produkte/

# Autor: Martin Vogel, Hochschule Bochum, 2017
# Lizenz: Namensnennung - Weitergabe unter gleichen Bedingungen 3.0 Deutschland
# (CC BY-SA 3.0 DE) https://creativecommons.org/licenses/by-sa/3.0/de/

import os

import utm
# utm ist ein Modul zur Koordinatenumrechnung
# von Tobias Bieniek <Tobias.Bieniek@gmx.de> (2012)
# https://pypi.python.org/pypi/utm

from tkinter import Tk
from tkinter.filedialog import askdirectory

print("Bitte wählen Sie den Ordner mit den ausgepackten xyz-Dateien aus:")
# Das Programm hat kein GUI-Fenster. Noch nicht.
Tk().withdraw() 
ordner = askdirectory(title="Bitte den gewünschten Ordner doppelklicken")
print("\nGewählter Ordner:\n",ordner)

print("\nGeben Sie einen Namen für die erzeugten Dateien an!")
print("Existierende Dateien mit diesem Namen und den Endungen .scr,")
print(".xyz und .stl werden überschrieben.\n")
name = input("Dateiname ohne Endung: ")

print("\nGeben Sie zwei gegenüberliegende Eckpunkte des zu modellierenden")
print("Rechtecks an! Das Eingabeformat ist (Breite, Länge) in Dezimalgrad,")
print("also beispielsweise 51.335757,7.479087 – Sie können die Koordinaten")
print("direkt aus Google Maps oder Ihrem GPS-Gerät übernehmen.\n")

lat1, lon1 = eval(input("Erstes Eckpunktkoordinatenpaar: "))
lat2, lon2 = eval(input("Gegenüberliegendes Koordinatenpaar: "))

# Die Eingabe oben war beliebig, wir brauchen aber gleich die Südwestecke
# unten links und die Nordostecke oben rechts:
ul_lat = min(lat1, lat2)
ul_lon = min(lon1, lon2)
or_lat = max(lat1, lat2)
or_lon = max(lon1, lon2)

# Ohne Tobias Bieniek wäre ich jetzt aufgeschmissen. Die Umrechnung der
# Geokoordinaten ist doch komplizierter als man naiverweise annimmt.
# Aus geographischem Längen- und Breitengrad wird ein Ostwert, ein Nordwert,
# eine Zonennummer und ein Zonenbuchstabe. Für NRW sind die beiden letzteren
# beispielsweise 32 und U.
e, n, zn, zb = utm.from_latlon(ul_lat, ul_lon)

# Die Zonennummer ist in den Dateien, die unter opengeodata.nrw.de
# heruntergeladen werden können, dem Ostwert vorangestellt.
ul_e = int("%i%i"%(zn,e))

# Der Nordwert wird übernommen.
ul_n = int(n)

# Das gleiche noch einmal für die Nordostecke unseres Rechtecks:
e, n, zn, zb = utm.from_latlon(or_lat, or_lon)
or_e = int("%i%i"%(zn,e))
or_n = int(n)

print("\nETRS89/UTM-Koordinaten des Rechtecks:")
print("%i,%i %i,%i in Zone %i%s" % (ul_e,ul_n,or_e,or_n,zn,zb))
print("Ausdehnung Ost-West: %i m" % (or_e-ul_e))
print("Ausdehnung Nord-Süd: %i m" % (or_n-ul_n))
print("Fläche: %i m²" % ((or_e-ul_e)*(or_n-ul_n)))

print("\nDie Auflösung der Daten beträgt 1 Meter in x- und y-Richtung, was")
print("zu extremen Dateigrößen und Verarbeitungszeiten der erzeugten Modelle")
print("führen kann.\n")
try:
    kl = int(input("Geben Sie einen ganzzahligen Wert für die Auflösung ein: "))
except:
    kl = 1

xmax = or_e - (or_e-ul_e) % kl
ymax = or_n - (or_n-ul_n) % kl

print("\nAbstand der neuen Punkte: %i m" % kl)
print("Daraus ergibt sich als obere rechte Ecke des Rechtecks:",xmax,ymax)
print("Ausdehnung Ost-West: %i m" % (xmax-ul_e))
print("Ausdehnung Nord-Süd: %i m" % (ymax-ul_n))
print("Fläche: %i m²" % ((xmax-ul_e)*(ymax-ul_n)))


# Liste der Dateien im zu durchsuchenden Ordner 
auszuwerten = os.listdir(ordner)

# Alle gefundenen Höhenwerte werden zunächst in ein Dictionary
# geschrieben, aus dem sie für die einzelnen Dateien wieder
# ausgelesen werden.

D = {}

for f in auszuwerten:
    # Bei den NRW-Höhendaten kann man am Dateinamen erkennen, ob gesuchte
    # Punkte in ihnen enthalten sein könnten.
    # Der Dateiname gibt die untere linke Ecke einer 2000x2000-m²-Kachel an.
    # Beispiel für einen Dateinamen: dgm1_32368_5700_2_nw.xyz
    # Hier sind die Ostwerte 32368000 bis 32369999 und die Nordwerte
    # 5700000 bis n5701999 enthalten.
    fs = f.split("_")

    try:
        fx = float(fs[1])
        fy = float(fs[2])
    except:
        fx = fy = -1

    if ul_e-2000 <= fx*1000 <= or_e and ul_n-2000 <= fy*1000 <= or_n :
        with open(ordner+"/"+f) as dgm:
            print("Verarbeite Datei",f,"…")
            for zeile in dgm:
                # Beispiel für eine Zeile aus Bochum:
                # 32372000.00 5706000.00   61.32
                try:
                    x,y,h = zeile.split()
                    x = int(float(x))
                    y = int(float(y))
                    if ul_e <= x <= or_e and ul_n <= y <= or_n:
                        D[x,y] = h
                except:
                    print("Abbruch, falsches Format:",zeile)
                    break

print()
print(len(D),"Punkte des gesuchten Rechtecks gefunden. (%.2f%%)" %
      (len(D)/(or_e-ul_e+1)/(or_n-ul_n+1)*100))
print(len(D)//kl**2,"Punkte werden für Ausgabe verwendet.")

ausname = name+".xyz"
print("\nSchreibe xyz-Ausgabedatei:",ausname)

with open(ausname,"w") as aus:
    for x in range(ul_e,or_e+1,kl):
        for y in range(ul_n,or_n+1,kl):
            aus.write("%i %i %s\n"%(x,y,D[(x,y)]))

ausname = name+".scr"
print("\nSchreibe CAD-Skriptdatei mit Prismenfeld:",ausname)

# Eigentlich wäre es universeller, anstelle der deutschsprachigen
# Bezeichner englischsprachige Bezeichner mit vorangestelltem Unterstrich
# zu verwenden. Dummerweise kommt BricsCAD damit nicht immer zurecht.
with open(ausname,"w") as aus:
    # Koordinateneingabe hat Vorrang vor Objektfang
    aus.write("OSnapCoord 1\n")
    # Ansicht schräg von Südwesten:
    aus.write("APunkt -1,-2,2\n")
    # Drahtmodell (Geschwindigkeit!)
    aus.write("-Vis A Drahtmodell\n")
    # Hier ist der interessante Bereich
    aus.write("Zoom F 0,0 %i,%i\n" % (or_e-ul_e, or_n-ul_n))

    for x in range(ul_e,or_e+1,kl):
        for y in range(ul_n,or_n+1,kl):
            aus.write("Quader %i,%i,0 %i,%i,%s\n" % (
                x-ul_e, y-ul_n, 
                x-ul_e+kl, y-ul_n+kl, D[(x,y)]))

ausname = name+".stl"
print("\nSchreibe STL-Datei für 3D-Druck:",ausname)

# Die Flächen umhüllen einen Quader, der unten auf Höhe null aufliegt
# und oben durch die Geländeoberfläche abgeschnitten wird. In Küstennähe
# sollte man sich etwas intelligenteres einfallen lassen.
# … im Hochgebirge bei kleinen Grundstücken auch :-)

with open(ausname,"w") as aus:
    aus.write("solid "+ausname+"\n")
    
    # Geländeoberfläche
    for x in range(ul_e,xmax+1-kl,kl):
      xu = x - ul_e
      for y in range(ul_n,ymax+1-kl,kl):
        yu = y - ul_n
        h1 = D[x,y]
        h2 = D[x+kl,y]
        h3 = D[x,y+kl]
        h4 = D[x+kl,y+kl]
        aus.write("facet normal 0 0 1\n"
                  "outer loop\n"
                  "vertex %i %i %s\n"%(xu, yu, h1) +
                  "vertex %i %i %s\n"%(xu+kl, yu+kl, h4) +
                  "vertex %i %i %s\n"%(xu+kl, yu, h2) +
                  "endloop\n"
                  "endfacet\n")
        aus.write("facet normal 0 0 1\n"
                  "outer loop\n"
                  "vertex %i %i %s\n"%(xu, yu, h1) +
                  "vertex %i %i %s\n"%(xu, yu+kl, h3) +
                  "vertex %i %i %s\n"%(xu+kl, yu+kl, h4) +
                  "endloop\n"
                  "endfacet\n")
        # Was für ein Aufwand für zwei popelige Dreiecke!
        
    # Unterseite
    for x in range(ul_e,xmax+1-kl,kl):
      xu = x - ul_e
      for y in range(ul_n,ymax+1-kl,kl):
        yu = y - ul_n
        aus.write("facet normal 0 0 -1\n"
                  "outer loop\n"
                  "vertex %i %i 0\n"%(xu, yu) +
                  "vertex %i %i 0\n"%(xu+kl, yu) +
                  "vertex %i %i 0\n"%(xu+kl, yu+kl) +
                  "endloop\n"
                  "endfacet\n")
        aus.write("facet normal 0 0 -1\n"
                  "outer loop\n"
                  "vertex %i %i 0\n"%(xu, yu) +
                  "vertex %i %i 0\n"%(xu+kl, yu+kl) +
                  "vertex %i %i 0\n"%(xu, yu+kl) +
                  "endloop\n"
                  "endfacet\n")
        
    # Linke Wand
    for y in range(ul_n,ymax+1-kl,kl):
        yu = y-ul_n
        h1 = D[ul_e,y]
        h2 = D[ul_e,y+kl]
        aus.write("facet normal -1 0 0\n"
                  "outer loop\n"
                  "vertex 0 %i 0\n"%(yu) +
                  "vertex 0 %i %s\n"%(yu+kl, h2) +
                  "vertex 0 %i %s\n"%(yu, h1) +
                  "endloop\n"
                  "endfacet\n")
        aus.write("facet normal -1 0 0\n"
                  "outer loop\n"
                  "vertex 0 %i 0\n"%(yu) +
                  "vertex 0 %i 0\n"%(yu+kl) +
                  "vertex 0 %i %s\n"%(yu+kl, h2) +
                  "endloop\n"
                  "endfacet\n")
        
    # Rechte Wand
    for y in range(ul_n,ymax+1-kl,kl):
        xu = xmax-ul_e
        yu = y-ul_n
        h1 = D[xmax,y]
        h2 = D[xmax,y+kl]
        aus.write("facet normal 1 0 0\n"
                  "outer loop\n"
                  "vertex %i %i 0\n" % (xu, yu) +
                  "vertex %i %i %s\n" % (xu, yu, h1) +
                  "vertex %i %i %s\n" % (xu, yu+kl, h2) +
                  "endloop\n"
                  "endfacet\n")
        aus.write("facet normal 1 0 0\n"
                  "outer loop\n"
                  "vertex %i %i 0\n" % (xu, yu) +
                  "vertex %i %i %s\n" % (xu, yu+kl, h2) +
                  "vertex %i %i 0\n" % (xu, yu+kl) +
                  "endloop\n"
                  "endfacet\n")
        
    # Vordere Wand
    for x in range(ul_e,xmax+1-kl,kl):
        h1 = D[x,ul_n]
        h2 = D[x+kl,ul_n]
        xu = x - ul_e
        aus.write("facet normal 0 -1 0\n"
                  "outer loop\n"
                  "vertex %i 0 0\n" % (xu) +
                  "vertex %i 0 %s\n" % (xu, h1) +
                  "vertex %i 0 %s\n" % (xu+kl, h2) +
                  "endloop\n"
                  "endfacet\n")
        aus.write("facet normal 0 -1 0\n"
                  "outer loop\n"
                  "vertex %i 0 0\n" % (xu) +
                  "vertex %i 0 %s\n" % (xu+kl, h2) +
                  "vertex %i 0 0\n" % (xu+kl) +
                  "endloop\n"
                  "endfacet\n")
        
    # Hintere Wand
    for x in range(ul_e,xmax+1-kl,kl):
        xu = x - ul_e
        yu = ymax - ul_n
        h1 = D[x,ymax]
        h2 = D[x+kl,ymax]
        aus.write("facet normal 0 1 0\n"
                  "outer loop\n"
                  "vertex %i %i 0\n"%(xu, yu) +
                  "vertex %i %i %s\n"%(xu+kl, yu, h2) +
                  "vertex %i %i %s\n"%(xu, yu, h1) +
                  "endloop\n"
                  "endfacet\n")
        aus.write("facet normal 0 1 0\n"
                  "outer loop\n"
                  "vertex %i %i 0\n"%(xu, yu) +
                  "vertex %i %i 0\n"%(xu+kl, yu) +
                  "vertex %i %i %s\n"%(xu+kl, yu, h2) +
                  "endloop\n"
                  "endfacet\n")


print("\nFertig. Die Ausgabedateien können nun weiterverarbeitet werden.")

--
Dipl.-Ing. Martin Vogel
Leiter des Bauforums

Heute schon programmiert? Einführung in Python 3 (PDF)

antworten
 



gesamter Thread:

zurück zum Forum
  Mix-Ansicht
Offenes Forum Bauingenieurwesen | Kontakt | Impressum
8419 Postings in 4025 Threads, 1093 registrierte User, 35 User online (0 reg., 35 Gäste)
powered by my little forum  RSS-Feed  ^
map | new