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

Martin Vogel ⌂ @, Dortmund / Bochum, Donnerstag, 16.03.2017, 15:14 (vor 1002 Tagen)

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

[Update: Das Programm wurde aktualisiert, siehe Antwort auf diesen Beitrag!]

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]

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

Heute kostenlos: Pythonbuch (PDF)

Tags:
GIS, NRW, DGM, Geländemodell, Geodaten, Python 3

Das Einhornkänguruh im Landschaftsschutzgebiet

Martin Vogel ⌂ @, Dortmund / Bochum, Freitag, 17.03.2017, 13:33 (vor 1001 Tagen) @ Martin Vogel

Durch Zufall entdeckt: Würde man das Landschaftsschutzgebiet hinter der BO fluten, hätte der neue See die Form eines Einhornkänguruhs :-)

[image]

Das Geländemodell im Bild hat eine Auflösung von 10 Metern, daher sind die Dreiecke, aus denen die Oberfläche geformt ist, noch recht gut erkennbar.

… kommt zufällig gerade jemand günstig an 220 Meter Spundwand?

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

Heute kostenlos: Pythonbuch (PDF)

Aktualisierte Version des Geländemodellierers

Martin Vogel ⌂ @, Dortmund / Bochum, Dienstag, 27.03.2018, 13:24 (vor 626 Tagen) @ Martin Vogel

Das Python-Programm zur Erstellung von 3D-CAD-Geländemodellen wurde nach der ersten Veröffentlichung noch um einige Funktionen erweitert. So kann man nun beispielsweise die Höhenauflösung, die standardmäßig 1 cm beträgt, auf einen gröberen Wert einstellen, um architektonische Höhenschichtenmodelle zu erstellen.

[image]
Der Kaisberg in Hagen-Vorhalle als 10-Meter-Quader-Modell
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)

Die STL-Dateien für den 3D-Druck werden nun auch im Binärformat geschrieben, sind dadurch nur noch halb so groß und werden schneller verarbeitet.
Die Unterseite des Modells ist nun nicht mehr fest auf Meeresspiegelhöhe eingestellt, sondern liegt (auf glatte 10 Meter gerundet) etwas unterhalb der geringsten Geländehöhe. Wussten Sie, dass es Punkte in NRW gibt, die fast 300 Meter unterhalb des Meerespiegels liegen? Werfen Sie mal einen Blick in den Tagebau Hambacher Forst, wo einmal der Ort Etzweiler lag: -298,99 m!
Für die Weiterverarbeitung mit CAD-Programmen gibt es nun auch einen direkten DXF-Export. Leider sind in der DXF-Datei aus skurrilen Lizenzierungsgründen keine 3D-Volumenkörper enthalten, sondern nur die umhüllende Fläche. Nach dem Laden sind daher noch zwei Bearbeitungsschritte nötig. Zuerst müssen mit dem AutoCAD/BricsCAD-Befehl REGION alle 3D-Flächen in Regionen umgewandelt werden. Sobald sich das CAD-Programm davon erholt hat, können Sie mit dem Befehl DMHEFTEN den von den Regionen umhüllten Volumenkörper erstellen lassen. Kochen Sie sich in der Zeit, die dazu benötigt wird, ruhig einen Tee. Oder zwei.
In manchen ländlichen Gegenden ist es gar nicht so einfach, herauszufinden, in welchen ZIP-Dateien des NRW-Geodatenservers die benötigten Geländekacheln zu finden sind. Das Programm ermittelt nun anhand eines Kataloges die Namen der maximal noch herunterzuladenden Archivdateien. Da sich die abgedeckten Bereiche überschneiden, genügen manchmal auch weniger Downloads.

Download des Python-Programms: Digitales Geländemodell.ZIP

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

Heute kostenlos: Pythonbuch (PDF)

Tags:
NRW, DGM, Geländemodell, Python 3

NRW-Geländedaten als TXT-Datei für BricsCAD 19

Martin Vogel ⌂ @, Dortmund / Bochum, Mittwoch, 22.05.2019, 10:44 (vor 205 Tagen) @ Martin Vogel

Das bekannte Pythonprogramm zur Geländemodellierung hat nach langer Zeit noch einmal eine Aktualisierung erfahren und erzeugt nun auch TXT-Dateien für den Geländemodellierungsbefehl SITE von BricsCAD 19. Zusammen mit dem ebenfalls neu hinzugekommenen Befehl GRADING lassen sich so auf einfachste Weise Böschungen für Aufschüttung und Aushub modellieren.

[image]

Download: ZIP-Datei

Eine Version für hochaufgelöste Geländedaten eines beliebigen Grundstücks in NRW mit im Schnitt 4 Punkten pro Quadratmeter wäre für dieses Ausgabeformat ebenfalls möglich. Wer Interesse daran hat, möge sich hier melden.

Hier der Quelltext der aktuellen Version:

  1. #!/usr/bin/env python3
  2.  
  3. import os
  4. import sys
  5. import webbrowser
  6. import struct
  7. from tkinter import Tk
  8. from tkinter.filedialog import askdirectory
  9. from math import sqrt, sin, cos, tan, radians, degrees, floor
  10.  
  11. print("""
  12. Dieses Python3-Skript erstellt Höhenmodelle für CAD und 3D-Druck.
  13.  
  14. Mit zwei Dezimalgradkoordinaten geben Sie die Eckpunkte einer rechteckigen
  15. Geländeoberfläche an, für die dieses Programm diverse 3D-Dateien erzeugt:
  16.  
  17. - ein Höhenmodell im STL-Format für 3D-Drucker,
  18. - eine DXF-Datei mit 3D-Flächen, die in Regionen umgewandelt und zu einem
  19. Volumenmodell geheftet werden können,
  20. - eine auf die Auswahl reduzierte XYZ-Datei und
  21. - diverse AutoCAD- bzw. BricsCAD-Befehlsskripte (SCR).
  22. """)
  23.  
  24. # Die CAD-Skripte enthalten Befehle zur Erzeugung von Prismenfeldern aus
  25. # 3D-Körpern, sodass Volumenbefehle angewendet werden können, sowie Befehle
  26. # zur Erzeugung von verschiedenen Oberflächenmodellen.
  27.  
  28. # Das aus Drei- und Vierecken zusammengesetzte Oberflächenmodell ist im
  29. # Ergebnis mit dem Inhalt der DXF-Datei identisch und lässt sich durch
  30. # Anwendung der beiden Befehle "Region" und "Heften" ebenfalls in ein
  31. # Volumenmodell umwandeln.
  32.  
  33. # Das Programm wurde für die öffentlichen Höhendaten für digitale
  34. # Geländemodelle in NRW geschrieben, lässt sich aber mit geringen Anpassungen
  35. # auch für andere DGM-Daten nutzen.
  36.  
  37. # Datenquelle:
  38. # http://www.bezreg-koeln.nrw.de/
  39. # brk_internet/geobasis/hoehenmodelle/gelaendemodelle/
  40. # bzw. https://www.opengeodata.nrw.de/produkte/
  41.  
  42. # Autor: Dipl.-Ing. Martin Vogel, Hochschule Bochum, 2017–2019
  43.  
  44. # Lizenz: Namensnennung - Weitergabe unter gleichen Bedingungen 3.0 Deutschland
  45. # (CC BY-SA 3.0 DE) https://creativecommons.org/licenses/by-sa/3.0/de/
  46.  
  47. # Version 13 vom 10. Mai 2019
  48. # XYZ-Daten als TXT-Datei für Bricscad-BIM-Geländeoberflächen
  49.  
  50. # Version 12 vom 3. Juni 2018
  51. # Anzeige und PDF-Export eines Höhendiagramms mittels Matplotlib.
  52.  
  53. # Version 11 vom 18. Juni 2017
  54. # Zu fehlenden Kacheln wird nun angegeben, in welchen Ortsdateien diese
  55. # zu finden sind. Dazu mussten alle 396 Dateien einmal heruntergeladen und
  56. # diese 269 Gigabyte ZIP-Dateien katalogisiert werden. -> Gelaendekatalog.csv
  57.  
  58. # Version 10 vom 15. Juni 2017
  59. # UTM-Koordinatenumrechnung jetzt ohne externes Modul und nach vereinfachtem
  60. # Algorithmus.
  61.  
  62. # Version 9 vom 13. Juni 2017
  63. # DXF-Export der Dreiecks- und Vierecksflächen aus Version 8.
  64.  
  65. # Version 8 vom 11. Juni 2017
  66. # Oberflächenmodell für CAD-Skript aus einzelnen Dreiecksflächen, wie in der
  67. # STL-Datei, Seitenflächen und Boden jedoch aus Vierecken.
  68.  
  69. # Version 7 vom 9. Juni 2017
  70. # Jetzt mit 3D-Netz, viel schneller als Volumenkörper im CAD-Programm.
  71. # Leider kann BricsCAD aus dieser Oberfläche keinen Volumenkörper heften.
  72.  
  73. # Version 6 vom 21. März 2017
  74. # Jetzt mit Log-Datei zur Garantie der Wiederholbarkeit.
  75.  
  76. # Version 5 vom 19. März 2017
  77. # Jetzt mit variabler Höhe der Unterseite. Es gibt in NRW tatsächlich Flächen,
  78. # die weit unterhalb des Meeresspiegels liegen. (Tagebau Etzweiler: -298,90m!)
  79.  
  80. # Version 4 vom 18. März 2017
  81. # Jetzt mit binärem STL-Export. Die Dateien sind nur noch halb so groß
  82. # und werden viel schneller gelesen als ASCII-STL-Dateien.
  83.  
  84. # Version 3 vom 17. März 2017
  85. # Einige mögliche Fehleingaben werden jetzt freundlich abgefangen.
  86. # Erweiterung um CAD-Skript für Dreiecksprismen, um den Minecraft-Effekt
  87. # loszuwerden.
  88.  
  89. # Version 2 vom 16. März 2017
  90. # Jetzt mit Eingabe der Höhenauflösung.
  91. # Außerdem werden nicht mehr alle gefundenen Punkte im Dictionary
  92. # abgelegt, sondern nur noch die später auch benötigten.
  93.  
  94. # Ideen für Erweiterungen:
  95. # - Grafische Oberfläche (mit Kontrolle der Kachelabdeckung) einbauen.
  96. # - Automatischer Download fehlender Kacheln.
  97. # - In der DXF-Datei noch Werte für Zoom, Blickwinkel, Geokoordinaten
  98. # usw. eintragen.
  99. # - Klären, warum FreeCAD nur eine einzige 3D-Fläche der DXF-Datei darstellt.
  100.  
  101. # Aufräumen:
  102. # - Alte C-Formatstrings durch f-Strings ersetzen.
  103.  
  104. def sysexit():
  105. input("\nProgramm wird beendet. [Enter]")
  106. sys.exit()
  107.  
  108. def log(s, sichtbar=True):
  109. "Bildschirmmeldung mit Protokollierung"
  110. with open(name+".log","a") as logfile:
  111. logfile.write(f"{s}\n")
  112. if sichtbar:
  113. print(f"\n{s}")
  114.  
  115. def utm(Bg,Lg):
  116. "Umrechnung von Breitengrad und Längengrad in UTM-Koordinaten"
  117. # http://www.ottmarlabonde.de/L1/UTMBeispielRechnung.htm
  118. # Literatur: A. Schödlbauer,
  119. # Rechenformeln und Rechenbeispiele zur Landesvermessung,
  120. # Teil 2, Herbert WichmannVerlag Karlsruhe
  121.  
  122. L = radians(Lg)
  123. B = radians(Bg)
  124. tB = tan(B)
  125. cB = cos(B)
  126.  
  127. # Halbachsen des WGS84-Ellipsoids:
  128. # a = 6378137.0
  129. # b = 6356752.314
  130.  
  131. # Radiusreduzierung
  132. mH = 0.9996
  133.  
  134. # c = a**2/b
  135. c = 6399593.626005325
  136.  
  137. # eq = (a**2-b**2)/b**2
  138. eq = 0.006739496819936062
  139.  
  140. # E0 = c*(1-3/4*eq+45/64*eq**2-175/256*eq**3+11025/16384*eq**4)
  141. E0 = 6367449.145759811
  142.  
  143. # E2 = c*(-3/8*eq+15/32*eq**2-525/1024*eq**3+2205/4096*eq**4)
  144. E2 = -16038.508797800609
  145.  
  146. # E4 = c*(15/256*eq**2-105/1024*eq**3+2205/16384*eq**4)
  147. E4 = 16.83262765753934
  148.  
  149. # E6 = c*(-35/3072*eq**3+315/12288*eq**4)
  150. E6 = -0.021980907677118407
  151.  
  152. LL = E0*B + E2*sin(2*B) + E4*sin(4*B) + E6*sin(6*B)
  153. x0 = mH*LL
  154.  
  155. LhL = 3 + 6*floor(degrees(L)/6)
  156. Zone = 30 + (3+LhL)/6
  157.  
  158. DL = L - radians(LhL)
  159.  
  160. etaq = eq*cB**2
  161. Nq = c/sqrt(1+etaq)
  162.  
  163. x2 = mH/2*Nq*tB*cB**2
  164. x4 = mH/24*Nq*tB*cB**4*(5-tB**2+9*etaq)
  165. x6 = mH/720*Nq*tB*cB**6*(61-58*tB**2+tB**4)
  166.  
  167. N = x0 + x2*DL**2 + x4*DL**4 + x6*DL**6
  168.  
  169. y1 = mH*Nq*cB
  170. y3 = mH/6*Nq*cB**3*(1-tB**2+etaq)
  171. y5 = mH/120*Nq*cB**5*(5-18*tB**2+tB**4+etaq*(14-58*tB**2))
  172.  
  173. E = y1*DL + y3*DL**3 + y5*DL**5 + 500000
  174.  
  175. return N, E, Zone
  176.  
  177.  
  178. # Schritt 1: Wie soll das Kind heißen?
  179. print("\nGeben Sie einen Basisnamen für die erzeugten Dateien an!\n"
  180. "Existierende Dateien mit diesem Namen und den Endungen .dxf, .scr,\n"
  181. ".xyz, .pdf und .stl werden ohne Warnung überschrieben.\n")
  182.  
  183. name = input("Dateiname ohne Endung: ")
  184. log("Basisname: "+name, sichtbar=False)
  185.  
  186. # Schritt 2: Wo sind die Geodaten?
  187. print("Bitte wählen Sie den Ordner mit den ausgepackten XYZ-Dateien aus:")
  188. # Für den Dialog müssen wir Tk laden.
  189. Fenster = Tk()
  190. # Wir brauchen das Hauptfenster hier nicht.
  191. Fenster.withdraw()
  192. # Dialogfenster
  193. ordner = askdirectory(title="Bitte den gewünschten Ordner doppelklicken")
  194.  
  195. if not ordner:
  196. print("\nKein Ordner ausgewählt.")
  197. sysexit()
  198.  
  199. log(f"Gewählter Ordner: {ordner}")
  200.  
  201. # Schritt 3: Welches Gebiet wollen wir modellieren?
  202. print("\nGeben Sie zwei gegenüberliegende Eckpunkte des zu modellierenden\n"
  203. "Rechtecks an! Das Eingabeformat ist (Breite, Länge) in Dezimalgrad,\n"
  204. "also beispielsweise 51.335757,7.479087 – Sie können die Koordinaten\n"
  205. "direkt aus Google Maps oder Ihrem GPS-Gerät übernehmen.\n")
  206.  
  207. try:
  208. lat1, lon1 = eval(input("Erstes Eckpunktkoordinatenpaar: "))
  209. lat2, lon2 = eval(input("Gegenüberliegendes Koordinatenpaar: "))
  210. except:
  211. print("Es wurden keine zwei mit einem Komma getrennte Zahlen erkannt.")
  212. sysexit()
  213.  
  214. # Die Eingabe oben war beliebig, wir brauchen aber gleich die Südwestecke
  215. # unten links und die Nordostecke oben rechts:
  216. ul_lat = min(lat1, lat2)
  217. ul_lon = min(lon1, lon2)
  218. or_lat = max(lat1, lat2)
  219. or_lon = max(lon1, lon2)
  220. log(f"Geokoordinaten: {ul_lat},{ul_lon} {or_lat},{or_lon}",
  221. sichtbar=False)
  222.  
  223. # Umrechnung der Dezimalgradkoordinaten ins UTM-System; aus Längen- und
  224. # Breitengrad werden ein Nordwert, ein Ostwert und eine Zonennummer.
  225. n, e, zn = utm(ul_lat, ul_lon)
  226.  
  227. # Die Zonennummer ist in den Dateien, die unter opengeodata.nrw.de
  228. # heruntergeladen werden können, dem Ostwert vorangestellt.
  229. ul_e = int("%i%i"%(zn,e))
  230.  
  231. # Der Nordwert wird übernommen.
  232. ul_n = int(n)
  233.  
  234. # Das gleiche noch einmal für die Nordostecke unseres Rechtecks.
  235. n, e, zn = utm(or_lat, or_lon)
  236. or_e = int("%i%i"%(zn,e))
  237. or_n = int(n)
  238.  
  239. log("\nLage und Größe des ausgewählten Rechtecks:")
  240. log("UTM-Koordinaten: %i,%i %i,%i in Zone %i" % (ul_e,ul_n,or_e,or_n,zn))
  241. log("Ausdehnung Ost-West: %i m" % (or_e-ul_e))
  242. log("Ausdehnung Nord-Süd: %i m" % (or_n-ul_n))
  243. log("Fläche: %i m²" % ((or_e-ul_e)*(or_n-ul_n)))
  244.  
  245. # Liste der zu öffnenden XYZ-Dateien zusammenbauen:
  246. # Bei den NRW-Höhendaten kann man am Dateinamen erkennen, ob gesuchte
  247. # Punkte in ihnen enthalten sein könnten.
  248. # Der Dateiname gibt die untere linke Ecke einer 2000x2000-m²-Kachel an.
  249. # Beispiel für einen Dateinamen: dgm1_32368_5700_2_nw.xyz
  250. # Hier sind die Ostwerte 32368000 bis 32369999 und die Nordwerte
  251. # 5700000 bis n5701999 enthalten.
  252.  
  253. print("\nUntersuche Vollständigkeit der Höhendaten…")
  254. xyz_Liste = []
  255. e_min = 2 * (ul_e // 2000)
  256. e_max = 2 * (or_e // 2000)
  257. n_min = 2 * (ul_n // 2000)
  258. n_max = 2 * (or_n // 2000)
  259. for e in range(e_min,e_max+1,2):
  260. for n in range(n_min,n_max+1,2):
  261. xyz_Liste.append("dgm1_%i_%i_2_nw.xyz" % (e,n))
  262.  
  263. # Sind alle Dateien vorhanden?
  264. fehlende_zip=[]
  265. for xyz_Name in xyz_Liste:
  266. if not os.path.isfile(ordner+"/"+xyz_Name):
  267. print("Die XYZ-Datei %s fehlt!" % xyz_Name)
  268. # print("Sie finden Sie in folgenden ZIP-Archiven:")
  269. with open("Gelaendekatalog.csv") as csv:
  270. for Zeile in csv:
  271. if xyz_Name in Zeile[1:]:
  272. zip_name = Zeile.split()[0]
  273. # print(zip_name)
  274. if zip_name not in fehlende_zip:
  275. fehlende_zip.append(zip_name)
  276.  
  277. if fehlende_zip:
  278. print("\nBitte laden Sie zuerst die fehlenden ZIP-Archive herunter und\n"
  279. "entpacken Sie die darin enthaltenen XYZ-Dateien in den Ordner\n"
  280. "%s.\n" % ordner)
  281. print("In diesen Archiven können Sie die fehlenden Kacheln finden:\n")
  282. print("\n".join(fehlende_zip))
  283. print("\nDie Downloadseite\n"
  284. "https://www.opengeodata.nrw.de/produkte/geobasis/dgm/dgm1/\n"
  285. "wird nun im Webbrowser aufgerufen …")
  286. webbrowser.open("https://www.opengeodata.nrw.de/"
  287. "produkte/geobasis/dgm/dgm1/")
  288. sysexit()
  289.  
  290. print("… alle benötigten Dateien sind vorhanden.")
  291.  
  292. # Alle Dateien sind vorhanden, jetzt kümmern wir uns um die Modelldetails:
  293. print("\nDie horizontale Auflösung der Daten beträgt einen Meter, was bei\n"
  294. "größeren Flächen zu extremen Dateigrößen und Verarbeitungszeiten der\n"
  295. "erzeugten Modelle führen kann. Tipp: beginnen Sie mit ungefähr 1000\n"
  296. "Punkten auf der gesamten Fläche und erhöhen Sie die Auflösung in\n"
  297. "weiteren Durchläufen schrittweise.\n")
  298.  
  299. # Vorschlag für die horizontale Auflösung ermitteln
  300. kl0 = max(1,int(((or_e-ul_e)*(or_n-ul_n)/1000)**0.5))
  301. p0 = ((or_e-ul_e)//kl0) * ((or_n-ul_n)//kl0)
  302. print("Bei %i m Auflösung würden Sie beispielsweise %i Punkte erhalten.\n" % (
  303. kl0, p0))
  304.  
  305. try:
  306. kl = int(input("Geben Sie einen ganzzahligen Wert ein [m]: "))
  307. except:
  308. kl = 1
  309.  
  310. print("\nDie vertikale Auflösung der Daten beträgt einen Zentimeter. Das ist\n"
  311. "normalerweise in Ordnung. Für einen Höhenschichteneffekt wie bei\n"
  312. "architektonischen Geländemodellen kann dieser Wert auch auf 100 oder\n"
  313. "mehr geändert werden, je nach gewünschter Effektstärke.\n"
  314. "Nebeneffekt: Das CAD-Volumenmodell benötigt bei größeren Werten\n"
  315. "deutlich weniger Speicherplatz.\n")
  316.  
  317. try:
  318. kh = int(input("Geben Sie einen ganzzahligen Wert ein [cm]: "))
  319. except:
  320. kh = 1
  321.  
  322. xmax = or_e - (or_e-ul_e) % kl
  323. ymax = or_n - (or_n-ul_n) % kl
  324.  
  325. log("Horizontale Auflösung [m]: %i" % kl, sichtbar=False)
  326. log("Vertikale Auflösung [cm]: %i" % kh, sichtbar=False)
  327. log("\nAbstand der neuen Punkte: %i m" % kl)
  328. log("\nAusdehnung Ost-West: %i m" % (xmax-ul_e))
  329. log("Ausdehnung Nord-Süd: %i m" % (ymax-ul_n))
  330. fqm = (xmax-ul_e)*(ymax-ul_n)
  331. log("Fläche: %i m² bzw. %.3f km²" % (fqm,fqm/1e6))
  332.  
  333. # Alle gefundenen Höhenwerte werden zunächst in ein Dictionary
  334. # geschrieben, aus dem sie für die einzelnen Dateien wieder
  335. # ausgelesen werden.
  336.  
  337. D = {}
  338.  
  339. # Die Höhe der Unterseite ist nicht null, sondern orientiert sich
  340. # am tatsächlichen Gelände.
  341.  
  342. minh = float("inf")
  343. maxh = -float("inf")
  344.  
  345. # Schleife über alle zu verwendenden Eingabedateien
  346. for dateiname in xyz_Liste:
  347. with open(ordner+"/"+dateiname) as dgm:
  348. log("Verwende XYZ-Datei %s" % dateiname)
  349. for zeile in dgm:
  350. # Beispiel für eine Zeile aus Bochum:
  351. # 32372000.00 5706000.00 61.32
  352. try:
  353. x,y,h = zeile.split()
  354. x = int(float(x))
  355. y = int(float(y))
  356. h = float(h)
  357. except ValueError:
  358. print("Abbruch, falsches Format:",zeile)
  359. break
  360. # Koordinaten der Zeile im gesuchten Rechteck?
  361. if ul_e <= x <= or_e and ul_n <= y <= or_n:
  362. # Liegt der Punkt im ausgedünnten Horizontalraster?
  363. if not (x - ul_e) % kl and not (y - ul_n) % kl:
  364. # Punkt wird bei gewählter Auflösung verwendet
  365. if kh != 1:
  366. # Höhenwerte werden gerundet
  367. h = round(round(h*100/kh)*kh/100, 2)
  368. D[x,y] = h
  369. minh = min(h, minh)
  370. maxh = max(h, maxh)
  371.  
  372. log(f"Größte gefundene Höhe: {maxh:.2f} Meter")
  373. log(f"Kleinste gefundene Höhe: {minh:.2f} Meter")
  374. minhs = minh # "echten" minh-Wert merken
  375. # Unterkante des Modells absenken:
  376. minh = 10 * floor(minh/10) - 10
  377. log(f"Setze Unterkante auf {minh:.2f} Meter.")
  378. log(f"Neue Modellhöhe: {maxh-minh:.2f} Meter")
  379.  
  380. # Diagramm anzeigen
  381. print("\nHöhendiagramm wird erzeugt.")
  382. try:
  383. import matplotlib.pyplot as plt
  384.  
  385. # Keine Interaktion, Diagramm nur anzeigen …
  386. plt.rcParams['toolbar'] = 'None'
  387. # Einheitlicher Maßstab für x und y
  388. plt.axis('equal')
  389.  
  390. # Für das Diagramm wird die untere linke Ecke auf (0,0)
  391. # gesetzt. Wer kann schon was mit UTM-Koordinaten anfangen?
  392.  
  393. # Liste der x-Werte
  394. xi = list(range(0,or_e-ul_e+1,kl))
  395. # Liste der y-Werte
  396. yi = list(range(0,or_n-ul_n+1,kl))
  397. # Matrix der Höhenwerte für alle x-y-Paare
  398. zi = [[D[(x+ul_e,y+ul_n)] for x in xi] for y in yi]
  399.  
  400. # Anzahl der Höhenlinien: etwa 10 (7 bis 14)
  401. nh = (maxh-minhs) * 100
  402. while nh > 70: nh = int(nh/10)
  403. if nh > 35: nh = int(nh/5)
  404. if nh > 14: nh = int(nh/2)
  405. # print("Zeichne",nh,"Höhenlinien.")
  406.  
  407. # Höhenlinien
  408. plt.contour(xi, yi, zi, nh, linewidths=1, colors="k")
  409. # farbige Oberfläche
  410. plt.pcolormesh(xi, yi, zi, cmap = plt.get_cmap('terrain'))
  411. # Legende
  412. plt.colorbar()
  413. # Überschrift
  414. plt.title(f"Ausschnittgröße {xmax-ul_e}×{ymax-ul_n} m\n"
  415. f"(0,0) bei {ul_lat}° Nord, {ul_lon}° Ost")
  416. plt.gcf().canvas.set_window_title(f'Höhendiagramm {name}')
  417. # anzeigen und zurück zum Programm …
  418. plt.show(block=False)
  419. # Bild als PDF speichern
  420. ausname = name+".pdf"
  421. log(f"Schreibe PDF-Ausgabedatei: {ausname}")
  422. plt.savefig(ausname, bbox_inches='tight')
  423.  
  424. except:
  425. print("\nFehler: Diagramm kann nicht angezeigt werden."
  426. " Ist matplotlib nicht installiert?\n")
  427. plt=None
  428.  
  429. # Alle Punkte als simple XYZ-Datei sichern
  430.  
  431. ausname = name+".xyz"
  432. log(f"Schreibe XYZ-Ausgabedatei: {ausname}")
  433.  
  434. with open(ausname,"w") as aus:
  435. for x in range(ul_e,or_e+1,kl):
  436. for y in range(ul_n,or_n+1,kl):
  437. aus.write(f"{x} {y} {D[(x,y)]:.2f}\n")
  438.  
  439. # XYZ-Daten für BricsCAD mit Komma trennen und zum Koordinatenursprung
  440. # verschieben.
  441.  
  442. ausname = name+".txt"
  443. log(f"Schreibe BricsCAD-BIM-Geländedatei: {ausname}")
  444.  
  445. with open(ausname,"w") as aus:
  446. for x in range(ul_e,or_e+1,kl):
  447. for y in range(ul_n,or_n+1,kl):
  448. aus.write(f"{x-ul_e},{y-ul_n},{D[(x,y)]:.2f}\n")
  449.  
  450. # DXF-Export
  451.  
  452. ausname = name+".dxf"
  453. log(f"Schreibe DXF-Datei mit 3D-Flächen: {ausname}")
  454.  
  455. # Weil 3D-Solids in einem obskuren „Geheimformat“ gespeichert werden, wird
  456. # hier nur die umhüllende Fläche erzeugt. Diese muss im CAD-Programm mittels
  457. # der Befehle REGION und HEFTEN in einen SOLID umgewandelt werden, um einen
  458. # richtigen Geländekörper zu erhalten.
  459.  
  460. with open(ausname,"w") as aus:
  461. aus.write("0\nSECTION\n2\nHEADER\n"
  462. "9\n$ACADVER\n1\nAC1006\n"
  463. "9\n$INSBASE\n10\n0.0\n20\n0.0\n30\n0.0\n"
  464. "9\n$INSUNITS\n70\n6\n"
  465. "9\n$EXTMIN\n10\n0.0\n20\n0.0\n"
  466. "9\n$EXTMAX\n10\n%f\n20\n%f\n" % (xmax-ul_e, ymax-ul_n) +
  467. "9\n$LIMMIN\n10\n0.0\n20\n0.0\n"
  468. "9\n$LIMMAX\n10\n%f\n20\n%f\n" % (xmax-ul_e, ymax-ul_n) +
  469. "0\nENDSEC\n"
  470. "0\nSECTION\n2\nENTITIES\n")
  471.  
  472. # Geländeoberfläche
  473. for x in range(ul_e,xmax+1-kl,kl):
  474. xu = x - ul_e
  475. for y in range(ul_n,ymax+1-kl,kl):
  476. yu = y - ul_n
  477. h1 = D[x,y]
  478. h2 = D[x+kl,y]
  479. h3 = D[x,y+kl]
  480. h4 = D[x+kl,y+kl]
  481. aus.write("0\n3DFACE\n"
  482. "10\n%i\n20\n%i\n30\n%s\n"%(xu, yu, h1) +
  483. "11\n%i\n21\n%i\n31\n%s\n"%(xu+kl, yu+kl, h4) +
  484. "12\n%i\n22\n%i\n32\n%s\n"%(xu+kl, yu, h2) +
  485. "13\n%i\n23\n%i\n33\n%s\n"%(xu+kl, yu, h2))
  486. aus.write("0\n3DFACE\n"
  487. "10\n%i\n20\n%i\n30\n%s\n"%(xu, yu, h1) +
  488. "11\n%i\n21\n%i\n31\n%s\n"%(xu, yu+kl, h3) +
  489. "12\n%i\n22\n%i\n32\n%s\n"%(xu+kl, yu+kl, h4) +
  490. "13\n%i\n23\n%i\n33\n%s\n"%(xu+kl, yu+kl, h4))
  491.  
  492.  
  493. # Unterseite
  494. aus.write("0\n3DFACE\n"
  495. "10\n%i\n20\n%i\n30\n%.2f\n"%(0,0, minh) +
  496. "11\n%i\n21\n%i\n31\n%.2f\n"%(xmax-ul_e, 0, minh) +
  497. "12\n%i\n22\n%i\n32\n%.2f\n"%(xmax-ul_e, ymax-ul_n, minh) +
  498. "13\n%i\n23\n%i\n33\n%.2f\n"%(0, ymax-ul_n, minh))
  499.  
  500. # Linke Wand
  501. for y in range(ul_n,ymax+1-kl,kl):
  502. yu = y-ul_n
  503. h1 = D[ul_e,y]
  504. h2 = D[ul_e,y+kl]
  505. aus.write("0\n3DFACE\n"
  506. "10\n0\n20\n%i\n30\n%.2f\n"%(yu, minh) +
  507. "11\n0\n21\n%i\n31\n%.2f\n"%(yu+kl, minh) +
  508. "12\n0\n22\n%i\n32\n%.2f\n"%(yu+kl, h2) +
  509. "13\n0\n23\n%i\n33\n%.2f\n"%(yu, h1))
  510.  
  511. # Rechte Wand
  512. for y in range(ul_n,ymax+1-kl,kl):
  513. xu = xmax-ul_e
  514. yu = y-ul_n
  515. h1 = D[xmax,y]
  516. h2 = D[xmax,y+kl]
  517. aus.write("0\n3DFACE\n"
  518. "10\n%i\n20\n%i\n30\n%.2f\n" % (xu, yu, minh) +
  519. "11\n%i\n21\n%i\n31\n%.2f\n" % (xu, yu, h1) +
  520. "12\n%i\n22\n%i\n32\n%.2f\n" % (xu, yu+kl, h2) +
  521. "13\n%i\n23\n%i\n33\n%.2f\n" % (xu, yu+kl, minh))
  522.  
  523. # Vordere Wand
  524. for x in range(ul_e,xmax+1-kl,kl):
  525. h1 = D[x,ul_n]
  526. h2 = D[x+kl,ul_n]
  527. xu = x - ul_e
  528. aus.write("0\n3DFACE\n"
  529. "10\n%i\n20\n0\n30\n%.2f\n" % (xu, minh) +
  530. "11\n%i\n21\n0\n31\n%.2f\n" % (xu, h1) +
  531. "12\n%i\n22\n0\n32\n%.2f\n" % (xu+kl, h2) +
  532. "13\n%i\n23\n0\n33\n%.2f\n" % (xu+kl, minh))
  533.  
  534. # Hintere Wand
  535. for x in range(ul_e,xmax+1-kl,kl):
  536. xu = x - ul_e
  537. yu = ymax - ul_n
  538. h1 = D[x,ymax]
  539. h2 = D[x+kl,ymax]
  540. aus.write("0\n3DFACE\n"
  541. "10\n%i\n20\n%i\n30\n%.2f\n"%(xu, yu, minh) +
  542. "11\n%i\n21\n%i\n31\n%.2f\n"%(xu+kl, yu, minh) +
  543. "12\n%i\n22\n%i\n32\n%.2f\n"%(xu+kl, yu, h2) +
  544. "13\n%i\n23\n%i\n33\n%.2f\n"%(xu, yu, h1))
  545.  
  546. aus.write("0\nENDSEC\n0\nEOF\n")
  547.  
  548.  
  549.  
  550. # Diverse Skripte für AutoCAD/BricsCAD
  551.  
  552. def scr_intro():
  553. "Nimmt Grundeinstellungen des Zeichnungseditors vor"
  554. # Rückgängigmachen ausschalten
  555. # Koordinateneingabe hat Vorrang vor Objektfang,
  556. # Einheit Meter,
  557. # Ansicht schräg von Südwesten,
  558. # Visueller Stil Drahtmodell (Geschwindigkeit!),
  559. # Zoom auf interessanten Quader
  560. return ("Zurück S K\n"
  561. "OSnapCoord 1\n"
  562. "_InsUnits 6\n"
  563. "APunkt -1,-2,2\n"
  564. "-Vis A Drahtmodell\n"
  565. "Zoom F 0,0,%f %i,%i,%f\n" % (minh, or_e-ul_e, or_n-ul_n, maxh))
  566.  
  567. # Einstellungen nach Skriptende:
  568. def scr_exit():
  569. # Rückgängigmachen einschalten
  570. # Zoom auf Grenzen
  571. return ("Zurück A\n"
  572. "Zoom G\n")
  573.  
  574.  
  575. ausname = name+".quadratprismen.scr"
  576. log(f"Schreibe CAD-Skriptdatei mit Quadratprismenfeld: {ausname}")
  577.  
  578. # Eigentlich wäre es universeller, anstelle der deutschsprachigen
  579. # Bezeichner englischsprachige Bezeichner mit vorangestelltem Unterstrich
  580. # zu verwenden. Dummerweise kommt BricsCAD damit nicht immer zurecht.
  581. with open(ausname,"w") as aus:
  582. aus.write(scr_intro())
  583.  
  584. for x in range(ul_e,or_e+1,kl):
  585. for y in range(ul_n,or_n+1,kl):
  586. aus.write("Quader %i,%i,%.2f %i,%i,%.2f\n" % (
  587. x-ul_e, y-ul_n, minh,
  588. x-ul_e+kl, y-ul_n+kl, D[(x,y)]))
  589.  
  590. aus.write(scr_exit())
  591.  
  592. ausname = name+".dreiecksprismen.scr"
  593. log(f"Schreibe CAD-Skriptdatei mit Dreiecksprismenfeld: {ausname}")
  594.  
  595. with open(ausname,"w") as aus:
  596. aus.write(scr_intro())
  597.  
  598. for x in range(ul_e,xmax+1-kl,kl):
  599. for y in range(ul_n,ymax+1-kl,kl):
  600. h1 = float(D[x,y])
  601. x1, y1 = x-ul_e, y-ul_n
  602. h2 = float(D[x+kl,y])
  603. x2, y2 = x1+kl, y1
  604. h3 = float(D[x+kl,y+kl])
  605. x3, y3 = x2, y1+kl
  606. h4 = float(D[x,y+kl])
  607. x4, y4 = x1, y3
  608. hmax = max(h1, h2, h3, h4)
  609. # Dreieck auf Nullebene zeichnen …
  610. aus.write("3DPoly %i,%i,%.2f %i,%i,%.2f %i,%i,%.2f S\n" % (
  611. x1,y1,minh, x2,y2,minh, x3,y3,minh))
  612. # … bis zum höchsten Punkt hochziehen …
  613. aus.write("_extrude L %f\n" % (hmax-minh))
  614. # … und oben schräg abschneiden.
  615. aus.write("Kappen L %i,%i,%f %i,%i,%f %i,%i,%f %i,%i,%f\n" %
  616. (x1,y1,h1, x2,y2,h2, x3,y3,h3, x1,y1,minh))
  617. # Und nochmal für das zweite Dreieck:
  618. aus.write("3DPoly %i,%i,%.2f %i,%i,%.2f %i,%i,%.2f S\n" % (
  619. x1,y1,minh, x4,y4,minh, x3,y3,minh))
  620. aus.write("_extrude L %f\n" % hmax)
  621. aus.write("Kappen L %i,%i,%f %i,%i,%f %i,%i,%f %i,%i,%f\n" %
  622. (x1,y1,h1, x4,y4,h4, x3,y3,h3, x1,y1,minh))
  623.  
  624. aus.write(scr_exit())
  625.  
  626. ausname = name+".mesh.scr"
  627. log(f"Schreibe CAD-Skriptdatei mit 3D-Netz (Mesh): {ausname}")
  628.  
  629. # Beim Quadratnetz hat BricsCAD die Einschränkung, dass es maximal
  630. # 256×256 Knoten haben darf. Daher muss hier ausgedünnt werden.
  631.  
  632. with open(ausname,"w") as aus:
  633. aus.write(scr_intro())
  634. n = 1
  635. klm = kl
  636. while (or_e-ul_e+1)//klm > 256 or (or_n-ul_n+1)//klm > 256:
  637. n += 1
  638. klm = kl * n
  639. if n>1:
  640. log("Verwende nur jeden %i. Punkt pro Richtung." % n)
  641.  
  642. # Gelände
  643. aus.write("3dnetz %i %i\n" % ((or_e-ul_e)//klm+1,(or_n-ul_n)//klm+1))
  644. for x in range(ul_e,or_e+1,klm):
  645. for y in range(ul_n,or_n+1,klm):
  646. aus.write("%i,%i,%f\n" % (x-ul_e,y-ul_n,D[(x,y)]))
  647. # Seitenflächen
  648. y = ul_n
  649. aus.write("3dnetz %i 2\n" % ((or_e-ul_e)//klm+1))
  650. for x in range(ul_e,or_e+1,klm):
  651. aus.write("%i,%i,%f\n" % (x-ul_e,y-ul_n,minh))
  652. aus.write("%i,%i,%f\n" % (x-ul_e,y-ul_n,D[(x,y)]))
  653. y = range(ul_n,or_n+1,klm)[-1]
  654. aus.write("3dnetz %i 2\n" % ((or_e-ul_e)//klm+1))
  655. for x in range(ul_e,or_e+1,klm):
  656. aus.write("%i,%i,%f\n" % (x-ul_e,y-ul_n,minh))
  657. aus.write("%i,%i,%f\n" % (x-ul_e,y-ul_n,D[(x,y)]))
  658. x = ul_e
  659. aus.write("3dnetz %i 2\n" % ((or_n-ul_n)//klm+1))
  660. for y in range(ul_n,or_n+1,klm):
  661. aus.write("%i,%i,%f\n" % (x-ul_e,y-ul_n,minh))
  662. aus.write("%i,%i,%f\n" % (x-ul_e,y-ul_n,D[(x,y)]))
  663. x = range(ul_e,or_e+1,klm)[-1]
  664. aus.write("3dnetz %i 2\n" % ((or_n-ul_n)//klm+1))
  665. for y in range(ul_n,or_n+1,klm):
  666. aus.write("%i,%i,%f\n" % (x-ul_e,y-ul_n,minh))
  667. aus.write("%i,%i,%f\n" % (x-ul_e,y-ul_n,D[(x,y)]))
  668. # Boden
  669. x_max = range(ul_e,or_e+1,klm)[-1]-ul_e
  670. y_max = range(ul_n,or_n+1,klm)[-1]-ul_n
  671. aus.write("3dfläche\n")
  672. aus.write("0,0,%f " % (minh))
  673. aus.write("0,%i,%f " % (y_max, minh))
  674. aus.write("%i,%i,%f " % (x_max, y_max, minh))
  675. aus.write("%i,0,%f \n" % (x_max, minh))
  676.  
  677. aus.write(scr_exit())
  678.  
  679. ausname = name+".3dflächen.scr"
  680. log(f"Schreibe CAD-Scriptdatei mit 3D-Flächen: {ausname}")
  681.  
  682. # Das 3D-Netz oben lässt sich in BricsCAD blöderweise nicht zu einem
  683. # Solid umformen. Hier wird daher nun derselbe Algorithmus verwendet,
  684. # der auch für STL-Dateien zum Einsatz kommt.
  685. # Für die Seiten und den Boden werden Vierecke verwendet.
  686.  
  687. with open(ausname,"w") as aus:
  688. aus.write(scr_intro())
  689.  
  690. # Geländeoberfläche
  691. for x in range(ul_e,xmax+1-kl,kl):
  692. xu = x - ul_e
  693. for y in range(ul_n,ymax+1-kl,kl):
  694. yu = y - ul_n
  695. h1 = D[x,y]
  696. h2 = D[x+kl,y]
  697. h3 = D[x,y+kl]
  698. h4 = D[x+kl,y+kl]
  699. aus.write("3dfläche\n"
  700. "%i,%i,%s\n"%(xu, yu, h1) +
  701. "%i,%i,%s\n"%(xu+kl, yu+kl, h4) +
  702. "%i,%i,%s\n\n\n"%(xu+kl, yu, h2))
  703. aus.write("3dfläche\n"
  704. "%i,%i,%s\n"%(xu, yu, h1) +
  705. "%i,%i,%s\n"%(xu, yu+kl, h3) +
  706. "%i,%i,%s\n\n\n"%(xu+kl, yu+kl, h4))
  707.  
  708.  
  709. # Unterseite
  710. aus.write("3dfläche\n"
  711. "%i,%i,%.2f\n"%(0,0, minh) +
  712. "%i,%i,%.2f\n"%(xmax-ul_e, 0, minh) +
  713. "%i,%i,%.2f\n"%(xmax-ul_e, ymax-ul_n, minh) +
  714. "%i,%i,%.2f\n\n"%(0, ymax-ul_n, minh))
  715.  
  716. # Linke Wand
  717. for y in range(ul_n,ymax+1-kl,kl):
  718. yu = y-ul_n
  719. h1 = D[ul_e,y]
  720. h2 = D[ul_e,y+kl]
  721. aus.write("3dfläche\n"
  722. "0,%i,%.2f\n"%(yu, minh) +
  723. "0,%i,%.2f\n"%(yu+kl, minh) +
  724. "0,%i,%.2f\n"%(yu+kl, h2) +
  725. "0,%i,%.2f\n\n"%(yu, h1))
  726.  
  727. # Rechte Wand
  728. for y in range(ul_n,ymax+1-kl,kl):
  729. xu = xmax-ul_e
  730. yu = y-ul_n
  731. h1 = D[xmax,y]
  732. h2 = D[xmax,y+kl]
  733. aus.write("3dfläche\n"
  734. "%i,%i,%.2f\n" % (xu, yu, minh) +
  735. "%i,%i,%.2f\n" % (xu, yu, h1) +
  736. "%i,%i,%.2f\n" % (xu, yu+kl, h2) +
  737. "%i,%i,%.2f\n\n" % (xu, yu+kl, minh))
  738.  
  739. # Vordere Wand
  740. for x in range(ul_e,xmax+1-kl,kl):
  741. h1 = D[x,ul_n]
  742. h2 = D[x+kl,ul_n]
  743. xu = x - ul_e
  744. aus.write("3dfläche\n"
  745. "%i,0,%.2f\n" % (xu, minh) +
  746. "%i,0,%.2f\n" % (xu, h1) +
  747. "%i,0,%.2f\n" % (xu+kl, h2) +
  748. "%i,0,%.2f\n\n" % (xu+kl, minh))
  749.  
  750. # Hintere Wand
  751. for x in range(ul_e,xmax+1-kl,kl):
  752. xu = x - ul_e
  753. yu = ymax - ul_n
  754. h1 = D[x,ymax]
  755. h2 = D[x+kl,ymax]
  756. aus.write("3dfläche\n"
  757. "%i,%i,%.2f\n"%(xu, yu, minh) +
  758. "%i,%i,%.2f\n"%(xu+kl, yu, minh) +
  759. "%i,%i,%.2f\n"%(xu+kl, yu, h2) +
  760. "%i,%i,%.2f\n\n"%(xu, yu, h1))
  761.  
  762. aus.write(scr_exit())
  763.  
  764. ausname = name+".ascii.stl"
  765. log(f"Schreibe STL-Datei für 3D-Druck (ASCII): {ausname}")
  766.  
  767. # Die Flächen umhüllen einen Quader, der unten auf Höhe minh aufliegt
  768. # und oben durch die Geländeoberfläche abgeschnitten wird.
  769.  
  770. # Für die „Flächennormalen“ wird hier immer eine der Achsenrichtungen
  771. # eingesetzt. Für die Geländeoberfläche ist das beispielsweise die
  772. # nach oben gerichtete z-Achse (0, 0, 1).
  773. # Der genaue Winkel ist auch gar nicht relevant. Hauptsache,
  774. # die „Normale“ zeigt irgendwie aus dem umhüllten Körper heraus
  775. # und nicht in ihn hinein.
  776.  
  777. with open(ausname,"w") as aus:
  778. aus.write("solid "+ausname+"\n")
  779.  
  780. # Geländeoberfläche
  781. for x in range(ul_e,xmax+1-kl,kl):
  782. xu = x - ul_e
  783. for y in range(ul_n,ymax+1-kl,kl):
  784. yu = y - ul_n
  785. h1 = D[x,y]
  786. h2 = D[x+kl,y]
  787. h3 = D[x,y+kl]
  788. h4 = D[x+kl,y+kl]
  789. aus.write("facet normal 0 0 1\n"
  790. "outer loop\n"
  791. "vertex %i %i %s\n"%(xu, yu, h1) +
  792. "vertex %i %i %s\n"%(xu+kl, yu+kl, h4) +
  793. "vertex %i %i %s\n"%(xu+kl, yu, h2) +
  794. "endloop\n"
  795. "endfacet\n")
  796. aus.write("facet normal 0 0 1\n"
  797. "outer loop\n"
  798. "vertex %i %i %s\n"%(xu, yu, h1) +
  799. "vertex %i %i %s\n"%(xu, yu+kl, h3) +
  800. "vertex %i %i %s\n"%(xu+kl, yu+kl, h4) +
  801. "endloop\n"
  802. "endfacet\n")
  803. # Was für ein Aufwand für zwei popelige Dreiecke!
  804.  
  805. # Die Unterseite und die Seitenflächen wiederholen im Moment
  806. # die Unterteilung der Oberseite. Dadurch hat die STL-Datei fast
  807. # doppelt so viele Dreiecke wie eigentlich nur nötig wären.
  808. # Durch geschickte Aufteilung könnte die Unterseite mit nur zwei
  809. # Dreiecken realisiert werden, wobei die Seitenflächen mit rund der
  810. # Hälfte der Dreiecke auskommen könnten. Bei großen Geländesteigungen
  811. # gibt es da aber ein paar Herausforderungen an den Algorithmus.
  812.  
  813. # Unterseite
  814. for x in range(ul_e,xmax+1-kl,kl):
  815. xu = x - ul_e
  816. for y in range(ul_n,ymax+1-kl,kl):
  817. yu = y - ul_n
  818. aus.write("facet normal 0 0 -1\n"
  819. "outer loop\n"
  820. "vertex %i %i %.2f\n"%(xu, yu, minh) +
  821. "vertex %i %i %.2f\n"%(xu+kl, yu, minh) +
  822. "vertex %i %i %.2f\n"%(xu+kl, yu+kl, minh) +
  823. "endloop\n"
  824. "endfacet\n")
  825. aus.write("facet normal 0 0 -1\n"
  826. "outer loop\n"
  827. "vertex %i %i %.2f\n"%(xu, yu, minh) +
  828. "vertex %i %i %.2f\n"%(xu+kl, yu+kl, minh) +
  829. "vertex %i %i %.2f\n"%(xu, yu+kl, minh) +
  830. "endloop\n"
  831. "endfacet\n")
  832.  
  833. # Linke Wand
  834. for y in range(ul_n,ymax+1-kl,kl):
  835. yu = y-ul_n
  836. h1 = D[ul_e,y]
  837. h2 = D[ul_e,y+kl]
  838. aus.write("facet normal -1 0 0\n"
  839. "outer loop\n"
  840. "vertex 0 %i %.2f\n"%(yu, minh) +
  841. "vertex 0 %i %.2f\n"%(yu+kl, h2) +
  842. "vertex 0 %i %.2f\n"%(yu, h1) +
  843. "endloop\n"
  844. "endfacet\n")
  845. aus.write("facet normal -1 0 0\n"
  846. "outer loop\n"
  847. "vertex 0 %i %.2f\n"%(yu, minh) +
  848. "vertex 0 %i %.2f\n"%(yu+kl, minh) +
  849. "vertex 0 %i %.2f\n"%(yu+kl, h2) +
  850. "endloop\n"
  851. "endfacet\n")
  852.  
  853. # Rechte Wand
  854. for y in range(ul_n,ymax+1-kl,kl):
  855. xu = xmax-ul_e
  856. yu = y-ul_n
  857. h1 = D[xmax,y]
  858. h2 = D[xmax,y+kl]
  859. aus.write("facet normal 1 0 0\n"
  860. "outer loop\n"
  861. "vertex %i %i %.2f\n" % (xu, yu, minh) +
  862. "vertex %i %i %.2f\n" % (xu, yu, h1) +
  863. "vertex %i %i %.2f\n" % (xu, yu+kl, h2) +
  864. "endloop\n"
  865. "endfacet\n")
  866. aus.write("facet normal 1 0 0\n"
  867. "outer loop\n"
  868. "vertex %i %i %.2f\n" % (xu, yu, minh) +
  869. "vertex %i %i %.2f\n" % (xu, yu+kl, h2) +
  870. "vertex %i %i %.2f\n" % (xu, yu+kl, minh) +
  871. "endloop\n"
  872. "endfacet\n")
  873.  
  874. # Vordere Wand
  875. for x in range(ul_e,xmax+1-kl,kl):
  876. h1 = D[x,ul_n]
  877. h2 = D[x+kl,ul_n]
  878. xu = x - ul_e
  879. aus.write("facet normal 0 -1 0\n"
  880. "outer loop\n"
  881. "vertex %i 0 %.2f\n" % (xu, minh) +
  882. "vertex %i 0 %.2f\n" % (xu, h1) +
  883. "vertex %i 0 %.2f\n" % (xu+kl, h2) +
  884. "endloop\n"
  885. "endfacet\n")
  886. aus.write("facet normal 0 -1 0\n"
  887. "outer loop\n"
  888. "vertex %i 0 %.2f\n" % (xu, minh) +
  889. "vertex %i 0 %.2f\n" % (xu+kl, h2) +
  890. "vertex %i 0 %.2f\n" % (xu+kl, minh) +
  891. "endloop\n"
  892. "endfacet\n")
  893.  
  894. # Hintere Wand
  895. for x in range(ul_e,xmax+1-kl,kl):
  896. xu = x - ul_e
  897. yu = ymax - ul_n
  898. h1 = D[x,ymax]
  899. h2 = D[x+kl,ymax]
  900. aus.write("facet normal 0 1 0\n"
  901. "outer loop\n"
  902. "vertex %i %i %.2f\n"%(xu, yu, minh) +
  903. "vertex %i %i %.2f\n"%(xu+kl, yu, h2) +
  904. "vertex %i %i %.2f\n"%(xu, yu, h1) +
  905. "endloop\n"
  906. "endfacet\n")
  907. aus.write("facet normal 0 1 0\n"
  908. "outer loop\n"
  909. "vertex %i %i %.2f\n"%(xu, yu, minh) +
  910. "vertex %i %i %.2f\n"%(xu+kl, yu, minh) +
  911. "vertex %i %i %.2f\n"%(xu+kl, yu, h2) +
  912. "endloop\n"
  913. "endfacet\n")
  914.  
  915.  
  916. ### Binäre STL-Datei ###
  917. #
  918. # https://de.wikipedia.org/wiki/STL-Schnittstelle
  919. #
  920. # Die struct-Codes unten bedeuten:
  921. # "<" niedrigwertiges Byte voran ("little-endian")
  922. # "I" vorzeichenlose Ganzzahl (4 Byte)
  923. # "12f" 12 Gleitkommawerte (je 4 Byte)
  924. # "h" kurze Ganzzahl (2 Byte)
  925.  
  926. ausname = name+".binär.stl"
  927. log(f"Schreibe STL-Datei für 3D-Druck (binär): {ausname}")
  928.  
  929. with open(ausname,"wb") as aus:
  930. # 80 Bytes ungenutzter Header
  931. aus.write(b'\0' * 80)
  932.  
  933. # Wie viele Dreiecke hat das Modell insgesamt?
  934. nx = (xmax - ul_e) // kl
  935. ny = (ymax - ul_n) // kl
  936. ngesamt = 4 * nx*ny + 4 * nx + 4 * ny
  937. aus.write(struct.pack('<I', ngesamt))
  938.  
  939. # Geländeoberfläche
  940. for x in range(ul_e,xmax+1-kl,kl):
  941. xu = x - ul_e
  942. for y in range(ul_n,ymax+1-kl,kl):
  943. yu = y - ul_n
  944. h1 = D[x,y]
  945. h2 = D[x+kl,y]
  946. h3 = D[x,y+kl]
  947. h4 = D[x+kl,y+kl]
  948. aus.write(struct.pack("<12fh",
  949. 0, 0, 1,
  950. xu, yu, h1,
  951. xu+kl, yu+kl, h4,
  952. xu+kl, yu, h2,
  953. 0))
  954. aus.write(struct.pack("<12fh",
  955. 0, 0, 1,
  956. xu, yu, h1,
  957. xu, yu+kl, h3,
  958. xu+kl, yu+kl, h4,
  959. 0))
  960.  
  961. # Unterseite
  962. for x in range(ul_e,xmax+1-kl,kl):
  963. xu = x - ul_e
  964. for y in range(ul_n,ymax+1-kl,kl):
  965. yu = y - ul_n
  966. aus.write(struct.pack("<12fh",
  967. 0, 0, -1,
  968. xu, yu, minh,
  969. xu+kl, yu+kl, minh,
  970. xu+kl, yu, minh,
  971. 0))
  972. aus.write(struct.pack("<12fh",
  973. 0, 0, -1,
  974. xu, yu, minh,
  975. xu, yu+kl, minh,
  976. xu+kl, yu+kl, minh,
  977. 0))
  978.  
  979. # Linke Wand
  980. for y in range(ul_n,ymax+1-kl,kl):
  981. yu = y-ul_n
  982. h1 = D[ul_e,y]
  983. h2 = D[ul_e,y+kl]
  984. aus.write(struct.pack("<12fh",
  985. -1, 0, 0,
  986. 0, yu, minh,
  987. 0, yu+kl, h2,
  988. 0, yu, h1,
  989. 0))
  990. aus.write(struct.pack("<12fh",
  991. -1, 0, 0,
  992. 0, yu, minh,
  993. 0, yu+kl, minh,
  994. 0, yu+kl, h2,
  995. 0))
  996.  
  997. # Rechte Wand
  998. xu = xmax-ul_e
  999. for y in range(ul_n,ymax+1-kl,kl):
  1000. yu = y-ul_n
  1001. h1 = D[xmax,y]
  1002. h2 = D[xmax,y+kl]
  1003. aus.write(struct.pack("<12fh",
  1004. 1, 0, 0,
  1005. xu, yu, minh,
  1006. xu, yu, h1,
  1007. xu, yu+kl, h2,
  1008. 0))
  1009. aus.write(struct.pack("<12fh",
  1010. 1, 0, 0,
  1011. xu, yu, minh,
  1012. xu, yu+kl, h2,
  1013. xu, yu+kl, minh,
  1014. 0))
  1015.  
  1016. # Vordere Wand
  1017. for x in range(ul_e,xmax+1-kl,kl):
  1018. h1 = D[x,ul_n]
  1019. h2 = D[x+kl,ul_n]
  1020. xu = x - ul_e
  1021. aus.write(struct.pack("<12fh",
  1022. 0, -1, 0,
  1023. xu, 0, minh,
  1024. xu, 0, h1,
  1025. xu+kl, 0, h2,
  1026. 0))
  1027. aus.write(struct.pack("<12fh",
  1028. 0, -1, 0,
  1029. xu, 0, minh,
  1030. xu+kl, 0, h2,
  1031. xu+kl, 0, minh,
  1032. 0))
  1033.  
  1034. # Hintere Wand
  1035. yu = ymax - ul_n
  1036. for x in range(ul_e,xmax+1-kl,kl):
  1037. xu = x - ul_e
  1038. h1 = D[x,ymax]
  1039. h2 = D[x+kl,ymax]
  1040. aus.write(struct.pack("<12fh",
  1041. 0, 1, 0,
  1042. xu, yu, minh,
  1043. xu+kl, yu, h2,
  1044. xu, yu, h1,
  1045. 0))
  1046. aus.write(struct.pack("<12fh",
  1047. 0, 1, 0,
  1048. xu, yu, minh,
  1049. xu+kl, yu, minh,
  1050. xu+kl, yu, h2,
  1051. 0))
  1052.  
  1053. log("Programmlauf erfolgreich beendet.\n\n")
  1054. print("Die Ausgabedateien können nun weiterverarbeitet werden.")
  1055.  
  1056. input("\nProgramm schließen mit [Enter]")
  1057. if plt:
  1058. plt.close()

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

Heute kostenlos: Pythonbuch (PDF)

Tags:
BIM, Site, Grading, Geländemodellierung, Aufschüttung, Auftrag, Aushub, Böschungswinkel

Farbige DXF-Höhenmodelle

Martin Vogel ⌂ @, Dortmund / Bochum, Samstag, 22.06.2019, 16:35 (vor 174 Tagen) @ Martin Vogel

Das in Python 3 geschriebene Geländemodellierungsprogramm für NRW hat noch ein kleines Extra bekommen: Die 3D-Flächen in der vom Programm aus den digitalen Geländedaten erzeugten DXF-Datei werden nun entsprechend ihrer relativen Höhe im Modell eingefärbt. Leider ist die mögliche Farbauswahl durch die praxisfernen „Indexfarben“ des DWG-Standards ziemlich beschränkt, das Ergebnis ist aber vielleicht dennoch ganz brauchbar.

Download: ZIP-Datei

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

Kurzanleitung:
ZIP-Datei entpacken und das Programm „Geländemodell.py“ starten. Gegebenenfalls muss vorher Python 3 heruntergeladen und installiert werden.
Das Programm möchte zuerst wissen, wie die erzeugten Dateien heißen sollen (zum Beispiel "Halde Gotthelf"), anschließend ist anzugeben, wo die schon heruntergeladenen oder noch herunterzuladenden DGM- und XYZ-Daten des NRW-Geodatenservers liegen (wenn Sie nicht wissen, was das bedeutet, drücken Sie einfach die Eingabetaste) und schließlich müssen nur noch die Eckpunkte des zu modellierenden Bereichs als Geokoordinaten eingegeben werden. Für das Modell im Bild waren das beispielsweise die beiden Koordinaten 51.468775,7.447726 und 51.472981,7.454811.

Die Geokoordinaten eines Geländepunktes erfahren Sie beispielsweise, indem sie diesen in Google Maps mit der Maus rechtsklicken und im dann erscheinenden Kontextmenü den Eintrag „was ist hier“ auswählen.

Bei einem Meter horizontaler Auflösung wird die DXF-Datei des hier abgebildeten Modells knapp 60 MB groß. Durch Verringerung der Auflösung auf 10 Meter reduziert sich die Dateigröße auf ein Hundertstel, dafür wird das aus Dreiecksflächen zusammengesetzte Modell gröber.

Ausser einem DXF-Modell werden auch STL-Dateien für den 3D-Druck erzeugt, sodass Sie das Geländemodell gleich körperlich herstellen können.
[image]
Etwas Geduld ist bei der heutigen 3D-Druckertechnik allerdings vonnöten. Mit 15 Stunden Druckzeit für ein 20 Zentimeter breites Modell müssen Sie schon rechnen.

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

Heute kostenlos: Pythonbuch (PDF)

Tags:
GIS, AutoCAD, BIM, BricsCAD, Geländemodellierung, DXF

RSS-Feed dieser Diskussion
powered by my little forum