remote-sensing-analysis

name	remote-sensing-analysis
description	Remote sensing pipeline with optical/SAR processing and ML classification. Use for satellite imagery analysis.
context	fork
globs	["**/*.tif","**/*.tiff","**/*.nc","**/scripts/rs_*.py"]
alwaysApply	false

Skill: Remote Sensing Analysis

Purpose

This skill equips any agent with expert technical knowledge in complete remote sensing, including:

• Optical image processing (multispectral, hyperspectral, panchromatic)
• SAR processing (radiometric calibration, speckle filtering, InSAR, polarimetric decomposition)
• Spectral index calculation and band math
• Radiometric and atmospheric correction
• Change detection (LandTrendr, CCDC, BFAST, neural networks)
• Time series: smoothing, gap-filling, phenological analysis
• ML/DL classification (Random Forest, U-Net, SegFormer, OBIA)
• Object detection in satellite/drone imagery
• Pansharpening and sensor fusion
• Photogrammetry and point cloud processing (LiDAR/SfM)
• LULC product analysis (MapBiomas, CGLS, ESRI, GLAD, ESA WorldCover)
• Spatial statistics and inter-product agreement metrics
• Accuracy assessment (Olofsson 2014) and spatial statistics
• Scientific literature search in indexed journals

---

1. Raster Image Processing

1.1 Reading, Reprojecting and Aligning

import rasterio
import numpy as np
from rasterio.warp import calculate_default_transform, reproject, Resampling
from pathlib import Path

def read_raster_safe(path: Path, band: int = 1) -> tuple[np.ndarray, dict]:
    """Safe raster reading with nodata handling."""
    with rasterio.open(path) as src:
        data = src.read(band)
        profile = src.profile.copy()
        nodata = src.nodata
        if nodata is not None:
            data = np.where(data == nodata, np.nan, data.astype(float))
    return data, profile

def reproject_match(src_path: Path, ref_path: Path,
                    resampling: str = 'nearest') -> tuple[np.ndarray, dict]:
    """
    Reproject raster to the CRS and grid of a reference raster.
    IMPORTANT: classifications → use 'nearest' (not 'bilinear').
    Reference: Foody (2002) RSE 80(1):185
    """
    from rasterio.enums import Resampling as RS
    method = getattr(RS, resampling)
    with rasterio.open(ref_path) as ref:
        ref_crs, ref_transform = ref.crs, ref.transform
        ref_height, ref_width = ref.height, ref.width
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, ref_crs, ref_width, ref_height,
            left=ref.bounds.left, bottom=ref.bounds.bottom,
            right=ref.bounds.right, top=ref.bounds.top)
        kwargs = src.meta.copy()
        kwargs.update({'crs': ref_crs, 'transform': transform,
                       'width': ref_width, 'height': ref_height})
        data = np.empty((src.count, ref_height, ref_width), dtype=src.dtypes[0])
        reproject(source=rasterio.band(src, list(range(1, src.count + 1))),
                  destination=data, src_transform=src.transform,
                  src_crs=src.crs, dst_transform=transform,
                  dst_crs=ref_crs, resampling=method)
    return data, kwargs

1.2 CRS Recommendations

CRS_RECOMMENDATIONS = {
    'global_analysis':    'ESRI:54012',   # Eckert IV (equal-area)
    'south_america':      'EPSG:4087',    # World Equidistant Cylindrical
    'brazil_national':    'EPSG:4674',    # SIRGAS 2000
    'landsat_default':    'EPSG:4326',    # WGS84 geographic
    'sentinel_utm_z22s':  'EPSG:32722',   # UTM Zone 22S
}
# ⚠️ ALWAYS verify CRS before spatial operations

---

2. Radiometric and Atmospheric Correction

2.1 Pipeline DN → Surface Reflectance

DN (Digital Number)
  ↓ radiometric calibration
Radiance (W/m²/sr/μm)
  ↓ divide by solar irradiance + cos(SZA)
TOA Reflectance (%)
  ↓ atmospheric correction (6S, Sen2Cor, LEDAPS, LaSRC)
Surface Reflectance (%)
  ↓ topographic correction (optional)
Normalised Reflectance

2.2 Atmospheric Correction Tools

Sensor	Tool	Method
Sentinel-2	Sen2Cor	LUT (Look-up Table)
Landsat 8/9 OLI	LaSRC	6SV radiative transfer
Landsat 5/7 TM/ETM+	LEDAPS	6S
Multi-sensor	Py6S	6S Python wrapper
Aquatic	ACOLITE	Dark Spectrum Fitting
Simple / relative	DOS (Dark Object Subtraction)	Empirical

# Example: TOA reflectance for Landsat 8 OLI
def dn_to_toa_reflectance(dn: np.ndarray, M_rho: float, A_rho: float,
                           sun_elevation_deg: float) -> np.ndarray:
    """
    Convert DN → TOA Reflectance (Landsat Collection 2).
    M_rho, A_rho: coefficients from MTL file (REFLECTANCE_MULT/ADD_BAND_x)
    Reference: USGS Landsat Collection 2 Science Product Guide
    """
    toa = M_rho * dn + A_rho
    sun_elevation_rad = np.deg2rad(sun_elevation_deg)
    return toa / np.sin(sun_elevation_rad)

---

3. Spectral Indices and Band Math

3.1 Core Indices

def compute_index(band_a: np.ndarray, band_b: np.ndarray,
                  formula: str = 'ndvi') -> np.ndarray:
    """Compute normalised spectral indices."""
    eps = 1e-8
    if formula == 'ndvi':     # (NIR - Red) / (NIR + Red)
        return (band_a - band_b) / (band_a + band_b + eps)
    elif formula == 'ndwi':   # (Green - NIR) / (Green + NIR)  Gao 1996
        return (band_b - band_a) / (band_b + band_a + eps)
    elif formula == 'nbr':    # (NIR - SWIR2) / (NIR + SWIR2)
        return (band_a - band_b) / (band_a + band_b + eps)

# Quick reference by category:
SPECTRAL_INDICES = {
    # Vegetation
    'NDVI':   '(NIR - Red) / (NIR + Red)',          # Tucker 1979
    'EVI':    '2.5*(NIR-Red)/(NIR+6*Red-7.5*Blue+1)', # Liu & Huete 1995
    'SAVI':   '(NIR-Red)/(NIR+Red+L)*(1+L)',        # Huete 1988, L=0.5
    'MSAVI':  '(2*NIR+1 - sqrt((2*NIR+1)²-8*(NIR-Red)))/2', # Qi 1994
    'NDRE':   '(RedEdge - Red) / (RedEdge + Red)',  # Gitelson 1994
    'LAI':    'from NDVI via empirical model',
    # Water
    'NDWI':   '(Green - NIR) / (Green + NIR)',      # McFeeters 1996
    'MNDWI':  '(Green - SWIR1) / (Green + SWIR1)', # Xu 2006
    'AWEI':   '4*(Green-SWIR1)-(0.25*NIR+2.75*SWIR2)', # Feyisa 2014
    # Urban / Bare Soil
    'NDBI':   '(SWIR1 - NIR) / (SWIR1 + NIR)',     # Zha 2003
    'BSI':    '((SWIR1+Red)-(NIR+Blue))/((SWIR1+Red)+(NIR+Blue))', # Rikimaru 2002
    # Fire
    'NBR':    '(NIR - SWIR2) / (NIR + SWIR2)',      # Key & Benson 2006
    'dNBR':   'pre_NBR - post_NBR',
    'RBR':    'dNBR / (pre_NBR + 1.001)',           # Parks 2014
    # SAR
    'VH/VV':  'VH / VV (log scale)',
    'RVI':    '4*VH / (VV + VH)',                   # Kim 2012
}

---

4. SAR Processing

4.1 Preprocessing Pipeline

SAR Image (SLC / GRD)
  ↓ Radiometric calibration (sigma0 / beta0 / gamma0)
  ↓ Speckle filtering (Lee / Refined Lee / Gamma-MAP)
  ↓ Terrain correction (RTC — Range-Doppler / Ellipsoid)
  ↓ Convert to dB: 10 * log10(sigma0)
  ↓ Geocoding (target EPSG)
ARD Product (Analysis-Ready SAR)

4.2 Speckle Filtering

from scipy.ndimage import uniform_filter
import numpy as np

def lee_filter(img: np.ndarray, size: int = 7) -> np.ndarray:
    """
    Lee filter for speckle reduction in SAR images.
    Reference: Lee (1980) IEEE Trans. Pattern Anal. Mach. Intell. 2(2):165
    """
    img_mean = uniform_filter(img, size)
    img_sq_mean = uniform_filter(img**2, size)
    img_variance = img_sq_mean - img_mean**2
    overall_variance = np.var(img)
    img_weights = img_variance / (img_variance + overall_variance)
    return img_mean + img_weights * (img - img_mean)

4.3 Polarimetric Decomposition (reference table)

Decomposition	Components	Typical Application
Freeman-Durden	Volume + Double-bounce + Odd	Forest, LULC
Touzi	Symmetric + Anti-symmetric	Wetlands
Cloude-Pottier	Entropy H, Anisotropy A, Alpha α	OBIA, classification
Yamaguchi (4-comp)	+ Helix	Urban areas

---

5. Change Detection

5.1 Main Approaches

Method	Type	Data	Reference
Image Differencing	Pixel	Optical	Singh 1989
CVA (Change Vector Analysis)	Pixel	Multi-band	Malila 1980
MAD / IR-MAD	Statistical	Multi-temporal	Nielsen 2007
LandTrendr	Temporal	Annual Landsat	Kennedy 2010
CCDC	Temporal	Dense Landsat	Zhu & Woodcock 2014
BFAST	Temporal	Any TS	Verbesselt 2010
BISE	Gaps	MODIS/NDVI	Lovell 2008
Siamese CNN	DL	Optical/SAR	Daudt 2018
ChangeFormer	DL	Optical	Bandara 2022

def image_differencing(img_t1: np.ndarray, img_t2: np.ndarray,
                        threshold: float = None) -> np.ndarray:
    """
    Change detection by simple image differencing.
    threshold: if None, uses mean + 2*std (automatic).
    """
    diff = np.abs(img_t2.astype(float) - img_t1.astype(float))
    if threshold is None:
        threshold = diff.mean() + 2 * diff.std()
    return diff > threshold

---

6. Time Series Analysis

6.1 Smoothing and Gap-filling

from scipy.signal import savgol_filter
import numpy as np

def savitzky_golay_smooth(ts: np.ndarray, window: int = 7,
                           polyorder: int = 3) -> np.ndarray:
    """
    Savitzky-Golay smoothing for vegetation time series.
    Preserves peaks and troughs better than moving averages.
    Reference: Chen et al. (2004) RSE 91(3-4):332
    """
    return savgol_filter(ts, window_length=window, polyorder=polyorder)

def harmonic_regression(t: np.ndarray, y: np.ndarray,
                         n_harmonics: int = 3) -> np.ndarray:
    """
    Harmonic regression for seasonality modelling (phenology).
    Reference: Zhu & Woodcock (2014) RSE 144:152
    """
    X = [np.ones_like(t), t]
    for k in range(1, n_harmonics + 1):
        X.append(np.cos(2 * np.pi * k * t / 365))
        X.append(np.sin(2 * np.pi * k * t / 365))
    X = np.column_stack(X)
    coeffs, _, _, _ = np.linalg.lstsq(X, y, rcond=None)
    return X @ coeffs

6.2 Temporal Frequency Metrics

def temporal_frequency(stack: np.ndarray, class_value: int,
                        nodata: int = 0, axis: int = 0) -> np.ndarray:
    """
    Compute the relative frequency of a class in a time series.
    Reference: Defries et al. (2004); Hansen et al. (2013) Science
    """
    valid_mask = stack != nodata
    class_mask = stack == class_value
    return np.where(
        valid_mask.sum(axis=axis) > 0,
        class_mask.sum(axis=axis) / valid_mask.sum(axis=axis),
        np.nan)

def temporal_stability_index(stack: np.ndarray, nodata: int = 0,
                               axis: int = 0) -> np.ndarray:
    """
    Temporal stability index (dominance of the most frequent class).
    Reference: Pérez-Hoyos et al. (2017) Remote Sensing 9(1):2
    """
    from scipy.stats import mode
    dominant_class, dominant_count = mode(stack, axis=axis)
    valid_count = (stack != nodata).sum(axis=axis)
    return np.where(valid_count > 0, dominant_count / valid_count, np.nan)

---

7. ML/DL Classification

7.1 Spatial Validation (No Data Leakage)

# ⚠️ NEVER use random split with spatial data — spatial autocorrelation
# will inflate accuracy. Use block spatial cross-validation.

from sklearn.model_selection import KFold
import geopandas as gpd

def spatial_block_cv(gdf: gpd.GeoDataFrame, n_splits: int = 5,
                     block_size_km: float = 50.0):
    """
    Cross-validation with spatial blocks to avoid data leakage.
    Reference: Roberts et al. (2017) Ecography 40(8):913
    
    Alternative: use scikit-spatial-models or spatial-cross-val.
    """
    # Generates a block grid over the dataset extent
    # and creates folds based on blocks (not samples)
    raise NotImplementedError("See implementation by Roberts et al. 2017")

7.2 Classification Quality Checklist

preprocessing:
  - [ ] Proper normalisation for RS (divide by 10000 for Landsat/S2 SR)
  - [ ] Balanced classes or class_weight handling
  - [ ] Spatially consistent augmentation (flips, rotations)
  - [ ] Bands in the correct order and documented

training:
  - [ ] Spatial cross-validation (not random split)
  - [ ] Per-class metrics: IoU, F1, precision, recall
  - [ ] Confusion matrix with % and absolutes
  - [ ] Learning curves plotted

validation:
  - [ ] Independent samples (not used in fit)
  - [ ] Stratified by class + region
  - [ ] Area-adjusted accuracy (Olofsson et al. 2014)

---

8. LULC Ensemble Methods

def majority_vote_ensemble(products: list[np.ndarray],
                            nodata: int = 0) -> np.ndarray:
    """
    Simple majority voting among multiple products.
    Reference: Herold et al. (2011) RSE 115(6):1559
    """
    stack = np.stack(products, axis=0)
    from scipy.stats import mode
    result, _ = mode(stack, axis=0)
    return result.squeeze()

def weighted_ensemble(products: list[np.ndarray],
                      weights: list[float]) -> np.ndarray:
    """
    Weighted ensemble by product confidence (one-hot + weighted sum).
    Reference: Fritz et al. (2011) ERL 6(4):044005
    """
    raise NotImplementedError("Requires one-hot encoding of categories")

---

9. Inter-Product Agreement Metrics

import numpy as np
from sklearn.metrics import cohen_kappa_score

def overall_accuracy(y_true, y_pred):
    """Proportion of agreeing pixels."""
    return np.mean(y_true == y_pred)

def cohens_kappa(y_true, y_pred):
    """Chance-corrected kappa — Cohen (1960)."""
    return cohen_kappa_score(y_true, y_pred)

def dice_coefficient(mask_a: np.ndarray, mask_b: np.ndarray) -> float:
    """Dice/Sørensen coefficient for binary classes."""
    intersection = np.logical_and(mask_a, mask_b).sum()
    return 2 * intersection / (mask_a.sum() + mask_b.sum() + 1e-8)

---

10. Accuracy Assessment (Olofsson 2014)

def area_adjusted_accuracy(confusion_matrix: np.ndarray,
                            mapped_areas: np.ndarray,
                            sample_counts: np.ndarray) -> dict:
    """
    Area-adjusted accuracy and area estimates (design-based).
    
    Reference: Olofsson et al. (2014) RSE 148:42
    REQUIRED for land use/cover change studies.
    
    Args:
        confusion_matrix: (n_classes × n_classes) sample counts
        mapped_areas: mapped area per class (same units)
        sample_counts: samples collected per stratum
    
    Returns:
        dict with adjusted OA, UA, PA, F1 + confidence intervals
    """
    # Implementation following equations 1-5 of Olofsson et al. (2014)
    raise NotImplementedError("Use `area` package or reference implementation")

# Olofsson 2014 CHECKLIST:
# [ ] Stratified sampling (strata = mapped classes)
# [ ] Sample size per stratum documented
# [ ] Confusion matrix in area proportions (not raw counts)
# [ ] Uncertainty (95% CI) reported for OA, UA, PA
# [ ] Estimated area with standard error and CI

---

11. LULC Products — Reference Tables

11.1 Common Class Hierarchy

Class	MapBiomas	CGLS	ESRI	GLAD
Forest	1	111,112,113	2	1
Savanna / Cerrado	3	121,122	—	3
Agriculture	18,39,20,21	40	5	6
Grassland	15	30	7	4
Urban	24,25	50	7	7
Water	33	80	1	8
Bare Soil	22,23	200	8	—

11.2 Resolution and Temporal Coverage

Product	Spatial Resolution	Temporal Resolution	Coverage
MapBiomas	30m (Landsat)	Annual (1985–present)	Brazil + South America
CGLS-LC100	100m	Annual (2015–2019)	Global
ESRI Land Cover	10m (Sentinel-2)	Annual (2017–present)	Global
GLAD ARD	30m (Landsat)	Annual (2000–present)	Global
ESA WorldCover	10m (Sentinel-1+2)	2020, 2021	Global

---

12. Quality Checklist for RS Analysis

Preprocessing

• CRS verified and consistent across all products
• Spatial resolution aligned (reprojection with appropriate method)
• Nodata handled explicitly (do not assume 0 = nodata)
• Geographic extent equal for all products
• Atmospheric correction documented and consistent with sensor

Agreement Analysis / Validation

• Reclassification dictionary documented and versioned
• Area-adjusted accuracy (Olofsson 2014)
• Spatial cross-validation (not random)
• 95% CI for all reported metrics
• Sample N reported with metrics

Time Series

• Cloud/shadow mask applied before smoothing
• Smoothing method justified (Savitzky-Golay vs Whittaker vs harmonic)
• Gap-filling documented
• Seasonality visually inspected

ML/DL

• Spatial cross-validation (not random split)
• Normalisation consistent between training and inference
• Class balancing documented
• Overfitting checked (train vs val curves)

---

13. Data Sources and APIs

# Scientific literature
semantic_scholar_api: "https://api.semanticscholar.org/graph/v1/paper/search?query={TERM}&fields=title,authors,year,abstract,citationCount"
crossref_api:         "https://api.crossref.org/works?query={TERM}&filter=type:journal-article,from-pub-date:2018"
mdpi_remote_sensing:  "https://www.mdpi.com/search?q={TERM}&journal=remotesensing"
eartharxiv:           "https://eartharxiv.org/search/?q={TERM}"

# Product documentation
mapbiomas_github:     "https://github.com/mapbiomas"
copernicus_land:      "https://land.copernicus.eu/global"
glad_umd:             "https://glad.umd.edu/dataset"
esa_worldcover:       "https://esa-worldcover.org"

# Data and catalogs
earthengine_catalog:  "https://developers.google.com/earth-engine/datasets"
stac_planetary:       "https://planetarycomputer.microsoft.com/catalog"
stac_element84:       "https://earth-search.aws.element84.com/v1"
usgs_earthexplorer:   "https://earthexplorer.usgs.gov"
copernicus_hub:       "https://scihub.copernicus.eu"

---

STAC MCP Server Integration

Discover and Query Satellite Data via MCP

The STAC (SpatioTemporal Asset Catalogs) protocol can be accessed through MCP servers, enabling programmatic satellite data discovery.

# Query STAC catalog via MCP tools
stac_query = {
    "collections": ["sentinel-2-l2a", "landsat-c2-l2"],
    "bbox": [-64.5, -12.5, -64.3, -12.3],  # Rondônia, Brazil
    "datetime": "2025-01-01/2025-12-31",
    "max_cloud_cover": 20
}

# Via MCP STAC server
items = await stac_search(stac_query)

# Returns list of items with: 
# - scene_id, datetime, cloud_cover, geometry
# - assets: B02, B03, B04, B08, B11, SCL, etc.
# - eo:bands with common names, center wavelength

Earth Engine Integration via MCP

# Google Earth Engine computation via MCP
gee_task = {
    "collection": "LANDSAT/LC08/C02/T1_L2",
    "region": {
        "type": "Point",
        "coordinates": [-64.4, -12.4]
    },
    "reducers": ["median"],
    "bands": ["SR_B4", "SR_B5"],  # Red, NIR
    "start_date": "2024-01-01",
    "end_date": "2024-12-31",
    "scale": 30
}

# Compute NDVI directly in Earth Engine
ndvi_result = await ee_compute_ndvi(gee_task)
# Returns computed raster data as GeoTIFF reference

STAC Browser Pattern for Data Discovery

# MCP-based STAC browser facade
class STACBrowser:
    """Discover satellite data across multiple catalogs via MCP."""
    
    CATALOGS = {
        "earthsearch": "https://earth-search.aws.element84.com/v1",
        "planetary_computer": "https://planetarycomputer.microsoft.com/api/stac/v1",
        "copernicus": "https://catalogue.dataspace.copernicus.eu/stac",
    }
    
    async def search_all(self, aoi: dict, date_range: tuple[str, str]) -> list[dict]:
        """Search across all catalogs in parallel."""
        results = await asyncio.gather(*[
            self._search_catalog(catalog, aoi, date_range)
            for catalog in self.CATALOGS
        ], return_exceptions=True)
        return self._merge_results(results)

Automated Literature Review with LLM

Paper Discovery & Synthesis

class LiteratureReviewAgent:
    """Automated literature discovery and synthesis for remote sensing."""
    
    @staticmethod
    async def search_papers(topic: str, max_results: int = 20) -> list[dict]:
        """Search across academic sources for relevant papers."""
        # Uses Semantic Scholar API
        # Uses arXiv API  
        # Uses CrossRef API
        # All accessed via MCP search tools
        pass
    
    @staticmethod
    async def synthesize_findings(papers: list[dict], focus: str) -> str:
        """Synthesize multiple papers into a coherent summary focused on a topic."""
        # LLM reads abstracts and full text where available
        # Produces structured synthesis: methods, results, limitations
        # Cross-references findings with current codebase patterns
        pass

Methodology Recommendation Engine

class MethodologyRecommender:
    """Recommend remote sensing methodologies based on problem description."""
    
    PROBLEM_PROMPT = """Given this remote sensing problem:
    
Problem: {problem_description}
Available data: {data_description}
Study area: {study_area}
    
Research and recommend:
1. Best classification/regression approach for this context
2. Preprocessing steps required
3. Accuracy assessment methodology (Olofsson 2014 stratified)
4. Inter-product agreement strategy if using multiple LULC products
5. Relevant papers and their key findings
"""
    # Uses LLM to search literature and synthesize recommendations
    # Cross-validates against known best practices in the codebase