Python Script to Split Large Datasets for Processing

#!/usr/bin/env python3
"""
Geospatial Dataset Partitioner for UAV Photogrammetry
Splits large imagery collections into spatially contiguous, buffered chunks
with explicit validation thresholds and memory-aware routing.
"""

import argparse
import json
import logging
import math
import os
from pathlib import Path
from typing import Optional, Dict, List, Tuple

import exifread
from shapely.geometry import Point, box
from shapely.affinity import translate

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")

def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(description="Geospatial partitioner for UAV imagery datasets.")
    parser.add_argument("--input-dir", type=Path, required=True, help="Root directory containing raw imagery.")
    parser.add_argument("--output-manifest", type=Path, required=True, help="JSON manifest path for chunk definitions.")
    parser.add_argument("--tile-size-m", type=float, default=1000.0, help="Target tile dimension in meters (default: 1000).")
    parser.add_argument("--overlap-pct", type=float, default=0.20, help="Overlap buffer as fraction of tile size (0.15-0.30).")
    parser.add_argument("--min-images", type=int, default=50, help="Minimum images per chunk before fallback routing.")
    parser.add_argument("--max-memory-mb", type=int, default=4096, help="Hard memory cap per chunk (MB).")
    parser.add_argument("--crs-tolerance-deg", type=float, default=0.0001, help="Max allowed CRS drift in decimal degrees.")
    return parser.parse_args()

def extract_gps_decimal(path: Path) -> Optional[Tuple[float, float]]:
    """Extract decimal degrees from EXIF. Returns None if tags are missing or malformed."""
    try:
        with open(path, "rb") as f:
            tags = exifread.process_file(f, details=False)
        
        lat_key, lon_key = "GPS GPSLatitude", "GPS GPSLongitude"
        if lat_key not in tags or lon_key not in tags:
            return None

        def dms_to_dec(dms_tag, ref_tag_name):
            vals = [float(v.num) / float(v.den) for v in dms_tag.values]
            dec = vals[0] + vals[1] / 60.0 + vals[2] / 3600.0
            ref = str(tags.get(ref_tag_name, "N"))
            return -dec if ref in ["S", "W"] else dec

        lat = dms_to_dec(tags[lat_key], "GPS GPSLatitudeRef")
        lon = dms_to_dec(tags[lon_key], "GPS GPSLongitudeRef")
        return lat, lon
    except Exception as e:
        logging.warning(f"EXIF read failed for {path.name}: {e}")
        return None

def build_grid(min_x: float, min_y: float, max_x: float, max_y: float, 
               tile_size: float, overlap_pct: float) -> List[Dict]:
    """Construct buffered spatial grid covering the dataset extent."""
    step = tile_size * (1.0 - overlap_pct)
    cols = math.ceil((max_x - min_x) / step) + 1
    rows = math.ceil((max_y - min_y) / step) + 1
    
    grid = []
    for r in range(rows):
        for c in range(cols):
            x0 = min_x + c * step
            y0 = min_y + r * step
            geom = box(x0, y0, x0 + tile_size, y0 + tile_size)
            grid.append({
                "tile_id": f"tile_{r:03d}_{c:03d}",
                "bounds": [x0, y0, x0 + tile_size, y0 + tile_size],
                "geometry": geom
            })
    return grid

def validate_chunk(chunk_images: List[Path], max_mem_mb: int, min_imgs: int) -> bool:
    """Enforce memory and density thresholds."""
    if len(chunk_images) < min_imgs:
        return False
    # Conservative memory estimate: ~100 MB per image for descriptor extraction
    estimated_mb = len(chunk_images) * 100
    return estimated_mb <= max_mem_mb

def main():
    args = parse_args()
    
    # Validation: enforce parameter bounds
    if not (0.15 <= args.overlap_pct <= 0.30):
        raise ValueError("--overlap-pct must be between 0.15 and 0.30")
    if args.tile_size_m < 200.0:
        raise ValueError("--tile-size-m must be >= 200m to maintain feature density")
    
    logging.info(f"Scanning {args.input_dir} for imagery...")
    image_paths = sorted(args.input_dir.glob("*.JPG")) + sorted(args.input_dir.glob("*.jpg")) + \
                  sorted(args.input_dir.glob("*.TIF")) + sorted(args.input_dir.glob("*.tif"))
    
    if not image_paths:
        raise FileNotFoundError("No supported imagery found in input directory.")

    # Extract coordinates & validate CRS tolerance
    coords = []
    valid_images = []
    for p in image_paths:
        gps = extract_gps_decimal(p)
        if gps:
            coords.append(gps)
            valid_images.append(p)
        else:
            logging.warning(f"Skipping {p.name}: missing or invalid GPS metadata")

    if len(valid_images) == 0:
        raise RuntimeError("No images with valid GPS coordinates. Cannot partition.")

    lats, lons = zip(*coords)
    # Simple UTM approximation for metric grid (production pipelines should use pyproj)
    # For this script, we treat decimal degrees as projected meters scaled by ~111,320 m/deg
    SCALE = 111320.0
    min_x, min_y = min(lons) * SCALE, min(lats) * SCALE
    max_x, max_y = max(lons) * SCALE, max(lats) * SCALE

    logging.info(f"Dataset extent: {len(valid_images)} valid images across "
                 f"{(max_x-min_x)/1000:.1f} x {(max_y-min_y)/1000:.1f} km")

    grid = build_grid(min_x, min_y, max_x, max_y, args.tile_size_m, args.overlap_pct)
    
    # Assign images to tiles
    chunk_manifest: Dict[str, List[str]] = {g["tile_id"]: [] for g in grid}
    for img, (lat, lon) in zip(valid_images, coords):
        pt = Point(lon * SCALE, lat * SCALE)
        for g in grid:
            if g["geometry"].contains(pt):
                chunk_manifest[g["tile_id"]].append(str(img))

    # Validate & prune underpopulated/overloaded chunks
    final_manifest = {}
    orphaned = []
    for tile_id, images in chunk_manifest.items():
        if validate_chunk(images, args.max_memory_mb, args.min_images):
            final_manifest[tile_id] = images
        else:
            orphaned.extend(images)
            logging.info(f"Chunk {tile_id} failed validation ({len(images)} images). Routing to fallback.")

    # Fallback routing: merge orphans into nearest valid chunk
    if orphaned:
        logging.warning(f"{len(orphaned)} images routed to fallback. Appending to largest valid chunk.")
        if final_manifest:
            largest_tile = max(final_manifest, key=lambda k: len(final_manifest[k]))
            final_manifest[largest_tile].extend(orphaned)
        else:
            raise RuntimeError("All chunks failed validation thresholds. Reduce --tile-size-m or --min-images.")

    # Write manifest
    args.output_manifest.parent.mkdir(parents=True, exist_ok=True)
    with open(args.output_manifest, "w") as f:
        json.dump({
            "metadata": {
                "tile_size_m": args.tile_size_m,
                "overlap_pct": args.overlap_pct,
                "total_images_processed": len(valid_images),
                "chunks_generated": len(final_manifest)
            },
            "chunks": final_manifest
        }, f, indent=2)

    logging.info(f"Manifest written to {args.output_manifest}")

if __name__ == "__main__":
    main()

python split_dataset.py \
  --input-dir /data/uav_survey_2024/raw \
  --output-manifest /data/uav_survey_2024/manifests/chunks_v1.json \
  --tile-size-m 800 \
  --overlap-pct 0.20 \
  --min-images 50 \
  --max-memory-mb 4096 \
  --crs-tolerance-deg 0.0001