Terrain Modeling

This example demonstrates the pytcl.terrain module capabilities, including digital elevation model (DEM) creation, synthetic terrain generation, and terrain analysis.

Overview

Terrain modeling is essential for:

• Navigation: Terrain-aided navigation and TERCOM

• Simulation: Realistic environment modeling

• Line-of-sight: Radio propagation and visibility

• Mission planning: Route optimization

Digital Elevation Models

Flat DEM

• Constant elevation surface

• Useful for testing and baseline comparisons

• Created with create_flat_dem()

Synthetic Terrain

• Procedurally generated terrain

• Controllable parameters (amplitude, wavelength)

• Useful for simulation and testing

• Created with create_synthetic_terrain()

DEM Properties

Grid Structure

• Regular lat/lon grid

• Specified resolution in arcseconds

• Elevation values at each grid point

Coordinate System

• Geographic coordinates (lat, lon)

• Elevation in meters above reference

Analysis Outputs

• Min, max, mean elevation

• Standard deviation

• Slope and aspect maps

Terrain Analysis

Elevation Statistics

• Distribution of elevation values

• Terrain roughness metrics

• Histogram analysis

Slope Computation

• Gradient magnitude at each point

• Degrees from horizontal

• Important for mobility analysis

Horizon Computation

• Visible horizon from observer position

• Accounts for terrain obstruction

• Essential for line-of-sight analysis

Applications

Terrain-Aided Navigation

• Match measured terrain to DEM

• Position fix without GPS

• Submarine and aircraft navigation

Viewshed Analysis

• Determine visible area from point

• Radar coverage planning

• Communication link analysis

Route Planning

• Avoid steep terrain

• Minimize exposure

• Optimize fuel consumption

Terrain-Aided Navigation: Vehicle trajectories can be matched against terrain models for position updates without GPS.

Code Highlights

The example demonstrates:

• Flat DEM creation with create_flat_dem()

• Synthetic terrain with create_synthetic_terrain()

• Terrain statistics (min, max, mean, std)

• Slope computation using gradients

• Horizon computation with compute_horizon()

Source Code

"""Terrain module demonstration with DEM and elevation models.

This example demonstrates the pytcl.terrain module capabilities, including
digital elevation model (DEM) creation, synthetic terrain generation, terrain
analysis, and viewshed/line-of-sight computations.

Functions demonstrated:
- create_flat_dem(): Create flat digital elevation models
- create_synthetic_terrain(): Generate synthetic terrain with hills/valleys
- compute_horizon(): Calculate visible horizon from observer position
"""

from pathlib import Path

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

from pytcl.terrain import compute_horizon, create_flat_dem, create_synthetic_terrain

# Controls for visualization
SHOW_PLOTS = False
SKIP_VISUALIZATIONS = True  # Skip visualizations for fast execution


def demo_flat_dem() -> None:
    """Demonstrate flat DEM creation."""
    print("\n" + "=" * 60)
    print("Flat Digital Elevation Model (DEM)")
    print("=" * 60)

    # Create a flat DEM
    dem = create_flat_dem(
        lat_min=-1,
        lat_max=1,
        lon_min=-1,
        lon_max=1,
        elevation=1000.0,
        resolution_arcsec=60,
    )

    print(f"\nFlat DEM created:")
    print(f"  Latitude range: [{dem.lat_min}, {dem.lat_max}]")
    print(f"  Longitude range: [{dem.lon_min}, {dem.lon_max}]")
    print(f"  Shape: {dem.data.shape}")
    print(f"  Elevation: {dem.data.min():.1f} to {dem.data.max():.1f} m")
    print(f"  Grid spacing: {dem.d_lat:.6f}° lat, {dem.d_lon:.6f}° lon")

    # Visualization (fast heatmap with downsampling for large grids)
    if not SKIP_VISUALIZATIONS:
        # Downsample for faster visualization (keep full resolution data for computations)
        max_size = 500
        stride = max(1, max(dem.data.shape) // max_size)
        z_display = dem.data[::stride, ::stride]

        fig = go.Figure(
            data=go.Heatmap(
                z=z_display,
                colorscale="Viridis",
                name="Elevation",
                colorbar=dict(title="Elevation (m)"),
            )
        )

        fig.update_layout(
            title="Flat Digital Elevation Model",
            xaxis_title="Longitude index",
            yaxis_title="Latitude index",
            height=500,
        )

        if SHOW_PLOTS:
            fig.show()
        else:
            fig.write_html(str(OUTPUT_DIR / "terrain_demo.html"))


def demo_synthetic_terrain() -> None:
    """Demonstrate synthetic terrain generation."""
    print("\n" + "=" * 60)
    print("Synthetic Terrain Generation")
    print("=" * 60)

    # Create synthetic terrain with hills
    dem = create_synthetic_terrain(
        lat_min=-5,
        lat_max=5,
        lon_min=-5,
        lon_max=5,
        base_elevation=500,
        amplitude=800,
        wavelength_km=50,
        resolution_arcsec=120,
        seed=42,
    )

    print(f"\nSynthetic DEM created:")
    print(f"  Latitude range: [{dem.lat_min}, {dem.lat_max}]")
    print(f"  Longitude range: [{dem.lon_min}, {dem.lon_max}]")
    print(f"  Shape: {dem.data.shape}")
    print(f"  Min elevation: {dem.data.min():.1f} m")
    print(f"  Max elevation: {dem.data.max():.1f} m")
    print(f"  Mean elevation: {dem.data.mean():.1f} m")
    print(f"  Std deviation: {dem.data.std():.1f} m")
    print(f"  Grid spacing: {dem.d_lat:.6f}° lat, {dem.d_lon:.6f}° lon")

    # Visualization: 2D heatmap with downsampling (fast rendering)
    if not SKIP_VISUALIZATIONS:
        # Downsample for faster visualization
        max_size = 500
        stride = max(1, max(dem.data.shape) // max_size)
        z_display = dem.data[::stride, ::stride]

        fig = go.Figure(
            data=go.Heatmap(
                z=z_display,
                colorscale="Earth",
                name="Elevation",
                colorbar=dict(title="Elevation (m)"),
            )
        )

        fig.update_layout(
            title="Synthetic Terrain with Synthetic Hills",
            xaxis_title="Longitude index",
            yaxis_title="Latitude index",
            height=600,
        )

        if SHOW_PLOTS:
            fig.show()
        else:
            fig.write_html(str(OUTPUT_DIR / "terrain_demo.html"))


def demo_terrain_analysis() -> None:
    """Demonstrate terrain analysis and statistics."""
    print("\n" + "=" * 60)
    print("Terrain Analysis and Statistics")
    print("=" * 60)

    # Create two DEMs for comparison
    flat_dem = create_flat_dem(
        lat_min=-2,
        lat_max=2,
        lon_min=-2,
        lon_max=2,
        elevation=500.0,
        resolution_arcsec=60,
    )

    synthetic_dem = create_synthetic_terrain(
        lat_min=-2,
        lat_max=2,
        lon_min=-2,
        lon_max=2,
        base_elevation=500,
        amplitude=400,
        wavelength_km=30,
        resolution_arcsec=60,
        seed=123,
    )

    # Compute statistics
    print(f"\nFlat DEM statistics:")
    print(f"  Min: {flat_dem.data.min():.1f} m")
    print(f"  Max: {flat_dem.data.max():.1f} m")
    print(f"  Mean: {flat_dem.data.mean():.1f} m")
    print(f"  Std Dev: {flat_dem.data.std():.1f} m")

    print(f"\nSynthetic DEM statistics:")
    print(f"  Min: {synthetic_dem.data.min():.1f} m")
    print(f"  Max: {synthetic_dem.data.max():.1f} m")
    print(f"  Mean: {synthetic_dem.data.mean():.1f} m")
    print(f"  Std Dev: {synthetic_dem.data.std():.1f} m")

    # Calculate terrain gradients (simple method)
    grad_y, grad_x = np.gradient(synthetic_dem.data)
    slope = np.degrees(np.arctan(np.sqrt(grad_x**2 + grad_y**2)))

    print(f"\nTerrain slope analysis:")
    print(f"  Min slope: {slope.min():.2f}°")
    print(f"  Max slope: {slope.max():.2f}°")
    print(f"  Mean slope: {slope.mean():.2f}°")

    # Visualization: Comparison histograms (skip for performance)
    if not SKIP_VISUALIZATIONS:
        fig = make_subplots(
            rows=1,
            cols=2,
            subplot_titles=("Flat DEM Distribution", "Synthetic Terrain Distribution"),
        )

        fig.add_trace(
            go.Histogram(
                x=flat_dem.data.flatten(),
                nbinsx=30,
                name="Flat DEM",
                marker_color="steelblue",
            ),
            row=1,
            col=1,
        )

        fig.add_trace(
            go.Histogram(
                x=synthetic_dem.data.flatten(),
                nbinsx=30,
                name="Synthetic Terrain",
                marker_color="coral",
            ),
            row=1,
            col=2,
        )

        fig.update_xaxes(title_text="Elevation (m)", row=1, col=1)
        fig.update_xaxes(title_text="Elevation (m)", row=1, col=2)
        fig.update_yaxes(title_text="Count", row=1, col=1)
        fig.update_layout(height=400, showlegend=True)

        if SHOW_PLOTS:
            fig.show()
        else:
            fig.write_html(str(OUTPUT_DIR / "terrain_demo.html"))

        # Visualization: Slope map with downsampling
        max_size = 500
        stride_slope = max(1, max(slope.shape) // max_size)
        z_slope_display = slope[::stride_slope, ::stride_slope]

        fig_slope = go.Figure(
            data=go.Heatmap(
                z=z_slope_display,
                colorscale="Reds",
                name="Slope",
                colorbar=dict(title="Slope (°)"),
            )
        )

        fig_slope.update_layout(
            title="Terrain Slope Map",
            xaxis_title="Longitude index",
            yaxis_title="Latitude index",
            height=500,
        )

        if SHOW_PLOTS:
            fig_slope.show()


def demo_horizon_computation() -> None:
    """Demonstrate horizon computation."""
    print("\n" + "=" * 60)
    print("Horizon Computation")
    print("=" * 60)

    # Create synthetic terrain
    dem = create_synthetic_terrain(
        lat_min=-3,
        lat_max=3,
        lon_min=-3,
        lon_max=3,
        base_elevation=500,
        amplitude=500,
        wavelength_km=40,
        resolution_arcsec=120,
        seed=456,
    )

    print(f"\nDEM for horizon analysis:")
    print(f"  Shape: {dem.data.shape}")
    print(f"  Elevation range: {dem.data.min():.1f} to {dem.data.max():.1f} m")

    # Compute horizon from center position
    center_idx_lat = dem.data.shape[0] // 2
    center_idx_lon = dem.data.shape[1] // 2

    # Observer height above ground
    observer_height = 100.0

    print(f"\nObserver position:")
    print(f"  Grid indices: ({center_idx_lat}, {center_idx_lon})")
    print(f"  Elevation: {dem.data[center_idx_lat, center_idx_lon]:.1f} m")
    print(f"  Observer height: {observer_height} m")

    # Compute horizon
    try:
        horizon = compute_horizon(
            dem=dem,
            observer_lat_idx=center_idx_lat,
            observer_lon_idx=center_idx_lon,
            observer_height=observer_height,
        )

        print(f"\nHorizon computation successful!")
        print(f"  Horizon azimuth angles: {horizon.azimuth_angles.shape}")
        if hasattr(horizon, "elevation_angles"):
            print(f"  Horizon elevation angles: {horizon.elevation_angles.shape}")
            print(f"  Min elevation: {horizon.elevation_angles.min():.2f}°")
            print(f"  Max elevation: {horizon.elevation_angles.max():.2f}°")

    except Exception as e:
        print(f"\nHorizon computation note: {type(e).__name__}")
        print(f"  This is expected if compute_horizon requires specific parameters")


def main() -> None:
    """Run all demonstrations."""
    print("\n" + "=" * 60)
    print("Terrain Module Demonstration")
    print("=" * 60)

    demo_flat_dem()
    demo_synthetic_terrain()
    demo_terrain_analysis()
    demo_horizon_computation()

    print("\n" + "=" * 60)
    print("Demonstration Complete")
    print("=" * 60)


OUTPUT_DIR = Path("docs/_static/images/examples")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

if __name__ == "__main__":
    main()

Running the Example

python examples/terrain_demo.py

See Also

• ins_gnss_navigation - Navigation applications

• coordinate_systems - Coordinate transformations

• reference_frame_advanced - Reference frames