Core Module

The core module provides the main Terrain class and essential functions for loading and processing DEM data.

Terrain Class

class src.terrain.core.Terrain(dem_data, dem_transform, dem_crs='EPSG:4326', cache_dir='terrain_cache', logger=None)

Bases: object

Core class for managing Digital Elevation Model (DEM) data and terrain operations.

Handles loading, transforming, and visualizing terrain data from raster sources. Supports coordinate reprojection, downsampling, color mapping, and 3D mesh generation for Blender visualization. Uses efficient caching to avoid recomputation of transforms.

Parameters:

• dem_data (np.ndarray)

• dem_transform (rasterio.Affine)

• dem_crs (str)

• cache_dir (str)

• logger (Optional[logging.Logger])

dem_shape

Shape of DEM array as (height, width).

Type:

tuple

dem_transform

Affine transform for geographic coordinates.

Type:

rasterio.Affine

data_layers

Dictionary of data layers (DEM, overlays, derived data).

Type:

dict

transforms

List of transform functions to apply.

Type:

list

vertices

Vertex positions for generated mesh.

Type:

np.ndarray

vertex_colors

RGBA colors for mesh vertices.

Type:

np.ndarray

Examples

>>> dem_data = np.random.rand(100, 100) * 1000
>>> transform = rasterio.Affine.identity()
>>> terrain = Terrain(dem_data, transform, dem_crs='EPSG:4326')
>>> terrain.apply_transforms()
>>> mesh = terrain.create_mesh(scale_factor=100.0)

Examples

Basic usage:

from src.terrain.core import Terrain, load_dem_files

dem, transform = load_dem_files('path/to/tiles')
terrain = Terrain(dem, transform)
terrain.apply_transforms()
mesh = terrain.create_mesh()

See Great Lakes Elevation Visualization for a complete example.

__init__(dem_data, dem_transform, dem_crs='EPSG:4326', cache_dir='terrain_cache', logger=None)

Initialize terrain from DEM data.

Parameters:

• dem_data (np.ndarray) – DEM array of shape (height, width) containing elevation values. Integer types are converted to float32. Must be 2D.

• dem_transform (rasterio.Affine) – Affine transform mapping pixel coordinates to geographic coordinates.

• dem_crs (str) – Coordinate reference system in EPSG format (default: ‘EPSG:4326’).

• cache_dir (str) – Directory for caching computations (default: ‘terrain_cache’).

• logger (logging.Logger, optional) – Logger instance for diagnostic output.

Raises:

• TypeError – If dem_data is not a numpy array or has unsupported dtype.

• ValueError – If dem_data is not 2D.

Return type:

None

visualize_dem(layer='dem', use_transformed=False, title=None, cmap='terrain', percentile_clip=True, clip_percentiles=(1, 99), max_pixels=500000, show_histogram=True)

Create diagnostic visualization of any terrain data layer.

Parameters:

• layer (str) – Name of data layer to visualize (default: ‘dem’)

• use_transformed (bool) – Whether to use transformed or original data (default: False)

• title (str) – Plot title (default: auto-generated based on layer)

• cmap (str) – Matplotlib colormap

• percentile_clip (bool) – Whether to clip extreme values

• clip_percentiles (tuple) – Tuple of (min, max) percentiles to clip (default: (1, 99))

• max_pixels (int) – Maximum number of pixels for subsampling

• show_histogram (bool) – Whether to show the histogram panel (default: True)

Return type:

None

add_transform(transform_func)

Add a transform function to the processing pipeline.

Parameters:

transform_func (callable) – Function that transforms DEM data. Should accept (dem_array: np.ndarray) and return transformed np.ndarray.

Returns:

Modifies internal transforms list in place.

Return type:

None

Examples

>>> terrain.add_transform(lambda dem: gaussian_filter(dem, sigma=2))
>>> terrain.apply_transforms()

add_data_layer(name, data, transform=None, crs=None, target_crs=None, target_layer=None, same_extent_as=None, resampling=Resampling.bilinear, nodata=None)

Add a data layer, optionally reprojecting to match another layer.

Stores data with geographic metadata (CRS and transform). Can automatically reproject and resample to match an existing layer’s grid for multi-layer analysis.

Parameters:

• name (str) – Unique name for this data layer (e.g., ‘dem’, ‘elevation’, ‘slope’).

• data (np.ndarray) – 2D array of data values, shape (height, width).

• transform (rasterio.Affine, optional) – Affine transform mapping pixel to geographic coords. Required unless same_extent_as is specified.

• crs (str, optional) – Coordinate reference system in EPSG format (e.g., ‘EPSG:4326’). Required unless same_extent_as is specified (inherits from reference layer).

• target_crs (str, optional) – Target CRS to reproject to. If None and target_layer specified, uses target layer’s CRS. If None and no target, uses input crs.

• target_layer (str, optional) – Name of existing layer to match grid and CRS. If specified, data is automatically reprojected and resampled to align.

• same_extent_as (str, optional) – Name of existing layer whose geographic extent this data covers. When specified, transform and CRS are automatically calculated from the reference layer’s bounds and the data’s shape. This is useful when score grids or overlays cover the same area as the DEM but at different resolutions. Implies target_layer if not specified.

• resampling (rasterio.enums.Resampling) – Resampling method for reprojection (default: Resampling.bilinear). See rasterio docs for options.

• nodata (float, optional) – No-data sentinel value used during reprojection. Source pixels with this value are treated as missing and won’t influence bilinear interpolation neighbours. Destination pixels with no source coverage are filled with this value. Defaults to np.nan for floating-point arrays and 0 for integer arrays.

Returns:

Modifies internal data_layers dictionary.

Return type:

None

Raises:

• KeyError – If target_layer or same_extent_as layer doesn’t exist.

• ValueError – If neither transform nor same_extent_as is provided.

Examples

>>> # Add elevation data with native CRS
>>> terrain.add_data_layer('dem', dem_array, transform, 'EPSG:4326')

>>> # Add overlay data, reproject to match DEM
>>> terrain.add_data_layer('landcover', lc_array, lc_transform, 'EPSG:3857',
...                        target_layer='dem')

>>> # Add score data that covers the same extent as DEM (automatic transform)
>>> terrain.add_data_layer('score', score_array, same_extent_as='dem')

>>> # Use nearest-neighbor for categorical data
>>> terrain.add_data_layer('zones', zone_array, zone_transform, 'EPSG:4326',
...                        target_layer='dem', resampling=Resampling.nearest)

get_bbox_wgs84(layer='dem')

Get bounding box in WGS84 coordinates (EPSG:4326).

Returns bbox in standard format used by OSM, web mapping APIs, etc. Handles reprojection from any source CRS back to WGS84.

Parameters:

layer (str) – Layer name to get bbox for (default: “dem”)

Returns:

Tuple of (south, west, north, east) in WGS84 degrees

Raises:

KeyError – If specified layer doesn’t exist

Return type:

tuple[float, float, float, float]

Example

>>> terrain = Terrain(dem_data, transform, dem_crs="EPSG:32617")  # UTM
>>> south, west, north, east = terrain.get_bbox_wgs84()
>>> print(f"Bounds: {south:.4f}°N to {north:.4f}°N, {west:.4f}°E to {east:.4f}°E")

compute_data_layer(name, source_layer, compute_func, transformed=False, cache_key=None)

Compute a new data layer from an existing one using a transformation function.

Allows creating derived layers (e.g., slope, aspect, hillshade) from existing data. Results are stored as new layer and optionally cached.

Parameters:

• name (str) – Name for the computed layer.

• source_layer (str) – Name of existing source layer to compute from.

• compute_func (Callable) – Function that accepts source array (np.ndarray) and returns computed array (np.ndarray). Can return same or different shape.

• transformed (bool) – If True, use already-transformed source data; if False, use original source data (default: False).

• cache_key (str, optional) – Custom cache identifier. If None, auto-generated from layer name and function name.

Returns:

The computed layer data array.

Return type:

np.ndarray

Raises:

• KeyError – If source_layer doesn’t exist.

• ValueError – If transformed=True but source hasn’t been transformed.

Examples

>>> # Compute slope from DEM using scipy
>>> from scipy.ndimage import sobel
>>> slope = terrain.compute_data_layer(
...     'slope', 'dem',
...     lambda dem: np.sqrt(sobel(dem, axis=0)**2 + sobel(dem, axis=1)**2)
... )

>>> # Compute hill-shade visualization
>>> from scipy.ndimage import gaussian_filter
>>> hillshade = terrain.compute_data_layer(
...     'hillshade', 'dem',
...     lambda dem: np.clip(gaussian_filter(dem, 2) * 0.5, 0, 1)
... )

>>> # Compute from transformed (downsampled) data
>>> downsampled_slope = terrain.compute_data_layer(
...     'slope_downsampled', 'dem',
...     lambda dem: np.gradient(dem)[0],
...     transformed=True
... )

apply_transforms(cache=False)

Apply all transforms to all data layers with optional caching.

Processes each data layer through the transform pipeline. Results are cached to avoid recomputation. Transforms are applied in order.

Parameters:

cache (bool) – Whether to cache results (default: False).

Returns:

Updates internal data_layers with ‘transformed_data’ for each layer.

Return type:

None

Examples

>>> terrain.add_transform(flip_raster(axis='horizontal'))
>>> terrain.apply_transforms(cache=True)
>>> dem_data = terrain.data_layers['dem']['transformed_data']

configure_for_target_vertices(target_vertices, method='average')

Configure downsampling to achieve approximately target_vertices.

This method calculates the appropriate zoom_factor to achieve a desired vertex count for mesh generation. It provides a more intuitive API than manually calculating zoom_factor from the original DEM shape.

Parameters:

• target_vertices (int) – Desired vertex count for final mesh (e.g., 500_000)

• method (str) – Downsampling method (default: “average”) - “average”: Area averaging - best for DEMs, no overshoot - “lanczos”: Lanczos resampling - sharp, minimal aliasing - “cubic”: Cubic spline interpolation - “bilinear”: Bilinear interpolation - safe fallback

Returns:

Calculated zoom_factor that was added to transforms

Raises:

ValueError – If target_vertices is invalid

Return type:

float

Example

terrain = Terrain(dem, transform) zoom = terrain.configure_for_target_vertices(500_000, method=”average”) print(f”Calculated zoom_factor: {zoom:.4f}”) terrain.apply_transforms() mesh = terrain.create_mesh(scale_factor=400.0)

detect_water_highres(slope_threshold=0.01, fill_holes=True, scale_factor=0.0001)

Detect water bodies on high-resolution DEM before downsampling.

This method properly handles water detection by: 1. Applying all NON-downsampling transforms to DEM (reproject, flip, etc.) 2. Detecting water on the high-resolution transformed DEM 3. Downsampling the water mask to match the final terrain resolution

This prevents false water detection that occurs when detecting on downsampled DEMs.

Parameters:

• slope_threshold (float) – Maximum slope magnitude to classify as water (default: 0.01)

• fill_holes (bool) – Whether to apply morphological hole filling (default: True)

• scale_factor (float) – Elevation scale factor to unscale before detection (default: 0.0001)

Returns:

Boolean water mask matching transformed_data shape

Return type:

np.ndarray

Raises:

• ValueError – If transforms haven’t been applied yet

• ImportError – If water detection module is not available

Example:

terrain = Terrain(dem, transform, dem_crs="EPSG:4326")
terrain.add_transform(reproject_raster(src_crs="EPSG:4326", dst_crs="EPSG:32617"))
terrain.add_transform(flip_raster(axis="horizontal"))
terrain.add_transform(scale_elevation(scale_factor=0.0001))
terrain.configure_for_target_vertices(target_vertices=1_000_000)
terrain.apply_transforms()  # Includes downsampling

# Detect water on high-res BEFORE it was downsampled
water_mask = terrain.detect_water_highres(slope_threshold=0.01)
# water_mask shape matches downsampled terrain

Note

This method requires that: - Transforms have been applied (apply_transforms() called) - The DEM layer exists in data_layers - Water detection module (src.terrain.water) is available

detect_water(slope_threshold=0.01, fill_holes=True)

Detect water bodies on the transformed (downsampled) DEM.

This is a fast, simple water detection method that works on the already- transformed DEM. It automatically handles elevation scale factor unscaling.

For higher accuracy at the cost of speed, use detect_water_highres() instead, which detects water on the full-resolution DEM before downsampling.

Parameters:

• slope_threshold (float) – Maximum slope magnitude to classify as water (default: 0.01). Water bodies have nearly zero slope.

• fill_holes (bool) – Whether to apply morphological hole filling (default: True)

Returns:

Boolean water mask matching transformed_data shape

Return type:

np.ndarray

Raises:

ValueError – If transforms haven’t been applied yet

Example

`python terrain = Terrain(dem, transform) terrain.apply_transforms() water_mask = terrain.detect_water() terrain.compute_colors(water_mask=water_mask) `

set_color_mapping(color_func, source_layers, *, color_kwargs=None, mask_func=None, mask_layers=None, mask_kwargs=None, mask_threshold=None)

Set up how to map data layers to colors (RGB) and optionally a mask/alpha channel.

Allows flexible color mapping by applying a function to one or more data layers. Optionally applies a separate mask function for transparency/alpha channel control. Color mapping is applied during mesh creation with compute_colors().

Parameters:

• color_func (Callable) – Function that accepts N data arrays (one per source_layers) and returns colored array of shape (H, W, 3) for RGB or (H, W, 4) for RGBA. Values should be in range [0, 1] for 8-bit output.

• source_layers (list[str]) – Names of data layers to pass to color_func, in order. E.g., [‘dem’] for single layer or [‘red’, ‘green’, ‘blue’] for composite.

• color_kwargs (dict, optional) – Additional keyword arguments passed to color_func.

• mask_func (Callable, optional) – Function producing alpha/mask values (0-1) for transparency. Takes layer arrays as input. If omitted, fully opaque.

• mask_layers (list[str] | str, optional) – Layer names for mask_func. If None, uses source_layers. Single string converted to list.

• mask_kwargs (dict, optional) – Additional keyword arguments for mask_func.

• mask_threshold (float, optional) – If mask_func is threshold-based, convenience parameter for threshold value (implementation-dependent).

Returns:

Modifies internal color mapping configuration.

Return type:

None

Raises:

ValueError – If source_layers or mask_layers refer to non-existent layers.

Examples

>>> # Single-layer elevation with viridis colormap
>>> from matplotlib.cm import viridis
>>> terrain.set_color_mapping(
...     lambda dem: viridis(dem / dem.max()),
...     ['dem']
... )

>>> # RGB composite from three layers
>>> terrain.set_color_mapping(
...     lambda r, g, b: np.stack([r, g, b], axis=-1),
...     ['red_band', 'green_band', 'blue_band']
... )

>>> # Elevation with water transparency mask
>>> terrain.set_color_mapping(
...     lambda dem: elevation_colormap(dem),
...     ['dem'],
...     mask_func=lambda dem: (dem > 0).astype(float),
...     mask_layers=['dem']
... )

>>> # Hillshade with elevation colors and slope transparency
>>> terrain.set_color_mapping(
...     lambda dem: dem_colors,
...     ['dem'],
...     mask_func=lambda dem: 1 - np.clip(np.gradient(dem)[0], 0, 1),
...     mask_layers=['dem']
... )

set_blended_color_mapping(base_colormap, base_source_layers, overlay_colormap, overlay_source_layers, overlay_mask, *, base_color_kwargs=None, overlay_color_kwargs=None)

Apply two different colormaps based on a spatial mask (hard transition).

Uses base colormap for most of terrain and overlay colormap for masked zones. This is useful for showing different data types in different regions, such as elevation colors for general terrain but suitability scores near parks.

Parameters:

• base_colormap (Callable[[ndarray], ndarray]) – Function mapping base layer(s) to RGB/RGBA. Takes N arrays (one per base_source_layers) and returns (H, W, 3) or (H, W, 4).

• base_source_layers (list[str]) – Layer names to pass to base_colormap, e.g., [‘dem’].

• overlay_colormap (Callable[[ndarray], ndarray]) – Function mapping overlay layer(s) to RGB/RGBA. Takes N arrays (one per overlay_source_layers) and returns (H, W, 3) or (H, W, 4).

• overlay_source_layers (list[str]) – Layer names to pass to overlay_colormap, e.g., [‘score’].

• overlay_mask (ndarray) – Boolean array of shape (num_vertices,) indicating where to use overlay colormap. True = overlay, False = base. Use compute_proximity_mask() to create this mask.

• base_color_kwargs (Dict[str, Any] | None) – Optional kwargs passed to base_colormap.

• overlay_color_kwargs (Dict[str, Any] | None) – Optional kwargs passed to overlay_colormap.

Returns:

Modifies internal color mapping to use blended approach.

Return type:

None

Raises:

• ValueError – If source layers don’t exist or overlay_mask has wrong shape.

• RuntimeError – If create_mesh() hasn’t been called yet (needed for mask validation).

Example

>>> from src.terrain.color_mapping import elevation_colormap
>>> # Compute proximity mask for park zones
>>> park_mask = terrain.compute_proximity_mask(
...     park_lons, park_lats, radius_meters=1000
... )
>>> # Set dual colormaps
>>> terrain.set_blended_color_mapping(
...     base_colormap=lambda elev: elevation_colormap(
...         elev, cmap_name="gist_earth"
...     ),
...     base_source_layers=["dem"],
...     overlay_colormap=lambda score: elevation_colormap(
...         score, cmap_name="cool", min_elev=0, max_elev=1
...     ),
...     overlay_source_layers=["score"],
...     overlay_mask=park_mask
... )
>>> terrain.compute_colors()  # Apply the blended mapping

set_multi_color_mapping(base_colormap, base_source_layers, overlays, *, base_color_kwargs=None)

Apply multiple data layers as overlays with different colormaps and priority.

This enables flexible data visualization where multiple geographic features (roads, trails, land use, power lines, etc.) can each be colored independently based on their own colormaps and source layers.

Parameters:

• base_colormap (Callable[[ndarray, ...], ndarray]) – Function mapping base layer(s) to RGB/RGBA. Takes N arrays (one per base_source_layers) and returns (H, W, 3) or (H, W, 4).

• base_source_layers (list[str]) – Layer names for base_colormap, e.g., [‘dem’].

• overlays (list[dict]) – List of overlay specifications. Each overlay dict contains colormap (function returning H,W,3/4), source_layers (list of layer names), colormap_kwargs (optional dict), threshold (value above which overlay applies, default 0.5), and priority (lower = higher priority).

• base_color_kwargs (Dict[str, Any] | None) – Optional kwargs passed to base_colormap.

Returns:

Modifies internal color mapping to use multi-overlay approach.

Return type:

None

Raises:

ValueError – If source layers don’t exist or overlay specs are invalid.

Example

>>> # Base elevation colors with roads and trails overlays
>>> terrain.set_multi_color_mapping(
...     base_colormap=lambda elev: elevation_colormap(elev, "michigan"),
...     base_source_layers=["dem"],
...     overlays=[
...         {
...             "colormap": lambda roads: colormap_roads(roads),
...             "source_layers": ["roads"],
...             "priority": 10,  # High priority roads show on top
...         },
...         {
...             "colormap": lambda trails: colormap_trails(trails),
...             "source_layers": ["trails"],
...             "priority": 20,  # Lower priority
...         },
...     ]
... )
>>> terrain.compute_colors()

compute_colors(water_mask=None)

Compute colors using color_func and optionally mask_func.

Supports three modes: - Standard: Single colormap applied to all vertices - Blended: Two colormaps blended based on proximity mask - Multi-overlay: Multiple overlays with different colormaps and priority

Parameters:

water_mask (np.ndarray, optional) – Boolean water mask in grid space (height × width). For blended mode, water pixels will be colored blue in the final vertex colors. For standard and multi-overlay modes, water detection is handled in create_mesh().

Returns:

RGBA color array.

Return type:

np.ndarray

create_mesh(base_depth=0.2, boundary_extension=True, scale_factor=100.0, height_scale=1.0, center_model=True, verbose=True, detect_water=False, water_slope_threshold=0.5, water_mask=None, two_tier_edge=True, edge_mid_depth=None, edge_base_material='clay', edge_blend_colors=False, smooth_boundary=False, smooth_boundary_window=5, use_catmull_rom=False, catmull_rom_subdivisions=10, use_rectangle_edges=True, use_fractional_edges=True, edge_sample_spacing=0.33)

Create a Blender mesh from transformed DEM data with both performance and control.

Generates vertices from DEM elevation values and faces for connectivity. Optionally creates boundary faces to close the mesh into a solid. Supports coordinate scaling and elevation scaling for visualization. Can optionally detect and apply water bodies to vertex alpha channel for water rendering.

Parameters:

• base_depth (float) – Positive depth offset below minimum surface elevation (default: 0.2). Creates a flat base plane at: min_surface_z - base_depth. Used when boundary_extension=True to create side faces. Positive values extend below surface, negative extend above.

• boundary_extension (bool) – Whether to create side faces around the terrain boundary to close the mesh (default: True). If False, creates open terrain surface.

• scale_factor (float) – Horizontal scale divisor for x/y coordinates (default: 100.0). Higher values produce smaller meshes. E.g., 100 means 100 DEM units = 1 Blender unit.

• height_scale (float) – Multiplier for elevation values (default: 1.0). Vertically exaggerates or reduces terrain features. Values > 1 exaggerate, < 1 flatten.

• center_model (bool) – Whether to center the model at origin (default: True). Centers XY coordinates but preserves absolute Z elevation values.

• verbose (bool) – Whether to log detailed progress information (default: True).

• detect_water (bool) – Whether to detect water bodies and apply to alpha channel (default: False). Uses slope-based detection on transformed DEM.

• water_slope_threshold (float) – Maximum slope magnitude to classify as water (default: 0.5). Only used if detect_water=True and water_mask is None.

• water_mask (np.ndarray) – Pre-computed boolean water mask (True=water, False=land). If provided, this mask is used instead of computing water detection. Allows water detection on unscaled DEM before elevation scaling transforms.

• two_tier_edge (bool) – Enable two-tier edge extrusion (default: True). Creates a small colored edge near the surface with a larger uniform base below.

• edge_mid_depth (float) – Positive depth offset below surface for middle tier (default: auto-calculated). If None, automatically set to base_depth * 0.25 (typically 0.05). Positive values extend below surface, negative extend above surface.

• edge_base_material (str | tuple) – Material for base layer (default: “clay”). Either a preset name (“clay”, “obsidian”, “chrome”, “plastic”, “gold”, “ivory”) or an RGB tuple (0-1 range).

• edge_blend_colors (bool) – Blend surface colors to mid tier in two-tier mode (default: False). If False, mid tier uses base_material color for sharp transition between mesh and edge.

• smooth_boundary (bool) – Apply smoothing to boundary points to eliminate stair-step edges (default: False). Useful for smoother mesh transitions when using two-tier edge.

• smooth_boundary_window (int) – Window size for boundary smoothing (default: 5). Larger values produce more smoothing. Only used when smooth_boundary=True.

• use_catmull_rom (bool) – Use Catmull-Rom curve fitting for smooth boundary geometry (default: False). When enabled, eliminates pixel-grid staircase pattern entirely by fitting smooth parametric curve through boundary points.

• catmull_rom_subdivisions (int) – Number of interpolated points per boundary segment when using Catmull-Rom curves (default: 10). Higher values = smoother curve but more vertices. Only used when use_catmull_rom=True.

• use_rectangle_edges (bool) – Use rectangle-edge sampling instead of morphological boundary detection (default: True). ~150x faster than morphological detection. Ideal for rectangular DEMs from raster sources. Uses original DEM shape to generate clean, regularly-sampled edge vertices.

• use_fractional_edges (bool) – Use fractional edge coordinates that preserve projection curvature (default: True). When enabled with use_rectangle_edges, surface tier aligns with mesh boundary (no gap) while mid/base tiers use fractional X,Y positions to follow smooth WGS84→UTM projection curves. Creates geographically accurate curved edges that connect seamlessly to terrain.

• edge_sample_spacing (float) – Pixel spacing for edge/skirt vertices (default: 0.33). Lower values = denser skirt (smoother but more memory). 0.33 = 3x denser than mesh edge. 1.0 = same density as mesh edge. 2.0 = half density (good for large meshes to avoid OOM).

Returns:

The created terrain mesh object, or None if creation failed.

Return type:

bpy.types.Object | None

Raises:

ValueError – If transformed DEM layer is not available (apply_transforms() not called).

Examples

# Default two-tier edge with smooth fractional edges and red clay base mesh = terrain.create_mesh(boundary_extension=True)

# Single-tier edge (backwards compatible) mesh = terrain.create_mesh(two_tier_edge=False)

# Two-tier with gold base material mesh = terrain.create_mesh(edge_base_material=”gold”)

# Two-tier with deeper edge (1 unit below min surface, 0.2 units below surface for mid) mesh = terrain.create_mesh(

base_depth=1.0, # Deep base (1 unit below min surface) edge_mid_depth=0.2 # Deeper colored edge (0.2 units below surface)

)

# Two-tier with custom RGB color mesh = terrain.create_mesh(edge_base_material=(0.6, 0.55, 0.5))

# Disable fractional edges for simple rectangular boundary mesh = terrain.create_mesh(use_fractional_edges=False)

geo_to_mesh_coords(lon, lat, elevation_offset=0.0, input_crs='EPSG:4326')

Convert geographic coordinates to Blender mesh coordinates.

Transforms lon/lat points through the same pipeline used for the terrain mesh, producing coordinates that align with the rendered terrain. Useful for placing markers, labels, or other objects at geographic locations on the terrain.

Parameters:

• lon (ndarray | float) – Longitude(s) or X coordinate(s). Can be a single float or array.

• lat (ndarray | float) – Latitude(s) or Y coordinate(s). Can be a single float or array (must match lon shape).

• elevation_offset (float) – Height above terrain surface in mesh units (default: 0.0). Positive values place points above the terrain.

• input_crs (str) – CRS of input coordinates (default: “EPSG:4326” for WGS84 lon/lat). Will be reprojected to match the transformed DEM’s CRS if different.

Returns:

Tuple of (x, y, z) arrays in Blender mesh coordinates. If single values were passed for lon/lat, returns single float values.

Raises:

• RuntimeError – If create_mesh() has not been called yet (model_params not set).

• ValueError – If the DEM layer has not been transformed yet.

Return type:

tuple[ndarray, ndarray, ndarray]

Example

>>> terrain = Terrain(dem, transform)
>>> terrain.add_transform(reproject_to_utm("EPSG:4326", "EPSG:32617"))
>>> terrain.apply_transforms()
>>> mesh = terrain.create_mesh()
>>> # Get mesh coords for a park at (-83.1, 42.4) in WGS84
>>> x, y, z = terrain.geo_to_mesh_coords(-83.1, 42.4, elevation_offset=0.1)
>>> # Create marker at (x, y, z) in Blender

compute_proximity_mask(lons, lats, radius_meters, input_crs='EPSG:4326', cluster_threshold_meters=None)

Create boolean mask for vertices within radius of geographic points.

Uses KDTree for efficient spatial queries to identify mesh vertices near specified geographic locations (e.g., parks, POIs). Optionally clusters nearby points first to create unified proximity zones.

Parameters:

• lons (ndarray) – Array of longitudes for points of interest.

• lats (ndarray) – Array of latitudes for points of interest (must match lons shape).

• radius_meters (float) – Radius in meters around each point/cluster to include.

• input_crs (str) – CRS of input coordinates (default: “EPSG:4326” for WGS84 lon/lat).

• cluster_threshold_meters (float | None) – Optional distance threshold for clustering nearby points using DBSCAN. Points within this distance are merged into single zones. If None, each point gets its own zone. Useful for merging nearby parks into continuous zones.

Returns:

Boolean array of shape (num_vertices,) where True indicates vertex is within radius of at least one point/cluster.

Raises:

• RuntimeError – If create_mesh() has not been called yet.

• ValueError – If lons and lats have different shapes.

Return type:

ndarray

Example

>>> # Create zones around parks
>>> park_lons = np.array([-83.1, -83.2, -83.15])
>>> park_lats = np.array([42.4, 42.5, 42.45])
>>> mask = terrain.compute_proximity_mask(
...     park_lons, park_lats,
...     radius_meters=1000,
...     cluster_threshold_meters=200  # Merge parks within 200m
... )
>>> # mask is True for vertices within 1km of park clusters

compute_proximity_mask_grid(lons, lats, radius_meters, input_crs='EPSG:4326', cluster_threshold_meters=None)

Create boolean mask for DEM grid cells within radius of geographic points.

Similar to compute_proximity_mask() but works directly on the transformed DEM grid WITHOUT requiring mesh creation. Returns a 2D boolean array matching the transformed DEM shape.

This is useful when you need to compute proximity zones before creating the mesh, avoiding the need for duplicate mesh creation.

Parameters:

• lons (ndarray) – Array of longitudes for points of interest.

• lats (ndarray) – Array of latitudes for points of interest (must match lons shape).

• radius_meters (float) – Radius in meters around each point/cluster to include.

• input_crs (str) – CRS of input coordinates (default: “EPSG:4326” for WGS84 lon/lat).

• cluster_threshold_meters (float | None) – Optional distance threshold for clustering nearby points using DBSCAN. If None, each point gets its own zone.

Returns:

Boolean array of shape (height, width) matching transformed DEM, where True indicates grid cell is within radius of at least one point/cluster.

Raises:

ValueError – If DEM has not been transformed yet or if lons/lats mismatch.

Return type:

ndarray

Example

>>> terrain = Terrain(dem, transform)
>>> terrain.add_transform(reproject_to_utm(...))
>>> terrain.apply_transforms()
>>> # Compute proximity BEFORE mesh creation
>>> park_mask = terrain.compute_proximity_mask_grid(
...     park_lons, park_lats, radius_meters=1000
... )
>>> # Use mask in color computations, then create mesh once

compute_ring_mask_grid(lons, lats, inner_radius_meters, outer_radius_meters, input_crs='EPSG:4326', cluster_threshold_meters=None)

Create boolean mask for ring (annulus) around geographic points.

Creates a ring-shaped mask where pixels are True if they are: - Within outer_radius_meters of a point, AND - Outside inner_radius_meters of all points

This is useful for creating visual outlines around areas of interest.

Parameters:

• lons (ndarray) – Array of longitudes for points of interest.

• lats (ndarray) – Array of latitudes for points of interest (must match lons shape).

• inner_radius_meters (float) – Inner radius of the ring in meters. Pixels closer than this to any point are excluded. Use 0 for filled circle.

• outer_radius_meters (float) – Outer radius of the ring in meters. Pixels farther than this from all points are excluded.

• input_crs (str) – CRS of input coordinates (default: “EPSG:4326” for WGS84 lon/lat).

• cluster_threshold_meters (float | None) – Optional distance threshold for clustering nearby points using DBSCAN. If None, each point gets its own ring.

Returns:

Boolean array of shape (height, width) matching transformed DEM, where True indicates grid cell is within the ring zone.

Raises:

ValueError – If inner_radius > outer_radius or if radii are negative.

Return type:

ndarray

Example

>>> # Create a 100m-wide ring at 2km from each park
>>> ring_mask = terrain.compute_ring_mask_grid(
...     park_lons, park_lats,
...     inner_radius_meters=1900,
...     outer_radius_meters=2000,
... )

apply_ring_color(ring_mask, ring_color=(0.1, 0.1, 0.1))

Apply a solid color to vertices within the ring mask.

This modifies self.colors in-place, setting the RGB values for vertices that fall within the ring mask to the specified color.

Parameters:

• ring_mask (ndarray) – 2D boolean array matching DEM shape. True = apply ring color.

• ring_color (tuple) – RGB tuple (0-1 range) for ring color. Default: dark gray.

Raises:

ValueError – If ring_mask shape doesn’t match DEM or colors not initialized.

Return type:

None

Loading Functions

src.terrain.core.load_dem_files(directory_path, pattern='*.hgt', recursive=False)

Load and merge DEM files from a directory into a single elevation dataset. Supports any raster format readable by rasterio (HGT, GeoTIFF, etc.).

Parameters:

• directory_path (str) – Path to directory containing DEM files

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

• recursive (bool) – Whether to search subdirectories recursively (default: False)

Returns:

(merged_dem, transform) where:

• merged_dem: numpy array containing the merged elevation data

• transform: affine transform mapping pixel to geographic coordinates

Return type:

tuple

Raises:

• ValueError – If no valid DEM files are found or directory doesn’t exist

• OSError – If directory access fails or file reading fails

• rasterio.errors.RasterioIOError – If there are issues reading the DEM files

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

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)
... )

Example:

dem, transform = load_filtered_hgt_files(
    'path/to/tiles',
    min_latitude=41  # Only load N41 and above
)

Color Mapping

See set_color_mapping() for single colormap usage.

See set_blended_color_mapping() for dual-colormap visualization, used in Combined Render: Full-Featured Example.

Data Layers

See add_data_layer() for adding georeferenced data layers, used in Snow Integration: Sledding Location Analysis for adding snow data.

See compute_proximity_mask() for creating proximity zones, used in Combined Render: Full-Featured Example for park proximity zones.