gemgis

GIS の接触点と DEM の標高

GemGIS を使用して、空間データから地質観測データを抽出します。サンプリングの前に、水平 CRS、垂直単位／測地基準系、ジオメトリタイプ、欠損データの扱いを確定してください。モデルの座標は、同じ投影座標系（単位：メートル）を使用する必要があります。水平座標を変換しても、標高やその測地基準系は変換されません。

接触点の頂点を抽出する

次の関数は、メートル単位の投影 CRS を使用し、標高の単位がメートルである単バンド DEM を受け取ります。dem は開いている Rasterio データセットであり、パスではありません。GemGIS はラインが複数のポイントになる際に属性を複製します。元のフィーチャインデックスを使って、それらの属性を上書きしないでください。

import gemgis as gg
from rasterio.transform import rowcol
import numpy as np
from pyproj import CRS

def contacts_at_dem(contacts, dem):
    crs = CRS.from_user_input(dem.crs)
    if (contacts.crs is None or not crs.is_projected or dem.count != 1
            or any(not np.isclose(a.unit_conversion_factor, 1.) for a in crs.axis_info[:2])):
        raise ValueError('Require a known vector CRS and a single-band metric DEM')
    if contacts.empty or not contacts.is_valid.all() or contacts.geometry.has_z.any():
        raise ValueError('Require valid nonempty 2D contact geometry')
    if 'formation' not in contacts or contacts['formation'].isna().any():
        raise ValueError('Formation labels are required')
    points = gg.vector.extract_xy(contacts.to_crs(dem.crs))
    xy = points[['X', 'Y']].to_numpy()
    rows, cols = np.asarray(rowcol(dem.transform, xy[:, 0], xy[:, 1]))
    if np.any((rows < 0) | (rows >= dem.height) | (cols < 0) | (cols >= dem.width)):
        raise ValueError('Vertices outside DEM')
    samples = np.ma.vstack(list(dem.sample(xy, indexes=1, masked=True)))[:, 0]
    if np.ma.getmaskarray(samples).any() or not np.isfinite(samples.data).all():
        raise ValueError('Missing DEM elevations')
    return gg.vector.extract_xyz(gdf=points, dem=dem)[['X', 'Y', 'Z', 'formation']]

with rasterio.open(dem_path) as dem: を使用して開き、データセットを開いている間にこの関数を呼び出します。ピクセルのサンプリングでは、その地点を含むセルを使用します。バイリニア補間ではありません。NoData または範囲外のサンプルを 0 で埋めないでください。

方位とヘルパー

傾斜角は 0～90 度です。出力の方位角は、DEM グリッドの北から時計回りに測った傾斜方向で、[0, 360) の範囲に収められます。走向の変換に (strike + 90) % 360 を使うのは、入力データが右手則に従うことが明示されている場合に限ります。指定された極性（+1／−1）を保持してください。このヘルパーが +1 を使用するのは、列が存在しない場合のみです。無効な傾斜角をクランプしないでください。方位がすでに DEM グリッドを基準としている場合は --azimuth-reference dem-grid を、地理的な北を基準としている場合は true-north を指定します。後者では、局所的に等角な投影法の場合のみ、局所子午線収差が適用されます。それ以外の投影法は拒否されます。磁北を基準とする方位には、まず根拠が明記された偏角補正が必要です。

準備用ヘルパー はこれらのルールを検証し、ベクターデータを DEM の CRS に再投影し、単位を ft と宣言された標高をメートルに変換して、新しい出力ディレクトリに CSV と spatial_metadata.json を書き込みます。--target-crs は DEM と一致する必要があることを示すアサーションです。必要に応じて、ラスターを別途リサンプリング／再投影してください。

python scripts/prepare_gempy_data.py --contacts contacts.gpkg \
  --orientations orientations.gpkg --azimuth-reference true-north --dem dem.tif --dem-z-unit m \
  --vertical-datum "documented source datum" --output-dir prepared

このコマンドは、インストール済みのスキルディレクトリを基準に解決してください。出力スキーマでは X／Y／Z と formation を使用し、方位データではさらに dip／azimuth／polarity を使用します。これらを現在のモデリングライブラリの入力 API に渡してください。すべての GemPy バージョンで同じキーワード引数が使えると仮定しないでください。

条件付きの操作

サンプリング済みプロファイルと方位の規約については、プロファイルと抽出を参照してください。範囲の順序、マスク、再投影、来歴については、CRS、クリッピング、ラスターの取り扱いを参照してください。既存の 3D ジオメトリにこのヘルパーを使う前に、測定済みの Z と DEM の Z のどちらを使用するかを明示的に選択してください。

検証の範囲

GemGIS 1.1.9、GeoPandas 1.1.4、Rasterio 1.4.4 を、実際の GeoTIFF／GeoPackage フィクスチャで検証しました。検証項目は、頂点数の不一致と連続していないインデックス、CRS 変換、フィート単位の標高、NoData、範囲外のポイント、方位、CLI による CSV の再読み込みです。合成標高で検証できるのはデータの取り扱いであり、実際の DEM の精度や垂直測地基準系の変換ではありません。

GemGIS 公式ソース、 Rasterio のサンプリング API、 2026-09-14 確認。