Kostenlose Digitale Geländemodelle aus ganz NRW für BricsCAD und AutoCAD (Allgemeines)

Martin Vogel ⌂ @, Dortmund / Bochum, Fri, 29.05.2020, 12:28 (vor 1390 Tagen)

Auf dem Geodatenserver der Bezirksregierung Köln hat es mehrere interessante Änderungen gegeben.

Zum einen sind die zu Gemeinden paketierten xyz-Höhendaten nun schon wieder in einem anderen Downloadverzeichnis zu finden (aktuell ist dieser Link gültig). Diese Daten wurden bisher von meinem Python-Programm zur Erstellung von Digitalen Geländemodellen beliebiger Grundstücke in NRW verwendet, das nun eigentlich erneut für diese URL anzupassen gewesen wäre.

Das geschieht aber nicht, denn zusätzlich zu den mehrere Gigabyte großen ZIP-Dateien gibt es nun alle 4-km²-Kacheln mit jeweils 4 Millionen zentimetergenauen NRW-Höhenwerten einzeln als nur noch knapp 20 MB kleine GNU-zip-Dateien. Das langwierige Herunterladen und Auspacken der riesigen ZIP-Archive entfällt endlich und das Programm lädt nun nur noch einzelne kleine Dateien herunter, falls es diese nicht ohnehin schon in einem früheren Programmlauf auf dem Rechner gespeichert hat.

Die dritte Änderung betrifft die Lizenz der Daten. Diese sind nun ohne jede Einschränkung nutzbar. Das Lizenzmodell heißt „Datenlizenz Deutschland - Zero - Version 2.0

[image]
Die Halden Hoppenbruch und Hoheward als TIN-Oberflächen in BricsCAD V20

Um das kostenlose Programm „Geländemodell.py“ auszuführen, benötigen Sie lediglich eine aktuelle Python-Version. Das Programm fragt Sie zuerst nach dem Verzeichnis für die herunterzuladenden Geländedaten und anschließend nach dem Basisnamen der anzulegenden Dateien. Angenommen, Sie möchten ein 3D-Modell der Halde Hoheward in Herten erstellen, so bietet sich als Basisname natürlich so etwas wie „Hoheward“ an. Als nächstes sind zwei Geokoordinaten einzugeben, die Sie zum Beispiel durch zweimaliges Anklicken eines Punktes in einer Google-Maps-Karte erhalten. Bei der Halde Hoheward sind das vielleicht die beiden Punkte „51.553891, 7.149416“ und „51.578594, 7.178255“.

Weitere Benutzereingaben sind gar nicht notwendig. Das Programm lädt die fehlenden Geländkacheln vom Geodatenserver herunter und baut daraus eine Anzahl verschiedener 3D-Dateien zusammen. Während die Dateien geschrieben werden, zeigt das Python-Programm einen Höhenplan des ausgewählten Geländes an.

[image]

Als Ausgabedateien erzeugt das Programm:

Download des Python-Programms:

Digitales Geländemodell V16.zip

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

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

Bücher:
CAD mit BricsCAD
Bauinformatik mit Python

Tags:
Python, CAD, GIS, NRW, AutoCAD, BIM, BricsCAD, Kostenlos, XYZ, Höhendaten, Höhenplan, Digitales Geländemodell


gesamter Thread:

 RSS-Feed dieser Diskussion

powered by my little forum