メッシュ後の地形#

地形サーフェスを使用して,3 D地形フォローメッシュを作成します.

メッシュに続く地形は,環境科学,例えば水文モデリング( Maxwell 2013ParFlow を参照)において一般的である.

この例では,特定の地形サーフェスに従う3 Dグリッド/メッシュを作成する簡単な方法を適用します.この例では,指定した数値標高モデル (DEM) が構造化されていることに注意してください(グリッド付きで3角形でない): これはDEMでは一般的です.

import numpy as np

import pyvista as pv
from pyvista import examples

グリッド地形サーフェス (DEM) をダウンロードします

dem = examples.download_crater_topo()
dem
HeaderData Arrays
ImageDataInformation
N Cells1677401
N Points1680000
X Bounds1.810e+06, 1.831e+06
Y Bounds5.640e+06, 5.658e+06
Z Bounds0.000e+00, 0.000e+00
Dimensions1400, 1200, 1
Spacing1.500e+01, 1.500e+01, 0.000e+00
N Arrays1
NameFieldTypeN CompMinMax
scalar1of1Pointsfloat6417.339e+022.787e+03


ここで,対象となる領域をサブサンプルして抽出し,この例を単純にしましょう(ロードしたばかりのDEMもかなり大きいです.).ロードしたDEMは pyvista.ImageData メッシュなので, pyvista.ImageData.extract_subset() filter: フィルタを使用できます.

subset = dem.extract_subset((500, 900, 400, 800, 0, 0), (5, 5, 1))
subset.plot(cpos="xy")
terrain mesh

メッシュに沿って地表の関心領域ができたので,DEMの3 Dサーフェスを作成します.

terrain = subset.warp_by_scalar()
terrain
HeaderData Arrays
StructuredGridInformation
N Cells6400
N Points6561
X Bounds1.818e+06, 1.824e+06
Y Bounds5.646e+06, 5.652e+06
Z Bounds1.441e+03, 2.769e+03
Dimensions81, 81, 1
N Arrays1
NameFieldTypeN CompMinMax
scalar1of1Pointsfloat6411.441e+032.769e+03


terrain.plot()
terrain mesh

これで,地形の3 D構造サーフェスが作成されました.構造化されたサーフェスを3 Dメッシュに拡張して,グリッドに沿った地形を形成できるようになりました.これを行うには,まずz方向(これらは地表から始まり)のセル間隔を設定します.次に,地表メッシュのXYZ構造化座標を繰り返し,各ZレベルをZセル間隔だけ下げます.これらの構造化された座標が得られたら, pyvista.StructuredGrid を作成できます.

z_cells = np.array([25] * 5 + [35] * 3 + [50] * 2 + [75, 100])

xx = np.repeat(terrain.x, len(z_cells), axis=-1)
yy = np.repeat(terrain.y, len(z_cells), axis=-1)
zz = np.repeat(terrain.z, len(z_cells), axis=-1) - np.cumsum(z_cells).reshape((1, 1, -1))

mesh = pv.StructuredGrid(xx, yy, zz)
mesh["Elevation"] = zz.ravel(order="F")
mesh
HeaderData Arrays
StructuredGridInformation
N Cells70400
N Points78732
X Bounds1.818e+06, 1.824e+06
Y Bounds5.646e+06, 5.652e+06
Z Bounds9.364e+02, 2.744e+03
Dimensions81, 81, 12
N Arrays1
NameFieldTypeN CompMinMax
ElevationPointsfloat6419.364e+022.744e+03


cpos = [
    (1826736.796308761, 5655837.275274233, 4676.8405505181745),
    (1821066.1790519988, 5649248.765538796, 943.0995128226014),
    (-0.2797856225380979, -0.27966946337594883, 0.9184252809434081),
]

mesh.plot(show_edges=True, lighting=False, cpos=cpos)
terrain mesh

Total running time of the script: (0 minutes 2.107 seconds)

Sphinx-Galleryによるギャラリー