[PYTHON] Exportieren Sie 3D-Daten aus QGIS
Einführung
Exportieren wir die GIS-Daten mithilfe der in QGIS integrierten Python-Konsole in eine OBJ-Datei.
- Für OBJ-Dateien legen Sie bitte den [Wikipedia-Link] vor (https://ja.wikipedia.org/wiki/Wavefront_.obj). Es ist ein sehr einfaches Format, das Sie als Text exportieren können. Es ist also einfach zu handhaben, wenn Sie Ihr eigenes exportieren möchten.
Ausgeschnitten
Schneiden Sie zuerst den Zielbereich aus. Wenn Sie die vorherigen Daten so verwenden, wie sie sind, sind sie zu groß und das Echtzeit-Rendering kann sie nicht verarbeiten. Versuchen Sie daher, sie auf einen geeigneten Bereich auszuschneiden. Dieses Mal werde ich die Gegend um Roppongi Hills ausschneiden. Wählen Sie zunächst "Ebene" - "Ebene erstellen" - "Neue temporäre Kratzschicht" aus dem Menü und wählen Sie im Dialogfeld ein Polygon mit einem geeigneten Namen und Geometrietyp aus. Stellen Sie sicher, dass sich die neu hinzugefügte Ebene im Bearbeitungsmodus befindet, und wählen Sie im Menü "Bearbeiten" - "Rechteck hinzufügen" - "Bereichsbereichsrechteck hinzufügen". Klicken Sie mit der linken Maustaste, um die Eckpunkte diagonal zueinander auszuwählen, und klicken Sie mit der rechten Maustaste, um das Hinzufügen abzuschließen.
In der Abbildung ist die Zeichnung des Rechtecks zum leichteren Verständnis halbtransparent.
Wählen Sie als Nächstes "Vektor" - "Raumberechnungswerkzeug" - "Zuschneiden", wählen Sie die Ebene der Gebäudedaten für die Eingabe und die Ebene mit dem rechteckigen Bereich aus, den Sie gerade für die Überlagerung erstellt haben, und führen Sie sie aus.
Der Bereich, den Sie exportieren möchten, ist jetzt beschnitten.
Finden Sie übrigens den Schwerpunkt dieses rechteckigen Bereichs. "Vektor" - "Geometrie-Werkzeug" - "Schwerpunkt". Sie können den Breiten- und Längengrad des Schwerpunkts ermitteln, indem Sie mit dem Tool zum Anzeigen von Funktionsinformationen auf den erstellten Punkt klicken.
QGIS-Makro? Python-Konsole
Exportieren wir diese Daten nun mithilfe der Python-Konsole in eine 3D-OBJ-Datei. Die Python-Konsole kann über das Menü aufgerufen werden, indem Sie "Plugins" - "Python-Konsole" auswählen.
Transformation koordinieren
Derzeit arbeitet das Koordinatensystem mit WGS84 (wahrscheinlich, wenn sich nichts geändert hat). Wenn Sie anfangen, die Details des Koordinatensystems zu erklären, wird es eine enorme Menge sein, also werde ich es hier weglassen. In WGS84 werden die Koordinaten durch Breiten- und Längengrade dargestellt. Konvertieren Sie von hier aus in das dreidimensionale orthogonale XYZ-Koordinatensystem. Die folgende Methode ist eine Koordinatenkonvertierungsmethode, die unter Bezugnahme auf das Programm im Buch World Survey System and Coordinate Conversion erstellt wurde. Die orthogonale Koordinatenumwandlung ist gemäß der Darstellung der Erde in WGS84 so genau wie möglich. Beachten Sie, dass diese Berechnung mit mindestens doppelter Genauigkeit durchgeführt werden muss. Da der Gleitkommawert von Python doppelt sein sollte, müssen Sie sich nicht zu viele Sorgen machen. Da die Vector3-Klasse von Unity usw. numerische Werte in float enthält, müssen Sie beim Portieren dieses Programms vorsichtig sein. Wenn Sie häufig mit geografischen Koordinaten arbeiten, ist es besser, mit Doppelkoordinaten umzugehen, damit Sie keine zusätzlichen Fehler machen.
def BLH2XYZ(b,l,h):
# WGS84
a = 6378137.0
f = 1.0 / 298.257223563
b = math.radians(b)
l = math.radians(l)
e2 = f * (2.0 - f)
N = a / math.sqrt(1.0 - e2 * math.pow(math.sin(b), 2.0))
X = (N + h) * math.cos(b) * math.cos(l);
Y = (N + h) * math.cos(b) * math.sin(l);
Z = (N * (1 - e2) + h) * math.sin(b);
return (X,Y,Z)
OBJ-Export
Erstellen Sie dann ein OBJ-Exportprogramm, das die oben beschriebene Koordinatenkonvertierungsmethode enthält.
def BLH2XYZ(b,l,h):
# WGS84
a = 6378137.0
f = 1.0 / 298.257223563
b = math.radians(b)
l = math.radians(l)
e2 = f * (2.0 - f)
N = a / math.sqrt(1.0 - e2 * math.pow(math.sin(b), 2.0))
X = (N + h) * math.cos(b) * math.cos(l);
Y = (N + h) * math.cos(b) * math.sin(l);
Z = (N * (1 - e2) + h) * math.sin(b);
return (X,Y,Z)
layer = iface.activeLayer()
features = layer.getFeatures()
i = 0
j = 0
fp = open("test.obj",mode='w')
fp.write("g "+str(i)+"\n")
#* Wechseln Sie zum Schwerpunkt des zuvor erfassten rechteckigen Bereichs
cx = 139.72957
cy = 35.66021
ox,oy,oz = BLH2XYZ(cy,cx,0)
for feature in features:
# print(feature.id())
mp = feature.geometry().asMultiPolygon()
try:
height = int(feature['height'])*2
if height < 1:
height = 2
except:
height = 5
for polygon in mp:
for polyline in polygon:
i=i+1
prev_p_index = j
for point in polyline:
x,y,z = BLH2XYZ(point.y(),point.x(),0)
x2,y2,z2 = BLH2XYZ(point.y(),point.x(),height)
x = x - ox
y = y - oy
z = z - oz
x2 = x2 - ox
y2 = y2 - oy
z2 = z2 - oz
s = math.radians(-cx)
rx = x * math.cos(s) - y * math.sin(s)
ry = x * math.sin(s) + y * math.cos(s)
s = math.radians(-cy)
rxx = rx * math.cos(s) - z * math.sin(s)
rz = rx * math.sin(s) + z * math.cos(s)
s = math.radians(-cx)
rx2 = x2 * math.cos(s) - y2 * math.sin(s)
ry2 = x2 * math.sin(s) + y2 * math.cos(s)
s = math.radians(-cy)
rxx2 = rx2 * math.cos(s) - z2 * math.sin(s)
rz2 = rx2 * math.sin(s) + z2 * math.cos(s)
fp.write("v "+str(ry)+" "+str(rxx)+" "+str(rz)+"\n")
fp.write("v "+str(ry2)+" "+str(rxx2)+" "+str(rz2)+"\n")
j=j+1
current = j - prev_p_index
for num in range(current):
p1 = (2*num+1)
p2 = (2*num+2)
p3 = (2*num+3)
p4 = (2*num+4)
if p3 > current*2:
p3 = p3 - current*2
if p4 > current*2:
p4 = p4 - current*2
p1 = p1 + prev_p_index*2
p2 = p2 + prev_p_index*2
p3 = p3 + prev_p_index*2
p4 = p4 + prev_p_index*2
fp.write("f "+str(p1)+" "+str(p2)+" "+str(p3)+" "+str(p1)+"\n")
fp.write("f "+str(p4)+" "+str(p3)+" "+str(p2)+" "+str(p4)+"\n")
fp.close()
print("Done")
In diesem Programm ist die Verarbeitung die Koordinatenkonvertierung von den geografischen Koordinaten in die orthogonalen Koordinaten, die Bewegung des Ursprungs, die Drehung des Ganzen und das Schreiben in OBJ.
Ergebnis
Recommended Posts