Data Loading Module

Functions for loading geographic data from various sources.

DEM Loading

src.terrain.data_loading.load_filtered_hgt_files(dem_dir, min_latitude=None, max_latitude=None, min_longitude=None, max_longitude=None, bbox=None, pattern='*.hgt')[source]

Load SRTM HGT files filtered by latitude/longitude range.

Works globally with standard SRTM naming convention. Filters files before loading to reduce memory usage for large DEM directories.

Parameters:

  • dem_dir (Path | str) – Directory containing HGT files

  • min_latitude (int) – Southern bound (e.g., -45 for S45, 42 for N42)

  • max_latitude (int) – Northern bound (e.g., 60 for N60)

  • min_longitude (int) – Western bound (e.g., -120 for W120)

  • max_longitude (int) – Eastern bound (e.g., 30 for E30)

  • bbox (tuple) – Bounding box as (west, south, east, north) tuple. If provided, overrides individual min/max parameters. Uses standard GIS convention: (min_lon, min_lat, max_lon, max_lat).

  • pattern (str) – File pattern to match (default: *.hgt)

Returns:

Tuple of (merged_dem, transform)

Raises:

ValueError – If no matching files found after filtering

Return type:

tuple[ndarray, Affine]

Example

>>> # Load only tiles in Michigan area using individual params
>>> dem, transform = load_filtered_hgt_files(
...     "/path/to/srtm",
...     min_latitude=41, max_latitude=47,
...     min_longitude=-90, max_longitude=-82
... )

>>> # Same area using bbox (west, south, east, north)
>>> dem, transform = load_filtered_hgt_files(
...     "/path/to/srtm",
...     bbox=(-90, 41, -82, 47)
... )

>>> # Alps region (Switzerland/Austria)
>>> dem, transform = load_filtered_hgt_files(
...     "/path/to/srtm",
...     bbox=(5, 45, 15, 48)
... )

Load SRTM HGT files with latitude filtering.

src.terrain.data_loading.save_score_grid(file_path, data, transform=None, data_key='data', **metadata)[source]

Save georeferenced raster data to an NPZ file.

Creates an NPZ file compatible with load_score_grid(). The transform is stored as a 6-element array that can be reconstructed as an Affine.

Parameters:

  • file_path (Path) – Output path for .npz file

  • data (ndarray) – 2D numpy array with raster data

  • transform (Affine | None) – Optional Affine transform for georeferencing

  • data_key (str) – Key name for the data array (default: “data”)

  • **metadata – Additional key=value pairs to store in the file

Returns:

Path to the saved file

Return type:

Path

Example

>>> from rasterio import Affine
>>> scores = compute_sledding_scores(dem)
>>> transform = Affine.translation(-83.5, 42.5) * Affine.scale(0.01, -0.01)
>>> save_score_grid("scores.npz", scores, transform, crs="EPSG:4326")

>>> # Load it back
>>> loaded_scores, loaded_transform = load_score_grid("scores.npz")
src.terrain.data_loading.find_score_file(name, search_dirs=None, subdirs=None)[source]

Search for a score file in common locations.

Useful for finding pre-computed score files that may be in various locations depending on how the pipeline was run.

Parameters:

  • name (str) – Base filename to search for (e.g., “sledding_scores.npz”)

  • search_dirs (list[Path]) – List of directories to search. Defaults to common locations.

  • subdirs (list[str]) – Subdirectories to check within each search_dir (e.g., [“sledding”, “xc_skiing”])

Returns:

Path to found file, or None if not found

Return type:

Path | None

Example

>>> path = find_score_file("sledding_scores.npz",
...                         search_dirs=[Path("docs/images"), Path("output")],
...                         subdirs=["sledding", ""])
>>> if path:
...     scores, transform = load_score_grid(path)

Score Data

src.terrain.scoring.compute_sledding_score(snow_depth, slope, coverage_months, roughness)[source]

Compute overall sledding suitability score.

Combines multiple factors using a multiplicative model where: 1. Deal breakers → immediate zero score 2. Base score = snow_score × slope_score × coverage_score 3. Final score = base_score × synergy_bonus

The multiplicative approach ensures that poor performance in any one factor significantly reduces the overall score, while synergies can boost exceptional combinations.

Parameters:

  • snow_depth – Snow depth in inches (scalar or array)

  • slope – Terrain slope in degrees (scalar or array)

  • coverage_months – Months of snow coverage (scalar or array)

  • roughness – Elevation std dev in meters (scalar or array) - physical terrain roughness

Returns:

Score(s) in range [0, ~1.5], same shape as inputs (Can exceed 1.0 due to synergy bonuses)

Return type:

ndarray

Example

>>> # Perfect conditions
>>> compute_sledding_score(12.0, 9.0, 4.0, 2.0)
1.12  # High score with bonuses
>>> # Deal breaker (too steep)
>>> compute_sledding_score(12.0, 45.0, 4.0, 2.0)
0.0  # Too steep (>40°)
>>> # Deal breaker (too rough)
>>> compute_sledding_score(12.0, 10.0, 4.0, 8.0)
0.0  # Too rough (>6m)

Used in Snow Integration: Sledding Location Analysis.

Example:

score = compute_sledding_score(
    terrain,
    depth_weight=0.4,
    coverage_weight=0.3,
    slope_weight=0.3
)

Roads

src.terrain.roads.get_roads_tiled(bbox, road_types=None, tile_size=2.0, retry_count=1, retry_delay=2.0, force_refresh=False)[source]

Fetch roads for large areas by tiling and merging results.

Works globally - anywhere OpenStreetMap has road coverage. Automatically tiles large bounding boxes to respect Overpass API limits and handles retries for failed tiles.

Parameters:

  • bbox (Tuple[float, float, float, float]) – (south, west, north, east) in WGS84 - works for any location worldwide

  • road_types (List[str] | None) – OSM highway tags to include. Default: [‘motorway’, ‘trunk’, ‘primary’]

  • tile_size (float) – Size of tiles in degrees (default: 2.0° recommended for Overpass API)

  • retry_count (int) – Number of retries for failed tiles (default: 1)

  • retry_delay (float) – Delay in seconds before retrying (default: 2.0)

  • force_refresh (bool) – Force fresh fetch, skip cache (default: False)

Returns:

Merged GeoJSON FeatureCollection with all road segments from all tiles, or empty FeatureCollection if all fetches fail.

Return type:

Dict[str, Any] | None

Example

>>> # Large area covering multiple US states
>>> bbox = (37.0, -95.0, 45.0, -85.0)
>>> roads = get_roads_tiled(bbox, road_types=["motorway", "trunk"])
>>> print(f"Found {len(roads['features'])} road segments")

>>> # Works globally - Alps region
>>> alps_bbox = (45.5, 6.0, 47.5, 14.0)
>>> roads = get_roads_tiled(alps_bbox, tile_size=1.5)

Fetch roads from OpenStreetMap via Overpass API.

src.terrain.roads.add_roads_layer(terrain, roads_geojson, bbox, resolution=30.0, road_width_pixels=3)[source]

Add roads as a data layer to terrain with automatic coordinate alignment.

This is the high-level API for road integration. Roads are rasterized to a grid with proper geographic metadata and added as a data layer. The library’s data layer pipeline ensures proper alignment even if terrain is downsampled or reprojected.

To color roads, use the multi-overlay color mapping system:

roads_geojson = get_roads(bbox)
terrain.add_roads_layer(terrain, roads_geojson, bbox, road_width_pixels=3)
terrain.set_multi_color_mapping(
    base_colormap=lambda dem: elevation_colormap(dem, 'michigan'),
    base_source_layers=['dem'],
    overlays=[{
        'colormap': road_colormap,
        'source_layers': ['roads'],
        'priority': 10,
    }],
)
terrain.compute_colors()
Parameters:

  • terrain – Terrain object (must have DEM data layer)

  • roads_geojson (Dict[str, Any]) – GeoJSON FeatureCollection with road LineStrings

  • bbox (Tuple[float, float, float, float]) – Bounding box as (south, west, north, east) in WGS84 degrees

  • resolution (float) – Pixel size in meters for rasterization (default: 30.0)

  • road_width_pixels (int) – Width of roads in raster pixels (default: 3). Higher values make roads more visible. At 30m resolution, 3 pixels ≈ 90m visual width.

Raises:

ValueError – If terrain missing DEM data layer

Return type:

None

Add roads as a terrain data layer. Used in Combined Render: Full-Featured Example.

src.terrain.roads.rasterize_roads_to_layer(roads_geojson, bbox, resolution=30.0, road_width_pixels=3)[source]

Rasterize GeoJSON roads to a layer grid with proper geographic transform.

Converts vector road data (GeoJSON LineStrings) to a raster grid where each pixel represents road presence/type. The result includes an Affine transform in WGS84 (EPSG:4326) coordinates for proper geographic alignment.

This is the key function that treats roads as a data layer - the output can be added to terrain via add_data_layer() and will automatically align to the DEM through reprojection and resampling.

Parameters:

  • roads_geojson (Dict[str, Any]) – GeoJSON FeatureCollection with road LineStrings

  • bbox (Tuple[float, float, float, float]) – Bounding box as (south, west, north, east) in WGS84 degrees

  • resolution (float) – Pixel size in meters (default: 30.0). At Detroit latitude, ~30m/pixel gives good detail without excessive memory use.

  • road_width_pixels (int) – Width of roads in raster pixels (default: 3). Draws roads with thickness instead of 1-pixel lines. Use 1 for thin roads, 3-5 for more visible roads. At 30m resolution, 3 pixels ≈ 90m visual width.

Returns:

  • road_grid: 2D array of uint8, values 0=no road, 1-4=road type (motorway > trunk > primary > secondary)

  • road_transform: rasterio.Affine transform mapping pixels to WGS84 coordinates

Return type:

Tuple of

src.terrain.roads.smooth_road_vertices(vertices, road_mask, y_valid, x_valid, smoothing_radius=2)[source]

Smooth Z coordinates of mesh vertices that are on roads.

Reconstructs a 2D Z grid from vertex positions, applies 2D Gaussian smoothing in grid space, then writes smoothed values back to road vertices only. Non-road vertices are unchanged.

This ensures smoothing follows the spatial layout of the road rather than the arbitrary vertex index order.

Parameters:

  • vertices (ndarray) – Mesh vertex positions (N, 3) array with [x, y, z] coords

  • road_mask (ndarray) – 2D array (H, W) where >0.5 indicates road pixels

  • y_valid (ndarray) – Array (N,) of y indices mapping vertices to road_mask rows

  • x_valid (ndarray) – Array (N,) of x indices mapping vertices to road_mask columns

  • smoothing_radius (int) – Gaussian smoothing sigma (default: 2, use 0 to disable)

Returns:

Modified vertices array with smoothed Z values on roads. X and Y coordinates are never modified.

Return type:

ndarray

src.terrain.roads.offset_road_vertices(vertices, road_mask, y_valid, x_valid, offset=0.0)[source]

Offset Z coordinates of mesh vertices that are on roads by a fixed amount.

A simpler alternative to smooth_road_vertices. Raises or lowers all road vertices by a constant offset, making roads visually distinct from terrain.

Parameters:

  • vertices (ndarray) – Mesh vertex positions (N, 3) array with [x, y, z] coords

  • road_mask (ndarray) – 2D array (H, W) where >0.5 indicates road pixels

  • y_valid (ndarray) – Array (N,) of y indices mapping vertices to road_mask rows

  • x_valid (ndarray) – Array (N,) of x indices mapping vertices to road_mask columns

  • offset (float) – Z offset to apply to road vertices. Positive = raise, negative = lower. Default: 0.0 (no change)

Returns:

Modified vertices array with offset Z values on roads. X and Y coordinates are never modified.

Return type:

ndarray

src.terrain.roads.smooth_road_mask(road_mask, sigma=1.0)[source]

Apply Gaussian blur to road mask for anti-aliased edges (GPU-accelerated).

The Bresenham line algorithm creates stair-step (aliased) edges. Applying Gaussian smoothing creates soft anti-aliased boundaries that render more smoothly, especially after the mask goes through resampling.

Uses PyTorch GPU acceleration when available (6x speedup on CUDA).

Parameters:

  • road_mask (ndarray) – 2D array of road values (0=no road, >0=road)

  • sigma (float) – Gaussian blur sigma in pixels (default: 1.0). Higher values = softer edges. Typical range: 0.5-2.0. - 0.5: Minimal softening - 1.0: Standard anti-aliasing (recommended) - 2.0: Very soft/blurry edges

Returns:

Smoothed road mask as float32 array. Values are now continuous (not binary) and may need thresholding if binary mask is needed.

Return type:

ndarray

Anti-alias road edges.

Diagnostics

src.terrain.diagnostics.generate_rgb_histogram(image_path, output_path)[source]

Generate and save an RGB histogram of a rendered image.

Creates a figure with histograms for each color channel (R, G, B) overlaid on the same axes with transparency. Useful for analyzing color balance and distribution in rendered outputs.

Parameters:

  • image_path (Path) – Path to the rendered image (PNG, JPEG, etc.)

  • output_path (Path) – Path to save the histogram image

Returns:

Path to saved histogram image, or None if failed

Return type:

Path | None

src.terrain.diagnostics.generate_luminance_histogram(image_path, output_path)[source]

Generate and save a luminance (B&W) histogram of a rendered image.

Shows distribution of brightness values with annotations for pure black and pure white pixel counts. Useful for checking exposure and clipping in rendered outputs, especially for print preparation.

Parameters:

  • image_path (Path) – Path to the rendered image (PNG, JPEG, etc.)

  • output_path (Path) – Path to save the histogram image

Returns:

Path to saved histogram image, or None if failed

Return type:

Path | None

src.terrain.diagnostics.plot_wavelet_diagnostics(original, denoised, output_path, title_prefix='Wavelet Denoising', nodata_value=nan, profile_row=None, cmap='terrain')[source]

Generate diagnostic plots showing wavelet denoising effects.

Creates a multi-panel figure showing: - Original DEM - Denoised DEM - Difference (noise removed) - Cross-section profile comparison

Parameters:

  • original (ndarray) – Original DEM before denoising

  • denoised (ndarray) – DEM after wavelet denoising

  • output_path (Path) – Path to save the diagnostic plot

  • title_prefix (str) – Prefix for plot titles

  • nodata_value (float) – Value treated as no data

  • profile_row (int | None) – Row index for cross-section (default: middle row)

  • cmap (str) – Colormap for elevation visualization

Returns:

Path to saved diagnostic plot

Return type:

Path

src.terrain.diagnostics.generate_upscale_diagnostics(original, upscaled, output_dir, prefix='score_upscale', scale=4, method='unknown', nodata_value=nan, cmap='viridis')[source]

Generate upscale diagnostic plot.

Parameters:

  • original (ndarray) – Original score grid before upscaling

  • upscaled (ndarray) – Score grid after upscaling

  • output_dir (Path) – Directory to save diagnostic plot

  • prefix (str) – Filename prefix

  • scale (int) – Upscaling factor used

  • method (str) – Upscaling method name

  • nodata_value (float) – Value treated as no data

  • cmap (str) – Colormap for score visualization

Returns:

Path to saved diagnostic plot

Return type:

Path