GitHub - da0bi/pysole: Field-data-based, physically-informed bedrock interpolation & 3D ray-based migration for geophysical signal traveltimes. PySole transforms limited survey points into robust, physically-constrained 3D bedrock models.

Physically-Informed Bedrock Interpolation & 3D Ray-Based Migration for Sparse Geophysical Datasets

PySole is designed to reconstruct the thickness distribution and basal topography (sole) of glaciers, landslides, and other gravity-driven, viscous flow phenomena. The highly simplified Shallow Ice Approximation describes gravity-driven, viscous flow by its fundamental inverse relation between depth of the creeping body and its surface slope, where greater depths correspond to gentler surface slopes and vice versa. PySole leverages real field data to adapt this basic viscous flow characteristic into simple, but powerful thickness predictors. It is specifically engineered for sparse geophysical datasets where 3D wavefield migration to image the bedrock is impossible due to insufficient spatial sampling. PySole transforms limited survey points into robust, physically-constrained 3D bedrock models.

Key Features

Baseline Workflow & Methodology

Installation

JSON Configuration File

JSON Configuration Parameter Reference

Technical & Methodological Notes
    1. Supported DEM Formats
    2. Parsing of Rock Outcrop Input Files
    3. Parsing of Survey Data CSV Files
    4. Wavenumbers and Wavelengths
    5. Variogram Binning with Minimum Pair Threshold
    6. High-Performance Dual Kriging Vector Engine
    7. Universal Kriging Drift Models
    8. Kriging Uncertainty Derivation
    9. Memory Guard and Kriging Uncertainty
    10. Spatial Smoothing of the Calculated DEMs

Package Architecture

Command Line Interface

Application Programming Interface

Real-World Example: Wurtenkees Glacier

Citation & References

Key Features

JSON Configuration & Terminal CLI Driven: All processing workflow options can be defined in a pysole.json file and executed via Python Application Programming Interface (API) or directly from the terminal using the pysole command line interface (CLI).
Multi-Core Parallel Acceleration: PySole supports multi-core CPU parallelization across computationally intensive processing steps.
Automated High-Resolution Diagnostic Plots: Automatically generates and exports diagnostic figures for each key processing milestone.
5 Supported Digital Elevation Model (DEM) Formats: Seamless loading and exporting of GeoTIFFs (.tif), ESRI ASCII Grids (.asc, .txt), 3-column or 2D CSV grids (.csv), NumPy binary arrays (.npy), and in-memory 2D arrays (np.ndarray) — see Supported DEM Formats.
Strict Coordinate Reference System & Spatial Alignment Verification: Verifies CRS and spatial bounds across all input layers (DEM, boundary outline, survey points, and optional mass balance / velocity rasters), halting with an explicit error on misalignment.
Flexible Survey Data Types: Accepts one- or two-way signal traveltimes as well as direct thickness/depth measurements (automatically skipping 3D migration when depth data are provided).
💡DEM Surface Slope Smoothing💡: The degree of DEM surface slope smoothing is crucial when estimating ice thickness with the Shallow Ice Approximation (SIA), which assumes a constant basal shear stress. By relaxing this rigid baseline constraint, Binder et al. (2009) derived an objective, field-data-driven optimization criterion for the surface slope smoothing process, which is implemented in PySole. The optimal degree of surface slope smoothing is derived by enforcing minimum spatial variance in basal shear stress as the optimization criterion:

$$\min_{k_{\text{c}}} \text{Var}_{xy}(\tau_{\text{b}})$$

In shallow ice dynamics, basal shear stress is given by:

$$\tau_{\text{b}} = \rho_{\text{ice}} \, g \, D \sin(\alpha)$$

where ice density, ρ_ice, and gravitational acceleration, g, are assumed to be constant. Thus, only the product of ice depth and surface slope, P = D sin(α), is evaluated during the optimization process. The surface DEM is smoothed by Gaussian low-pass filtering in the 2D Fast Fourier Transform (FFT) frequency domain, defined by the spatial cutoff wavenumber k_c [rad/m] or cutoff wavelength λ_c [m] — see Wavenumbers and Wavelengths. The spatial variance of τ_b is quantified via variogram analysis, with an optional interactive mode to refine spectral cutoffs and variogram parameters. The resulting optimal surface slope field directly informs pre- and post-migration interpolation, 3D ray-based migration, surface curvature estimation, and flow velocity filtering. Both the FFT low-pass filtering and product variogram evaluations are multi-threaded across CPU cores.
3D Ray-Based Migration: Optional 3D ray-based migration introduced by Binder et al. (2009) for sparse geophysical signal traveltimes, using the optimally smoothed surface slope field for numerically stable ray displacement vectors.
Kriging Interpolation: Native, parallelized 2D Dual Kriging engine supporting Ordinary and Universal Kriging (plus Regression Kriging via PyKrige) — see Dual Kriging Vector Engine and Baseline Workflow & Methodology. Optional Kriging estimation uncertainty fields are evaluated alongside predicted grids when calculate_uncertainty: true — see Kriging Uncertainty Derivation.
Interactive Drift Analyzer: Automated diagnostic engine that evaluates and ranks up to 120 candidate models built from 8 physical and spatial drift primitives ("sia", "z_dem", "curvature_dem", "mass_balance", "velocity", "velocity_sia", "linear_xy", "quadratic_xy") plus Ordinary Kriging, with support for multi-model comparison suites, per-model diagnostic plots, and cross-model comparison tables — see Drift Analyzer Guide.
ML Hole Filling & Geomorphological Margin Blending: Optional scikit-learn Random Forest regression gap filling and geomorphological margin blending to smoothly taper bedrock elevations into surrounding surface DEM terrain.
Final DEMs Spatial Smoothing: Optional post-processing spatial smoothing ("gaussian", "median", "fft_lowpass") for the calculated thickness and bedrock DEMs — see Spatial Smoothing of the Calculated DEMs.
Survey Planner: Automatically generates physically-informed longitudinal flowline and transverse cross-profile survey layouts for unprobed glaciers based on a synthetic SIA thickness model, exporting field-ready GPX and GeoJSON tracks — see Survey Planner Guide.

Baseline Workflow & Methodology

Figure 1: End-to-end baseline workflow of the PySole solver dual-path processing pipeline. Click diagram to view in high resolution.

1. DEM Loading & Metadata

Grid spacing (dx, dy) and spatial bounds are automatically extracted from GeoTIFF and ESRI ASCII DEM metadata headers (and autodetected for 3-column .csv tables). If target dx and dy pixel resolutions are specified, 2D bilinear grid resampling is performed automatically. For headerless DEM formats (2D .csv matrices, .npy, or np.ndarray), spatial parameters (origin, crs, dx, dy) are provided under spatial_parameters — see Supported DEM Formats.

Caution

All input datasets (DEM, survey data, outline geometries, and optional mass balance / velocity rasters) must share a projected, metric coordinate reference system ("crs", e.g. UTM). Geographic coordinates (latitude/longitude in degrees) will cause invalid distance, surface slope, and basal shear stress calculations.

2. Pre-Migration Traveltime Interpolation

Across spatial wavenumber cutoffs k_c, the point products of traveltime observations and corresponding low-pass filtered surface slopes, P_T,i = T_i sin(α_smoothed,i), are evaluated. Once the minimum-variance optimization criterion is satisfied, the optimally smoothed surface slope field, sin(α_opt(x,y)) (clamped at slope_floor_deg, default 5°, to prevent division singularities), is deployed in the subsequent Kriging interpolation:

BSS-derived Product "P" (Default Baseline): Interpolates P_T,i using Ordinary Kriging (constant local mean without external drift) with optional zero-traveltime boundary conditions (T = 0 s), then reconstructs the continuous traveltime field T(x,y) via slope division:

$$T(x,y) = \frac{P_{\text{T}}(x,y)}{\sin(\alpha_{\text{opt}}(x,y))}$$
Direct "T" Baseline: Directly interpolates signal traveltimes T_i using Universal Kriging with the physically-informed ["sia"] drift model by default (or any candidate selected via the Drift Analyzer).

Caution

Avoid interpolating a BSS-derived product target ("P") with Universal Kriging using ["sia"] or ["velocity_sia"]! Dividing P by sin(α_opt) after fitting an inverse-slope drift introduces an unphysical 1/sin²(α) double-scaling artifact at low slopes. Use "sia" and "velocity_sia" only with direct targets ("T" or "D").

3. 3D Ray-Based Migration

Solves the Eikonal equation to relocate subsurface reflection points along 3D ray paths, improving the imaging of steep valley flanks and basal overdeepenings. An optional interactive mode allows testing different signal propagation velocities and inspecting horizontal survey point displacements.

4. Post-Migration Surface Slope Optimization

Re-evaluates the spectral optimization criterion on the migrated depths, P_D,i = D_i sin(α_smoothed,i), across wavenumber cutoffs k_c to determine the post-migration optimal surface slope field sin(α_opt(x,y)).

5. Final Depth Interpolation & Uncertainty Display

Reconstructs continuous ice thickness D(x,y) and bedrock elevation Z_bed(x,y) using either:

BSS-derived Product "P" (Default Baseline): Interpolates P_D,i = D_i sin(α_opt,i) via Ordinary Kriging (with optional D = 0 m boundary conditions) and divides by sin(α_opt(x,y)):

$$D(x,y) = \frac{P_{\text{D}}(x,y)}{\sin(\alpha_{\text{opt}}(x,y))}$$
Direct "D" Baseline: Directly interpolates (migrated) depths D_i using Universal Kriging (defaulting to ["sia"], or evaluated across up to 120 drift models with the Drift Analyzer).

When calculate_uncertainty: true, Kriging estimation variance is converted to standard error (σ_D = √σ²) to map spatial thickness uncertainty and total volume uncertainty — see Kriging Uncertainty Derivation.

Installation

Direct Installation from GitHub (Recommended)

You can install PySole directly from GitHub using pip:

pip install git+https://github.com/da0bi/pysole.git

PyPI Installation

🚨 NOT AVAILABLE YET 🚨

pip install pysole

Local / Development Installation

To clone and install PySole directly from source in editable mode:

git clone https://github.com/da0bi/pysole.git
cd pysole
pip install -e .

JSON Configuration File

All execution options can be fully defined in a single JSON configuration file, which by default is named pysole.json. Whenever an optional parameter key or entire section is omitted from the JSON file, PySole automatically populates it with its default value and continues without error. This is, however, not true for the survey_data_path key - if this is left empty, or totally omitted, the Survey Planner is automatically started - see here for details. Unknown keys in the configuration file never stop a PySole run - they are reported in a single warning and ignored.

{
    "inputs": {
        "dem_path": null,
        "outline_path": null,
        "survey_data_path": null,
        "survey_data_type": "one_way_traveltime",
        "survey_profile_column": "profile_id",
        "mass_balance_path": null,
        "velocity_path": null,
        "ice_density": 900.0,
        "g": 9.81,
        "n_cores": -1,
        "log_level": "INFO",
        "show_progress": true
    },
    "spatial_parameters": {
        "dx": null,
        "dy": null,
        "bounds": null,
        "origin": null,
        "crs": null
    },
    "migration_parameters": {
        "perform_migration": true,
        "velocity": 0.16,
        "interactive_migration": false,
        "horizontal_displacement_visualization": "arrows"
    },
    "optimization_parameters": {
        "fft_filter_metric": "wavenumber",
        "kc_max": null,
        "kc_min": null,
        "lambda_min": null,
        "lambda_max": null,
        "n_steps": null,
        "nrbins": null,
        "slope_floor_deg": 5.0,
        "interactive_optimization": false
    },
    "kriging_parameters": {
        "engine": "native",
        "max_memory_fraction": 0.5,
        "pre_migration": {
            "interpolation_target": "P",
            "method": "ordinary",
            "drift_analyzer": false,
            "drift_terms": null,
            "variogram_model": "spherical",
            "include_zero_boundary_condition": true,
            "calculate_uncertainty": false
        },
        "post_migration": {
            "interpolation_target": "P",
            "method": "ordinary",
            "drift_analyzer": false,
            "drift_terms": null,
            "variogram_model": "spherical",
            "include_zero_boundary_condition": true,
            "calculate_uncertainty": false
        }
    },
    "finalization_parameters": {
        "interactive_finalization": false,
        "random_forest_gap_filling": false,
        "apply_margin_blend": false,
        "min_gap_dist": 50.0,
        "smooth_bedrock": false,
        "smoothing_method": "gaussian",
        "smoothing_sigma": 1.5,
        "smoothing_kernel_size": 3,
        "smoothing_kc_cutoff": null
    },
    "outputs": {
        "output_dir": null,
        "output_format": "tif",
        "output_prefix": "",
        "plots_dir": "figures",
        "save_bedrock_grid": true,
        "save_traveltime_grid": false,
        "save_traveltime_uncertainty": false,
        "save_migrated_points": false,
        "save_thickness_grid": false,
        "save_thickness_uncertainty": false,
        "save_basal_shear_stress": false,
        "save_basal_shear_stress_uncertainty": false
    }
}

JSON Configuration Parameter Reference

Section	Parameter	Type	Default	Description
`inputs`	`dem_path`	`str`	`null`	(Required) File path to the surface Digital Elevation Model (`.asc`, `.tif`, `.csv`, `.npy`).
	`outline_path`	`str`	`null`	File path to creeping body / glacier boundary polygon (`.shp`, `.geojson`, `.gpkg`, `.csv`). If `null`, domain is derived from non-NaN DEM pixels.
	`survey_data_path`	`str`	`null`	File path to signal traveltime or thickness observations CSV `[ X, Y, (Z), survey_value, (profile_id)]`. Columns in brackets are optional - see here for CSV parsing details. Column `profile_id` (default `inputs.survey_profile_column`) is considered for the Drift Analyzer's Leave-One-Profile-Out (LOPO) cross-validation. See Drift Analyzer Guide for full details. If `null`, `PySole` automatically starts the interactive `Survey Planner` and saves the results to the file paths defined in the `outputs` section. See `here` for full details.
	`survey_data_type`	`str`	`"one_way_traveltime"`	Observation data type: `"one_way_traveltime"`, `"two_way_traveltime"`, or `"depth"` (skips 3D migration).
	`survey_profile_column`	`str`	`"profile_id"`	Column name in the imported survey dataset CSV (`inputs.survey_data_path`) specifying survey line / profile IDs for the Drift Analyzer's Leave-One-Profile-Out (LOPO) cross-validation (default: `"profile_id"`). Set to `null` to disable profile extraction when using headerless CSVs.
	`mass_balance_path`	`str`	`null`	Optional file path to a mean surface mass balance raster (`.tif`, `.asc`, `.csv`, `.npy`). Automatically co-registered to the DEM grid (kept unsmoothed to preserve local debris/avalanche features). Enables the `"mass_balance"` drift term and serves as an additional feature in Random Forest gap filling.
	`velocity_path`	`str`	`null`	Optional file path to a mean surface flow velocity magnitude raster (`.tif`, `.asc`, `.csv`, `.npy`). Automatically co-registered to the DEM grid, zero-clamped outside the glacier outline, and spectrally low-pass filtered at optimal k_c. Enables the `"velocity"` and `"velocity_sia"` drift terms and serves as an additional feature in Random Forest gap filling.
	`ice_density`	`float`	`900.0`	Density of the creeping medium in kg/m³ (`900.0` kg/m³ for temperate glacier ice by default). Used to calculate basal shear stress τ_b.
	`g`	`float`	`9.81`	Gravitational acceleration constant in m/s² (`9.81` m/s²). Used to calculate basal shear stress τ_b.
	`n_cores`	`int`	`-1`	Number of CPU cores applied across all parallelized processes (`-1` for all available cores).
	`log_level`	`str`	`"INFO"`	Package logging level verbosity: `"INFO"` (default), `"DEBUG"`, `"WARNING"`, `"ERROR"`, or `"CRITICAL"`. Appends timestamped logs to `pysole.log`.
	`show_progress`	`bool`	`true`	If `true` (default), displays terminal progress bars during heavy processing steps. Set automatically to `false` if running `PySole` in batch processing scripts.
`spatial_parameters`	`dx`	`float`	`null`	Target grid resolution along X in meters. If defined, automatically resamples the DEM grid. If `null`, native resolution is kept.
	`dy`	`float`	`null`	Target grid resolution along Y in meters. If defined, automatically resamples the DEM grid. If `null`, native resolution is kept.
	`bounds`	`list[float]`	`null`	Optional spatial bounding box `[minx, miny, maxx, maxy]`. Leave `null` by default. Use only to manually override invalid or missing spatial bounds in DEM raster headers (`.asc`, `.tif`). `bounds` are automatically calculated from `origin`, `dx`, `dy`, and grid dimensions for headerless DEMs (`.csv`, `.npy`).
	`origin`	`list[float]`	`null`	Lower-left coordinate origin `[xll, yll]` in projected metric units for headerless DEM formats (`.csv`, `.npy`, `np.ndarray`). If `null`, defaults to `(0.0, 0.0)`.
	`crs`	`str` / `int`	`null`	Coordinate Reference System (optional, e.g. `"EPSG:32633"`, `32633`, or PROJ string). When provided, validates metric projection, checks spatial alignment with vector outlines, and embeds EPSG metadata into output GeoTIFFs.
`migration_parameters`	`perform_migration`	`bool`	`true`	If `true`, performs 3D ray-based migration on signal traveltimes. If `false`, migration is skipped.
	`velocity`	`float`	`0.16`	Signal propagation velocity (default value of `0.16` m/ns is characteristic for radar wave propagation in temperate ice).
	`interactive_migration`	`bool`	`false`	If `true`, enables interactive velocity testing with visual migrated depths and horizontal displacement vector plots.
	`horizontal_displacement_visualization`	`str`	`"arrows"`	Horizontal displacement visualization mode: `"arrows"` (default, quiver displacement vectors overlaid on traveltime grid) or `"grad_traveltime_grid"` (displays traveltime field without arrows, and creates a companion figure `##_<output_prefix>_3d_migration_hor_displacements.png` showing traveltime gradient grid and original and migrated survey points).
`optimization_parameters`	`fft_filter_metric`	`str`	`"wavenumber"`	Frequency spectrum filter metric: `"wavenumber"` [rad/m] or `"wavelength"` [m].
	`kc_max`	`float`	`null`	Maximum corner frequency cutoff k_c,max in [rad/m]. If `null`, defaults to grid Nyquist wavenumber k_Nyquist. For full details, see here.
	`kc_min`	`float`	`null`	Minimum corner frequency cutoff k_c,min in [rad/m]. If `null`, defaults to half of the DEM extent. For full details, see here.
	`lambda_min`	`float`	`null`	Minimum spatial cutoff wavelength λ_min in [m]. If `null`, falls back directly to `kc_max`. For full details, see here.
	`lambda_max`	`float`	`null`	Maximum spatial cutoff wavelength λ_max in [m]. If `null`, falls back directly to `kc_min`. For full details, see here.
	`n_steps`	`int`	`null`	Number of evaluation steps across k_c spectrum using logarithmic/geometric spacing (`np.geomspace`). If `null`, dynamically calculated from DEM grid dimensions and constrained by [10, 50]. For full details, see here.
	`nrbins`	`int`	`null`	Number of variogram lag distance bins. If `null` (default), dynamically calculated to receive ~30 point pairs per bin. The number of bins is constrained by a maximum of 30 and an absolute minimum of 3 bins. For full details, see here.
	`slope_floor_deg`	`float`	`5.0`	Minimum surface slope angle threshold in degrees [°] enforced during surface slope optimization to prevent numerical division singularities.
	`interactive_optimization`	`bool`	`false`	If `true`, enables interactive CLI prompt to inspect BSS variance curve and adjust corner frequency spectrum parameters (`kc_min`, `kc_max`, `lambda_min`, `lambda_max`, `n_steps`), lag distance bin count (`nrbins`), and correlation range (`a_range`).
`kriging_parameters`	`engine`	`str`	`"native"`	Kriging calculation engine: `"native"` (default, high-performance Dual Kriging solver) or `"pykrige"` (uses external `PyKrige` package - optional dependency in `pyproject.toml`). Both engines support `"ordinary"` and `"universal"` with every drift term, however, it is recommended to use the numerically optimized `"native"` engine. `"regression"` is only available through `PyKrige` (if `engine` is `"native"` automatically switches to `"pykrige"` and posts a log message). See the engine support matrix below.
	`max_memory_fraction`	`float`	`0.5`	Fraction of the available RAM (`0 < f ≤ 0.9`) that the native Kriging solver may use. If the estimate exceeds it, the worker threads are reduced first; the Kriging variance is skipped only as a last resort. Also available as `Solver(max_memory_fraction=...)`. See here for details.
	`pre_migration`	`dict`	Sub-section	Configuration for pre-migration traveltime field, T(x,y), interpolation.
	`pre_migration.interpolation_target`	`str`	`"P"`	Pre-migration interpolation target: `"P"` for BSS-derived products (P = T_i · sin α_{opt, i}, default) or `"T"` for direct signal traveltimes (T_i).
	`pre_migration.method`	`str`	by target	Kriging approach: `"ordinary"` (constant mean, default for `"P"`), `"universal"` (external drift, default for `"T"`), or `"regression"` (only available through `PyKrige`). `"ordinary"` overrides `pre_migration.drift_analyzer` and `pre_migration.drift_terms` (both are ignored with a warning). See How `drift_terms` is interpreted.
	`pre_migration.drift_analyzer`	`bool`	`false`	If `true`, the Drift Analyzer cross-validates candidate drift models for `method: "universal"` and prompts for selection (batch mode selects the top-ranked model). When `drift_terms` is `null`, evaluates the dynamic default suite (24, 48, 60, or 120 candidates depending on `mass_balance_path` and `velocity_path`). Ignored for `"ordinary"`, skipped with a warning for `"regression"`, and auto-disabled if `drift_terms` defines a single candidate.
	`pre_migration.drift_terms`	`list` / `null`	`null`	`null` or a flat/nested list of primitives (`"sia"`, `"z_dem"`, `"curvature_dem"`, `"mass_balance"`, `"velocity"`, `"velocity_sia"`, `"linear_xy"`, `"quadratic_xy"`). See Universal Kriging Drift Models and How `drift_terms` is interpreted.
	`pre_migration.variogram_model`	`str`	`"spherical"`	Theoretical variogram model (`"spherical"`, `"exponential"`, `"gaussian"`, `"linear"`). `"linear"` is only supported under `engine: "pykrige"` (falls back to `"spherical"` under `"native"`).
	`pre_migration.include_zero_boundary_condition`	`bool`	`true`	If `true` (default), includes zero traveltime boundary points (T = 0 s) along the perimeter and rock outcrop margin outline(s).
	`pre_migration.calculate_uncertainty`	`bool`	`false`	If `true`, computes the Kriging estimation variance during pre-migration Kriging. Defaults to `false` to save computational time.
	`post_migration`	`dict`	Sub-section	Configuration for final bedrock depth, D(x,y), interpolation.
	`post_migration.interpolation_target`	`str`	`"P"`	Post-migration interpolation target: `"P"` for BSS-derived products (P = D_i · sin α_{opt, i}, default) or `"D"` for direct (migrated) depths (D_i).
	`post_migration.method`	`str`	by target	Kriging approach: `"ordinary"` (constant mean, default for `"P"`), `"universal"` (external drift, default for `"D"`), or `"regression"` (only available through `PyKrige`). `"ordinary"` overrides `post_migration.drift_analyzer` and `post_migration.drift_terms` (both are ignored with a warning). See How `drift_terms` is interpreted.
	`post_migration.drift_analyzer`	`bool`	`false`	If `true`, the Drift Analyzer cross-validates candidate drift models for `method: "universal"` and prompts for selection (batch mode selects the top-ranked model). When `drift_terms` is `null`, evaluates the dynamic default suite (24, 48, 60, or 120 candidates depending on `mass_balance_path` and `velocity_path`). Ignored for `"ordinary"`, skipped with a warning for `"regression"`, and auto-disabled if `drift_terms` defines a single candidate.
	`post_migration.drift_terms`	`list` / `null`	`null`	`null` or a flat/nested list of primitives (`"sia"`, `"z_dem"`, `"curvature_dem"`, `"mass_balance"`, `"velocity"`, `"velocity_sia"`, `"linear_xy"`, `"quadratic_xy"`). See Universal Kriging Drift Models and How `drift_terms` is interpreted.
	`post_migration.variogram_model`	`str`	`"spherical"`	Theoretical variogram model (`"spherical"`, `"exponential"`, `"gaussian"`, `"linear"`). `"linear"` is only supported under `engine: "pykrige"` (falls back to `"spherical"` under `"native"`).
	`post_migration.include_zero_boundary_condition`	`bool`	`true`	If `true` (default), includes zero thickness boundary points (D = 0 m) along the perimeter and rock outcrop margin outline(s).
	`post_migration.calculate_uncertainty`	`bool`	`false`	If `true`, computes the Kriging estimation variance during post-migration Kriging. Defaults to `false` to save computational time.
`finalization_parameters`	`interactive_finalization`	`bool`	`false`	If `true`, enables interactive terminal prompts to inspect and adjust finalization parameters (smoothing, RF gap filling, margin blending) with real-time figure 07 refresh before finalizing.
	`random_forest_gap_filling`	`bool`	`false`	If `true`, applies Random Forest machine learning gap filling across unmeasured interior regions.
	`apply_margin_blend`	`bool`	`false`	If `true`, applies geomorphological margin blending to seamlessly transition calculated bedrock elevation to surrounding surface DEM terrain.
	`min_gap_dist`	`float`	`50.0`	Minimum gap distance in meters inside which bedrock is smoothly tapered and blended into surface DEM terrain.
	`smooth_bedrock`	`bool`	`false`	If `true`, applies spatial DEM post-processing smoothing directly to the depth field D(x,y) to eliminate high-frequency slope-division noise.
	`smoothing_method`	`str`	`"gaussian"`	Final bedrock DEM smoothing algorithm choice: `"gaussian"` (default), `"median"`, or `"fft_lowpass"`.
	`smoothing_sigma`	`float`	`1.5`	Smoothing strength (radius in pixels) for `"gaussian"` filtering. Higher values produce smoother bedrock terrain.
	`smoothing_kernel_size`	`int`	`3`	Window kernel size (k × k) for `"median"` filtering (must be an odd integer). Higher values produce smoother bedrock terrain.
	`smoothing_kc_cutoff`	`float`	`null`	Corner frequency cutoff wavenumber (k_c,smooth) for `"fft_lowpass"`. If `null`, defaults to k_c,opt from the surface slope smoothing optimization. Lower values produce smoother bedrock terrain.
`outputs`	`output_dir`	`str`	`null`	General workspace directory for all output files. `"."` defines `output_dir` as the parent directory of the `pysole.json` configuration file. If `null` (default), automatically creates a `pysole` folder inside the parent directory of `survey_data_path`. All relative output paths are resolved relative to `output_dir`. Absolute output paths override `output_dir`.
	`output_format`	`str` / `list[str]`	`"tif"`	Desired export format(s): `"tif"`, `"asc"`, `"csv"`, `"npy"`, a list of formats (e.g. `["tif", "asc", "csv"]`), or `"all"` to export all four formats.
	`output_prefix`	`str`	`""`	Optional filename prefix or absolute filepath prefix prepended to export filenames. When empty (default `""`), exports use clean standard names with applied drift model slug (e.g. `bedrock_<model_slug>.<ext>`, `thickness_<model_slug>.<ext>`) and diagnostic figures follow `##_<name>.png`. When specified, exports become `<output_prefix>_<name>_<model_slug>.<ext>` and figures follow `##_<output_prefix>_<name>.png` (migrated points CSV is exported as `[<output_prefix>_]migrated_points.csv` without drift slug). Absolute filepath prefix overrides `output_dir` and diagnostic figures are written beside it in a `figures/` subfolder.
	`plots_dir`	`str`	`"figures"`	Diagnostic figures are always saved - `plots_dir` controls where. `null` is equivalent to `"figures"`. A relative path is resolved inside `output_dir`, or absolute `output_prefix`. `PySole` data products and figures always share the same location, except an absolute `plots_dir` is defined.
	`save_bedrock_grid`	`bool`	`true`	If `true` (default), exports the predicted bedrock elevation raster grid Z_bed(x,y) to `<output_prefix>_bedrock_<model_slug>.<ext>`.
	`save_traveltime_grid`	`bool`	`false`	If `true`, exports the pre-migration interpolated traveltime raster grid, T(x,y), masked by the creeping body outline to `<output_prefix>_traveltime_<model_slug>.<ext>`.
	`save_traveltime_uncertainty`	`bool`	`false`	If `true`, exports the Kriging traveltime uncertainty raster grid, σ_T(x,y), in traveltime units and masked by the creeping body outline to `<output_prefix>_traveltime_uncertainty_<model_slug>.<ext>`. Requires `pre_migration.calculate_uncertainty: true`. Unavailable (with a warning) if migration is skipped or the survey data are depths.
	`save_migrated_points`	`bool`	`false`	If `true`, exports the 3D ray-migrated survey points to `<output_prefix>_migrated_points.csv` (`x, y, z_surface, depth_migrated`). If migration is skipped, it is automatically set to `false` and a log message is printed.
	`save_thickness_grid`	`bool`	`false`	If `true`, exports the final thickness raster grid, D(x,y), masked by the creeping body outline to `<output_prefix>_thickness_<model_slug>.<ext>`.
	`save_thickness_uncertainty`	`bool`	`false`	If `true`, exports the Kriging thickness uncertainty raster grid, σ_D(x,y), masked by the creeping body outline to `<output_prefix>_thickness_uncertainty_<model_slug>.<ext>`. Requires `post_migration.calculate_uncertainty: true`. Skipped with a warning if `calculate_uncertainty` is `false` or if the memory guard had to skip the variance.
	`save_basal_shear_stress`	`bool`	`false`	If `true`, exports the final basal shear stress raster grid, τ_b(x,y), masked by the creeping body outline to `<output_prefix>_basal_shear_stress_<model_slug>.<ext>`.
	`save_basal_shear_stress_uncertainty`	`bool`	`false`	If `true`, exports the basal shear stress Kriging uncertainty raster grid, σ_{τ_b}(x,y), masked by the creeping body outline to `<output_prefix>_basal_shear_stress_uncertainty_<model_slug>.<ext>`. Requires `post_migration.calculate_uncertainty: true`. Skipped with a warning if `calculate_uncertainty` is `false` or if the memory guard had to skip the variance.

Technical & Methodological Notes

1. Supported DEM Formats

PySole supports 5 distinct DEM formats:

GeoTIFF (.tif, .tiff, .geotiff): Recommended. Automatically extracts spatial bounds, CRS, pixel resolution, and nodata values using rasterio.
ESRI ASCII Grid (.asc, .txt): Standard 6-line header GIS raster format (ncols, nrows, xllcorner, yllcorner, cellsize, nodata_value). Explicitly defining spatial_parameters.crs is recommended.
CSV Grid (.csv): Supports both 2D elevation matrices ([Z_nx × Z_ny], metadata from spatial_parameters) and 3-column XYZ grid tables ([X, Y, Z], spacing and bounds inferred automatically).
NumPy Binary Array (.npy): Fast 2D binary array format; spatial metadata (origin, crs, dx, dy) defined via spatial_parameters.
NumPy 2D Array (np.ndarray): In-memory array passed directly into the Solver constructor (dem=dem_grid); spatial metadata defined via spatial_parameters.

2. Parsing of Rock Outcrop Input Files

Vector polygon files (.shp, .geojson, .gpkg) or CSV outline files containing interior rings (separated by NaN rows) are automatically parsed as polygon holes. PySole treats pixels inside rock outcrop holes as exposed bedrock (Z = Z_Surface). When boundary conditions are enabled, the perimeters of interior rings are included with T = 0 s or D = 0 m.

For rock outcrop holes to be detected from a Shapefile (.shp):

Topology: Stored as interior rings (polygon.interiors) within a Polygon / MultiPolygon feature (e.g. QGIS "Add Ring" or ArcGIS "Construct Hole").
Winding Order & Validity: Standard OGC orientation (Clockwise exterior, Counter-Clockwise interior rings) without self-intersections (validate_and_extract_polygons() auto-repairs minor topology issues on load).
CRS Alignment: Must match the DEM raster projection.

3. Parsing of Survey Data CSV Files

PySole ingests point survey datasets (inputs.survey_data_path) containing scattered soundings of two-way traveltime, one-way traveltime, or direct ice thickness. Internally, soundings are mapped into a standardized array:

$$\begin{bmatrix} X & Y & Z_{\text{surface}} & \text{value} & [\text{profile\_id}] \end{bmatrix}$$

Missing Surface Elevation Z_surface (3- vs. 4-Column Inputs): For 3-column inputs [X, Y, value], PySole automatically samples Z_surface from the surface DEM (inputs.dem_path) using 2D bilinear interpolation. For 4-column inputs [X, Y, Z, value], the explicit surface elevation Z is preserved directly.
With Header Row (e.g. profile_id,x,y,z,twt or x,y,twt,profile_id): PySole extracts the profile column matching inputs.survey_profile_column (default: "profile_id") by name while the remaining columns have to follow the standard [X,Y,(Z),value] sequence. This allows the profile column to be placed anywhere. This enables the Leave-One-Profile-Out cross-validation in the Drift Analyzer. If Z_surface is omitted (e.g. x,y,twt,profile_id), PySole samples Z_surface from the DEM automatically.
Without Header Row (Headerless Numeric CSV): Read positionally via np.loadtxt. Columns must strictly follow standard positional ordering: [X, Y, value], [X, Y, Z_surface, value], or [X, Y, Z_surface, value, profile_id]. Setting "survey_profile_column": null bypasses header searching and suppresses warnings. Non-standard ordering (e.g. placing profile IDs in Column 0, or omitting Z in a 4-column headerless file) causes positional misinterpretation or spatial DEM bounds errors.

Tip

Rule of Thumb for Survey CSVs:

Headerless CSV: Follow strict positional ordering [X, Y, value], [X, Y, Z_surface, value], or [X, Y, Z_surface, value, profile_id]. Set "survey_profile_column": null.
Header Present: The profile column can appear in any position. If Z_surface is omitted or profile_id is in a non-standard position, always include a header row and set "survey_profile_column" to match the header name.

4. Wavenumbers and Wavelengths

In PySole, 2D spatial Gaussian low-pass smoothing operates in the physical 2D spatial frequency domain with wavenumber components k_x, k_y [rad/m] and magnitude k = √(k_x² + k_y²):

$$k_x = \frac{2\pi \cdot \text{fftfreq}(N_x)}{dx}, \quad k_y = \frac{2\pi \cdot \text{fftfreq}(M_y)}{dy}, \quad k_c = \frac{2\pi}{\lambda_c} \Longleftrightarrow \lambda_c = \frac{2\pi}{k_c}$$

where dx, dy [m] are grid cell resolutions and N_x, M_y are grid dimensions. For a DEM with dx = dy = 5.0 m, the Nyquist wavelength is λ_Nyquist = 2 · min(|dx|, |dy|) = 10.0 m (k_Nyquist = π / 5.0 ≈ 0.6283 rad/m):

Cutoff Wavenumber k_c [rad/m]	Spatial Wavelength λ_c [m]	Surface Feature Scale
k_Nyquist ≈ 0.6283 rad/m	10.0 m	Nyquist limit (2 · min(
0.3142 rad/m	20.0 m	Fine spatial smoothing (filters features < 20 m)
0.1257 rad/m	50.0 m	Medium-fine spatial smoothing
0.0628 rad/m	100.0 m	Medium spatial smoothing
0.0314 rad/m	200.0 m	Broad spatial smoothing
0.0100 rad/m	628.3 m	Very broad regional smoothing

Table 1: Physical Reference Conversion Table.

Derivation of FFT Filter Parameter Defaults

When kc_max, kc_min, lambda_min, lambda_max, or n_steps are left as null in pysole.json, PySole derives defaults from the grid resolution and maximum domain extent L_max = max(N_x · |dx|, M_y · |dy|):

Maximum Wavenumber kc_max / Minimum Wavelength lambda_min: Defaults to the Nyquist limit k_c,max = k_Nyquist = π / min(|dx|, |dy|) (λ_min = 2 · min(|dx|, |dy|)). Custom lambda_min converts via k_c,max = 2π / λ_min.
Minimum Wavenumber kc_min / Maximum Wavelength lambda_max: Defaults to the half-domain limit k_c,min = 4π / L_max (λ_max = L_max / 2). Custom lambda_max converts via k_c,min = 2π / λ_max.
Evaluation Step Count n_steps: Counts the integer Fourier modes across [k_c,min, k_c,max] clamped to [10, 50] and spaced geometrically (np.geomspace):
$$n_{\text{steps}} = \min\left(50,\ \max\left(10,\ \left\lfloor \frac{(k_{\text{c,max}} - k_{\text{c,min}}) \cdot L_{\text{max}}}{2\pi} \right\rfloor\right)\right)$$

After each surface slope optimization pass, both k_c,opt [rad/m] and λ_c,opt [m] are logged to the console and pysole.log.

5. Variogram Binning with Minimum Pair Threshold

Experimental variogram lag distance bins are calculated from pairwise Euclidean distances across N survey points. When nrbins is null (default), PySole sets the bin count from the number of point pairs within the variogram range h ≤ h_max (half of the maximum pairwise distance) to ensure at least 30 point pairs per bin (Webster & Oliver, 2007), bounded to [3, 30]:

$$ \text{nrbins} = \min\left(30,\ \max\left(3,\ \left\lfloor \frac{N_{\text{pairs}}(h \le h_{\text{max}})}{30} \right\rfloor\right)\right) $$

If a user-specified nrbins yields fewer than 30 average pairs per bin, a warning is logged while honoring the explicit choice (enforcing a floor of 3 bins). For large surveys (N > 5000), lag statistics are accumulated in chunks using N_pairs = N(N − 1)/2.

6. High-Performance Dual Kriging Vector Engine

PySole features a native geostatistical engine based on Dual Kriging (Matheron, 1981). Unlike Primal Kriging (which solves a separate linear system at every target grid node), Dual Kriging solves the global (N + L) × (N + L) linear system once via LU decomposition (scipy.linalg.lu_factor):

$$\mathbf{K} \mathbf{w}_z = \mathbf{z}_{\text{aug}}, \qquad Z_{\text{grid}}(x,y) = \mathbf{w}_z \cdot \mathbf{k}_{\text{grid}} = \mathbf{w}_{\text{sample}} \cdot \boldsymbol{\Gamma}_{\text{grid}} + \mathbf{w}_{\text{drift}} \cdot \mathbf{F}_{\text{grid}}$$

where K is the augmented sample variogram and drift matrix, w_z = [w_sample^T, w_drift^T]^T is the solved global dual weight vector, and k_grid = [Γ_grid^T, F_grid^T]^T combines the 1D sample-to-grid cross-variogram vector Γ_grid and the drift basis vector F_grid evaluated at (x,y). Point prediction across m grid nodes reduces to an O(N · m) BLAS dot product (~180x speedup over loop-based solvers, interpolating 300,000+ nodes in under 50 ms). The production solver and the Drift Analyzer share the same unified KrigingSystem.

Thread-Safe Level-3 BLAS Variance Engine

While point prediction requires only an O(N · m) vector dot product, evaluating Kriging estimation variance σ²(x,y) = C(0) − k_grid^T K⁻¹ k_grid scales as O(N² · m) (disabled by default via calculate_uncertainty: false). When enabled, PySole pre-computes the explicit inverse K⁻¹ once via lu_solve(lu_piv, np.eye(N_aug)) prior to grid chunking, and worker threads (ThreadPoolExecutor) evaluate chunk variance weights via lock-free Level-3 BLAS matrix products (W_sub = K⁻¹ K_rhs,sub) while threadpoolctl prevents BLAS thread oversubscription.

Adaptive Diagonal Tikhonov Regularization & Scale Invariance

To prevent ill-conditioning when soundings are densely spaced along GPR profiles, PySole applies adaptive diagonal Tikhonov regularization strictly to the N × N sample observation block K_ii += λ_reg (i = 1, …, N), scaled by the empirical variogram sill S:

$$\lambda_{\text{reg}} = \epsilon_{\text{reg}} \cdot S = 10^{-6} \cdot S$$

Because both Γ and λ_regI scale proportionally with S while the drift Lagrange blocks remain unperturbed (F^T w_sample = 0), predictions are strictly invariant to linear rescaling of the input data.

7. Universal Kriging Drift Models

The Drift Analyzer can be enabled for pre- and/or post-migration Universal Kriging (drift_analyzer: true) to automatically evaluate, rank, and compare candidate drift models using spatial cross-validation metrics — see the Drift Analyzer Guide for full details. Because statistical rankings should always be paired with glaciological judgment, the governing assumptions of all primitive and multi-drift models are summarized below.

Single Drift Models (8 Primitives + Ordinary Kriging)

Shallow Ice Approximation Physical Drift Model (["sia"]): Rearranging basal shear stress τ_b = ρ_ice g D sin(α) for ice depth D yields an inverse relationship with the optimally smoothed surface slope sin(α_opt(x,y)) (clamped at slope_floor_deg, default 5.0°):

$$U(x,y) = \sin(\alpha_{\text{opt}}(x,y))^{-1}$$

Default drift model for direct targets ("T" or "D"). Produces a terrain-conforming physical background trend across unmeasured regions without requiring assumptions about absolute τ_b. Use only with direct targets ("T" or "D") to avoid 1/sin²(α) double-scaling on "P".
Surface Elevation Drift Model (["z_dem"]): Uses the unsmoothed DEM surface elevation Z_dem(x,y):

$$U(x,y) = Z_{\text{dem}}(x,y)$$

Captures elevation-dependent thickness patterns (e.g. thicker ice in lower valley basins/confluences or thinner ice on steep high-altitude ridges and summits).
Surface Curvature Drift Model (["curvature_dem"]): Uses the 2D Laplacian surface curvature C_{k_c}(x,y) derived from the optimally smoothed DEM surface:

$$C_{k_c}(x,y) = \nabla^2 Z_{\text{smooth}, k_c}(x,y) = \frac{\partial^2 Z}{\partial x^2} + \frac{\partial^2 Z}{\partial y^2}$$

Captures morphometric terrain curvature—predicting smaller depths over convex ridges (∇²Z < 0) and larger depths in concave bowls and troughs (∇²Z > 0).
Surface Mass Balance Drift Model (["mass_balance"]): Uses a mean surface mass balance raster ḃ(x,y) (configured via inputs.mass_balance_path or Solver(..., mass_balance=...)), co-registered to the DEM grid and kept unsmoothed to preserve localized accumulation and ablation anomalies (supraglacial debris cover, avalanche cones, wind drift, topographic shading):

$$U(x,y) = \dot{b}(x,y)$$

Captures mass-conservation forcing where localized mass turnover anomalies drive ice thickness variations independent of elevation alone.
Surface Flow Velocity Magnitude Drift Model (["velocity"]): Uses a mean surface flow velocity magnitude raster u_{s,smooth,k_c}(x,y) = |u_s(x,y)|_{k_c} (configured via inputs.velocity_path or Solver(..., velocity=...)), zero-clamped outside the glacier outline (~outline_mask -> 0) and spectrally low-pass filtered at the optimal cutoff wavenumber k_c:

$$U(x,y) = u_{s,\text{smooth}, k_c}(x,y)$$

Differentiates fast-flowing, deep main glacier trunks from stagnant, shallow margins or dead-ice zones that share similar surface slopes.
Glen's Flow Law Coupled SIA–Velocity Drift Model (["velocity_sia"]): Couples the spectrally smoothed surface velocity u_{s,smooth,k_c}(x,y) and the floored optimal surface slope sin(α_opt(x,y)) via Glen's flow law (n = 3) for internal ice creep (u_d = [2A / (n + 1)] (ρ_ice g sin α)ⁿ Dⁿ⁺¹ ⇒ D ∝ u_s^1/4 / (sin α)^3/4):

$$U(x,y) = \frac{\left(\max\left(u_{s,\text{smooth}, k_c}(x,y),\, 0\right)\right)^{1/4}}{\left(\sin(\alpha_{\text{opt}}(x,y))\right)^{3/4}}$$

Provides a physically coupled creep-thickness covariate without requiring a known creep rate factor A. Because "velocity_sia" already couples "velocity" and "sia", it is never combined with "sia" or "velocity" in the same multi-drift model. Like "sia", use "velocity_sia" only with direct targets ("T" or "D").
Linear Surface Drift Model (["linear_xy"]): Fits a 1st-order planar coordinate trend across normalized map axes (X, Y):

$$U(x,y) = a_1 X + a_2 Y$$

Captures a constant regional spatial gradient across the domain.
Quadratic Surface Drift Model (["quadratic_xy"]): Fits a 2nd-order parabolic surface across normalized map axes (X, Y):

$$U(x,y) = a_1 X + a_2 Y + a_3 X^2 + a_4 Y^2 + a_5 X Y$$

Captures regional spatial bends, basin bowls, ice cap domes, or saddle geometry across map space.
No Drift ([] - Ordinary Kriging): Assumes a constant local spatial mean without external drift covariates (equivalent to method: "ordinary"). Inside a Drift Analyzer candidate suite, include [] as a nested entry (e.g. drift_terms: [[], "sia", "z_dem"]) to evaluate Ordinary Kriging alongside Universal Kriging candidates. An empty list [] at the top level is rejected (use null instead, or set method: "ordinary").

Multi-Drift Model Structure

PySole combines the physical raster drifts ("sia", "z_dem", "curvature_dem", "mass_balance", "velocity", "velocity_sia") with each other and with the spatial polynomial coordinate trends ("linear_xy", "quadratic_xy") into an augmented multi-drift Universal Kriging system.

With all 8 primitives active, the combinatorial space is governed by three independent groups after applying the two collinearity pruning rules:

Group 1 - SIA & Velocity ("sia", "velocity", "velocity_sia" — 5 valid states):

[] (no SIA or velocity term)
["sia"]
["velocity"]
["sia", "velocity"]
["velocity_sia"] (never combined with "sia" or "velocity" in the same model because it already couples both)

Group 2 - Independent scalar rasters ("z_dem", "curvature_dem", "mass_balance" — 2³ = 8 valid states):

From [] (no scalar raster term) up to ["z_dem", "curvature_dem", "mass_balance"] (each of the 3 rasters independently included or excluded).

Group 3 - Polynomial coordinate trends ("linear_xy", "quadratic_xy" — 3 valid states):

[] (no polynomial term), ["linear_xy"], or ["quadratic_xy"] (never combined together because "quadratic_xy" already contains the linear [x, y] columns).

In total, 5 × 8 × 3 = 120 candidate models are available: selecting [] across all three groups yields Ordinary Kriging ([]), while the remaining 119 combinations are Universal Kriging drift models (8 single-drift + 111 multi-drift).

Breakdown of the 120 Candidates

Category	Count	Details
Ordinary Kriging	1	`[]` (no drift terms: Group 1 = `[]`, Group 2 = `[]`, Group 3 = `[]`)
Single Drift Models	8	6 single-raster (`["sia"]`, `["z_dem"]`, `["curvature_dem"]`, `["mass_balance"]`, `["velocity"]`, `["velocity_sia"]` + 2 single-polynomial (`["linear_xy"]`, `["quadratic_xy"]`)
Multi-Raster Drift Models (no polynomial)	33	5 × 8 − 1 − 6 = 33 combinations of 2 to 5 raster terms: 5 × 8 = 40 raster states (Group 1 × Group 2, with Group 3 = `[]`), minus 1 empty state (`[]`, Ordinary Kriging) and minus the 6 single-raster models above)
Multi-Drift Models with `"linear_xy"`	39	5 × 8 − 1 = 39 models (6 single-raster + 33 multi-raster paired with `"linear_xy"`): 5 × 8 = 40 raster states paired with `"linear_xy"`, minus 1 where the raster state is `[]` (which is the single-drift model `["linear_xy"]` above)
Multi-Drift Models with `"quadratic_xy"`	39	5 × 8 − 1 = 39 models (6 single-raster + 33 multi-raster paired with `"quadratic_xy"`): 5 × 8 = 40 raster states paired with `"quadratic_xy"`, minus 1 where the raster state is `[]` (which is the single-drift model `["quadratic_xy"]` above)
Total Universal Kriging Drift Models	119	8 Single + 111 Multi Drift Models (33 + 39 + 39 = 111)
Total Evaluated by `DriftAnalyzer`	120	119 Universal Kriging Drift Models + Ordinary Kriging

Dynamic Suite Size by Configured Input Rasters (`drift_terms: null`, `drift_analyzer: true`)

Configured Input Rasters	Active Primitives	Single Drift	Multi Drift	Ordinary	Total Suite Size
DEM only (default)	`sia`, `z_dem`, `curvature_dem`, `linear_xy`, `quadratic_xy` (5)	5	18	1	24 candidates
DEM + `mass_balance_path`	Core 5 + `mass_balance` (6)	6	41	1	48 candidates
DEM + `velocity_path`	Core 5 + `velocity`, `velocity_sia` (7)	7	52	1	60 candidates
DEM + `mass_balance_path` + `velocity_path`	All 8 primitives	8	111	1	120 candidates

Note

For the complete catalog of all 111 Multi-Drift Models, their mathematical formulations, parameter counts (n_drift), and glaciological interpretations, see the Drift Analyzer & Multi-Drift Guide (docs/drift_analyzer.md).

How `drift_terms` is Interpreted

The Kriging parameters method and drift_analyzer determine how drift_terms is resolved:

`method`	`drift_analyzer`	Meaning of `drift_terms`
`"ordinary"`	`any`	Ignored (`ordinary` overrides the analyzer and the terms; one warning if either is set).
`"universal"`	`false`	Flat Drift Model List: All listed terms (e.g. `["sia", "z_dem"]`) are applied together as ONE multi-drift model (a single nested candidate `[["sia", "z_dem"]]` unwraps to `["sia", "z_dem"]`). Non-Flat Drift Model List (e.g. `[["sia", "z_dem"], ["sia"]]`): Executed as a Non-Interactive Multi-Model Run using the default spatial CV exclusion buffer radius.
`"universal"`	`true`	Flat Drift Model List: Evaluates an interactive combinatorial candidate suite of all valid non-redundant subsets of the listed terms (e.g. `["sia", "z_dem"]` tests `["sia"]`, `["z_dem"]`, and `["sia", "z_dem"]`; nesting `[]` adds Ordinary Kriging). Default (`null`): Dynamically evaluates 24 candidates (DEM only), 48 candidates (DEM + `mass_balance_path`), 60 candidates (DEM + `velocity_path`), or 120 candidates (DEM + `mass_balance_path` + `velocity_path`). Non-Flat Drift Model List (e.g. `[["sia", "z_dem"], ["sia"], []]`): Evaluates the exact listed models without combinatorial expansion.
`"regression"`	`false`	Random Forest feature covariates: Predictor features used to train a Random Forest trend model before Ordinary Kriging of residuals (only available through `PyKrige`; automatically switches engine if `engine: "native"`). Default (`null`): Uses `["sia", "z_dem", "curvature_dem", "quadratic_xy"]` plus `"mass_balance"` and `"velocity"` when configured.
`"regression"`	`true`	Analyzer skipped with warning: Logs a warning, skips the analyzer, and runs Regression Kriging with `drift_terms`.

Candidate Suite & Single-Model Rules

Single-model auto-disable: When drift_analyzer: true resolves to a single candidate (e.g. ["sia"] or [["sia", "z_dem"]]), PySole logs an INFO notice, sets drift_analyzer = false, and runs Universal Kriging directly.
Collinearity pruning: In an analyzer suite (drift_analyzer: true), candidates never combine linear_xy with quadratic_xy or velocity_sia with sia / velocity within the same model. In a single model (drift_analyzer: false), redundant "linear_xy" or "sia"/"velocity" terms are automatically dropped with a warning in favor of "quadratic_xy" or "velocity_sia".
Top-level empty list: An empty list [] at the top level is rejected (use null instead).
Product target ("P") safeguard: For interpolation_target: "P", the Drift Analyzer evaluates every candidate containing "sia" or "velocity_sia" but flags it with a warning and sorts it below non-penalized models due to 1/sin²(α) double-scaling.
Multi-Model Comparison Suite: Selecting multiple models interactively (e.g. 1-3, 5) or passing an explicit multi-model list exports all rasters to <output_dir> with <model_slug> suffixes, routes diagnostic plots to figures/<model_slug>/, and saves drift_model_comparison.{txt,csv,md} in <output_dir> — see Multi-Model Comparison Suite.
engine: "pykrige" performance note: PyKrige's UniversalKriging computes an internal leave-one-out statistic on instantiation that scales as O(N³–N⁴) in the number of survey points. Use engine: "native" for Universal Kriging on datasets with more than a few hundred points.

8. Kriging Uncertainty Derivation

When calculate_uncertainty: true is enabled, PySole evaluates the 2D Kriging estimation variance field σ²(x,y) and converts it into the Kriging Standard Error σ(x,y) (±1.00 σ corresponds to the 68.3% Gaussian confidence interval; ±1.96 σ to 95%):

Direct Target Interpolation (interpolation_target: "D" or "T"): The estimation variance field (σ_D²(x,y) [m²] or σ_T²(x,y) [s²]) is evaluated directly:

$$\sigma_{\text{D or T}}(x,y) = \sqrt{\sigma_{\text{D or T}}^2(x,y)} \quad [\pm\,\text{m or s}]$$
BSS-derived Product Target Interpolation (interpolation_target: "P"): The product variance σ_P²(x,y) is propagated through the optimal surface slope field sin(α_opt(x,y)):

$$\sigma_{\text{D or T}}(x,y) = \sqrt{\frac{\sigma_{\text{P}}^2(x,y)}{\sin^2(\alpha_{\text{opt}}(x,y))}} \quad [\pm\,\text{m or s}]$$
Total Ice Volume Uncertainty (σ_V): Across N valid glacier cells of area ΔA = |dx · dy| (A_total = N · ΔA), Kriging interpolation errors are spatially autocorrelated. Following Binder et al. (2009), total volume uncertainty σ_V is integrated from the domain-average thickness standard error σ̄_D:

$$\sigma_V = A_{\text{total}} \cdot \bar{\sigma}_{\text{D}} = \sum_{i \in \text{glacier}} \sigma_{\text{D}, i} \cdot \Delta A \quad \Longleftrightarrow \quad \frac{\sigma_V}{V} = \frac{\bar{\sigma}_{\text{D}}}{\bar{D}}$$

9. Memory Guard and Kriging Uncertainty

When calculate_uncertainty: true, post-migration standard error (σ_D) is plotted in diagnostic figures and exported when save_*_uncertainty: true (save_traveltime_uncertainty exports σ_T when pre-migration uncertainty is enabled). Because Kriging standard error depends only on sample geometry and the variogram, it represents a lower bound on total uncertainty.

Before allocating the (N + n_drift)² inverse matrix and per-thread grid chunks, PySole's memory guard compares estimated peak RAM against kriging_parameters.max_memory_fraction (default 0.5, maximum 0.9) of available RAM (respecting Linux cgroup container quotas in Docker/Kubernetes/HPC):

Worker thread count is reduced first (preserving exact numerical results).
Only if a single thread with variance still exceeds the budget is variance skipped (logging a warning and displaying "Uncertainty not available" in diagnostic plots).
If available RAM cannot be determined, the guard is disabled.

10. Spatial Smoothing of the Calculated Depth and Bedrock DEMs

When smooth_bedrock: true, PySole applies the spatial smoothing operator S ("gaussian", "median", or "fft_lowpass") directly to the ice depth field D(x,y) before subtracting from the surface DEM:

$$D_{\text{smooth}}(x,y) = S(D(x,y)), \qquad Z_{\text{bed}}(x,y) = Z_{\text{surface}}(x,y) - D_{\text{smooth}}(x,y)$$

Smoothing D(x,y) directly prevents high-frequency surface DEM roughness residuals Z_surface − S(Z_surface) from imprinting onto the ice thickness map.

Package Architecture

Figure 2: Overview of PySole package architecture, class structure, sub-engine modules, and API methods. High-level Pipeline and Solver methods are color-coded by their underlying sub-engine domain. Click diagram to view in high resolution.

Command Line Interface

PySole provides a convenient terminal entry point for executing processing workflows directly from your command line without writing Python scripts.

1. Convenient Execution & `$PATH` Setup

When you install PySole (pip install . or pip install -e .), pip automatically creates the pysole executable binary script in your Python environment's binary directory (bin/ on Linux/macOS or Scripts\ on Windows).

Virtual Environment (Recommended): When your virtual environment is activated (source .venv/bin/activate or conda activate myenv), its bin/ directory is automatically prepended to your $PATH. You can immediately execute pysole from any working directory:
```
pysole pysole.json
```
User Installation without Virtual Environment (pip install --user .): If installed without a virtual environment using --user, Python places entry point scripts in ~/.local/bin (Linux/macOS) or %APPDATA%\Python\Scripts (Windows). If pysole is not recognized by your terminal, add ~/.local/bin to your $PATH variable by placing the following line in your ~/.bashrc:
```
export PATH="$HOME/.local/bin:$PATH"
```

2. Available CLI Commands

Generate a Default pysole.json Template:
```
pysole --init
```
Creates a clean template pysole.json configuration file in your current working directory.
Run PySole with Default pysole.json Configuration:
```
pysole
```
Automatically loads and executes pysole.json in the current working directory.
Run PySole with a Custom Configuration File:
```
pysole path/to/custom_config.json
```
Executes the full pipeline defined in custom_config.json.
Run in Non-Interactive Batch Mode:
```
pysole pysole.json --batch
```
Runs PySole in non-interactive batch mode, automatically disabling all interactive prompts and terminal dialogs. Built-in mode for HPC cluster jobs and automated background pipelines.
Run with Verbose / Debug Logging:
```
pysole pysole.json --verbose
```
Enables DEBUG level logging verbosity for detailed computational diagnostics (-v, --verbose and --debug are equivalent). Acts as a temporary CLI runtime override taking precedence over the PYSOLE_LOG_LEVEL environment variable and the log_level defined in pysole.json.
Enable the Interactive Drift Analyzer:
```
pysole pysole.json --drift-analyzer
```
Runs the Drift Analyzer helper before each Universal Kriging pass as an on-the-fly override across all passes (equivalent to enabling drift_analyzer: true on both Kriging passes). Omitting the flag preserves the per-pass drift_analyzer settings configured in pysole.json. Ignored (with a warning) for Ordinary Kriging.
```
pysole pysole.json --drift-analyzer --profile-col line_id
```
Specifies a custom column header name in the imported survey data CSV file (inputs.survey_data_path) that identifies individual survey profiles or lines. It enables Leave-One-Profile-Out (LOPO) cross-validation in the Drift Analyzer and overrides the default inputs.survey_profile_column ("profile_id").

See here for full details and CLI options.

Survey Planner (no survey data required):

# Interactive mode: prompts for omitted parameters (tau_0, lambda_c, max-km)
pysole plan-survey --dem dem.tif --outline outline.geojson

# Non-interactive mode: fully defined flags (or pass --batch)
pysole plan-survey --dem dem.tif --outline outline.geojson --lambda-c 200 --tau 100 --max-km 5 --out-dir results --format tif

Computes a synthetic SIA thickness model and proposes optimal GPR/seismic survey tracks within a length budget. If parameters are omitted, PySole interactively prompts with defaults in the terminal (or applies defaults in batch mode).

See here for full details and CLI options.

Show the Installed Version:
```
pysole --version
```
Prints PySole <version> and exits (-V is equivalent). python -m pysole ... is also supported.
Display CLI Help & Usage Options:
```
pysole --help
```
Prints the complete list of available CLI flags, arguments, and usage descriptions and exits (-h is equivalent).

Error handling: a missing configuration file, invalid JSON, or an invalid configuration value (e.g. a non-boolean value such as "maybe" for a boolean parameter) is reported as a single pysole: error: ... line on stderr with exit status 1. Add -v / --debug to get the full traceback. Boolean parameters accept JSON true/false (also 0/1 and the strings "true"/"false", "yes"/"no", "on"/"off", case-insensitive); null falls back to the default.

Application Programming Interface

PySole provides a dual-layer Python API: a High-Level Solver Orchestrator for streamlined pipeline execution, and Decoupled Specialized Sub-Engines for custom advanced research workflows.

1. High-Level One-Liner Execution

Run the complete pipeline from a pysole.json configuration file in a single line:

import pysole

# Execute complete workflow using default 'pysole.json'
final_bedrock = pysole.run_from_config()

# Or execute with a custom configuration file and log level override
final_bedrock = pysole.run_from_config("custom_config.json", log_level="DEBUG")

# Inspect spatial metadata and export final bedrock raster
print(final_bedrock.shape, final_bedrock.bounds)
final_bedrock.save("final_bedrock.tif")

2. High-Level `Solver` Orchestrator API

Initialize the central Solver facade to steer the pipeline through processing milestones. When instantiated, Solver creates and manages underlying specialized sub-engine instances (self.migrator, self.bss_optimizer, self.kriging_engine, self.finalizer):

import pysole

# 1. Initialize High-Level Solver Orchestrator
model = pysole.Solver(
    dem="surface_dem.asc",
    outline="creeping_body.shp",
    survey_data_type="one_way_traveltime",
    # mass_balance="mass_balance_decadal.tif",  # Optional surface mass balance raster
    # velocity="velocity_decadal.tif",          # Optional surface flow velocity raster
    pre_interpolation_target="T",
    pre_kriging_method="universal",
    pre_drift_terms=["sia"],
    post_interpolation_target="D",
    post_kriging_method="universal",
    post_drift_terms=["sia"],
    perform_migration=True,
)

# Inspect underlying specialized sub-engine instances managed by Solver
print(model.geometry)        # GridGeometry spatial container instance
print(model.migrator)        # EikonalMigrator sub-engine instance
print(model.bss_optimizer)   # BSSOptimizer sub-engine instance
print(model.kriging_engine)  # KrigingEngine sub-engine instance
print(model.finalizer)       # BedrockFinalizer sub-engine instance

# Option A: Run complete end-to-end pipeline in a single call
final_bedrock = model.run_pipeline(survey_data_path="sparse_survey.csv")

# Option B: Steer workflow step-by-step through individual milestones:
# 2. Migrate sparse GPR/Seismic traveltimes (delegates to model.migrator)
bedrock_pts = model.migrate_eikonal(
    travel_times="sparse_survey.csv",
    velocity=0.16,
)

# 3. Iterative BSS variance optimization (delegates to model.bss_optimizer)
model.optimize_bss(kc_max=0.3, kc_min=0.01, n_steps=20)

# 4. Primary Kriging spatial interpolation (delegates to model.kriging_engine)
kriged_bedrock, kriged_variance = model.calculate_bedrock(
    method="universal",
)

# 5. Finalize bedrock (delegates to model.finalizer for gap filling & margin blending)
final_bedrock = model.finalize_bedrock(
    interactive=False,
    smooth_bedrock=True,
    smoothing_method="gaussian",
)

# 6. Export predicted bedrock elevation raster
final_bedrock.save("final_bedrock.tif")

3. Decoupled Sub-Engines API (Standalone Custom Workflows)

Advanced users and researchers can bypass the Solver orchestrator entirely and instantiate specialized sub-engines directly for modular, custom research scripts:

import pysole
from pysole import (
    load_dem,
    load_outline,
    load_survey_points,
    GridGeometry,
    EikonalMigrator,
    BSSOptimizer,
    KrigingEngine,
    BedrockFinalizer,
    PipelineExporter,
)

# 1. Load raster inputs and construct spatial GridGeometry container
dem_grid, meta = load_dem("surface_dem.asc")
outline_mask = load_outline("creeping_body.shp", dem_grid=dem_grid, meta=meta)
geometry = GridGeometry.create(dem_grid.shape, dx=meta["dx"], dy=meta["dy"], bounds=meta["bounds"])
survey_pts = load_survey_points("sparse_survey.csv", bounds=geometry.bounds, dem_grid=dem_grid)

# 2. Standalone 3D Ray Migration Sub-Engine (pysole.migration.EikonalMigrator)
migrator = EikonalMigrator(dem=dem_grid, geometry=geometry, outline_mask=outline_mask)
# mig_res = migrator.migrate(travel_time_grid=tt_grid, survey_points=survey_pts, velocity=0.16)

# 3. Standalone Basal Shear Stress Optimization Sub-Engine (pysole.variogram.BSSOptimizer)
bss_optimizer = BSSOptimizer(dem=dem_grid, geometry=geometry)
opt_res = bss_optimizer.optimize(survey_points=survey_pts, kc_max=0.3, kc_min=0.01, n_steps=20)
print(f"Optimal Corner Frequency: {opt_res.optimal_kc:.4f} rad/m")

# 4. Standalone Dual Kriging Vector Sub-Engine (pysole.interpolation.KrigingEngine)
krig_engine = KrigingEngine(dem=dem_grid, geometry=geometry, outline_mask=outline_mask)
krig_res = krig_engine.interpolate(sample_points=survey_pts, method="universal", drift_terms=["z_dem"], variogram_model="spherical")
# Returns typed KrigingResult dataclass containing krig_res.bedrock_grid and krig_res.variance_grid

# 5. Standalone Bedrock Finalizer Sub-Engine (pysole.interpolation.BedrockFinalizer)
finalizer = BedrockFinalizer(dem=dem_grid, geometry=geometry, outline_mask=outline_mask)
rf_filled = finalizer.fill_holes(krig_res.bedrock_grid)
blended_bedrock = finalizer.blend_margin(rf_filled, min_gap_dist=50.0)

# 6. Standalone Disk Exporter Sub-Engine (pysole.pipeline.PipelineExporter)
# exporter = PipelineExporter(model)
# exporter.export_raster(blended_bedrock, suffix="final_bedrock", name="Bedrock Elevation")

Real-World Example: Wurtenkees Glacier

You can run a complete real-world demonstration on the Wurtenkees Glacier dataset (Hohe Tauern, Eastern Alps, Austria) using the provided configuration file examples/wuk/pysole_wuk.json.

An executable example script is provided in examples/wuk/run_wuk_example.py:

python examples/wuk/run_wuk_example.py

The GPR dataset and DEM inputs are sourced from the Master's thesis by Binder (2011, written in German and available from ResearchGate).

Citation & References

If you use PySole in your research, please cite the underlying methodology introduced by Binder et al. (2009):

BibTeX

@article{binder2009determination,
  title        = {Determination of total ice volume and ice-thickness distribution of two glaciers in the Hohe Tauern region, Eastern Alps, from GPR data},
  author       = {Binder, Daniel and Br{\"u}ckl, Ewald and Roch, Karl-Heinz and Behm, Michael and Sch{\"o}ner, Wolfgang and Hynek, Bernhard},
  journal      = {Annals of Glaciology},
  volume       = {50},
  number       = {51},
  pages        = {71--79},
  year         = {2009},
  publisher    = {Cambridge University Press},
  doi          = {10.3189/172756409789097522}
}

Full Reference List (APA)

Binder, D., Brückl, E., Roch, K.H., Behm, M., Schöner, W., & Hynek, B. (2009). Determination of total ice volume and ice-thickness distribution of two glaciers in the Hohe Tauern region, Eastern Alps, from GPR data. Annals of Glaciology, 50(51), 71–79. doi:10.3189/172756409789097522
Binder, D. (2011). Bestimmung der Eismächtigkeitsverteilung dreier Gletscher der Hohen Tauern auf Basis von Ground Penetrating Radar (GPR) Daten (Master's thesis, Vienna University of Technology, Vienna, Austria). Available from ResearchGate.
Matheron, G. (1981). Splines and kriging: Their formal equivalence. In D. F. Merriam (Ed.), Down-to-Earth statistics: Solutions looking for geological problems (Vol. 8, pp. 77–95). Syracuse University. (Syracuse University Geology Contribution No. 8).
Webster, R., & Oliver, M. A. (2007). Geostatistics for Environmental Scientists (2nd ed.). John Wiley & Sons. doi:10.1002/9780470517277

Name		Name	Last commit message	Last commit date
.github/workflows		.github/workflows
docs		docs
examples		examples
images		images
src/pysole		src/pysole
tests		tests
.gitignore		.gitignore
.python-version		.python-version
CHANGELOG.md		CHANGELOG.md
CODE_OF_CONDUCT.md		CODE_OF_CONDUCT.md
LICENSE		LICENSE
MANIFEST.in		MANIFEST.in
README.md		README.md
pyproject.toml		pyproject.toml
pysole.json		pysole.json
uv.lock		uv.lock

Latest commit

History

Folders and files

Repository files navigation

Table of Contents