Computational Methods - Black-Lights/planetscope-py GitHub Wiki

Computational Methods

PlanetScope-py implements three enhanced computational methods for spatial density analysis, each optimized for different scenarios with coordinate system fixes and performance improvements built-in.

Enhanced Features

  • Fixed coordinate system handling - Proper north-to-south orientation with corrected transforms
  • Rasterization as optimized default - Changed from vector overlay to rasterization for better performance
  • Built-in coordinate fixes - Automatic correction of mirrored/flipped outputs
  • Enhanced chunk processing - Improved mosaicking and memory management
  • Robust error handling - PROJ/CRS compatibility with multiple fallbacks
  • Increased performance thresholds - Optimized method selection logic

Method Comparison Overview

Method Speed Memory Precision Coordinate Fixes Best For
Rasterization Fastest Low High ✓ Built-in Most use cases, default choice
Vector Overlay Moderate High Highest ✓ Built-in Complex geometries, precision-critical
Adaptive Grid Variable Lowest Good ✓ Built-in Large areas, memory constraints

1. Enhanced Rasterization Method

Enhanced Features

The enhanced rasterization method includes critical coordinate system fixes:

  • Corrected transform orientation - Proper north-to-south pixel ordering
  • Fixed northwest origin - Geographic alignment with (0,0) at northwest corner
  • Automatic coordinate validation - Built-in checks for geographic correctness
  • Enhanced performance - Optimized as new default method

How It Works

# Enhanced rasterization with coordinate system fixes
def _calculate_enhanced_rasterization_density():
    # CRITICAL COORDINATE SYSTEM FIX
    if config.coordinate_system_fixes:
        # Create corrected transform with proper orientation
        pixel_width = (bounds[2] - bounds[0]) / width    # Positive (west to east)
        pixel_height = -(bounds[3] - bounds[1]) / height # Negative (north to south)
        
        # Transform: start at northwest corner (proper geographic alignment)
        transform = Affine(
            pixel_width, 0.0, bounds[0],      # X: west to east
            0.0, pixel_height, bounds[3],     # Y: north to south (negative height)
            0.0, 0.0, 1.0
        )

Enhanced Configuration

from planetscope_py import SpatialDensityEngine, DensityConfig, DensityMethod

# Enhanced rasterization (now default)
config = DensityConfig(
    method=DensityMethod.RASTERIZATION,       # Enhanced default method
    resolution=30.0,                          # Optimized default resolution
    coordinate_system_fixes=True,             # Enable coordinate fixes (default)
    chunk_size_km=200.0,                      # Increased chunk size
    max_memory_gb=16.0,                       # Increased memory limit
    parallel_workers=4                        # Parallel processing
)

engine = SpatialDensityEngine(config)
result = engine.calculate_density(scenes, roi)

print(f"Method used: {result.method_used.value}")
print(f"Coordinate fixes applied: {result.coordinate_system_corrected}")
print(f"Computation time: {result.computation_time:.3f}s")

Enhanced Performance Characteristics

# Enhanced performance with coordinate fixes (Milan dataset: 43 scenes, 355 km²)
Resolution    Grid Size      Time      Coordinate Fixes    Performance
100m         266×133        0.03s     ✓ Applied          500,000+ cells/sec
30m          888×444        0.06s     ✓ Applied          2,600,000+ cells/sec  
10m          2,665×1,333    0.15s     ✓ Applied          23,700,000+ cells/sec
3m           8,870×4,435    0.45s     ✓ Applied          87,000,000+ cells/sec

When to Use Enhanced Rasterization

  • Default choice for all applications (enhanced default)
  • High-resolution analysis with coordinate accuracy
  • Regular and irregular geometries
  • Performance-critical applications with coordinate fixes
  • Large number of scenes (>100 scenes) with proper orientation

2. Enhanced Vector Overlay Method

Enhanced Features

  • Coordinate system validation - Built-in geographic correctness checks
  • Improved memory management - Better handling of complex geometries
  • Enhanced spatial indexing - Optimized intersection calculations
  • Robust error handling - Graceful degradation for complex cases

Enhanced Configuration

# Enhanced vector overlay with coordinate fixes
config = DensityConfig(
    method=DensityMethod.VECTOR_OVERLAY,
    resolution=100.0,                         # Moderate resolution for performance
    coordinate_system_fixes=True,             # Coordinate validation enabled
    parallel_workers=2,                       # Reduced for memory control
    max_memory_gb=16.0,                       # Increased memory allowance
    chunk_size_km=25.0                        # Smaller chunks for precision
)

engine = SpatialDensityEngine(config)
result = engine.calculate_density(scenes, roi)

print(f"Precision: Exact geometric calculations with coordinate fixes")
print(f"Memory usage optimized: {result.grid_info}")

Enhanced Performance

# Enhanced performance with coordinate fixes (Milan dataset: 43 scenes, 355 km²)
Resolution    Grid Size      Time      Coordinate Fixes    Intersections/sec
100m         267×134        53s       ✓ Applied           29,000
50m          533×267        203s      ✓ Applied           35,000
30m          889×445        >300s     ✓ Applied           33,000

When to Use Enhanced Vector Overlay

  • Precision-critical analysis with coordinate accuracy
  • Complex irregular geometries requiring exact calculations
  • Small to medium datasets (<500 scenes) with enhanced accuracy
  • Research applications requiring geographic precision
  • Quality validation with coordinate system verification

3. Enhanced Adaptive Grid Method

Enhanced Features

  • Coordinate-aware refinement - Hierarchical processing with proper orientation
  • Enhanced memory efficiency - Optimized for large areas with coordinate fixes
  • Smart refinement criteria - Geographic pattern recognition
  • Robust chunking - Proper coordinate handling across chunks

Enhanced Configuration

# Enhanced adaptive grid with coordinate fixes
config = DensityConfig(
    method=DensityMethod.ADAPTIVE_GRID,
    resolution=100.0,                         # Target resolution
    coordinate_system_fixes=True,             # Geographic correctness
    max_memory_gb=8.0,                        # Memory-efficient operation
    parallel_workers=4,                       # Balanced processing
    chunk_size_km=50.0                        # Optimized chunk size
)

engine = SpatialDensityEngine(config)
result = engine.calculate_density(scenes, large_roi)

print(f"Adaptive refinement with coordinate fixes applied")
print(f"Memory efficiency: {result.grid_info}")

Enhanced Hierarchical Processing

# Enhanced adaptive grid output with coordinate information
{
    'method_used': 'adaptive_grid',
    'coordinate_system_corrected': True,
    'grid_info': {
        'width': 133,
        'height': 66,
        'output_resolution': 200.0,
        'max_level': 2,
        'total_cells': 8978,
        'coordinate_fixes_applied': True,
        'level_statistics': {
            'level_0': {'cell_count': 2278, 'resolution_m': 400.0},
            'level_1': {'cell_count': 2278, 'resolution_m': 200.0},
            'level_2': {'cell_count': 9112, 'resolution_m': 100.0}
        }
    }
}

Enhanced Automatic Method Selection

Updated Selection Logic

The enhanced AUTO method selection favors rasterization and includes coordinate fixes:

def enhanced_select_method(scenes, roi, resolution, config):
    """Enhanced method selection with coordinate fixes and performance optimization."""
    
    # Calculate characteristics
    roi_area_km2 = calculate_area(roi)
    scene_count = len(scenes)
    estimated_raster_mb = estimate_memory(roi, resolution)
    
    # Enhanced selection logic (updated thresholds)
    if estimated_raster_mb > config.max_memory_gb * 1024 * 0.7:
        method = "adaptive_grid"    # Large memory requirement
        logger.info("Very large raster detected, using adaptive grid method")
    elif scene_count > 2000:       # Increased threshold from 1000
        method = "rasterization"    # Many scenes (enhanced default)
        logger.info("Many scenes detected, using rasterization method")
    else:
        method = "rasterization"    # Enhanced: Default to rasterization
        logger.info("Standard dataset, using rasterization method (enhanced default)")
    
    return method

Enhanced Selection Examples

# Enhanced selection with coordinate fixes enabled
config = DensityConfig(coordinate_system_fixes=True)

# Small ROI, few scenes → Enhanced Rasterization (changed from Vector Overlay)
roi = box(9.1, 45.4, 9.2, 45.5)  # 100 km²
scenes = 25
# Selected: rasterization (enhanced default)

# Medium ROI, many scenes → Enhanced Rasterization  
roi = box(9.0, 45.0, 9.5, 45.8)  # 2000 km²
scenes = 150
# Selected: rasterization (with coordinate fixes)

# Large ROI, any scenes → Enhanced Adaptive Grid
roi = box(8.0, 44.0, 12.0, 47.0)  # 50,000 km²
scenes = 500
# Selected: adaptive_grid (with coordinate fixes)

Enhanced Performance Benchmarks

Speed Comparison with Coordinate Fixes

# Enhanced benchmark results (Milan dataset: 43 scenes, 355 km²)
Method              100m Grid    50m Grid     30m Grid     Coordinate Fixes
Enhanced Raster    0.03s        0.06s        0.09s        ✓ Applied
Enhanced Vector    53s          203s         >300s        ✓ Applied
Enhanced Adaptive  9s           15s          25s          ✓ Applied

Accuracy Validation

All enhanced methods produce geographically correct results:

# Enhanced density statistics with coordinate validation
Method                Min   Max   Mean   Std    Coordinate Accuracy
Enhanced Raster      12    21    16.2   1.8    ✓ Validated
Enhanced Vector      13    21    16.3   1.7    ✓ Validated
Enhanced Adaptive    12    21    16.2   1.8    ✓ Validated

Enhanced Memory Usage

# Enhanced memory usage with coordinate fixes (Milan ROI)
Resolution    Enhanced Raster    Enhanced Vector    Enhanced Adaptive
100m         0.1 MB             2.5 MB             0.8 MB
30m          14 MB              45 MB              3.5 MB
10m          150 MB             500 MB             25 MB
3m           1.5 GB             5.5 GB             180 MB

Enhanced Method Selection Guide

Updated Decision Matrix

1. Need coordinate system accuracy?
   ALWAYS → All methods include coordinate fixes ✓

2. Is your ROI > 1000 km²?
   YES → Use Enhanced Adaptive Grid
   NO → Continue to 3

3. Do you have > 2000 scenes? (increased threshold)
   YES → Use Enhanced Rasterization
   NO → Continue to 4

4. Need maximum precision with coordinate accuracy?
   YES → Use Enhanced Vector Overlay
   NO → Continue to 5

5. Default recommendation?
   → Use Enhanced Rasterization (new default)

Enhanced Quick Reference

Scenario Enhanced Method Configuration Coordinate Fixes
Quick analysis Enhanced Rasterization resolution=100, method=RASTERIZATION ✓ Automatic
High precision Enhanced Vector Overlay resolution=30, method=VECTOR_OVERLAY ✓ Automatic
Large areas Enhanced Adaptive Grid resolution=100, method=ADAPTIVE_GRID ✓ Automatic
High resolution Enhanced Rasterization resolution=3, method=RASTERIZATION ✓ Automatic
Default workflow Enhanced Auto resolution=30, method=AUTO ✓ Automatic

Enhanced Configuration Reference

DensityConfig Enhanced Parameters

config = DensityConfig(
    resolution=30.0,                          # Grid resolution (default: 30.0, optimized)
    method=DensityMethod.RASTERIZATION,       # Enhanced default method
    chunk_size_km=200.0,                      # Increased chunk size (default: 200.0)
    max_memory_gb=16.0,                       # Increased memory limit (default: 16.0)
    parallel_workers=4,                       # Parallel processing (default: 4)
    no_data_value=-9999.0,                    # NoData value (default: -9999.0)
    coordinate_system_fixes=True,             # Enable coordinate fixes (default: True)
    force_single_chunk=False,                 # Force single chunk (default: False)
    validate_geometries=True                  # Validate inputs (default: True)
)

Enhanced Method Options

Method Value Enhanced Features Default Use
DensityMethod.RASTERIZATION "rasterization" Coordinate fixes, optimized performance ✓ New default
DensityMethod.VECTOR_OVERLAY "vector_overlay" Enhanced precision, coordinate validation Precision-critical
DensityMethod.ADAPTIVE_GRID "adaptive_grid" Memory-efficient, coordinate-aware refinement Large areas
DensityMethod.AUTO "auto" Enhanced selection logic, coordinate fixes Smart default

Enhanced Integration Examples

One-Line Analysis with Coordinate Fixes

from planetscope_py import quick_planet_analysis

# Enhanced analysis with automatic coordinate fixes
result = quick_planet_analysis(
    milan_roi, "last_month",
    resolution=30.0,                          # Optimized default
    method="rasterization"                    # Enhanced default method
)

# Validate coordinate fixes applied
if result['summary'].get('coordinate_system_corrected', False):
    print("✓ Enhanced coordinate system fixes applied automatically")
    print(f"✓ Method: {result['density_result'].method_used.value}")
    print(f"✓ Performance: {result['summary']['computation_time_s']:.2f}s")

Enhanced Method Comparison

# Compare enhanced methods with coordinate fixes
methods = ['rasterization', 'vector_overlay', 'adaptive_grid']
comparison_results = {}

for method in methods:
    print(f"Testing enhanced {method} method...")
    
    config = DensityConfig(
        method=getattr(DensityMethod, method.upper()),
        resolution=100.0,
        coordinate_system_fixes=True          # Enable for all methods
    )
    
    engine = SpatialDensityEngine(config)
    result = engine.calculate_density(scenes, roi)
    
    comparison_results[method] = {
        'computation_time': result.computation_time,
        'coordinate_fixes': result.coordinate_system_corrected,
        'method_used': result.method_used.value,
        'geographic_accuracy': 'validated' if result.coordinate_system_corrected else 'unvalidated'
    }

# Display enhanced comparison
print("\nEnhanced Method Comparison:")
print("Method           | Time    | Coord Fixes | Geographic Accuracy")
print("-" * 60)
for method, stats in comparison_results.items():
    print(f"{method:15} | {stats['computation_time']:6.2f}s | {stats['coordinate_fixes']!s:10} | {stats['geographic_accuracy']}")

Enhanced Troubleshooting

Coordinate System Issues (Now Fixed)

# Problem: Mirrored or flipped output maps
# Solution: Enhanced coordinate fixes (automatic in v4.0+)

config = DensityConfig(
    coordinate_system_fixes=True              # Enabled by default
)

result = engine.calculate_density(scenes, roi)

# Verify fixes applied
if result.coordinate_system_corrected:
    print("✓ Coordinate system issues automatically resolved")
    print("✓ Geographic alignment validated")
    print("✓ North-to-south orientation corrected")

Enhanced Performance Issues

# Enhanced performance optimization
def optimize_method_selection(roi, scenes, target_time_seconds=60):
    """Enhanced method selection for performance targets."""
    
    roi_area_km2 = calculate_area_km2(roi)
    scene_count = len(scenes)
    
    if target_time_seconds < 10:               # Very fast
        return DensityConfig(
            method=DensityMethod.RASTERIZATION,
            resolution=100.0,
            coordinate_system_fixes=True       # Always enabled
        )
    elif target_time_seconds < 60:             # Fast
        return DensityConfig(
            method=DensityMethod.RASTERIZATION,
            resolution=30.0,
            coordinate_system_fixes=True
        )
    else:                                      # Precision over speed
        return DensityConfig(
            method=DensityMethod.VECTOR_OVERLAY,
            resolution=50.0,
            coordinate_system_fixes=True
        )

# Usage
optimal_config = optimize_method_selection(roi, scenes, target_time_seconds=30)

Next Steps