[PYTHON] Exporter des données 3D depuis QGIS
introduction
Exportons les données SIG dans un fichier OBJ en utilisant la console Python intégrée à QGIS.
- Pour les fichiers OBJ, veuillez présenter le [lien Wikipédia](fichier https://ja.wikipedia.org/wiki/Wavefront_.obj). C'est un format très simple que vous pouvez exporter sous forme de texte, il est donc facile à gérer si vous souhaitez exporter le vôtre.
Coupé
Commencez par découper la plage cible. Si vous utilisez les données précédentes telles quelles, elles seront trop volumineuses et le rendu en temps réel ne pourra pas les gérer, je vais donc les découper dans une plage appropriée. Cette fois, je vais découper la zone autour de Roppongi Hills. Tout d'abord, sélectionnez "Calque" - "Créer un calque" - "Nouveau calque temporaire" dans le menu, et dans la boîte de dialogue, sélectionnez un polygone avec un nom et un type de géométrie appropriés. Assurez-vous que le nouveau calque que vous avez ajouté est en mode édition, puis sélectionnez "Edition" - "Ajouter un rectangle" - "Ajouter un rectangle de plage de zone" dans le menu. Cliquez avec le bouton gauche pour sélectionner un sommet et un sommet diagonal, et cliquez avec le bouton droit pour terminer l'ajout.
Sur la figure, le dessin du rectangle est rendu semi-transparent pour une compréhension facile.
Ensuite, sélectionnez "Vecteur" - "Outil de calcul spatial" - "Recadrer", sélectionnez le calque des données du bâtiment comme entrée, puis sélectionnez le calque avec la zone rectangulaire que vous venez de créer comme superposition, puis exécutez-le.
La zone que vous souhaitez exporter est maintenant rognée.
Au fait, trouvez le centre de gravité de cette zone rectangulaire. "Vecteur" - "Outil de géométrie" - "Centre de gravité". Vous pouvez trouver la latitude et la longitude du centre de gravité en cliquant sur le point créé avec l'outil d'affichage des informations sur les entités.
Macro QGIS? Console Python
Maintenant, exportons ces données dans un fichier OBJ 3D à l'aide de la console Python. La console Python peut être appelée à partir du menu en sélectionnant "Plugins" - "Python Console".
Transformation de coordonnées
Actuellement, le système de coordonnées fonctionne avec WGS84 (probablement si rien n'a changé). Si vous commencez à expliquer les détails du système de coordonnées, ce sera une quantité énorme, donc je vais l'omettre ici. Dans WGS84, les coordonnées sont représentées par la latitude et la longitude. À partir de là, passez au système de coordonnées orthogonales tridimensionnelles XYZ. La méthode suivante est une méthode de conversion de coordonnées créée en se référant au programme dans le livre World Survey System and Coordinate Conversion. La conversion des coordonnées orthogonales est aussi précise que possible selon la représentation de la terre en WGS84. Notez que ce calcul doit être effectué avec une précision au moins double. La virgule flottante de Python doit être double, vous n'avez donc pas à vous en soucier trop, mais la classe Vector3 d'Unity, etc. a des nombres en float, alors soyez prudent lorsque vous portez ce programme. Lorsqu'il s'agit souvent de coordonnées géographiques, il est préférable de traiter les doubles afin de ne pas faire d'erreurs inutiles.
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)
Export OBJ
Ensuite, créez un programme d'exportation OBJ qui incorpore la méthode de conversion de coordonnées ci-dessus.
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")
#* Changement du centre de gravité de la zone rectangulaire acquise plus tôt
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")
Dans ce programme, le traitement consiste à convertir les coordonnées géographiques en coordonnées orthogonales, à déplacer l'origine, à faire pivoter le tout et à écrire dans OBJ.
résultat
Recommended Posts