| """Re-score a Blue Marble Clouds model on CloudSEN12+ (L1C, official high-quality test split: 975 images, 4 classes). |
| |
| pip install numpy onnxruntime rasterio huggingface_hub |
| python evaluate.py --model medium/model_fp32.onnx # or medium/model_int8.onnx |
| |
| Data: CloudSEN12+ test_509_high.mlstac, about 2.7 GB, downloaded once.""" |
| import numpy as np |
|
|
| |
| MEAN = np.array([2357.089, 2137.385, 2018.788, 2082.986, 2295.651, 2854.537, 3122.849, 3040.560, 3306.481, 1473.847, 506.070, 2472.825, 1838.929], np.float32) |
| STD = np.array([1624.683, 1675.806, 1557.708, 1833.702, 1823.738, 1733.977, 1732.131, 1679.732, 1727.26, 1024.687, 442.165, 1331.411, 1160.419], np.float32) |
|
|
|
|
| def normalize(dn, bands=None): |
| """dn: [C, H, W] digital numbers for the 13 bands, or for the band indices in `bands`. Missing bands are 0 after normalisation.""" |
| out = np.zeros((13, *dn.shape[1:]), np.float32) |
| idx = list(range(13)) if bands is None else list(bands) |
| out[idx] = (dn.astype(np.float32) - MEAN[idx, None, None]) / STD[idx, None, None] |
| return out |
|
|
|
|
| def session(path, threads=0): |
| import onnxruntime as ort |
| so = ort.SessionOptions() |
| if threads: so.intra_op_num_threads = threads |
| prov = [p for p in ("CUDAExecutionProvider", "CPUExecutionProvider") if p in ort.get_available_providers()] |
| return ort.InferenceSession(path, so, providers=prov) |
|
|
|
|
| def predict(sess, x): |
| """Whole image in one pass: pad H and W to multiples of 16 by edge replication, crop the logits back.""" |
| h, w = x.shape[-2:] |
| xp = np.pad(x, ((0, 0), (0, (-h) % 16), (0, (-w) % 16)), mode="edge")[None].astype(np.float32) |
| return sess.run(None, {sess.get_inputs()[0].name: xp})[0][0, :, :h, :w].argmax(0) |
|
|
|
|
| class MIoU: |
| def __init__(self, classes): |
| self.classes, n = classes, len(classes) |
| self.cm = np.zeros((n, n), np.int64) |
|
|
| def add(self, y, p): |
| n = len(self.classes); m = y >= 0 |
| self.cm += np.bincount(n * y[m].astype(np.int64) + p[m], minlength=n * n).reshape(n, n) |
|
|
| def report(self, title): |
| inter = np.diag(self.cm).astype(float); union = self.cm.sum(0) + self.cm.sum(1) - inter |
| iou = inter / np.maximum(union, 1) |
| print(f"{title}: mIoU {100 * iou.mean():.2f} (" + ", ".join(f"{c} {100 * v:.2f}" for c, v in zip(self.classes, iou)) + ")", flush=True) |
| return float(iou.mean()) |
|
|
|
|
| import argparse, json, os, struct |
| import numpy as np |
| from rasterio.io import MemoryFile |
|
|
|
|
| def fetch(data): |
| p = f"{data}/test_509_high.mlstac" |
| if os.path.exists(p): return p |
| from huggingface_hub import hf_hub_download |
| return hf_hub_download("isp-uv-es/CloudSEN12Plus", "test/test_509_high.mlstac", repo_type="dataset", local_dir=data) |
|
|
|
|
| def images(path): |
| """Container: b'#y', uint64 index length, JSON index {key: [offset, length]}, then one GeoTIFF per image |
| (bands 1-13: L1C B01-B12, band 14: expert label 0 clear, 1 thick cloud, 2 thin cloud, 3 cloud shadow).""" |
| with open(path, "rb") as f: |
| assert f.read(2) == b"#y" |
| n = struct.unpack("<Q", f.read(8))[0]; idx = json.loads(f.read(n)); off = 10 + n |
| for k in sorted(idx): |
| start, length = idx[k][:2] |
| f.seek(off + start) |
| with MemoryFile(f.read(length)) as mf, mf.open() as s: a = s.read() |
| yield a[:13].astype(np.uint16), a[13].astype(np.int64) |
|
|
|
|
| def main(): |
| ap = argparse.ArgumentParser(); ap.add_argument("--model", required=True); ap.add_argument("--data", default="data") |
| a = ap.parse_args(); os.makedirs(a.data, exist_ok=True) |
| sess = session(a.model); m = MIoU(["clear", "thick cloud", "thin cloud", "cloud shadow"]); n = 0 |
| for dn, y in images(fetch(a.data)): |
| m.add(y, predict(sess, normalize(dn))); n += 1 |
| m.report(f"CloudSEN12+ test ({n} images)") |
|
|
|
|
| if __name__ == "__main__": main() |
|
|
