Snow Analysis Module

SNODAS snow data processing and slope statistics computation.

This package provides tools for loading, processing, and analyzing snow data from NOAA’s SNODAS (Snow Data Assimilation System) and computing terrain slope statistics for scoring.

SNODAS Processing

Pipeline for loading and computing snow statistics:

  1. batch_process_snodas_data - Load and reproject SNODAS files

  2. calculate_snow_statistics - Compute aggregated seasonal statistics

  3. load_snodas_stats - High-level orchestrator with caching and fallback

SNODAS data processing utilities.

Functions for loading, processing, and computing statistics from NOAA SNODAS (Snow Data Assimilation System) snow depth data.

Pipeline steps: - batch_process_snodas_data: Load and reproject SNODAS files - calculate_snow_statistics: Compute aggregated snow statistics

High-level orchestration: - load_snodas_stats: Complete pipeline with tiling, caching, and mock fallback

src.snow.snodas.batch_process_snodas_data(snodas_dir, extent, target_shape, processed_dir='processed_snodas', max_workers=14)[source]

Step 1: Load and process SNODAS files.

Process all SNODAS .dat.gz files in snodas_dir, cropping to ‘extent’ & reprojecting to ~1/120 deg (EPSG:4326).

Parameters:

  • snodas_dir – Directory containing SNODAS .dat.gz files

  • extent (Tuple[float, float, float, float]) – (minx, miny, maxx, maxy) bounding box

  • target_shape (Tuple[int, int]) – (height, width) for output arrays

  • processed_dir (str) – Directory for processed .npz files

  • max_workers (int) – Number of parallel workers

Returns:

Dict with processed_files mapping dates to file paths

Return type:

Dict[str, Any]

src.snow.snodas.calculate_snow_statistics(input_data, snow_season_start_month=11, snow_season_end_month=4)[source]

Step 2: Compute snow statistics from loaded files.

Group processed SNODAS by winter season and compute aggregated stats.

Parameters:

  • input_data (Dict[str, Any]) – Dict with “processed_files” from batch_process_snodas_data

  • snow_season_start_month (int) – First month of snow season (default: November)

  • snow_season_end_month (int) – Last month of snow season (default: April)

Returns:

  • median_max_depth: Median of seasonal max depths

  • mean_snow_day_ratio: Average fraction of days with snow

  • interseason_cv: Year-to-year variability

  • mean_intraseason_cv: Within-winter variability

  • metadata: Processing metadata

  • failed_files: List of files that failed to process

Return type:

Dict with final statistics

src.snow.snodas.load_snodas_stats(terrain=None, snodas_dir=None, *, cache_dir=PosixPath('snodas_cache'), cache_name='snodas', tile_config=None, mock_data=False, mock_shape=(500, 500))[source]

Load SNODAS snow statistics with automatic tiling, caching, and mock fallback.

Orchestrates the standard SNODAS pipeline: 1. batch_process_snodas_data (load and reproject) 2. calculate_snow_statistics (aggregate seasonal stats)

Uses GriddedDataLoader for memory-safe tiled processing. Falls back to mock data when real data is unavailable or loading fails.

Parameters:

  • terrain – Terrain object providing extent and resolution context. Required for real data; can be None with mock_data=True.

  • snodas_dir (Path | None) – Directory containing SNODAS .dat.gz files.

  • cache_dir (Path) – Directory for pipeline caching.

  • cache_name (str) – Base name for cache files.

  • tile_config – Optional TiledDataConfig for memory-safe processing. If None, uses GriddedDataLoader defaults.

  • mock_data (bool) – If True, return mock data without attempting real loading.

  • mock_shape (tuple) – Shape for mock data arrays.

Returns:

  • median_max_depth: Median of seasonal max depths (mm)

  • mean_snow_day_ratio: Average fraction of days with snow

  • interseason_cv: Year-to-year variability

  • mean_intraseason_cv: Within-winter variability

Return type:

Dict with snow statistics arrays

Slope Statistics

Tiled slope computation at full DEM resolution with geographic transform-aware aggregation. Ensures cliff faces within a pixel are captured rather than hidden by downsampling.

Tiled slope statistics computation using geographic transforms.

Computes slope and related statistics at full DEM resolution, then aggregates per output pixel using exact geographic coordinate mapping.

This ensures that: 1. Cliff faces within a pixel are captured (not hidden by downsampling) 2. Downsampling respects actual geographic transforms (not uniform strides) 3. Memory usage stays bounded by processing in tiles with halos

class src.snow.slope_statistics.TiledSlopeConfig(target_tile_outputs=100, halo=1, max_tile_size=4096)[source]

Bases: object

Configuration for tiled slope computation.

Parameters:

  • target_tile_outputs (int)

  • halo (int)

  • max_tile_size (int)

target_tile_outputs: int = 100

Target number of output pixels per tile dimension.

halo: int = 1

Halo size (pixels) for gradient computation at tile boundaries.

max_tile_size: int = 4096

Maximum tile size in source pixels for memory safety.

__init__(target_tile_outputs=100, halo=1, max_tile_size=4096)
Parameters:

  • target_tile_outputs (int)

  • halo (int)

  • max_tile_size (int)

Return type:

None

class src.snow.slope_statistics.TileSpec(src_slice, core_slice, out_slice, row_stride, col_stride)[source]

Bases: object

Specification for a single tile during processing.

Parameters:
src_slice: Tuple[slice, slice]

Slice into source DEM (with halo).

core_slice: Tuple[slice, slice]

Slice within tile to extract core region (no halo).

out_slice: Tuple[slice, slice]

Slice into output arrays.

row_stride: int

Average source rows per output row in this tile.

col_stride: int

Average source columns per output column in this tile.

__init__(src_slice, core_slice, out_slice, row_stride, col_stride)
Parameters:
Return type:

None

class src.snow.slope_statistics.SlopeStatistics(slope_mean, slope_max, slope_min, slope_std, slope_p95, roughness, aspect_sin, aspect_cos)[source]

Bases: object

Aggregated slope statistics at output resolution.

Parameters:
slope_mean: ndarray

Mean slope in degrees for each output pixel.

slope_max: ndarray

Maximum slope (cliff detection).

slope_min: ndarray

Minimum slope (flat spots).

__init__(slope_mean, slope_max, slope_min, slope_std, slope_p95, roughness, aspect_sin, aspect_cos)
Parameters:
Return type:

None

slope_std: ndarray

Standard deviation of slopes (terrain consistency).

slope_p95: ndarray

95th percentile slope (robust cliff indicator).

roughness: ndarray

Elevation standard deviation (surface bumpiness).

aspect_sin: ndarray

sin(aspect) for vector averaging across output pixels.

aspect_cos: ndarray

cos(aspect) for vector averaging across output pixels.

property dominant_aspect: ndarray

Compute dominant aspect from vector-averaged sin/cos components.

Handles circular nature of aspect (0-360 degrees).

Returns:

Dominant aspect in degrees (0=North, 90=East, 180=South, 270=West)

property aspect_strength: ndarray

Compute aspect strength (consistency of slope direction).

Returns:

Strength from 0 (uniform directions) to 1 (all same direction)

src.snow.slope_statistics.compute_pixel_mapping(dem_shape, dem_transform, dem_crs, target_shape, target_transform, target_crs)[source]

Compute precise pixel-to-pixel mapping using geographic transforms.

For each output pixel (i, j): 1. Calculate its geographic bounds using target_transform 2. Transform to source CRS if needed 3. Use dem_transform inverse to find source pixel bounds 4. Return mapping: output_pixel_ij → [list of source pixels]

Parameters:

  • dem_shape (Tuple[int, int]) – (height, width) of source DEM

  • dem_transform (Affine) – Affine transform for source DEM

  • dem_crs (str) – CRS string for source DEM (e.g., ‘EPSG:4326’)

  • target_shape (Tuple[int, int]) – (height, width) of target output grid

  • target_transform (Affine) – Affine transform for target grid

  • target_crs (str) – CRS string for target grid

Returns:

  • ‘row_stride’: average source pixels per output row

  • ’col_stride’: average source pixels per output column

  • ’mapping’: detailed per-pixel mapping (for irregular transforms)

Return type:

Mapping structure containing

Handles:

  • Non-integer scaling ratios

  • Rotation, skew, reprojection

  • Edge effects and partial pixel overlaps

src.snow.slope_statistics.compute_tile_layout(dem_shape, target_shape, pixel_mapping, config=None)[source]

Compute adaptive tile size based on pixel mapping.

Strategy: 1. Estimate average source pixels per output pixel 2. Tile size = target_tile_outputs × stride (in source pixels) 3. Clamp to max_tile_size for memory safety 4. Ensure tile boundaries align with output pixel boundaries

Parameters:

  • dem_shape (Tuple[int, int]) – (height, width) of source DEM

  • target_shape (Tuple[int, int]) – (height, width) of target grid

  • pixel_mapping (Dict) – Result from compute_pixel_mapping()

  • config (TiledSlopeConfig | None) – TiledSlopeConfig (uses defaults if None)

Returns:

List of TileSpec for each tile to process

Return type:

List[TileSpec]

src.snow.slope_statistics.compute_tile_slopes(tile_data, core_slice, cell_size_m)[source]

Compute slope and aspect at full resolution for a tile.

Parameters:

  • tile_data (ndarray) – Tile elevation data (with halo for proper boundaries)

  • core_slice (Tuple[slice, slice]) – Slice to extract core region (without halo)

  • cell_size_m (float) – Cell size in meters

Returns:

  • slope_deg: Slope in degrees (core region only)

  • aspect_deg: Aspect in degrees (core region only)

Return type:

Tuple of

src.snow.slope_statistics.aggregate_by_geographic_mapping(slope_full_res, aspect_full_res, elevation_full_res, row_stride, col_stride, output_shape)[source]

Aggregate statistics using geographic pixel mapping.

For now uses uniform stride-based aggregation. This could be extended to handle irregular mappings from reprojection.

Parameters:

  • slope_full_res (ndarray) – Slope in degrees at full resolution

  • aspect_full_res (ndarray) – Aspect in degrees at full resolution

  • elevation_full_res (ndarray) – Elevation at full resolution

  • row_stride (int) – Source pixels per output row

  • col_stride (int) – Source pixels per output column

  • output_shape (Tuple[int, int]) – (height, width) of output grid

Returns:

  • ‘mean’: mean slope per output pixel

  • ’max’: max slope per output pixel

  • ’min’: min slope per output pixel

  • ’std’: std dev of slopes

  • ’p95’: 95th percentile slope

  • ’roughness’: elevation std dev

  • ’aspect_sin’: vector-averaged aspect sin component

  • ’aspect_cos’: vector-averaged aspect cos component

Return type:

Dictionary of statistics arrays

src.snow.slope_statistics.compute_tiled_slope_statistics(dem, dem_transform, dem_crs, target_shape, target_transform, target_crs, config=None)[source]

Compute slope statistics at target resolution using tiled processing.

Processes the full-resolution DEM in tiles to preserve cliff faces and extreme slopes that would be lost by downsampling first.

Parameters:

  • dem (ndarray) – Full-resolution DEM array (height x width, float32 preferred)

  • dem_transform (Affine) – Affine transform for DEM

  • dem_crs (str) – CRS for DEM (e.g., ‘EPSG:4326’)

  • target_shape (Tuple[int, int]) – Output shape (out_height, out_width)

  • target_transform (Affine) – Affine transform for output grid

  • target_crs (str) – CRS for output grid

  • config (TiledSlopeConfig | None) – Configuration parameters (uses defaults if None)

Returns:

SlopeStatistics dataclass with all computed statistics

Return type:

SlopeStatistics

See Also