Spectral extraction pipeline¶

1. Introduction & User Guide¶

Why is this useful?¶

In remote sensing and forestry (especially with hyperspectral or high-resolution lidar data), raster files are often significantly larger than available system RAM. A standard approach of "load everything, then process" leads to memory error crashes or system freezing.

Phytospatial solves this by offering an adaptive extraction engine that abstracts away the complexity of memory management. It provides four strategies:

1. In-Memory: Fastest, but requires the whole file to fit in RAM.
2. Tiled: Standard streaming (512x512 windows). Safe for any file size but slightly slower due to overhead.
3. Blocked: Optimized streaming that reads data using the file's native internal structure (chunks/strips), minimizing disk seeking.
4. Auto: The "production" mode. It statically analyzes your hardware and the file structure to pick the best strategy automatically.

This tutorial demonstrates how to extract spectral information from rasters (either as statistics or raw pixel data) for individual tree crowns.

2. Setup and Imports¶

First, we import the necessary modules. We will use phytospatial's tiered architecture:

• Tier 1 (IO/Layer): For direct data handling.
• Tier 2 (Resources): For safety checks.
• Tier 3 (Extract): The high-level orchestration engine.

In [ ]:

import logging
import sys
from pathlib import Path
import pandas as pd

# Phytospatial imports
from phytospatial.raster.io import load
from phytospatial.raster.resources import estimate_memory_safety, analyze_structure, ProcessingMode
from phytospatial.extract import extract_to_dataframe

# Configure logging to see phytospatial's decision making process
logging.basicConfig(
    level=logging.INFO, 
    format='%(asctime)s - %(levelname)s - %(message)s',
    force=True
)
log = logging.getLogger(__name__)

import logging import sys from pathlib import Path import pandas as pd # Phytospatial imports from phytospatial.raster.io import load from phytospatial.raster.resources import estimate_memory_safety, analyze_structure, ProcessingMode from phytospatial.extract import extract_to_dataframe # Configure logging to see phytospatial's decision making process logging.basicConfig( level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s', force=True ) log = logging.getLogger(__name__)

3. Data Configuration¶

Define your input paths here.

• Crowns: A vector file (Shapefile, GeoJSON, GeoPackage) containing polygons.
• Raster: A geospatial raster (TIFF, ENVI, etc.).

In [ ]:

# Update these paths to point to your actual data
DATA_DIR = Path("data")
CROWNS_FILE = DATA_DIR / "crowns.shp"
RASTER_FILE = DATA_DIR / "input_hdrs/mosaic_test.hdr" # or .tif
OUTPUT_DIR = DATA_DIR / "results"

OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

print(f"Crowns: {CROWNS_FILE}")
print(f"Raster: {RASTER_FILE}")

# Update these paths to point to your actual data DATA_DIR = Path("data") CROWNS_FILE = DATA_DIR / "crowns.shp" RASTER_FILE = DATA_DIR / "input_hdrs/mosaic_test.hdr" # or .tif OUTPUT_DIR = DATA_DIR / "results" OUTPUT_DIR.mkdir(parents=True, exist_ok=True) print(f"Crowns: {CROWNS_FILE}") print(f"Raster: {RASTER_FILE}")

4. Perform checks (Resource Analysis)¶

Before running any heavy code, it is best practice to perform a "Pre-flight Check." This uses estimate_memory_safety and analyze_structure to peek at the file metadata without loading the pixel data.

Why do this?

• Prevents crashes before they happen.
• Tells you if your file is "efficient" (tiled internally) or "inefficient" (striped/scanline), which informs your processing strategy.

In [ ]:

if not CROWNS_FILE.exists() or not RASTER_FILE.exists():
    print("⚠️ Error: One or more input files do not exist. Please check paths.")
else:
    print(f"--- Inspecting {RASTER_FILE.name} --- make")

    # Check Memory Safety
    estimate = estimate_memory_safety(RASTER_FILE)
    
    print(f"Raw Data Size:    {estimate.raw_bytes / (1024**3):.2f} GB")
    print(f"Required RAM:     {estimate.total_required_bytes / (1024**3):.2f} GB (includes safety overhead)")
    print(f"Available RAM:    {estimate.available_system_bytes / (1024**3):.2f} GB")
    print(f"Is Safe to Load?  {'✅ YES' if estimate.is_safe else '❌ NO'}")
    
    # Check File Structure
    struct = analyze_structure(RASTER_FILE)
    print(f"\nInternal Structure: {'Tiled' if struct.is_tiled else 'Striped/Other'}")
    print(f"Block Shape:        {struct.block_shape}")
    print(f"Efficiency Score:   {struct.efficiency_score:.2f} (1.0 is perfect)")
    
    print(f"\nRecommended Mode:   {estimate.recommendation.value.upper()}")

if not CROWNS_FILE.exists() or not RASTER_FILE.exists(): print("⚠️ Error: One or more input files do not exist. Please check paths.") else: print(f"--- Inspecting {RASTER_FILE.name} --- make") # Check Memory Safety estimate = estimate_memory_safety(RASTER_FILE) print(f"Raw Data Size: {estimate.raw_bytes / (1024**3):.2f} GB") print(f"Required RAM: {estimate.total_required_bytes / (1024**3):.2f} GB (includes safety overhead)") print(f"Available RAM: {estimate.available_system_bytes / (1024**3):.2f} GB") print(f"Is Safe to Load? {'✅ YES' if estimate.is_safe else '❌ NO'}") # Check File Structure struct = analyze_structure(RASTER_FILE) print(f"\nInternal Structure: {'Tiled' if struct.is_tiled else 'Striped/Other'}") print(f"Block Shape: {struct.block_shape}") print(f"Efficiency Score: {struct.efficiency_score:.2f} (1.0 is perfect)") print(f"\nRecommended Mode: {estimate.recommendation.value.upper()}")

5. Pipeline A: The "Happy Path" (In-Memory)¶

If the check says "YES", In-Memory is the fastest method. It loads the entire raster object into Python's heap.

Use case: Small study areas, individual tiles, or machines with massive RAM (64GB+).

In [ ]:

%%time

# Check safety first to avoid crashing the kernel
estimate = estimate_memory_safety(RASTER_FILE)

if estimate.is_safe:
    print("Loading raster into memory...")
    # Explicitly load the raster object first
    raster = load(RASTER_FILE)
    
    print("Extracting features...")
    df_memory = extract_to_dataframe(
        raster_input=raster,
        vector_input=CROWNS_FILE,
        threshold=0.001
    )
    
    output_path = OUTPUT_DIR / "stats_memory.csv"
    df_memory.to_csv(output_path, index=False)
    print(f"✓ Saved {len(df_memory)} rows to {output_path}")
    
else:
    print("⚠️ Skipping In-Memory pipeline: Raster is too large for this machine.")

%%time # Check safety first to avoid crashing the kernel estimate = estimate_memory_safety(RASTER_FILE) if estimate.is_safe: print("Loading raster into memory...") # Explicitly load the raster object first raster = load(RASTER_FILE) print("Extracting features...") df_memory = extract_to_dataframe( raster_input=raster, vector_input=CROWNS_FILE, threshold=0.001 ) output_path = OUTPUT_DIR / "stats_memory.csv" df_memory.to_csv(output_path, index=False) print(f"✓ Saved {len(df_memory)} rows to {output_path}") else: print("⚠️ Skipping In-Memory pipeline: Raster is too large for this machine.")

6. Pipeline B: Tiled Streaming (The Safety Net)¶

If your raster is massive (for example 50GB), you cannot load it. Tiled mode creates a virtual grid over the file and processes it one window at a time.

Use case: Massive mosaics or files with inefficient internal structures (scanline/striped TIFFs).

In [ ]:

%%time

print("Starting Tiled extraction (Virtual 512x512 grid)...")

# Notice we pass the file PATH, not the loaded object
df_tiled = extract_to_dataframe(
    raster_input=RASTER_FILE,
    vector_input=CROWNS_FILE,
    tile_mode="tiled",
    tile_size=512,  # Adjustable window size
    threshold=0.001
)

output_path = OUTPUT_DIR / "stats_tiled.csv"
    df_tiled.to_csv(output_path, index=False)
print(f"✓ Saved {len(df_tiled)} rows to {output_path}")

%%time print("Starting Tiled extraction (Virtual 512x512 grid)...") # Notice we pass the file PATH, not the loaded object df_tiled = extract_to_dataframe( raster_input=RASTER_FILE, vector_input=CROWNS_FILE, tile_mode="tiled", tile_size=512, # Adjustable window size threshold=0.001 ) output_path = OUTPUT_DIR / "stats_tiled.csv" df_tiled.to_csv(output_path, index=False) print(f"✓ Saved {len(df_tiled)} rows to {output_path}")

7. Pipeline C: Blocked Streaming (Native IO)¶

If your file is a Cloud Optimized GeoTIFF (COG) or uses internal tiling, Blocked mode is superior to Tiled. It reads data exactly as it sits on the disk, creating zero read-overhead.

Use case: Optimized formats (COGs) or ENVI files with good interleaving.

In [ ]:

%%time

# Only run this if the structure analysis suggested it's efficient
struct = analyze_structure(RASTER_FILE)

if struct.is_efficient:
    print(f"Starting Blocked extraction (Native {struct.block_shape})...")

    df_blocked = extract_to_dataframe(
        raster_input=RASTER_FILE,
        vector_input=CROWNS_FILE,
        tile_mode="blocked",
        threshold=0.001
    )
    
    output_path = OUTPUT_DIR / "stats_blocked.csv"
    df_blocked.to_csv(output_path, index=False)
    print(f"✓ Saved {len(df_blocked)} rows to {output_path}")

else:
    print("⚠️ Skipping Blocked pipeline: File structure is inefficient (likely striped).")

%%time # Only run this if the structure analysis suggested it's efficient struct = analyze_structure(RASTER_FILE) if struct.is_efficient: print(f"Starting Blocked extraction (Native {struct.block_shape})...") df_blocked = extract_to_dataframe( raster_input=RASTER_FILE, vector_input=CROWNS_FILE, tile_mode="blocked", threshold=0.001 ) output_path = OUTPUT_DIR / "stats_blocked.csv" df_blocked.to_csv(output_path, index=False) print(f"✓ Saved {len(df_blocked)} rows to {output_path}") else: print("⚠️ Skipping Blocked pipeline: File structure is inefficient (likely striped).")

8. Pipeline D: Auto Mode (Production Ready)¶

In a production script, you don't want to manually check if estimate.is_safe every time. Use tile_mode="auto".

How it works:

1. It checks your available RAM.
2. If it fits, it uses In-Memory (Fastest).
3. If not, it checks the file structure.
4. If the structure is good, it uses Blocked.
5. If the structure is bad, it falls back to Tiled.

In [ ]:

%%time

print("Starting Auto extraction (Letting the engine decide)...")

df_auto = extract_to_dataframe(
    raster_input=RASTER_FILE,
    vector_input=CROWNS_FILE,
    tile_mode="auto",  # The magic switch
    threshold=0.001
)

output_path = OUTPUT_DIR / "stats_auto.csv"
df_auto.to_csv(output_path, index=False)
print(f"✓ Saved {len(df_auto)} rows to {output_path}")

%%time print("Starting Auto extraction (Letting the engine decide)...") df_auto = extract_to_dataframe( raster_input=RASTER_FILE, vector_input=CROWNS_FILE, tile_mode="auto", # The magic switch threshold=0.001 ) output_path = OUTPUT_DIR / "stats_auto.csv" df_auto.to_csv(output_path, index=False) print(f"✓ Saved {len(df_auto)} rows to {output_path}")

9. Inspecting Results¶

Finally, let's look at the data we extracted. The dataframe will contain the Crown ID, Species, and the calculated statistics for every band in the raster.

In [ ]:

# Display the first few rows of the result
print(f"Extraction Shape: {df_auto.shape}")
df_auto.head()

# Display the first few rows of the result print(f"Extraction Shape: {df_auto.shape}") df_auto.head()