QGISで3D地形を表示する

概要：
ざっくりとした手順と解説：

概要：

目的：QGISのpythonコンソールからDEMを書き換える(2022年8月)．

準備するもの：特になし

広告　【Udemy】オンラインコース：PythonによるGISデータの分析と可視化

ざっくりとした手順と解説：

前回までで3D表示はできたのですが，DEM欠測値の処理が甘かったために3Dモデルが陥没しています．これを修正します．

「プラグイン」→「Pythonコンソール」を立ち上げます．

Pythonコンソールが開いたら，「エディタの表示」→「新しいエディタ」をクリックしてカラの編集画面を開き，コードを貼り付けます．コードを貼り付けたら，修正するDEM（ここでは「merge1」）をアクティブレイヤにして（レイヤ領域でクリックします）コードをRUNします．

ここでは，①DEMの欠測値=0.0をwaterLebel=22.0[m]で置き換える，②置き換えられたDEMを平面直角座標系に再投影する，③heightmap（0～1に正規化された標高ラスタ）を作成する，の三つの処理を行っています．

import numpy as np
from osgeo import gdal
from osgeo import osr

# PARAMETERS ################################################
minH = 20.0 # [m]
waterLevel = 22.0 # [m]
maxH = 30.0 # [m]
layname = "marge1_filled"
layname1 = "marge1_filled_JGD2011"
#############################################################

def find_layer(name):
    lay = ''
    for a in QgsProject.instance().mapLayers().values():
        if (a.name()==name):
            lay = a
            break
    return lay
def load_tif(tifname, layname, crs0):
    iface.addRasterLayer(tifname, layname)
    lay = iface.activeLayer()
    print("Loadeding... ", lay.name())
    lay.setCrs(crs0)
    
# Opening the active layer
lay = iface.activeLayer()
print("Opening ", lay.name())
uri = lay.dataProvider().dataSourceUri()
src = gdal.Open(uri, gdal.GA_ReadOnly)
crs0 = lay.crs()
h0 = src.GetRasterBand(1) # for single band ONLY
h1 = h0.ReadAsArray()
width = h1.shape[1]
height = h1.shape[0]
dim = h1.ndim
size = h1.size
print('width=',width,'height=',height)
geotransform = src.GetGeoTransform()
print("geotransform ", geotransform)
originY = geotransform[3]
originX = geotransform[0]
pixelW = geotransform[1]
pixelH = geotransform[5]

# fill nodata with "waterlevel"
h1 = np.where(h1<minH, waterLevel, h1)
h1 = np.where(h1>maxH, maxH, h1)

# save
dst = uri[:-4]+"_filled.tif"
print("saving", dst)
driver = gdal.GetDriverByName('GTiff')
width = int(width)
height = int(height)
dst_raster = driver.Create(dst, width, height, 1, gdal.GDT_Float32)
dst_raster.SetGeoTransform( (originX, pixelW, 0.0, originY, 0.0, pixelH) )
dst_band = dst_raster.GetRasterBand(1)
dst_band.WriteArray(h1)
dst_band.FlushCache()
dst_raster = None

# load to the project
load_tif(dst, layname, crs0)

# warp
dst1 = dst[:-4]+"_JGD2011.tif"
print("saving ", dst1)
rlayer = QgsProject.instance().mapLayersByName(layname)[0]
processing.run("gdal:warpreproject", 
                    {"INPUT":dst,
                    "SOURCE_CRS":"EPSG:4612",
                    "TARGET_CRS":"EPSG:6677",
                    "RESAMPLING":1,
                    "TARGET_RESOLUTION":5.0,
                    'OUTPUT': dst1
                    })

# load to the project
crs0 = QgsCoordinateReferenceSystem("EPSG:6677")
load_tif(dst1, layname1, crs0)
lay = iface.activeLayer()
uri = lay.dataProvider().dataSourceUri()
src = gdal.Open(uri, gdal.GA_ReadOnly)
h0 = src.GetRasterBand(1) # for single band ONLY
h1 = h0.ReadAsArray()
width = h1.shape[1]
height = h1.shape[0]
print('width=',width,'height=',height)
geotransform = src.GetGeoTransform()
print("geotransform ", geotransform)

# heightmap
h1 = np.where(h1<minH, waterLevel, h1) # fill "nodata"s generated by WARP
h2 = (h1 - minH)/(maxH-minH)
print("hmin, hmax = ",h2.min(),h2.max())
dst = uri[:-4]+"_heightmap.tif"
print("saving ", dst)
driver = gdal.GetDriverByName('GTiff')
width = int(width)
height = int(height)
dst_raster = driver.Create(dst, width, height, 1, gdal.GDT_Float32)
dst_raster.SetGeoTransform( (originX, pixelW, 0.0, originY, 0.0, pixelH) )
dst_band = dst_raster.GetRasterBand(1)
dst_band.WriteArray(h2)
dst_band.FlushCache()
dst_raster = None

print('finished')

5～11行：入力部分．waterLevelに適切な値（この場合は22[m]）を指定します．minH,maxHはきりのよい数字にします．laynameは欠測値を書き換えたDEMレイヤの名前，layname1は平面直角座標に再投影されたDEMレイヤの名前．
27行：アクティブレイヤを探しています．RUNの前に「merge1」をアクティブにしておく必要があります．
47-48行：DEMを処理しているのはたったの2行，その他は入出力操作にすぎません．

正常にRUNが終了すると，ふたつのラスタデータがプロジェクトに追加されています．申し遅れましたが，ここまでラスタレイヤはすべて緯度経度で管理されていました．ここで平面直角座標系に投影しなおしたのは，Unity，Three.jsやCADなどで扱いやすくするためです．0～1に正規化したheightmapを生成しておくのも同様の理由からです．