Topographische NRW-Karten mit Matplotlib (Software)

Martin Vogel ⌂ @, Dortmund / Bochum, Freitag, 15.06.2018, 13:18 (vor 101 Tagen)

Mein Python-3-Programm, das die Höhendaten des Geodatenservers NRW in dreidimensionale CAD-Geländemodelle übertragen kann (siehe 3D-Modell jedes beliebigen Grundstücks in NRW und Aktualisierte Version des Geländemodellierers) erzeugt nun auch eine farbige topograhische Karte mit Höhenlinien. Voraussetzung ist, dass das Modul „matplotlib“ installiert wurde.

Download: ZIP-Datei

Die Karte lässt sich in CAD-Programmen wie AutoCAD oder BricsCAD als Materialtextur auf die Modelloberfläche legen. Die Grafik zeigt links ein Orthofoto vom NRW-Geodatenserver und rechts die mit matplotlib erzeugte Karte. Darunter wurden beide Grafiken als Oberflächentextur des 3D-Geländemodells verwendet.

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

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

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

Tags:
Python, GIS, NRW, Maps, Karte

RSS-Feed dieser Diskussion
powered by my little forum