implement pop zonal stats fo outage areas (#512)

Author: iferencik

rapida/cli/connectivity.py

@@ -94,7 +94,7 @@ def parse_intervals(ctx, param, value):
 
 @click.option('--popvar', required=False, multiple=True,
               type=str, callback=validate_variables,
-              help=f"Open or more RAPIDA population variable to compute zonal stats for withing the connectivity zones"
+              help=f"One or more RAPIDA population variable to compute zonal stats for withing the connectivity zones"
               )
 
 @click.option(

rapida/cli/ntl.py

@@ -4,7 +4,7 @@
 from typing import Iterable
 import click
 import tempfile
-
+from rapida.cli.assess import get_variables_by_components
 from rapida.cli import RapidaCommandGroup
 from rapida.ntl.nasa.const import ARCHIVE, OPERATIONAL, PROCESSING_LEVEL_NAMES, PRODUCT_NAMES, NTL_FILENAME_PATTERN, ROUTES, COLLECTIONS
 from rapida.ntl.nasa.search import search as nasa_search
@@ -16,10 +16,25 @@
 from rich.table import Table
 from rapida.ntl.fetch import DELIVERABLES, fetch as fetch_ntl
 from rapida.ntl.outage import detect_outage
+from rapida.cli.assess import assess
+from rapida.project.project import Project
+from tempfile import TemporaryDirectory
+
 
 logger = logging.getLogger(__name__)
 
 
+
+def validate_variables(ctx, param, value):
+    """
+    click callback function to validate  polulation
+    """
+    valid_vars = get_variables_by_components(['population'])['population']
+    invalid = [v for v in value if v not in valid_vars]
+    if invalid:
+        raise click.BadParameter(f"Invalid variable{'s' if len(invalid) > 1 else ''}: {', '.join(invalid)} for population . Valid options: {', '.join(valid_vars)}")
+    return value
+
 class ProcessingLevelChoiceOption(click.Option):
     """
     Custom Click option that dynamically validates choices based on the value
@@ -499,6 +514,11 @@ async def fetch(ctx, bbox:tuple[float, float, float, float]=None, nominal_date:d
     )
 
 
+@click.option('--popvar', required=False, multiple=True,
+              type=str, callback=validate_variables,
+              help=f"One or more RAPIDA population variable to compute zonal stats for outages"
+              )
+
 @click.option('-ot', "--percentage_drop",
     type=int,
     default=50,
@@ -528,12 +548,13 @@ async def fetch(ctx, bbox:tuple[float, float, float, float]=None, nominal_date:d
 
 @click.pass_context
 async def detect(ctx, bbox:tuple[float, float, float, float]=None, nominal_date:datetime=None, deliverable:str=None,
-                 mask_clouds:bool=True, dst_dir:str=None, percentage_drop:int=None, display:bool=False):
+                 mask_clouds:bool=True, dst_dir:str=None, popvar:str|tuple[str]=None, percentage_drop:int=None, display:bool=False):
     progress = ctx.obj.get('progress')
     with progress:
         return await detect_outage(
             bbox=bbox, nominal_date=nominal_date, deliverable=deliverable, dst_dir=dst_dir,
-            progress=progress, mask_clouds=mask_clouds, percentage_drop=percentage_drop, display=display
+            progress=progress, mask_clouds=mask_clouds, percentage_drop=percentage_drop, display=display,
+            pop_vars=popvar
         )
 
 

rapida/ntl/fetch.py

@@ -150,7 +150,6 @@ async def fetch(bbox:tuple[numbers.Number]=None, nominal_date:datetime=None,
             # Break out of the 'routes' loop to stop searching for this processing_level
             break
 
-    print(downloaded_files)
     return downloaded_files
 
 

rapida/ntl/nasa/io.py

@@ -189,14 +189,15 @@ def extract_bb(image_files: list[str] = None, sds_name:str=None, return_gt=False
     ds = gdal.Translate(destName='', srcDS=vrt_path, **translate_options)
 
     gt = ds.GetGeoTransform()
+    srs = ds.GetSpatialRef()
     [gdal.Unlink(e) for e in to_unlink]
     band = ds.GetRasterBand(1)
     array = band.ReadAsArray()
 
 
     ds = None
     if return_gt:
-        return array, gt
+        return array, gt, srs
     return array
 
 async def extract(image_files: list[str] = None, sds_name:str=None, product:str=None, dst_dir:str=None,

rapida/ntl/noaa/cmask.py

@@ -101,7 +101,7 @@ def select_required_granules(sorted_granules: list, bbox: tuple, progress:Progre
                     return combo  # Exits immediately with the absolute minimum required set
 
         logger.warning("Exhausted all combinations. BBOX cannot be fully covered by available data.")
-        logger.warning(f"Selecting max of {combos}.")
+
 
         max_cov = max(combos)
         return combos[max_cov]

rapida/ntl/outage.py

@@ -12,15 +12,23 @@
 import numpy as np
 import os
 from rapida.ntl import utils
-from rapida.admin.util import bbox_to_geojson_polygon
+from rapida.util.bbox_param_type import get_best_semantic_label
+from osgeo import gdal, gdal_array, ogr
+from rapida.util.geo import polygonize_raster_mask
+from rapida.cli.assess import assess
+import click
+from rapida.project.project import Project
+from tempfile import TemporaryDirectory
+import geopandas as gpd
+import datetime
 logger = logging.getLogger('rapida')
-
+gdal.UseExceptions()
 
 
 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,
-        display: bool = False, progress: Progress = 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):
 
     logger.info(f'Fetching best imagery for {deliverable} {bbox}-{nominal_date} ')
     # with open(os.path.join('/tmp', 'bbox.geojson'), "w") as ff:
@@ -44,7 +52,7 @@ async def detect_outage(
     monthly_timestamp, monthly_files = next(iter(monthly_results.items()))
     monthly_product, monthly_image_files = next(iter(monthly_files.items()))
 
-    monthly_data, gt = extract_bb(image_files=monthly_image_files, sds_name=nasa_const.SUB_DATASETS['A3'],
+    monthly_data, gt, rsrs = extract_bb(image_files=monthly_image_files, sds_name=nasa_const.SUB_DATASETS['A3'],
                                     bbox=bbox, progress=progress, return_gt=True)
 
     monthly_data_label = f'{monthly_product}_{monthly_timestamp}'
@@ -165,11 +173,92 @@ async def detect_outage(
                 arrays[f'{daily_data_label}_ZSCORE'] = zscore
                 arrays[f'{daily_data_label}_OUTAGE'] = outage
 
-
-    file_name = utils.get_custom_bbox_label(bbox)
+    # --- 5. UNIFIED DISPLAY & EXPORT ---
+    #file_name = utils.get_custom_bbox_label(bbox)
+    file_name = get_best_semantic_label(bbox=bbox)
     outage_tif_path = os.path.join(dst_dir, f'{deliverable}_{file_name}.tif')
+    outage_gpkg_path = os.path.join(dst_dir, f'{deliverable}_{file_name}.gpkg')
     write_outage_tif(src_arrays=arrays, gt=gt, dst_path=outage_tif_path)
-    # --- 5. UNIFIED DISPLAY & EXPORT ---
+
+    if pop_vars:
+        # vectorize outage
+        gtiff_driver = gdal.GetDriverByName("GTiff")
+
+        # Ensure parent directory exists
+        os.makedirs(os.path.dirname(os.path.abspath(outage_gpkg_path)), exist_ok=True)
+
+        if not os.path.exists(outage_gpkg_path):
+            # 1. Ask GDAL to sniff out the right driver name based on the path structure/extension
+            v_driver = ogr.GetDriverByName('GPKG')
+
+            # 2. Use the driver to initialize the empty file header
+            with v_driver.CreateDataSource(outage_gpkg_path) as init_ds:
+                pass  # Successfully generated empty stub container
+        for k, arr in arrays.items():
+            if 'outage' in k.lower():
+                # Ensure data type is compatible (e.g., forcing int64/bool to int32/uint8 if needed)
+                if arr.dtype == bool:
+                    arr = arr.astype("uint8")
+                elif arr.dtype == "int64":
+                    arr = arr.astype("int32")  # GDAL standard doesn't always love 64-bit ints
+                gdal_dt = gdal_array.NumericTypeCodeToGDALTypeCode(arr.dtype)
+                rows, cols = arr.shape
+                # Generate a unique virtual memory path
+                vsimem_path = f"/vsimem/{k}.tif"
+
+
+                try:
+                    # Create, write, and safely close the dataset within the context manager
+                    with gtiff_driver.Create(vsimem_path, cols, rows, 1, gdal_dt) as ds:
+                        ds.SetGeoTransform(gt)
+                        ds.SetSpatialRef(rsrs)
+                        band = ds.GetRasterBand(1)
+                        band.WriteArray(arr)
+                        band.FlushCache()
+
+                    # The dataset is now closed, but the file exists in /vsimem/
+                    polygonize_raster_mask(
+                        raster_ds=vsimem_path,  # Passes the string path
+                        dst_dataset=outage_gpkg_path,
+                        dst_layer=k,
+                        geom_type=ogr.wkbPolygon
+                    )
+
+                    with TemporaryDirectory(dir=dst_dir, delete=True) as project_folder:
+                        project = Project(path=project_folder, polygons=outage_gpkg_path,
+                                          comment='temp project for outage zonal stats')
+                        with click.Context(assess) as ctx:
+                            ctx.ensure_object(dict)
+                            ctx.obj['progress'] = progress
+                            # 2. Use invoke. Do NOT pass 'ctx' manually here.
+                            # Click intercepts this and injects it as the first argument automatically.
+                            ctx.invoke(
+                                assess,
+                                components=('population',),
+                                variables=pop_vars,
+                                year=year,
+                                project=project.path,
+                                force=False
+                            )
+                            stat_gpkg_path = os.path.join(project_folder, 'data', f'{project.name}.gpkg')
+                            pop_stat_gdf = gpd.read_file(stat_gpkg_path, layer='stats.population')
+
+                            pop_stat_gdf = pop_stat_gdf.to_crs('EPSG:4326')
+
+                            pop_stat_gdf.to_file(
+                                filename=outage_gpkg_path,
+                                driver="gpkg",
+                                engine="pyogrio",
+                                mode="w",
+                                layer=k,
+                                promote_to_multi=True,
+                                index=False
+                            )
+
+
+                finally:
+                    # Unlink (delete) the virtual file to free up system memory
+                    gdal.Unlink(vsimem_path)
     if display:
         from rapida.ntl import vis
         vis.display2(data=arrays,

rapida/util/geo.py

@@ -146,15 +146,23 @@ def polygonize_raster_mask(raster_ds=None, band=1, dst_dataset=None, dst_layer=N
 
     with gdal.OpenEx(raster_ds, gdal.OF_RASTER|gdal.OF_READONLY) as rds:
         with gdal.OpenEx(dst_dataset, gdal.OF_VECTOR | gdal.OF_UPDATE) as vds:
-            logger.info(f'Polygonizing {raster_ds} ')
+            logger.info(f'Polygonizing {raster_ds} into {dst_dataset}:{dst_layer} ')
             mband = rds.GetRasterBand(band)
-            mask_lyr = vds.CreateLayer(dst_layer, geom_type=geom_type, srs=rds.GetSpatialRef())
+            layer_index = vds.GetLayerByName(dst_layer)
+            rsrs = rds.GetSpatialRef()
+            if layer_index is not None:
+                vds.DeleteLayer(dst_layer)
+            mask_lyr = vds.CreateLayer(dst_layer, geom_type=geom_type, srs=rsrs)
 
             r = gdal.Polygonize(srcBand=mband, maskBand=mband,outLayer=mask_lyr,iPixValField=-1)
 
             assert r == 0, f'Failed to polygonize {raster_ds}'
+
             for feature in mask_lyr:
+                if rsrs is None:continue
+                if rsrs.IsGeographic():continue
                 geom = feature.GetGeometryRef()
+
                 simplified_geom = geom.Simplify(constants.DEFAULT_MASK_POLYGONIZATION_SMOOTHING_BUFFER)  # Use SimplifyPreserveTopology(tolerance) if needed
                 smoothed_geom = simplified_geom.Buffer(constants.DEFAULT_MASK_POLYGONIZATION_SMOOTHING_BUFFER).Buffer(-constants.DEFAULT_MASK_POLYGONIZATION_SMOOTHING_BUFFER)
                 feature.SetGeometry(smoothed_geom)