|
| 1 | +--- |
| 2 | +tags: |
| 3 | +- remote_sensing |
| 4 | +- sentinel_streaming |
| 5 | +- streamable |
| 6 | +- innovation-summit-2025 |
| 7 | +--- |
| 8 | + |
| 9 | +Sentinel-2 STAC Quicklook Streaming |
| 10 | +================ |
| 11 | +ESIIL, 2024-06-01 |
| 12 | + |
| 13 | +This snippet streams a low-resolution overview of Sentinel-2 imagery |
| 14 | +using the STAC API and GDAL's `/vsicurl` interface so that no local |
| 15 | +GeoTIFF download is required. |
| 16 | + |
| 17 | +We used the following package versions when testing the code: |
| 18 | + |
| 19 | +- requests==2.31.0 |
| 20 | +- rasterio==1.3.9 |
| 21 | +- numpy==1.26.4 |
| 22 | +- pandas==2.2.2 |
| 23 | +- matplotlib==3.8.4 |
| 24 | + |
| 25 | +```python |
| 26 | +import json, os, io |
| 27 | +from typing import List, Dict, Any, Optional, Tuple |
| 28 | +import requests, rasterio |
| 29 | +import numpy as np |
| 30 | +import pandas as pd |
| 31 | +import matplotlib.pyplot as plt |
| 32 | + |
| 33 | + |
| 34 | +def _ensure_ok(r: requests.Response): |
| 35 | + try: |
| 36 | + r.raise_for_status() |
| 37 | + except requests.HTTPError as e: |
| 38 | + msg = "" |
| 39 | + try: |
| 40 | + msg = "\nServer says: " + json.dumps(r.json(), indent=2)[:1000] |
| 41 | + except Exception: |
| 42 | + msg = "\nServer says: " + (r.text[:1000] if r.text else "<no body>") |
| 43 | + raise requests.HTTPError(f"{e}{msg}") from None |
| 44 | + |
| 45 | + |
| 46 | +def search_stac_first_item_resilient( |
| 47 | + stac_api: str, |
| 48 | + collections: List[str], |
| 49 | + datetime_range: str, |
| 50 | + bbox: Optional[List[float]] = None, |
| 51 | + cloud_lt: Optional[float] = 20, |
| 52 | +) -> Dict[str, Any]: |
| 53 | + """ |
| 54 | + Robust STAC search: try 'query', then no filter, then CQL2-JSON. |
| 55 | + Returns the first Feature; raises if none. |
| 56 | + """ |
| 57 | + url = stac_api.rstrip("/") + "/search" |
| 58 | + base = {"collections": collections, "datetime": datetime_range, "limit": 1} |
| 59 | + if bbox is not None: |
| 60 | + base["bbox"] = bbox |
| 61 | + |
| 62 | + if cloud_lt is not None: |
| 63 | + payload1 = {**base, "query": {"eo:cloud_cover": {"lt": cloud_lt}}} |
| 64 | + r = requests.post(url, json=payload1, timeout=90, headers={"Content-Type": "application/json"}) |
| 65 | + if r.status_code < 400: |
| 66 | + _ensure_ok(r) |
| 67 | + feats = r.json().get("features", []) |
| 68 | + if feats: |
| 69 | + return feats[0] |
| 70 | + |
| 71 | + r = requests.post(url, json=base, timeout=90, headers={"Content-Type": "application/json"}) |
| 72 | + if r.status_code < 400: |
| 73 | + _ensure_ok(r) |
| 74 | + feats = r.json().get("features", []) |
| 75 | + if feats: |
| 76 | + return feats[0] |
| 77 | + |
| 78 | + if cloud_lt is not None: |
| 79 | + payload3 = { |
| 80 | + **base, |
| 81 | + "filter-lang": "cql2-json", |
| 82 | + "filter": { |
| 83 | + "op": "<", |
| 84 | + "args": [ |
| 85 | + {"property": "eo:cloud_cover"}, |
| 86 | + cloud_lt |
| 87 | + ], |
| 88 | + }, |
| 89 | + } |
| 90 | + r = requests.post(url, json=payload3, timeout=90, headers={"Content-Type": "application/json"}) |
| 91 | + _ensure_ok(r) |
| 92 | + feats = r.json().get("features", []) |
| 93 | + if feats: |
| 94 | + return feats[0] |
| 95 | + |
| 96 | + _ensure_ok(r) |
| 97 | + raise ValueError("Search succeeded but returned no features for the given parameters.") |
| 98 | + |
| 99 | + |
| 100 | +def choose_stac_asset(assets: Dict[str, Dict[str, Any]], asset_key: Optional[str] = None) -> Tuple[str, str]: |
| 101 | + if asset_key and asset_key in assets and assets[asset_key].get("href"): |
| 102 | + return asset_key, assets[asset_key]["href"] |
| 103 | + |
| 104 | + def score(k, a): |
| 105 | + href = (a.get("href") or "").lower() |
| 106 | + roles = [r.lower() for r in a.get("roles", [])] |
| 107 | + t = (a.get("type") or a.get("media_type") or "").lower() |
| 108 | + s = 0 |
| 109 | + if "cog" in json.dumps(a).lower() or href.endswith((".tif", ".tiff")): |
| 110 | + s += 5 |
| 111 | + if "data" in roles: |
| 112 | + s += 3 |
| 113 | + if "image/tiff" in t or "geotiff" in t: |
| 114 | + s += 2 |
| 115 | + if k.lower() in ("data", "analytic", "reflectance", "visual", "b04", "b08"): |
| 116 | + s += 1 |
| 117 | + return s |
| 118 | + |
| 119 | + if not assets: |
| 120 | + raise ValueError("STAC item has no assets.") |
| 121 | + key, a = max(assets.items(), key=lambda kv: score(*kv)) |
| 122 | + href = a.get("href") |
| 123 | + if not href: |
| 124 | + raise ValueError("Chosen asset has no href.") |
| 125 | + return key, href |
| 126 | + |
| 127 | + |
| 128 | +def stream_stac_quicklook( |
| 129 | + stac_api: str, |
| 130 | + collection: str, |
| 131 | + datetime_range: str, |
| 132 | + bbox: Optional[List[float]] = None, |
| 133 | + cloud_lt: Optional[float] = 20, |
| 134 | + asset_key: Optional[str] = None, |
| 135 | + overview_level: int = 4, |
| 136 | +) -> Tuple[np.ndarray, Dict[str, Any], plt.Figure]: |
| 137 | + item = search_stac_first_item_resilient(stac_api, [collection], datetime_range, bbox=bbox, cloud_lt=cloud_lt) |
| 138 | + |
| 139 | + item_self = next((l.get("href") for l in item.get("links", []) if l.get("rel") == "self"), None) |
| 140 | + if item_self is None: |
| 141 | + item_self = f"{stac_api.rstrip('/')}/collections/{item.get('collection')}/items/{item.get('id')}" |
| 142 | + r = requests.get(item_self, timeout=90) |
| 143 | + _ensure_ok(r) |
| 144 | + item_full = r.json() |
| 145 | + |
| 146 | + key, href = choose_stac_asset(item_full.get("assets", {}), asset_key=asset_key) |
| 147 | + vsi_href = f"/vsicurl/{href}" if href.startswith("http") else href |
| 148 | + |
| 149 | + gdal_env = { |
| 150 | + "AWS_NO_SIGN_REQUEST": os.environ.get("AWS_NO_SIGN_REQUEST", "YES"), |
| 151 | + "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": os.environ.get("CPL_VSIL_CURL_ALLOWED_EXTENSIONS", "tif,tiff") |
| 152 | + } |
| 153 | + with rasterio.Env(**gdal_env): |
| 154 | + with rasterio.open(vsi_href) as ds: |
| 155 | + if ds.overviews(1): |
| 156 | + olist = ds.overviews(1) |
| 157 | + idx = max(0, min(overview_level - 1, len(olist) - 1)) |
| 158 | + out_h = max(1, ds.height // olist[idx]) |
| 159 | + out_w = max(1, ds.width // olist[idx]) |
| 160 | + else: |
| 161 | + scale = 8 if max(ds.width, ds.height) > 2048 else 2 |
| 162 | + out_h = max(1, ds.height // scale) |
| 163 | + out_w = max(1, ds.width // scale) |
| 164 | + arr = ds.read(1, out_shape=(out_h, out_w)) |
| 165 | + meta = ds.meta.copy() |
| 166 | + meta.update({"chosen_asset_key": key, "href": href, "vsi_href": vsi_href, "item_id": item.get("id")}) |
| 167 | + |
| 168 | + fig = plt.figure(figsize=(6, 5)) |
| 169 | + plt.imshow(arr, interpolation="nearest") |
| 170 | + plt.axis("off") |
| 171 | + plt.title(f"{collection} • {key} • {item.get('id')}") |
| 172 | + plt.tight_layout() |
| 173 | + return arr, meta, fig |
| 174 | + |
| 175 | +``` |
| 176 | + |
| 177 | +To fetch a quicklook over Rocky Mountain National Park: |
| 178 | + |
| 179 | +```python |
| 180 | +bbox_rmnp = [-105.9, 40.1, -105.3, 40.6] # west, south, east, north |
| 181 | + |
| 182 | +arr, meta, fig = stream_stac_quicklook( |
| 183 | + stac_api="https://earth-search.aws.element84.com/v1", |
| 184 | + collection="sentinel-2-l2a", |
| 185 | + datetime_range="2024-08-01T00:00:00Z/2024-08-07T23:59:59Z", |
| 186 | + bbox=bbox_rmnp, |
| 187 | + cloud_lt=20, |
| 188 | + asset_key="visual" # or None |
| 189 | +) |
| 190 | +fig |
| 191 | +``` |
0 commit comments