Spaces:
Sleeping
Sleeping
| """ | |
| Curate real Delhi before/after pairs from free, openly-licensed Sentinel-2 | |
| imagery (Copernicus / AWS Open Data, via the Earth Search STAC API), and | |
| register them into docs/delhi_eval/manifest.json. | |
| Why: DDA (the client) hasn't handed off sample GeoTIFFs yet (see | |
| docs/IMPLEMENTATION_PLAN_DDA.md, "Blocked until DDA provides..."), and no | |
| imagery exists in library_sources/. This gets a *real*, freely-licensed Delhi | |
| dataset in place now so calibration/harness work isn't blocked on that | |
| handoff. Caveat: Sentinel-2 is 10m GSD — good for large land-use change | |
| (new colonies, vegetation loss, water bodies), NOT reliable for individual | |
| building/road-level detail. Treat this as a coarse-GSD complement to finer | |
| imagery (DDA's own GeoTIFFs, or drone/high-res satellite), not a replacement. | |
| GT masks are NOT generated here (diff-based "ground truth" would be circular | |
| and would bake in seasonal/illumination noise as if it were real change). | |
| Masks still need hand-labeling per docs/delhi_eval/README.md. | |
| Usage: | |
| python scripts/build_delhi_pairs_sentinel2.py --count 30 | |
| """ | |
| import argparse | |
| import json | |
| import os | |
| from pathlib import Path | |
| import numpy as np | |
| import rasterio | |
| import rasterio.warp | |
| from rasterio.windows import Window | |
| from rasterio.windows import transform as window_transform | |
| os.environ.setdefault("AWS_NO_SIGN_REQUEST", "YES") | |
| os.environ.setdefault("GDAL_DISABLE_READDIR_ON_OPEN", "EMPTY_DIR") | |
| ROOT = Path(__file__).resolve().parent.parent | |
| MANIFEST_PATH = ROOT / "docs" / "delhi_eval" / "manifest.json" | |
| LIBRARY_ROOT = ROOT / "library_sources" | |
| # Same MGRS tile (43RFM), near-zero cloud cover, same season (mid/late June) | |
| # ~7 years apart — season-matched so diffs reflect real structural change | |
| # rather than crop-calendar swings (Delhi's agri belt looks wildly different | |
| # in Nov vs Jul purely from monsoon/rabi-kharif cycles). Picked via STAC. | |
| BEFORE_URL = ("https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/" | |
| "43/R/FM/2019/6/S2A_43RFM_20190629_1_L2A/TCI.tif") | |
| AFTER_URL = ("https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/" | |
| "43/R/FM/2026/6/S2B_43RFM_20260617_0_L2A/TCI.tif") | |
| DATE_BEFORE = "2019-06-29" | |
| DATE_AFTER = "2026-06-17" | |
| # Delhi NCT bbox (west, south, east, north), clipped to this tile's footprint | |
| # (tile covers lon <= ~77.15, so far-east Delhi across the Yamuna is excluded). | |
| DELHI_BBOX_LONLAT = (76.84, 28.40, 77.15, 28.80) | |
| CELL_PX = 300 # ~3km x 3km at 10m GSD | |
| NODATA_FRACTION_LIMIT = 0.15 # skip cells that are >15% black fill | |
| def _lonlat_to_pixel_bbox(ds, bbox_lonlat): | |
| from rasterio.warp import transform_bounds | |
| west, south, east, north = transform_bounds("EPSG:4326", ds.crs, *bbox_lonlat) | |
| row_start, col_start = ds.index(west, north) | |
| row_stop, col_stop = ds.index(east, south) | |
| return max(0, row_start), min(ds.height, row_stop), max(0, col_start), min(ds.width, col_stop) | |
| def _read_cell(ds, row, col, size): | |
| window = Window(col, row, size, size) | |
| arr = ds.read([1, 2, 3], window=window) | |
| return np.transpose(arr, (1, 2, 0)) # HWC | |
| def main(): | |
| parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) | |
| parser.add_argument("--count", type=int, default=30) | |
| parser.add_argument("--cell-px", type=int, default=CELL_PX) | |
| args = parser.parse_args() | |
| print(f"Opening COGs (windowed HTTP reads, no full download)...\n before={BEFORE_URL}\n after={AFTER_URL}") | |
| with rasterio.open(BEFORE_URL) as before_ds, rasterio.open(AFTER_URL) as after_ds: | |
| assert before_ds.transform == after_ds.transform and before_ds.crs == after_ds.crs, \ | |
| "before/after grids are not pixel-aligned" | |
| row0, row1, col0, col1 = _lonlat_to_pixel_bbox(before_ds, DELHI_BBOX_LONLAT) | |
| size = args.cell_px | |
| candidates = [] | |
| rows = range(row0, row1 - size, size) | |
| cols = range(col0, col1 - size, size) | |
| print(f"Scanning {len(list(rows))}x{len(list(cols))} grid of {size}x{size}px cells for change signal...") | |
| for row in range(row0, row1 - size, size): | |
| for col in range(col0, col1 - size, size): | |
| before = _read_cell(before_ds, row, col, size) | |
| after = _read_cell(after_ds, row, col, size) | |
| nodata_frac = max( | |
| np.mean(np.all(before == 0, axis=-1)), | |
| np.mean(np.all(after == 0, axis=-1)), | |
| ) | |
| if nodata_frac > NODATA_FRACTION_LIMIT: | |
| continue | |
| diff_score = float(np.mean(np.abs(before.astype(np.int16) - after.astype(np.int16)))) | |
| green_before = before[..., 1].astype(np.int16) | |
| green_after = after[..., 1].astype(np.int16) | |
| veg_signal = float(np.mean(np.abs(green_before - green_after))) | |
| easting, northing = rasterio.transform.xy(before_ds.transform, row + size // 2, col + size // 2) | |
| lon, lat = rasterio.warp.transform(before_ds.crs, "EPSG:4326", [easting], [northing]) | |
| lon, lat = lon[0], lat[0] | |
| candidates.append({ | |
| "row": row, "col": col, "diff_score": diff_score, | |
| "veg_signal": veg_signal, "lon": lon, "lat": lat, | |
| "before": before, "after": after, | |
| }) | |
| print(f"{len(candidates)} valid (non-nodata) cells scanned.") | |
| candidates.sort(key=lambda c: -c["diff_score"]) | |
| selected = candidates[: args.count] | |
| manifest = json.loads(MANIFEST_PATH.read_text(encoding="utf-8")) if MANIFEST_PATH.exists() else {"pairs": []} | |
| existing_ids = {p["pair_id"] for p in manifest["pairs"]} | |
| next_num = len(manifest["pairs"]) + 1 | |
| before_dir = LIBRARY_ROOT / DATE_BEFORE[:4] | |
| after_dir = LIBRARY_ROOT / DATE_AFTER[:4] | |
| before_dir.mkdir(parents=True, exist_ok=True) | |
| after_dir.mkdir(parents=True, exist_ok=True) | |
| base_crs = before_ds.crs | |
| added = 0 | |
| for cand in selected: | |
| while f"delhi_{next_num:04d}" in existing_ids: | |
| next_num += 1 | |
| pair_id = f"delhi_{next_num:04d}" | |
| fname = f"sentinel2_delhi_r{cand['row']}_c{cand['col']}.tif" | |
| transform = window_transform(Window(cand["col"], cand["row"], size, size), before_ds.transform) | |
| profile = dict( | |
| driver="GTiff", width=size, height=size, count=3, dtype="uint8", | |
| crs=base_crs, transform=transform, compress="deflate", | |
| ) | |
| before_path = before_dir / fname | |
| after_path = after_dir / fname | |
| with rasterio.open(before_path, "w", **profile) as dst: | |
| dst.write(np.transpose(cand["before"], (2, 0, 1))) | |
| with rasterio.open(after_path, "w", **profile) as dst: | |
| dst.write(np.transpose(cand["after"], (2, 0, 1))) | |
| change_types = ["mixed_gsd"] | |
| if cand["veg_signal"] > cand["diff_score"] * 0.6: | |
| change_types.append("vegetation") | |
| manifest["pairs"].append({ | |
| "pair_id": pair_id, | |
| "before_path": str(before_path.relative_to(ROOT)), | |
| "after_path": str(after_path.relative_to(ROOT)), | |
| "date_before": DATE_BEFORE, | |
| "date_after": DATE_AFTER, | |
| "gsd": 10.0, | |
| "zone": f"Delhi ({cand['lat']:.4f}N, {cand['lon']:.4f}E)", | |
| "change_types": change_types, | |
| "gt_mask": None, | |
| "notes": "Sentinel-2 L2A TCI, 10m GSD — coarse resolution, NOT reliable for individual " | |
| "building/road detail. Real Delhi coordinates/dates. GT mask still needs hand-labeling.", | |
| }) | |
| existing_ids.add(pair_id) | |
| added += 1 | |
| print(f" {pair_id} ({cand['lat']:.4f}N, {cand['lon']:.4f}E) diff_score={cand['diff_score']:.1f}" | |
| f" tags={change_types}") | |
| MANIFEST_PATH.write_text(json.dumps(manifest, indent=2), encoding="utf-8") | |
| print(f"\nAdded {added} pair(s). Manifest now has {len(manifest['pairs'])} total.") | |
| print(f"Images written to {before_dir.relative_to(ROOT)}/ and {after_dir.relative_to(ROOT)}/") | |
| if __name__ == "__main__": | |
| main() | |