gemgis

GIS contacts and DEM elevations

Use GemGIS to extract geological observations from spatial data. Establish the horizontal CRS, vertical units/datum, geometry type and missing-data policy before sampling. Model coordinates must share a projected metre CRS; converting horizontal coordinates does not transform elevations or their datum.

Extract contact vertices

The following function accepts a single-band DEM already in a metre projected CRS, with elevations in metres. dem is an open Rasterio dataset, not a path. GemGIS repeats attributes when a line becomes multiple points: do not overwrite those attributes using the original feature index.

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']]

Open with with rasterio.open(dem_path) as dem: and call this function while the dataset is open. Pixel sampling uses the containing cell; it is not bilinear interpolation. Do not fill NoData or outside-extent samples with zero.