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