Spaces:
Paused
Paused
Download satquery_engine/services/radiometry.py from SM737/satquery-api: direct link, hf CLI and curl.
- Browser
- Download file 4.05 kB
-
https://huggingface.co/spaces/SM737/satquery-api/resolve/main/satquery_engine/services/radiometry.py
- Command line
-
hf download hf://spaces/SM737/satquery-api/satquery_engine/services/radiometry.py
-
curl -L -o radiometry.py https://huggingface.co/spaces/SM737/satquery-api/resolve/main/satquery_engine/services/radiometry.py
4.05 kB
| """One explicit RGB scientific-to-model contract, shared by every tile.""" | |
| import numpy as np | |
| from scipy import ndimage | |
| from rasterio.enums import Resampling | |
| from satquery_engine.services.bands import detect_band_map | |
| def rgb_indexes(src): | |
| mapping = detect_band_map(src) | |
| if all(k in mapping.indices for k in ("red","green","blue")): | |
| return [mapping.indices[k] for k in ("red","green","blue")] | |
| if src.count == 3 and not mapping.warnings and not any(src.descriptions): | |
| return [1,2,3] # documented ordinary RGB upload convention | |
| # FLAIR-style VHR chips have a fixed five-band layout (RGB, NIR, DSM), | |
| # uint8 storage, and 20 cm georeferencing but commonly omit descriptions. | |
| # Keep this deliberately narrow so arbitrary multispectral rasters are not | |
| # silently routed into RGB models. | |
| if ( | |
| src.count == 5 | |
| and all(dtype == "uint8" for dtype in src.dtypes) | |
| and not any(src.descriptions) | |
| and not mapping.warnings | |
| and 0.15 <= max(abs(src.res[0]), abs(src.res[1])) <= 0.25 | |
| ): | |
| return [1,2,3] | |
| raise ValueError("RGB analysis requires unambiguous red, green and blue bands.") | |
| def rgb_unit_data(src, *, window=None, out_shape=None, scene_valid=None): | |
| indexes = rgb_indexes(src) | |
| shape = (3,*out_shape) if out_shape else None | |
| masked = src.read(indexes, window=window, out_shape=shape, masked=True, resampling=Resampling.nearest).astype("float32") | |
| values = masked.filled(np.nan) | |
| good = np.all(np.isfinite(values) & ~np.ma.getmaskarray(masked),axis=0) | |
| scales = np.asarray([src.scales[i-1] for i in indexes],dtype="float32")[:,None,None] | |
| offsets = np.asarray([src.offsets[i-1] for i in indexes],dtype="float32")[:,None,None] | |
| values = values*scales+offsets | |
| byte = all(src.dtypes[i-1]=="uint8" for i in indexes) and np.all(scales==1) and np.all(offsets==0) | |
| # Registered rasters retain the explicit original encoding in metadata. | |
| encoding = src.tags().get("satquery_rgb_encoding") | |
| if byte or encoding == "uint8": | |
| values /= 255.0 | |
| method = "uint8 / 255" | |
| else: | |
| method = "reflectance using metadata scales and offsets" | |
| if good.any() and (np.min(values[:,good])<0 or np.max(values[:,good])>1.00001): | |
| raise ValueError("RGB radiometry is unsupported: use 8-bit RGB or metadata-scaled reflectance in 0-1. No per-tile scaling was guessed.") | |
| # Exact black, border-connected padding is NoData evidence, not water or | |
| # shadow. Interior dark surfaces and nonzero dark water remain valid. | |
| exact_black = np.all(values == 0, axis=0) & good | |
| if window is not None: | |
| # Tile edges are not image edges: an interior black roof or pond must | |
| # retain the same validity in every overlapping inference window. | |
| if scene_valid is None: | |
| _, scene_valid, _ = rgb_unit_data(src) | |
| if scene_valid.shape != (src.height, src.width): | |
| raise ValueError("RGB scene validity must match the native source grid.") | |
| rows = np.floor(window.row_off + (np.arange(good.shape[0]) + .5) * window.height / good.shape[0]).astype(int) | |
| cols = np.floor(window.col_off + (np.arange(good.shape[1]) + .5) * window.width / good.shape[1]).astype(int) | |
| good &= scene_valid[np.ix_(rows, cols)] | |
| elif exact_black.any(): | |
| labels, count = ndimage.label(exact_black) | |
| border_ids = set(labels[0]) | set(labels[-1]) | set(labels[:, 0]) | set(labels[:, -1]) | |
| border_ids.discard(0) | |
| if border_ids: | |
| padding = np.isin(labels, list(border_ids)) | |
| good[padding] = False | |
| values[:,~good] = 0 | |
| selection = "explicit band metadata" | |
| if not all(k in detect_band_map(src).indices for k in ("red", "green", "blue")): | |
| selection = "ordinary three-band RGB convention" if src.count == 3 else "FLAIR-style VHR RGB/NIR/DSM contract" | |
| return values.astype("float32"), good, {"indexes":indexes,"selection":selection,"method":method,"scales":scales.ravel().tolist(),"offsets":offsets.ravel().tolist()} | |