From 574060d769a94af5c2740c1ec18ba11951e87d95 Mon Sep 17 00:00:00 2001 From: Jin Igarashi Date: Sun, 26 Jul 2026 11:58:11 +0900 Subject: [PATCH 1/5] feat: implement NTL in assess command --- rapida/cli/assess.py | 45 +++- rapida/components/ntl/__init__.py | 341 +++++++++++++++++++++++++++- rapida/ntl/outage.py | 2 + rapida/util/bbox_param_type.py | 22 ++ rapida/util/download_remote_file.py | 14 ++ 5 files changed, 415 insertions(+), 9 deletions(-) diff --git a/rapida/cli/assess.py b/rapida/cli/assess.py index 5216883..622115f 100644 --- a/rapida/cli/assess.py +++ b/rapida/cli/assess.py @@ -91,6 +91,21 @@ def validate_variables(ctx, param, value): return value +def validate_popvars(ctx, param, value): + """ + click callback function to validate --popvar values against population component variables. + """ + if not value: + return value + valid_vars = get_variables_by_components(['population']).get('population', []) + invalid = [v for v in value if v not in valid_vars] + if invalid: + raise click.BadParameter( + f"Invalid popvar{'s' if len(invalid) > 1 else ''}: {', '.join(invalid)}. " + f"Valid options: {', '.join(valid_vars)}") + return value + + def validate_datetime_range(ctx, param, value): """ click callback function to validate --datetime value @@ -172,6 +187,23 @@ def build_variable_help(): @click.option("--outage-date", "outage_date", type=click.DateTime(formats=["%Y-%m-%d"]), required=False, help='The human experience of a specific night, local time zone matched to the center of bbox') +@click.option('--popvar', 'popvar', required=False, multiple=True, + type=str, callback=validate_popvars, + help="Optional. One or more population variables to compute affected-population zonal stats " + "over detected outages (ntl component only). If omitted, the population count-variables " + "already present in the project are auto-selected.") + +@click.option('-ot', '--percentage-drop', 'percentage_drop', type=int, default=50, show_default=True, + help="Outage threshold: minimum %% radiance drop to flag a pixel as an outage (ntl component only).") + +@click.option('-cm', '--cmask', 'mask_clouds', is_flag=True, default=False, + help="Enable cloud masking during NTL outage detection (ntl component only). Off by default to " + "avoid over-masking (NASA A2/water-vapor is frequently flagged cloudy).") + +@click.option('--bbox-buffer', 'bbox_buffer', type=float, default=0.0, show_default=True, + help="Optional. Enlarge the project bbox by this many meters before NTL outage detection " + "(NTL has ~500m pixels; helps small AOIs). 0 = no buffer (ntl component only).") + @click.option('--cloud-cover', '-cc', required=False, type=int, multiple=False, default=5, show_default=True,help=f"Optional. Minimum cloud cover rate to search items for landuse component.") @click.option('-p', '--project', @@ -182,7 +214,7 @@ def build_variable_help(): help=f'Force assess components. Downloaded data or computed data will be ignored and recomputed.') @click.pass_context -def assess(ctx, all=False, components=None, variables=None, year=None, datetime_range=None, outage_date=None, cloud_cover=None, project: str = None, force=False): +def assess(ctx, all=False, components=None, variables=None, year=None, datetime_range=None, outage_date=None, popvar=None, percentage_drop=50, mask_clouds=False, bbox_buffer=0.0, cloud_cover=None, project: str = None, force=False): """ Assess/evaluate a specific geospatial exposure components/variables @@ -210,6 +242,13 @@ def assess(ctx, all=False, components=None, variables=None, year=None, datetime rapida assess -c landuse -dt 2025-02-01/2025-05-31 -cc 10: Search Sentinel 2 item which is less than 10% of cloud cover from February to May 2025. + rapida assess -c ntl -v noaa_outage --outage-date 2026-07-20 -ot 40 --bbox-buffer 2000: NTL outage detection for a specific night. + + NTL notes: the 'ntl' component requires --outage-date. NOAA (noaa_outage) is usually more reliable than the NASA + archive (nasa_outage), whose bottom-of-atmosphere product is frequently flagged as cloudy. For small areas of + interest, enlarge the search window with --bbox-buffer (NTL pixels are ~500m) and tune the outage threshold with + -ot/--percentage-drop. To inspect input data quality for the same bbox and date, use `rapida ntl search`. + """ progress = ctx.obj.get('progress') if not is_rapida_initialized(): @@ -267,6 +306,10 @@ def assess(ctx, all=False, components=None, variables=None, year=None, datetime year=year, datetime_range=datetime_range, outage_date=outage_date, + pop_vars=popvar, + percentage_drop=percentage_drop, + mask_clouds=mask_clouds, + bbox_buffer=bbox_buffer, cloud_cover=cloud_cover, force=force) diff --git a/rapida/components/ntl/__init__.py b/rapida/components/ntl/__init__.py index 29224ad..c406f67 100644 --- a/rapida/components/ntl/__init__.py +++ b/rapida/components/ntl/__init__.py @@ -2,12 +2,64 @@ from rapida.core.variable import Variable from rapida.project.project import Project from rapida.session import Session +import asyncio +import datetime +import glob import logging import os +import re +import click logger = logging.getLogger('rapida') class NtlComponent(Component): + """ + Nighttime Lights (NTL) power-outage detection component for `rapida assess -c ntl`. + + Detects power outages by comparing a target night's VIIRS nighttime-lights radiance + against a recent monthly baseline, then summarizes the result per admin polygon in a + ``stats.ntl`` layer of the project geopackage. + + Variables (data streams), pass with ``-v``: + - ``noaa_outage`` : NOAA real-time VIIRS data. Usually the most reliable choice. + - ``nasa_nrt_outage`` : NASA Black Marble operational / near-real-time (LANCE). Last ~7 days only. + - ``nasa_outage`` : NASA Black Marble archived (LAADS). Bottom-of-atmosphere; frequently + flagged as cloudy, so results can be sparse. + If ``-v`` is omitted, all three variables are assessed. + + Options (defined on the assess command): + - ``--outage-date YYYY-MM-DD`` : REQUIRED. The target night to assess. + - ``--popvar NAME`` : Population variable(s) whose affected population is summed over + outages. If omitted, the population count-variables already + present in the project are auto-selected. + - ``-ot / --percentage-drop N``: Radiance-drop threshold (%) to flag a pixel as an outage. Default 50. + - ``-cm / --cmask`` : Enable cloud masking. Off by default to avoid over-masking. + - ``--bbox-buffer METERS`` : Enlarge the project bbox before detection (NTL pixels are ~500m, + which helps small areas of interest). Default 0 (no buffer). + + Credentials (read from the environment / .env): + - NASA streams require ``EARTHDATA_TOKEN``. + - NOAA stream requires ``SPACETRACK_USER`` and ``SPACETRACK_PASSWORD``. + + Outputs: + - ``/data/ntl/_.tif`` : multi-band outage raster. + - ``stats.ntl`` layer in the project geopackage with, per admin polygon: + ``_outage_pixels`` (outage cell count), ``_outage_pct`` (% of area in outage), + and ``__affected`` (population inside outage areas), reusing the population + rasters already in the project. + + Examples: + rapida assess -c ntl -v noaa_outage --outage-date 2026-07-20 + Detect outages for one night using NOAA data. + + rapida assess -c ntl -v noaa_outage --outage-date 2026-07-20 -ot 40 --bbox-buffer 2000 + Loosen the threshold and widen the search window for a small area of interest. + + rapida assess -c ntl -v nasa_nrt_outage --outage-date 2026-07-20 --popvar total + Also compute the affected total population per admin polygon. + + Tip: to inspect input data quality for the same bbox and date, use `rapida ntl search`. + """ def __call__(self, variables: list[str], **kwargs): if not variables: @@ -44,17 +96,290 @@ def __init__(self, **kwargs): output_filename = f"{self.name}.tif" self.local_path = os.path.join(os.path.dirname(geopackage_path), self.component, output_filename) - def download(self,force=False, **kwargs): - pass - def download(self, **kwargs): pass + def resolve(self, **kwargs): pass - def compute(self, **kwargs): - pass - def evaluate(self, **kwargs): - pass + def _check_credentials(self, deliverable): + """ + Fail fast with a clear message if the credentials required by the target + data stream are not set, instead of raising deep inside the pipeline. + """ + missing = [] + if 'NASA' in deliverable and not os.environ.get('EARTHDATA_TOKEN'): + missing.append('EARTHDATA_TOKEN') + if 'NOAA' in deliverable: + for var_name in ('SPACETRACK_USER', 'SPACETRACK_PASSWORD'): + if not os.environ.get(var_name): + missing.append(var_name) + if missing: + raise click.UsageError( + f"Missing credentials for {deliverable}: {', '.join(missing)}. " + f"Set them in your .env / environment." + ) + + def _run_async(self, coro): + """ + Run an async coroutine from the synchronous assess flow. Mirrors the + AsyncCommand pattern in rapida/cli/aclick.py: nest_asyncio is applied + globally at rapida.cli import time, so run_until_complete tolerates being + called from an already-async context. + """ + try: + loop = asyncio.get_event_loop() + except RuntimeError: + loop = asyncio.new_event_loop() + asyncio.set_event_loop(loop) + return loop.run_until_complete(coro) + + def compute(self, force=False, outage_date=None, progress=None, year=None, pop_vars=None, + percentage_drop=50, mask_clouds=False, bbox_buffer=0.0, **kwargs): + # Imported here (not at module top) to avoid a circular import: + # rapida.ntl.outage imports rapida.cli.assess. + from rapida.ntl.outage import detect_outage + from rapida.util.bbox_param_type import buffer_bbox + + if outage_date is None: + raise click.UsageError("The ntl component requires --outage-date (YYYY-MM-DD).") + + deliverable = self.name.upper() + self._check_credentials(deliverable) + + project = Project(path=os.getcwd()) + bbox = project.geobounds + if bbox_buffer: + bbox = buffer_bbox(bbox, bbox_buffer) + + dst_dir = os.path.dirname(self.local_path) + os.makedirs(dst_dir, exist_ok=True) + + # NOTE: pop_vars is intentionally NOT passed to detect_outage. detect_outage's own + # pop_vars path builds a throwaway Project from the outage polygons and re-downloads + # population into a separate gpkg. In the assess flow we instead reuse the population + # data already in the project and write a stats.ntl layer (see self.evaluate). + coro = detect_outage( + bbox=bbox, + nominal_date=outage_date, + deliverable=deliverable, + dst_dir=dst_dir, + mask_clouds=mask_clouds, + percentage_drop=percentage_drop, + pop_vars=None, + display=False, + progress=progress, + year=year or datetime.datetime.now().year, + ) + outage_tif = self._run_async(coro) + + if outage_tif and os.path.exists(outage_tif): + self.evaluate(outage_tif=outage_tif, pop_vars=pop_vars or None, year=year, progress=progress) + else: + logger.info(f'No outage output produced for {deliverable}; nothing to evaluate.') + return outage_tif + + def evaluate(self, outage_tif=None, pop_vars=None, year=None, progress=None, **kwargs): + """ + Compute per-admin-polygon NTL statistics from the outage raster and write them + into the project geopackage as a stats. (stats.ntl) layer. + + Statistics per polygon: + - _outage_pixels : count of outage cells + - _outage_pct : percentage of polygon area flagged as outage + - __affected : population inside outage areas (population raster x outage), + reusing the population rasters already in the project. + """ + import geopandas as gpd + from rapida.util import geo + from rapida.stats.raster_zonal_stats import zst + + project = Project(path=os.getcwd()) + ntl_dir = os.path.dirname(self.local_path) + + # When no pop_vars were requested, auto-select the population count-variables that + # already have data in this project, so affected population is computed out of the box. + if not pop_vars: + pop_vars = self._discover_population_vars(project, year) + if pop_vars: + logger.info( + f'Auto-selected population variables for affected-population stats: {", ".join(pop_vars)}' + ) + else: + logger.info( + 'No population variables found in the project; writing outage extent only. ' + 'Run `rapida assess -c population` first to include affected population.' + ) + + # 1. Collapse the OUTAGE band(s) into a single 0/1 mask in the source CRS (EPSG:4326). + outage_4326 = os.path.join(ntl_dir, f'{self.name}_outage_mask_4326.tif') + self._extract_outage_mask(outage_tif, outage_4326) + + # 2. Reproject/crop/align the outage mask to the project CRS for zonal stats. + outage_proj = os.path.join(ntl_dir, f'{self.name}_outage_mask.tif') + geo.import_raster( + source=outage_4326, dst=outage_proj, target_srs=project.target_srs, + crop_ds=project.geopackage_file_path, crop_layer_name=project.polygons_layer_name, + progress=progress, + ) + + # 3. Outage extent per polygon (sum = cell count, mean = fraction -> percentage). + src_rasters = [outage_proj, outage_proj] + vars_ops = [(f'{self.name}_outage_pixels', 'sum'), (f'{self.name}_outage_fraction', 'mean')] + + # 4. Affected population per requested pop_var, reusing project population rasters. + for pv in (pop_vars or ()): + pop_raster = self._find_population_raster(project, pv, year) + if not pop_raster: + logger.warning( + f'Population raster for "{pv}" not found under ' + f'{os.path.join(project.data_folder, "population", pv)}. ' + f'Run `rapida assess -c population -v {pv}` first. ' + f'Skipping affected-population for "{pv}".' + ) + continue + affected = self._compute_affected_population(pop_raster, outage_4326, ntl_dir, pv) + src_rasters.append(affected) + vars_ops.append((f'{self.name}_{pv}_affected', 'sum')) + + # 5. Zonal stats against the project polygons; merge into stats.ntl if it already exists. + dst_layer = f'stats.{self.component}' + lnames = gpd.list_layers(project.geopackage_file_path).name.tolist() + polygons_layer = dst_layer if dst_layer in lnames else project.polygons_layer_name + gdf = zst( + src_rasters=src_rasters, polygon_ds=project.geopackage_file_path, + polygon_layer=polygons_layer, vars_ops=vars_ops, progress=progress, + ) + + frac_col = f'{self.name}_outage_fraction' + if frac_col in gdf.columns: + gdf[f'{self.name}_outage_pct'] = gdf[frac_col] * 100.0 + gdf.drop(columns=[frac_col], inplace=True) + + self._write_stats_layer(project, gdf, dst_layer) + logger.info(f'Wrote NTL statistics to {project.geopackage_file_path}:{dst_layer}') + + def _extract_outage_mask(self, outage_tif, dst_4326): + """Collapse every band whose description contains 'OUTAGE' into a single 0/1 GeoTIFF.""" + import numpy as np + from osgeo import gdal + + ds = gdal.Open(outage_tif) + try: + gt = ds.GetGeoTransform() + srs = ds.GetSpatialRef() + xsize, ysize = ds.RasterXSize, ds.RasterYSize + combined = None + for i in range(1, ds.RasterCount + 1): + band = ds.GetRasterBand(i) + desc = band.GetDescription() or '' + if 'OUTAGE' in desc.upper(): + arr = band.ReadAsArray() + mask = np.nan_to_num(arr, nan=0.0) > 0 + combined = mask if combined is None else (combined | mask) + finally: + ds = None + + if combined is None: + raise click.UsageError(f'No OUTAGE band found in {outage_tif}') + + driver = gdal.GetDriverByName('GTiff') + dst = driver.Create(dst_4326, xsize, ysize, 1, gdal.GDT_Byte) + dst.SetGeoTransform(gt) + if srs is not None: + dst.SetSpatialRef(srs) + dst.GetRasterBand(1).WriteArray(combined.astype('uint8')) + dst.FlushCache() + dst = None + return dst_4326 + + def _discover_population_vars(self, project, year): + """ + Auto-discover population count-variables (operator 'sum') that already have a raster + in this project. Ratio variables (e.g. dependency) are excluded because summing + population x outage is only meaningful for additive counts. + """ + discovered = [] + try: + with Session() as session: + pop_defs = session.get_component('population') + except Exception as e: + logger.debug(f'Could not read population component definitions: {e}') + return discovered + for name, var_data in pop_defs.items(): + if var_data.get('operator') != 'sum': + continue + if self._find_population_raster(project, name, year): + discovered.append(name) + return discovered + + def _find_population_raster(self, project, pop_var, year): + """Locate the population raster already produced in the project for pop_var.""" + base = os.path.join(project.data_folder, 'population', pop_var) + if not os.path.isdir(base): + return None + if year: + candidate = os.path.join(base, f'{pop_var}_{year}.tif') + if os.path.exists(candidate): + return candidate + matches = [ + p for p in sorted(glob.glob(os.path.join(base, f'{pop_var}_*.tif'))) + if re.fullmatch(rf'{re.escape(pop_var)}_\d+\.tif', os.path.basename(p)) + ] + if matches: + return matches[-1] + fallback = os.path.join(base, f'{pop_var}.tif') + return fallback if os.path.exists(fallback) else None + + def _compute_affected_population(self, pop_raster, outage_4326, ntl_dir, pop_var): + """population x outage, aligned to the population grid, so a zonal sum yields affected people.""" + from osgeo import gdal + from osgeo_utils.gdal_calc import Calc + + pds = gdal.Open(pop_raster) + try: + gt = pds.GetGeoTransform() + width, height = pds.RasterXSize, pds.RasterYSize + psrs = pds.GetSpatialRef() + finally: + pds = None + + minx, maxy = gt[0], gt[3] + maxx = minx + gt[1] * width + miny = maxy + gt[5] * height + + outage_aligned = os.path.join(ntl_dir, f'{self.name}_{pop_var}_outage_on_popgrid.tif') + gdal.Warp( + outage_aligned, outage_4326, format='GTiff', dstSRS=psrs, + outputBounds=(minx, miny, maxx, maxy), width=width, height=height, + resampleAlg='near', dstNodata=0, + ) + + affected = os.path.join(ntl_dir, f'{self.name}_{pop_var}_affected.tif') + ds = Calc( + calc='a*b', a=pop_raster, b=outage_aligned, outfile=affected, + projectionCheck=True, format='GTiff', quiet=True, overwrite=True, + ) + ds = None + return affected + + def _write_stats_layer(self, project, gdf, dst_layer): + """Overwrite (or create) the stats layer in the project geopackage from a GeoDataFrame.""" + import io + from osgeo import gdal + from pyogrio import write_dataframe + + with io.BytesIO() as bio: + fpath = f'/vsimem/{dst_layer}.fgb' + write_dataframe(df=gdf, path=bio, layer=dst_layer, driver='FlatGeobuf') + gdal.FileFromMemBuffer(fpath, bio.getbuffer()) + bio.seek(0) + with gdal.OpenEx(fpath) as src: + options = gdal.VectorTranslateOptions( + format='GPKG', accessMode='overwrite', layerName=dst_layer, makeValid=True, + ) + gdal.VectorTranslate(destNameOrDestDS=project.geopackage_file_path, srcDS=src, options=options) + gdal.Unlink(fpath) + def __call__(self, *args, **kwargs): - print(self.name) \ No newline at end of file + return self.compute(**kwargs) diff --git a/rapida/ntl/outage.py b/rapida/ntl/outage.py index 35b54a1..2917d9c 100644 --- a/rapida/ntl/outage.py +++ b/rapida/ntl/outage.py @@ -263,3 +263,5 @@ async def detect_outage( from rapida.ntl import vis vis.display2(data=arrays, title=f'Outage inputs and results for {deliverable} at {bbox} on {nominal_date.date()}') + + return outage_tif_path diff --git a/rapida/util/bbox_param_type.py b/rapida/util/bbox_param_type.py index 7adb98f..cb9cb2a 100644 --- a/rapida/util/bbox_param_type.py +++ b/rapida/util/bbox_param_type.py @@ -1,10 +1,32 @@ import logging +import math import reverse_geocoder as rg import click logger = logging.getLogger(__name__) +def buffer_bbox(bbox: tuple[float, float, float, float], meters: float) -> tuple[float, float, float, float]: + """ + Enlarge a geographic (WGS84 lon/lat) bbox by a distance in meters on all sides. + + Meters are converted to degrees using an approximation at the bbox center latitude. + Useful for coarse rasters (e.g. NTL ~500m pixels) where a small AOI would otherwise + contain too few pixels. + + :param bbox: (min_lon, min_lat, max_lon, max_lat) in EPSG:4326 + :param meters: buffer distance in meters; 0/None returns the bbox unchanged + :return: the enlarged bbox + """ + if not meters: + return bbox + minlon, minlat, maxlon, maxlat = bbox + center_lat = (minlat + maxlat) / 2.0 + dlat = meters / 111320.0 + dlon = meters / (111320.0 * max(math.cos(math.radians(center_lat)), 1e-6)) + return (minlon - dlon, minlat - dlat, maxlon + dlon, maxlat + dlat) + + class BboxParamType(click.ParamType): name = "bbox" def convert(self, value, param, ctx): diff --git a/rapida/util/download_remote_file.py b/rapida/util/download_remote_file.py index b26155b..6f097fc 100644 --- a/rapida/util/download_remote_file.py +++ b/rapida/util/download_remote_file.py @@ -105,6 +105,20 @@ async def download_file(file_url=None, dst_file_path=None, remote_size_str = response.headers.get('Content-Length') remote_size = int(remote_size_str) if remote_size_str else None + # Auth/authorization failures (e.g. NASA Earthdata / LANCE NRT) are + # served as an HTML login page with HTTP 200, which would otherwise be + # saved with the data extension and later crash cryptically in GDAL. + # Detect it here so we fail clearly and never cache the bad file. + content_type = response.headers.get('Content-Type', '') + if 'text/html' in content_type.lower(): + raise Exception( + f'Authentication/authorization failed for {file_url}: server ' + f'returned an HTML page (Content-Type: {content_type}) instead of ' + f'the data file. For NASA NRT (LANCE) data, ensure the account behind ' + f'EARTHDATA_TOKEN has authorized the required application at ' + f'https://urs.earthdata.nasa.gov (Applications -> Authorized Apps).' + ) + if os.path.exists(dst_file_path): # Only compare sizes if remote_size is known if not force and remote_size is not None and os.path.getsize(dst_file_path) == remote_size: From 8e3ffffc361c328f8a2c07b9b30232132a100b6c Mon Sep 17 00:00:00 2001 From: Jin Igarashi Date: Sun, 26 Jul 2026 13:11:04 +0900 Subject: [PATCH 2/5] fix: exclude ntl from --all option --- rapida/cli/assess.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/rapida/cli/assess.py b/rapida/cli/assess.py index 622115f..5f1f8b7 100644 --- a/rapida/cli/assess.py +++ b/rapida/cli/assess.py @@ -169,7 +169,7 @@ def build_variable_help(): @click.command(short_help='assess/evaluate a specific geospatial exposure components/variables', no_args_is_help=True) @click.option( '--all', '-a', is_flag=True, default=False, - help="compute all components and variables if this option is set" + help="compute all auto-runnable components and variables (excludes landuse and ntl) if this option is set" ) @click.option( '--components', '-c', required=False, multiple=True, @@ -230,7 +230,7 @@ def assess(ctx, all=False, components=None, variables=None, year=None, datetime Usage: - rapida assess --all: assess all components + rapida assess --all: assess all auto-runnable components (excludes landuse and ntl) rapida assess -c rwi: assess RWI component only. @@ -276,7 +276,7 @@ def assess(ctx, all=False, components=None, variables=None, year=None, datetime target_components = components if len(components) == 0: if all: - target_components = set(filter(lambda x: x != "landuse", all_components)) + target_components = set(filter(lambda x: x not in {"landuse", "ntl"}, all_components)) else: logger.warning(f"At least one component is required. If you want to assess all components, use --all option") return From 5654cadcfdeed2083bca9c85a35a91f3098e5437 Mon Sep 17 00:00:00 2001 From: Jin Igarashi Date: Sun, 26 Jul 2026 20:27:49 +0900 Subject: [PATCH 3/5] fix: export matplotlib file --- rapida/components/ntl/__init__.py | 5 +++++ rapida/ntl/outage.py | 12 ++++++++---- rapida/ntl/vis.py | 22 ++++++++++++++++++++-- 3 files changed, 33 insertions(+), 6 deletions(-) diff --git a/rapida/components/ntl/__init__.py b/rapida/components/ntl/__init__.py index c406f67..3e2ea27 100644 --- a/rapida/components/ntl/__init__.py +++ b/rapida/components/ntl/__init__.py @@ -155,6 +155,10 @@ def compute(self, force=False, outage_date=None, progress=None, year=None, pop_v dst_dir = os.path.dirname(self.local_path) os.makedirs(dst_dir, exist_ok=True) + # Save the outage analysis visualization (same 6-panel plot as + # `ntl detect --display`) into the project's ntl folder. + display_path = os.path.join(dst_dir, f"{self.name}_outage.png") + # NOTE: pop_vars is intentionally NOT passed to detect_outage. detect_outage's own # pop_vars path builds a throwaway Project from the outage polygons and re-downloads # population into a separate gpkg. In the assess flow we instead reuse the population @@ -168,6 +172,7 @@ def compute(self, force=False, outage_date=None, progress=None, year=None, pop_v percentage_drop=percentage_drop, pop_vars=None, display=False, + display_path=display_path, progress=progress, year=year or datetime.datetime.now().year, ) diff --git a/rapida/ntl/outage.py b/rapida/ntl/outage.py index 2917d9c..06d2975 100644 --- a/rapida/ntl/outage.py +++ b/rapida/ntl/outage.py @@ -28,7 +28,7 @@ async def detect_outage( bbox: tuple[numbers.Number] = None, nominal_date: datetime = None, deliverable: str = None, dst_dir: str = None, mask_clouds:bool = True, percentage_drop:int = None, pop_vars:str|tuple[str]=None, - display: bool = False, progress: Progress = None, year=datetime.datetime.now().year): + display: bool = False, display_path: str = None, progress: Progress = None, year=datetime.datetime.now().year): logger.info(f'Fetching best imagery for {deliverable} {bbox}-{nominal_date} ') # with open(os.path.join('/tmp', 'bbox.geojson'), "w") as ff: @@ -259,9 +259,13 @@ async def detect_outage( finally: # Unlink (delete) the virtual file to free up system memory gdal.Unlink(vsimem_path) - if display: + if display or display_path: from rapida.ntl import vis - vis.display2(data=arrays, - title=f'Outage inputs and results for {deliverable} at {bbox} on {nominal_date.date()}') + vis.display2( + data=arrays, + title=f'Outage inputs and results for {deliverable} at {bbox} on {nominal_date.date()}', + save_path=display_path, # None for interactive-only detect + show=display, # keep plt.show() only when --display asked + ) return outage_tif_path diff --git a/rapida/ntl/vis.py b/rapida/ntl/vis.py index e40e70b..8f50978 100644 --- a/rapida/ntl/vis.py +++ b/rapida/ntl/vis.py @@ -114,11 +114,17 @@ def plot(array): plt.show() -def display2(data=dict(), interpolation='nearest', title='', max_discrete_vals=5): +def display2(data=dict(), interpolation='nearest', title='', max_discrete_vals=5, + save_path=None, show=True): """ Improved display function that maximizes screen real estate, ensures perfectly aligned subplots, and automatically detects discrete classification maps to build custom legends. + + @args + @save_path - if set, the figure is written to this path (PNG) instead of + only being shown interactively. Parent directory is created. + @show - if True, the figure is shown interactively via plt.show(). """ n = len(data) if n == 0: return @@ -190,4 +196,16 @@ def display2(data=dict(), interpolation='nearest', title='', max_discrete_vals=5 axes_flat[j].axis('off') fig.suptitle(title, fontsize=20) - plt.show() \ No newline at end of file + + if save_path: + import os + os.makedirs(os.path.dirname(save_path), exist_ok=True) + fig.savefig(save_path, dpi=150, bbox_inches='tight') + + if show: + plt.show() + + if save_path: + # Free the figure once it has been persisted to avoid leaking figures + # when many variables are rendered in a single run. + plt.close(fig) \ No newline at end of file From 95fd06e324642adcd427bb8269b613d92f471dac Mon Sep 17 00:00:00 2001 From: Jin Igarashi Date: Sun, 26 Jul 2026 21:22:44 +0900 Subject: [PATCH 4/5] fix: pass cloudmask param to fetch instead of hardcoding to pass True always --- rapida/ntl/fetch.py | 5 +++-- rapida/ntl/noaa/search.py | 6 +++++- rapida/ntl/outage.py | 2 +- 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/rapida/ntl/fetch.py b/rapida/ntl/fetch.py index 8e06818..b90bc9d 100644 --- a/rapida/ntl/fetch.py +++ b/rapida/ntl/fetch.py @@ -31,7 +31,8 @@ async def download_and_track(granule, dest_dir, prog_bar): return granule.timestamp, result_dict async def fetch(bbox:tuple[numbers.Number]=None, nominal_date:datetime=None, - deliverable:str=None, dst_dir:str=None, progress:Progress=None): + deliverable:str=None, dst_dir:str=None, progress:Progress=None, + mask_clouds:bool=False): """ Indentify and download the BEST available data suitable to detect outages. :param dst_dir: @@ -48,7 +49,7 @@ async def fetch(bbox:tuple[numbers.Number]=None, nominal_date:datetime=None, logger.info(f'Going to predict VIIRS satellite passes for {nominal_date.date()} over target area: {bbox}') granules = await async_search_granules( satellites=None, nominal_date=nominal_date, bbox=bbox, - cmask=True, progress=progress) + cmask=mask_clouds, progress=progress) if not granules: logger.info(f'No descending granules for found for {nominal_date.date()} over target area {bbox}') return diff --git a/rapida/ntl/noaa/search.py b/rapida/ntl/noaa/search.py index e13eb6e..6257e4d 100644 --- a/rapida/ntl/noaa/search.py +++ b/rapida/ntl/noaa/search.py @@ -50,6 +50,7 @@ class Granule: elevation:float cloud_cover = None pint:float = None + url:str = None @property def id(self): return f"{self.start_time:%Y%m%d%H%M%S}{self.start_time.microsecond // 100000}" # @@ -427,9 +428,12 @@ async def night_granules_async(self, bbox:Iterable[float]=None, nominal_date:dat current_granule.start_time = start_time logger.debug(f'Replacing granule {old_timestamp} with {current_granule.timestamp}') + # Always resolve the public URL so downstream consumers (e.g. + # select_required_granules) can read granule.url regardless of cmask. + url = public_url(file_path=file_path, satellite=self.satellite, source=source) + current_granule.url = url if cmask: # Use the unique URL as the key (Always unique) - url = public_url(file_path=file_path, satellite=self.satellite, source=source) selected_granules[url] = current_granule else: # Use the unique file_path as the key to prevent SNPP/N20/N21 overwrites diff --git a/rapida/ntl/outage.py b/rapida/ntl/outage.py index 06d2975..63e5fed 100644 --- a/rapida/ntl/outage.py +++ b/rapida/ntl/outage.py @@ -36,7 +36,7 @@ async def detect_outage( #fetch daily data, source independent # --- 2. FETCH DAILY TARGET DATA --- daily_results = await fetch(bbox=bbox, nominal_date=nominal_date, deliverable=deliverable, - progress=progress, dst_dir=dst_dir) + progress=progress, dst_dir=dst_dir, mask_clouds=mask_clouds) if not daily_results: logger.info(f'No imagery was found for {nominal_date:"%Y%m%d"} over {bbox} {deliverable.split("_")[0]}') logger.info(f'Consider adjusting source, date or the bounding box') From 2396ee416d61a7688a0404e17d03e087cb8ec2d8 Mon Sep 17 00:00:00 2001 From: Jin Igarashi Date: Sun, 26 Jul 2026 21:49:42 +0900 Subject: [PATCH 5/5] Snap the crop bounds outward to the source tile pixel grid --- rapida/ntl/nasa/io.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/rapida/ntl/nasa/io.py b/rapida/ntl/nasa/io.py index c3cd68d..9264e47 100644 --- a/rapida/ntl/nasa/io.py +++ b/rapida/ntl/nasa/io.py @@ -7,6 +7,7 @@ import h5py import secrets import asyncio +import math import os from rich.progress import Progress from rapida.ntl import cache @@ -134,6 +135,23 @@ def get_val(key): gdal.FileFromMemBuffer(tile_vrt_path, tile_xml.encode('utf-8')) tile_vrts.append(tile_vrt_path) + # Snap the crop bounds outward to the source tile pixel grid. gdal.BuildVRT with + # a fractional-pixel offset between the mosaic origin and the tile grid raises + # "Wrong values in SrcRect" when mosaicking adjacent tiles. Bounds that are + # already grid-aligned (e.g. integer degrees from `ntl detect -b`) are a no-op; + # sub-pixel reprojection noise from project.geobounds gets aligned to whole + # pixels. All tiles share the same global grid, so the last tile's origin + # (west/north) and pixel size (px_w/px_h) define the grid phase. + if bbox is not None: + minx, miny, maxx, maxy = bbox + ph = abs(px_h) + eps = 1e-6 # tolerance so already-aligned edges don't bump a whole pixel + minx = west + math.floor((minx - west) / px_w + eps) * px_w + maxx = west + math.ceil((maxx - west) / px_w - eps) * px_w + maxy = north - math.floor((north - maxy) / ph + eps) * ph + miny = north - math.ceil((north - miny) / ph - eps) * ph + bbox = (minx, miny, maxx, maxy) + # 3. Mosaic all virtual tiles into the final master VRT vrt_opts = gdal.BuildVRTOptions( outputBounds=bbox,