地形図#

これは テクスチャを適用する の例と非常に似ていますが,地形メッシュの上にGeoTIFFからの航空画像をプロットすることに焦点を当てています.

import matplotlib as mpl
import matplotlib.pyplot as plt

import pyvista as pv
from pyvista import examples

まずは標高データと地形図を読み込むことから始めましょう.

# Load the elevation data as a surface
elevation = examples.download_crater_topo().warp_by_scalar()
# Load the topographic map from a GeoTiff
topo_map = examples.download_crater_imagery()
topo_map = topo_map.flip_y()  # flip to align to our dataset

elevation
HeaderData Arrays
StructuredGridInformation
N Cells1677401
N Points1680000
X Bounds1.810e+06, 1.831e+06
Y Bounds5.640e+06, 5.658e+06
Z Bounds7.339e+02, 2.787e+03
Dimensions1400, 1200, 1
N Arrays1
NameFieldTypeN CompMinMax
scalar1of1Pointsfloat6417.339e+022.787e+03


先ほど読み込んだイメージを点検してみましょう.

mpl.rcParams['figure.dpi'] = 500
plt.imshow(topo_map.to_array())
topo map
<matplotlib.image.AxesImage object at 0x7f2183346030>

サーフェスメッシュ(ここでは pyvista.StructuredGrid を使います)として地形メッシュをロードし, pyvista.read_texture() を使用して:class:pyvista.Texture としてイメージをロードしたら,次のようにイメージをサーフェスメッシュにマッピングできます.

# Bounds of the aerial imagery - given to us
bounds = (1818000, 1824500, 5645000, 5652500, 0, 3000)
# Clip the elevation dataset to the map's extent
local = elevation.clip_box(bounds, invert=False)
# Apply texturing coordinates to associate the image to the surface
local.texture_map_to_plane(use_bounds=True, inplace=True)
HeaderData Arrays
UnstructuredGridInformation
N Cells436733
N Points222110
X Bounds1.818e+06, 1.825e+06
Y Bounds5.645e+06, 5.653e+06
Z Bounds1.381e+03, 2.787e+03
N Arrays2
NameFieldTypeN CompMinMax
scalar1of1Pointsfloat6411.381e+032.787e+03
Texture CoordinatesPointsfloat3220.000e+001.000e+00


表示します.イメージが予想どおりに調整されていることに注意してください.

local.plot(texture=topo_map, cpos="xy")
topo map

これが3 D遠近法です.

local.plot(texture=topo_map)
topo map

また,周辺領域を抽出し,テクスチャマッピングされた局所的な地形と外部領域をプロットすることによって,領域全体を表示することもできます.

# Extract surrounding region from elevation data
surrounding = elevation.clip_box(bounds, invert=True)

# Display with a shading technique
p = pv.Plotter()
p.add_mesh(local, texture=topo_map)
p.add_mesh(surrounding, color="white")
p.enable_eye_dome_lighting()
p.camera_position = [
    (1831100.0, 5642142.0, 8168.0),
    (1820841.0, 5648745.0, 1104.0),
    (-0.435, 0.248, 0.865),
]
p.show()
topo map

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

Sphinx-Galleryによるギャラリー