地形上の地質図

地形面上のGeoTIFFのテクスチャマッピング．

GeoTIFFから地形に画像や地図を重ね合わせるためには，テクスチャ（GeoTIFF）を貼り付けたいメッシュや地形の延長線上にテクスチャ座標（「テクスチャマッピング」）が一致する必要があるのです．

これは，GeoTIFF 自体の空間リファレンスを使用することで，イメージのない部分を切り取ることなく，テクスチャがドレープされているメッシュ全体を保持することができるからです．この例では，テクスチャが地形サーフェス上に伸びるように明示的に設定し，テクスチャ/GeoTIFF が地形サーフェスよりはるかに大きな範囲を持つようにしています．

原文はこちら： pyvista/pyvista-support#14

import os
import tempfile

import numpy as np
import pyvista as pv
import requests
from pyvista import examples

path = examples.download_file("topo_clean.vtk")
topo = pv.read(path)
topo

UnstructuredGrid	Information
N Cells	824278
N Points	413250
X Bounds	3.299e+05, 3.442e+05
Y Bounds	4.253e+06, 4.271e+06
Z Bounds	1.494e+03, 2.723e+03
N Arrays	0

GeoTIFF/テクスチャを読み込みます(ダウンロードに1分程度かかる場合があります) https://dl.dropbox.com/s/bp9j3fl3wbi0fld/downsampled_Geologic_map_on_air_photo.tif?dl=0

url = "https://dl.dropbox.com/s/bp9j3fl3wbi0fld/downsampled_Geologic_map_on_air_photo.tif?dl=0"

response = requests.get(url)  # noqa: S113
filename = os.path.join(tempfile.gettempdir(), "downsampled_Geologic_map_on_air_photo.tif")  # noqa: PTH118
open(filename, "wb").write(response.content)  # noqa: SIM115, PTH123

8175934

以下のブロックでは， get_gcps 関数を使用してラスターの Ground Control Point を取得することができますが，これは GDAL に依存します．このチュートリアルでは，ユーザーがGDALをインストールするのを避けるために，GCPをハードコードするつもりです．

def get_gcps(filename):
    """
    Helper function retrieves the Ground Control
    Points of a GeoTIFF. Note that this requires gdal.
    """
    import rasterio

    def get_point(gcp):
        return np.array([gcp.x, gcp.y, gcp.z])

    # Load a raster
    src = rasterio.open(filename)
    # Grab the Groung Control Points
    points = np.array([get_point(gcp) for gcp in src.gcps[0]])
    # Now Grab the three corners of their bounding box
    # -- This guarantees we grab the right points
    bounds = pv.PolyData(points).bounds
    origin = [bounds[0], bounds[2], bounds[4]]  # BOTTOM LEFT CORNER
    point_u = [bounds[1], bounds[2], bounds[4]]  # BOTTOM RIGHT CORNER
    point_v = [bounds[0], bounds[3], bounds[4]]  # TOP LEFT CORNER
    return origin, point_u, point_v

# Fetch the GCPs
# origin, point_u, point_v = get_gcps(filename)

# Hard code GCPs
origin = [310967.75148705335, 4238841.045453942, 0.0]
point_u = [358682.9364281533, 4238841.045453942, 0.0]
point_v = [310967.75148705335, 4276281.98755258, 0.0]

# Use the GCPs to map the texture coordinates onto the topography surface
topo.texture_map_to_plane(origin, point_u, point_v, inplace=True)

Header	Data Arrays
UnstructuredGrid Information N Cells 824278 N Points 413250 X Bounds 3.299e+05, 3.442e+05 Y Bounds 4.253e+06, 4.271e+06 Z Bounds 1.494e+03, 2.723e+03 N Arrays 1	Name Field Type N Comp Min Max Texture Coordinates Points float32 2 3.737e-01 8.576e-01

GCPをトポサーフェスに関連づけ，テクスチャ座標を表示する

p = pv.Plotter()
p.add_point_labels(
    np.array(
        [
            origin,
            point_u,
            point_v,
        ]
    ),
    ["Origin", "Point U", "Point V"],
    point_size=5,
)

p.add_mesh(topo)
p.show(cpos="xy")

スタティックなシーン

インタラクティブなシーン

GeoTIFFをPyVistaで Texture として読み込みます:

texture = pv.read_texture(filename)

# Now plot the topo surface with the texture draped over it
# And make window size large for a high-res screenshot
p = pv.Plotter(window_size=np.array([1024, 768]) * 3)
p.add_mesh(topo, texture=texture)
p.camera_position = [
    (337461.4124956896, 4257141.430658634, 2738.4956020899253),
    (339000.40935731295, 4260394.940646875, 1724.0720826501868),
    (0.10526647627366331, 0.2502863297360612, 0.962432190920575),
]
p.show()

スタティックなシーン

インタラクティブなシーン