tessl/pypi-srtm-py

Python parser for the Shuttle Radar Topography Mission elevation data

—

Pending

Quality

Pending

Does it follow best practices?

Impact

Pending

No eval scenarios have been run

Overview

Eval results

Files

Advanced Interpolation

Name: tessl/pypi-srtm-py
Author: tessl

Sophisticated algorithms for elevation estimation including Inverse Distance Weighted (IDW) interpolation and approximation methods for improved accuracy between data points.

Capabilities

Inverse Distance Weighted (IDW) Interpolation

Advanced interpolation method that estimates elevation values by weighting nearby data points based on their inverse distance, providing smoother and more accurate elevation estimates.

Note: These are private methods (indicated by the underscore prefix) intended for internal use, but they can be accessed directly if needed for advanced interpolation control.

class GeoElevationData:
    def _IDW(self, latitude: float, longitude: float, radius: float = 1) -> Optional[float]: ...

class GeoElevationFile:
    def _InverseDistanceWeighted(self, latitude: float, longitude: float, radius: float = 1) -> Optional[float]: ...

Parameters:

latitude (float): Target latitude coordinate
longitude (float): Target longitude coordinate
radius (float): Interpolation radius (1 or 2). Default: 1
- radius=1: 3x3 matrix of surrounding points
- radius=2: 5x5 matrix of surrounding points

Returns: Interpolated elevation value in meters, or None if insufficient data

Usage Example:

import srtm

elevation_data = srtm.get_data()

# Standard elevation lookup
standard_elevation = elevation_data.get_elevation(45.8566, 7.8566)

# IDW interpolation with 3x3 matrix (radius=1)
idw_elevation_r1 = elevation_data._IDW(45.8566, 7.8566, radius=1)

# IDW interpolation with 5x5 matrix (radius=2) 
idw_elevation_r2 = elevation_data._IDW(45.8566, 7.8566, radius=2)

print(f"Standard elevation: {standard_elevation}m")
print(f"IDW interpolation (r=1): {idw_elevation_r1}m")
print(f"IDW interpolation (r=2): {idw_elevation_r2}m")

# Calculate differences
if all(x is not None for x in [standard_elevation, idw_elevation_r1, idw_elevation_r2]):
    diff_r1 = abs(standard_elevation - idw_elevation_r1)
    diff_r2 = abs(standard_elevation - idw_elevation_r2)
    print(f"Difference (r=1): {diff_r1:.2f}m")
    print(f"Difference (r=2): {diff_r2:.2f}m")

Distance-Weighted Approximation

Simpler approximation method using distance-weighted averaging of nearby elevation points for smoother elevation profiles.

class GeoElevationFile:
    def approximation(self, latitude: float, longitude: float) -> Optional[float]: ...

Parameters:

latitude (float): Target latitude coordinate
longitude (float): Target longitude coordinate

Returns: Approximate elevation value in meters, or None if no data available

Usage Example:

import srtm

elevation_data = srtm.get_data()

# Get specific SRTM file
srtm_file = elevation_data.get_file(45.8566, 7.8566)

if srtm_file:
    # Direct elevation lookup
    direct_elevation = srtm_file.get_elevation(45.8566, 7.8566)
    
    # Distance-weighted approximation
    approx_elevation = srtm_file.approximation(45.8566, 7.8566)
    
    print(f"Direct elevation: {direct_elevation}m")
    print(f"Approximated elevation: {approx_elevation}m")
    
    if direct_elevation and approx_elevation:
        difference = abs(direct_elevation - approx_elevation)
        print(f"Approximation difference: {difference:.2f}m")

Comparative Interpolation Analysis

Compare different interpolation methods to understand their effects on elevation accuracy and smoothness.

Usage Example:

import srtm
import numpy as np

def compare_interpolation_methods(latitude, longitude):
    """Compare different interpolation methods for a given coordinate."""
    elevation_data = srtm.get_data()
    srtm_file = elevation_data.get_file(latitude, longitude)
    
    results = {}
    
    # Direct lookup (no interpolation)
    results['direct'] = elevation_data.get_elevation(latitude, longitude, approximate=False)
    
    # Built-in approximation
    results['approximate'] = elevation_data.get_elevation(latitude, longitude, approximate=True)
    
    # IDW with different radii
    results['idw_r1'] = elevation_data._IDW(latitude, longitude, radius=1)
    results['idw_r2'] = elevation_data._IDW(latitude, longitude, radius=2)
    
    # File-level approximation
    if srtm_file:
        results['file_approx'] = srtm_file.approximation(latitude, longitude)
        results['file_idw'] = srtm_file._InverseDistanceWeighted(latitude, longitude, radius=1)
    
    return results

# Test coordinates in mountainous area
test_coords = [
    (45.8566, 7.8566),  # Mont Blanc region
    (46.5197, 7.9553),  # Jungfrau region
    (45.9763, 7.6584)   # Matterhorn region
]

for lat, lon in test_coords:
    print(f"\nInterpolation comparison for ({lat}, {lon}):")
    results = compare_interpolation_methods(lat, lon)
    
    for method, elevation in results.items():
        if elevation is not None:
            print(f"  {method:12}: {elevation:7.2f}m")
        else:
            print(f"  {method:12}: No data")
    
    # Calculate variations
    valid_results = {k: v for k, v in results.items() if v is not None}
    if len(valid_results) > 1:
        elevations = list(valid_results.values())
        std_dev = np.std(elevations)
        range_val = max(elevations) - min(elevations)
        print(f"  Standard deviation: {std_dev:.2f}m")
        print(f"  Range: {range_val:.2f}m")

Batch Interpolation Processing

Apply interpolation methods to multiple coordinates efficiently for large-scale elevation analysis.

Usage Example:

import srtm
import csv
from typing import List, Tuple

def batch_interpolation(coordinates: List[Tuple[float, float]], method='idw_r1'):
    """Apply interpolation to a batch of coordinates."""
    elevation_data = srtm.get_data(batch_mode=True)
    
    results = []
    
    for i, (lat, lon) in enumerate(coordinates):
        if i % 100 == 0:
            print(f"Processing coordinate {i+1}/{len(coordinates)}")
        
        if method == 'direct':
            elevation = elevation_data.get_elevation(lat, lon, approximate=False)
        elif method == 'approximate':
            elevation = elevation_data.get_elevation(lat, lon, approximate=True)
        elif method == 'idw_r1':
            elevation = elevation_data._IDW(lat, lon, radius=1)
        elif method == 'idw_r2':
            elevation = elevation_data._IDW(lat, lon, radius=2)
        else:
            elevation = elevation_data.get_elevation(lat, lon)
        
        results.append({
            'latitude': lat,
            'longitude': lon,
            'elevation': elevation,
            'method': method
        })
    
    return results

# Generate test coordinates along a hiking trail
trail_coordinates = []
start_lat, start_lon = 45.8, 7.8
for i in range(50):
    lat = start_lat + (i * 0.001)  # Move north
    lon = start_lon + (i * 0.0005)  # Move slightly east
    trail_coordinates.append((lat, lon))

# Process with different methods
methods = ['direct', 'approximate', 'idw_r1', 'idw_r2']

for method in methods:
    print(f"\nProcessing with method: {method}")
    results = batch_interpolation(trail_coordinates, method)
    
    # Save results to CSV
    filename = f"elevation_trail_{method}.csv"
    with open(filename, 'w', newline='') as csvfile:
        fieldnames = ['latitude', 'longitude', 'elevation', 'method']
        writer = csv.DictWriter(csvfile, fieldnames=fieldnames)
        writer.writeheader()
        writer.writerows(results)
    
    # Calculate statistics
    valid_elevations = [r['elevation'] for r in results if r['elevation'] is not None]
    if valid_elevations:
        min_elev = min(valid_elevations)
        max_elev = max(valid_elevations)
        avg_elev = sum(valid_elevations) / len(valid_elevations)
        print(f"  Elevation range: {min_elev:.1f}m to {max_elev:.1f}m")
        print(f"  Average elevation: {avg_elev:.1f}m")
        print(f"  Valid points: {len(valid_elevations)}/{len(results)}")
    
    print(f"  Results saved to {filename}")

Interpolation Quality Assessment

Evaluate the quality and consistency of different interpolation methods.

Usage Example:

import srtm
import random
import statistics

def assess_interpolation_quality(region_bounds, sample_size=100):
    """Assess quality of interpolation methods over a region."""
    lat_min, lat_max, lon_min, lon_max = region_bounds
    elevation_data = srtm.get_data()
    
    # Generate random sample points
    sample_points = []
    for _ in range(sample_size):
        lat = random.uniform(lat_min, lat_max)
        lon = random.uniform(lon_min, lon_max)
        sample_points.append((lat, lon))
    
    # Test each interpolation method
    methods = {
        'direct': lambda lat, lon: elevation_data.get_elevation(lat, lon, approximate=False),
        'approximate': lambda lat, lon: elevation_data.get_elevation(lat, lon, approximate=True),
        'idw_r1': lambda lat, lon: elevation_data._IDW(lat, lon, radius=1),
        'idw_r2': lambda lat, lon: elevation_data._IDW(lat, lon, radius=2)
    }
    
    results = {method: [] for method in methods}
    
    for lat, lon in sample_points:
        for method_name, method_func in methods.items():
            elevation = method_func(lat, lon)
            if elevation is not None:
                results[method_name].append(elevation)
    
    # Analyze results
    print(f"Interpolation quality assessment over region:")
    print(f"  Bounds: {lat_min}-{lat_max}°N, {lon_min}-{lon_max}°E")
    print(f"  Sample size: {sample_size} points")
    print()
    
    for method, elevations in results.items():
        if elevations:
            mean_elev = statistics.mean(elevations)
            std_elev = statistics.stdev(elevations) if len(elevations) > 1 else 0
            min_elev = min(elevations)
            max_elev = max(elevations)
            coverage = len(elevations) / sample_size * 100
            
            print(f"{method:12}:")
            print(f"  Coverage: {coverage:5.1f}%")
            print(f"  Mean:     {mean_elev:7.1f}m")
            print(f"  Std Dev:  {std_elev:7.1f}m")
            print(f"  Range:    {min_elev:7.1f}m - {max_elev:7.1f}m")
            print()
    
    return results

# Assess quality over Alpine region
alpine_bounds = (45.5, 46.5, 7.0, 8.0)  # (lat_min, lat_max, lon_min, lon_max)
quality_results = assess_interpolation_quality(alpine_bounds, sample_size=200)

Algorithm Details

IDW Algorithm

The Inverse Distance Weighted interpolation uses the formula:

Z(x,y) = Σ(wi * zi) / Σ(wi)

Where:

wi = 1/di^p (weight based on inverse distance)
di is the distance from point (x,y) to known point i
zi is the elevation value at known point i
p is the power parameter (typically 2)

Performance Considerations

Computational Cost: IDW with larger radius requires more calculations
Accuracy Trade-off: Higher radius generally provides smoother results but may lose local detail
Memory Usage: Larger interpolation matrices require more memory
Data Availability: Interpolation quality depends on density of available SRTM data points

When to Use Each Method

Direct Lookup: Fastest, use when exact SRTM grid values are sufficient
Built-in Approximation: Good balance of speed and accuracy for most applications
IDW (radius=1): Better accuracy with moderate computational cost
IDW (radius=2): Smoothest results, highest computational cost
File Approximation: Use when working directly with GeoElevationFile objects