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.

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.

20191203_fig_cliprect.PNG

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.

20191203_fig_cliped.PNG

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

20191203_fig_3d.PNG

Recommended Posts

Exporter des données 3D depuis QGIS
Exporter les données CASTable
Extraction de données depuis S3
Obtenir les données structurelles de CHEMBLID
De l'installation d'Elasticsearch à la saisie des données
Exporter les données DB au format json
Python: exclure les balises des données html
Frappez les données du trésor de Python Pandas
Obtenir des données de Quandl en Python
Exporter mp4 de maya en utilisant ffmpeg
Extraire des données spécifiques d'un JSON complexe
Structure de données persistante créée à partir de zéro
[Bases de la science des données] Acquisition de données à partir de l'API
Obtenez des données de Twitter avec Tweepy
[Python scipy] Augmentation / réduction de l'échelle des données 2D