Satellite Remote Sensing

import planetary_computer
import pystac_client
import rioxarray

# Search Sentinel-2 imagery via Microsoft Planetary Computer
catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

# Search for cloud-free imagery over a region
search = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=[11.0, 46.0, 12.0, 47.0],  # Tyrol, Austria
    datetime="2025-06-01/2025-08-31",
    query={"eo:cloud_cover": {"lt": 10}},
)

items = search.item_collection()
print(f"Found {len(items)} scenes with <10% cloud cover")

# Load a specific band as xarray DataArray
item = items[0]
red = rioxarray.open_rasterio(item.assets["B04"].href)
nir = rioxarray.open_rasterio(item.assets["B08"].href)

# Cloud masking for Sentinel-2 using the SCL band
import numpy as np

def mask_clouds_sentinel2(scl_band: np.ndarray) -> np.ndarray:
    """
    Create cloud mask from Sentinel-2 Scene Classification Layer.
    SCL values: 0=no_data, 1=saturated, 2=dark_area, 3=cloud_shadow,
    4=vegetation, 5=bare_soil, 6=water, 7=unclassified,
    8=cloud_medium, 9=cloud_high, 10=cirrus, 11=snow
    """
    cloud_classes = {0, 1, 3, 8, 9, 10}
    mask = np.isin(scl_band, list(cloud_classes))
    return mask  # True where clouds/invalid

import rasterio
from rasterio.merge import merge
from rasterio.warp import calculate_default_transform, reproject, Resampling

def reproject_raster(src_path: str, dst_path: str, dst_crs: str = "EPSG:4326"):
    """Reproject a raster to a target coordinate reference system."""
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds
        )
        kwargs = src.meta.copy()
        kwargs.update({
            "crs": dst_crs,
            "transform": transform,
            "width": width,
            "height": height,
        })
        with rasterio.open(dst_path, "w", **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.bilinear,
                )

def compute_indices(red: np.ndarray, nir: np.ndarray,
                    green: np.ndarray, swir: np.ndarray) -> dict:
    """
    Compute common spectral indices from surface reflectance bands.
    All inputs should be float arrays with values in [0, 1].
    """
    eps = 1e-10  # avoid division by zero
    ndvi = (nir - red) / (nir + red + eps)
    ndwi = (green - nir) / (green + nir + eps)
    nbr = (nir - swir) / (nir + swir + eps)
    evi = 2.5 * (nir - red) / (nir + 6 * red - 7.5 * 0.0001 + 1 + eps)
    savi = 1.5 * (nir - red) / (nir + red + 0.5 + eps)

    return {
        "NDVI": ndvi,   # vegetation vigor [-1, 1]
        "NDWI": ndwi,   # water bodies [-1, 1]
        "NBR": nbr,     # burn severity [-1, 1]
        "EVI": evi,     # enhanced vegetation
        "SAVI": savi,   # soil-adjusted vegetation
    }

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score

# Stack bands into feature array: (n_pixels, n_bands)
# training_labels: land cover classes from ground truth polygons
bands = np.stack([blue, green, red, nir, swir1, swir2, ndvi, ndwi], axis=-1)
n_rows, n_cols, n_bands = bands.shape
X = bands.reshape(-1, n_bands)

# Train Random Forest classifier
rf = RandomForestClassifier(n_estimators=200, max_depth=20, n_jobs=-1)
scores = cross_val_score(rf, X_train, y_train, cv=5, scoring="f1_macro")
print(f"5-fold F1: {scores.mean():.3f} +/- {scores.std():.3f}")

rf.fit(X_train, y_train)
classification = rf.predict(X).reshape(n_rows, n_cols)