conn bar and no barr (#510)

* conn bar and no barr

* move valid vars in connectivity

Author: iferencik

rapida/cli/assess.py

@@ -147,7 +147,7 @@ def build_variable_help():
         if vars_list:
             vars_str = ", ".join(vars_list)
             parts.append(f"{comp} ({vars_str})")
-        if comp == 'population': print(len(vars_list))
+
     return "The variable/s to be assessed. Will be filtered by selected components. Available variables per component:\n\n" + "\n\n".join(parts)
 
 
@@ -211,7 +211,7 @@ 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.
 
     """
-
+    progress = ctx.obj.get('progress')
     if not is_rapida_initialized():
         return
 
@@ -231,7 +231,7 @@ def assess(ctx, all=False, components=None,  variables=None, year=None, datetime
         sys.exit(0)
 
     logger.info(f'Current project/folder: {prj.path}')
-    with Progress(disable=False, console=None) as progress:
+    with progress:
         with Session() as session:
             all_components = session.get_components()
             target_components = components

rapida/cli/connectivity.py

@@ -1,4 +1,4 @@
-from typing import Union
+from typing import Union, Iterable
 
 import click
 import logging
@@ -7,10 +7,21 @@
 from rapida.connectivity import run_connectivity_analysis
 from rapida.cli.aclick import AsyncCommand
 from rapida.connectivity.isochrone import MODE_MAP
-
+from rapida.cli.assess import get_variables_by_components
 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
+
 def parse_intervals(ctx, param, value):
     """Parses a comma-separated string of numbers into a list of integers."""
     if not value:
@@ -46,6 +57,20 @@ def parse_intervals(ctx, param, value):
     help="Comma-separated time intervals in minutes for the catchment areas."
 )
 
+
+@click.option(
+    '-sd', '--sites-dataset',
+    type=click.Path(exists=True, file_okay=True, dir_okay=False, readable=True),
+    default=None,
+    help="Path to an OGR-supported vector data source (e.g., GPKG, Shapefile) containing sites."
+)
+@click.option(
+    '-sl', '--sites-layer',
+    type=str,
+    default="0",
+    help="Name or index of the layer to use from sites dataset. Defaults to the first layer (layer 0)."
+)
+
 @click.option(
     '-bd', '--barriers-dataset',
     type=click.Path(exists=True, file_okay=True, dir_okay=False, readable=True),
@@ -59,12 +84,19 @@ def parse_intervals(ctx, param, value):
     help="Name or index of the layer to use from barriers dataset. Defaults to the first layer (layer 0)."
 )
 
+
 @click.option('-bb', "--barriers-buffer",
     type=int,
     default=5,
     required=False,
     help="The value in meters to used to buffer the geometries in barriers/dataset/layer in case the barriers are lines"
 )
+
+@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"
+              )
+
 @click.option(
     "--dst-dir",
     "-d",           # Short option
@@ -84,12 +116,15 @@ def parse_intervals(ctx, param, value):
 @click.pass_context
 async def connectivity(ctx, bbox:tuple[float, float, float, float]=None, travel_mode:str=None,
                        time_intervals:list[int] =None, dst_dir:str=None,
-                       barriers_dataset:str=None, barriers_layer:str=None, barriers_buffer:int=None
+                       barriers_dataset:str=None, barriers_layer:str=None, barriers_buffer:int=None,
+                       sites_dataset:str=None, sites_layer:str=None,popvar:str|tuple[str]=None
     ):
-    logger.info(f'Running connectivity analysis ')
+    logger.info(f'Running connectivity analysis')
     progress = ctx.obj.get('progress')
     with progress:
         return await run_connectivity_analysis(
             bbox=bbox, dst_dir=dst_dir, travel_mode=travel_mode, time_intervals=time_intervals,
-            barriers_dataset=barriers_dataset, barriers_layer=barriers_layer, barriers_buffer=barriers_buffer, progress=progress
+            barriers_dataset=barriers_dataset, barriers_layer=barriers_layer, barriers_buffer=barriers_buffer,
+            sites_dataset=sites_dataset, sites_layer=sites_layer, pop_vars=popvar,
+            progress=progress
         )

rapida/connectivity/__init__.py

@@ -2,27 +2,66 @@
 import os.path
 from rapida.util.bbox_param_type import get_best_semantic_label
 from rich.progress import Progress
-from rapida.connectivity.io import prepare_osm_pbf,extract_health_sites, extract_origins_from_geojson
+from rapida.connectivity.io import prepare_osm_pbf,extract_health_sites, extract_origins_from_geojson, extract_origins
 from rapida.connectivity.graph import compile_valhalla_graph
 from rapida.connectivity.isochrone import connectivity_areas
+# from rapida.cli.assess import assess
+# import click
+# from rapida.project.project import Project
+# from tempfile import TemporaryDirectory
 
 
 
 async def run_connectivity_analysis(
         bbox:tuple[float, float, float, float]=None, travel_mode:str=None, time_intervals:list[int] =None,
-        dst_dir:str=None, barriers_dataset:str=None, barriers_layer:str=None, barriers_buffer:int=None, progress:Progress=None
+        dst_dir:str=None, barriers_dataset:str=None, barriers_layer:str=None, barriers_buffer:int=None,
+        sites_dataset:str=None, sites_layer:str=None,pop_vars:str|tuple[str]=None,
+        progress:Progress=None
     ):
     bbox_label = get_best_semantic_label(bbox=bbox)
     dest_dir = os.path.join(dst_dir, bbox_label)
     bbox_pbf = await prepare_osm_pbf(bbox=bbox, dst_dir=dest_dir, progress=progress)
-    health_sites = await extract_health_sites(pbf_path=bbox_pbf, dst_dir=dest_dir, progress=progress)
+    if sites_dataset is None:
+        sites = await extract_health_sites(pbf_path=bbox_pbf, dst_dir=dest_dir, progress=progress)
+    else:
+        sites = sites_dataset
+
     dag_tar_path = await compile_valhalla_graph(pbf_path=bbox_pbf,dst_dir=dest_dir, progress=progress)
-    origins = extract_origins_from_geojson(geojson_path=health_sites)
+    origins = extract_origins(sites_dataset=sites, src_layer=sites_layer)
+
+
     results = await connectivity_areas(
+        tar_path=dag_tar_path, origins=origins, travel_mode=travel_mode, intervals_minutes=time_intervals)
+    isochrones_path = os.path.join(dest_dir,  'isochrones.geojson')
+    with open(isochrones_path, "w") as f:
+        json.dump(results, f, indent=2)
+
+    # with TemporaryDirectory(dir=dest_dir, delete=False) as project_folder:
+    #     project = Project(path=project_folder, polygons=isochrones_path, comment='temp project for conn isochrones')
+    #
+    #
+    #     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.
+    #         await ctx.invoke(
+    #             assess,
+    #             components=('population',),
+    #             variables=pop_vars,
+    #             year=2026,
+    #             project=project.path,
+    #             force=False
+    #         )
+
+
+    if barriers_dataset is not None:
+        barrier_results = await connectivity_areas(
         tar_path=dag_tar_path, origins=origins, travel_mode=travel_mode, intervals_minutes=time_intervals,
         barriers_dataset=barriers_dataset, barriers_layer=barriers_layer, barriers_buffer=barriers_buffer
-                                       )
-    with open(os.path.join(dest_dir, 'isochrones.geojson'), "w") as f:
-        json.dump(results, f, indent=2)
+                             )
+        with open(os.path.join(dest_dir, 'isochrones_with_barriers.geojson'), "w") as f:
+            json.dump(barrier_results, f, indent=2)
+
 
     return

rapida/connectivity/io.py

@@ -11,7 +11,7 @@
 from rich.progress import Progress
 import json
 from shapely.geometry import shape, mapping
-from osgeo import gdal
+from osgeo import gdal, ogr, osr
 from shapely.wkb import loads as load_wkb
 from shapely.ops import orient, transform
 import numpy as np
@@ -58,6 +58,7 @@ async def prepare_osm_pbf(bbox: tuple[float, float, float, float], dst_dir: str
 
         # Perform quick spatial overlay check
         gdf = gpd.GeoDataFrame.from_features(features, crs="EPSG:4326")
+
         intersecting = gdf[gdf.intersects(bbox_geom)]
 
         if intersecting.empty:
@@ -214,6 +215,67 @@ def extract_origins_from_geojson(geojson_path: str) -> list[tuple[float, float]]
     return origins
 
 
+def extract_origins(sites_dataset: str=None, src_layer: str = None) -> list[tuple[float, float]]:
+    """
+    Extracts a list of (longitude, latitude) tuples from a spatial file using OGR.
+    Handles reprojection to WGS84 (EPSG:4326) if the source is not in lat/lon.
+    """
+    # Open the dataset
+    with gdal.OpenEx(sites_dataset, gdal.OF_VECTOR) as src_ds:
+
+        try:
+            layer = src_ds.GetLayer(int(src_layer))
+        except ValueError:
+            layer = src_ds.GetLayerByName(str(src_layer))
+
+        if layer is None:
+            raise ValueError(f"Layer '{src_layer}' could not be found in the dataset {sites_dataset}.")
+
+
+
+        # Set up coordinate transformation to WGS84 (EPSG:4326)
+        source_srs = layer.GetSpatialRef()
+        target_srs = osr.SpatialReference()
+        target_srs.ImportFromEPSG(4326)
+
+        # OGR 3+ strict axis mapping strategy (ensures Longitude/Latitude order)
+        target_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+
+        transform = None
+        if source_srs and not source_srs.IsSame(target_srs):
+            transform = osr.CoordinateTransformation(source_srs, target_srs)
+
+        origins = []
+        layer.ResetReading()
+        # Iterate through features
+        for feature in layer:
+            # If you still need the filter: if feature.GetField("osm_id") != 80: continue
+
+            geom = feature.GetGeometryRef()
+            if geom is not None:
+                # Clone geometry to avoid modifying original layer data during transform
+                geom_clone = geom.Clone()
+
+                # Reproject if necessary
+                if transform:
+                    geom_clone.Transform(transform)
+
+                # Helper function to recursively extract points from nested collections
+                def extract_points(g):
+                    name = g.GetGeometryName()
+                    if name == "POINT":
+                        origins.append((g.GetX(), g.GetY()))
+                    elif name in ("MULTIPOINT", "GEOMETRYCOLLECTION"):
+                        for i in range(g.GetGeometryCount()):
+                            sub_geom = g.GetGeometryRef(i)
+                            if sub_geom is not None:
+                                extract_points(sub_geom)
+
+                extract_points(geom_clone)
+
+        return origins
+
+
 def read_barriers_grid(src_path: str, src_layer: str = None, barriers_buffer:float=None) -> list:
     """Reads a vector source and cuts features into micro-tiles to stay under Valhalla's limit."""
     if not src_path:

rapida/project/project.py

@@ -16,7 +16,7 @@
 from geopandas import GeoDataFrame
 from osgeo import gdal, ogr, osr
 from azure.storage.fileshare import ShareClient
-
+import reverse_geocoder as rg
 from rapida import constants
 from rapida.az.blobstorage import check_blob_exists, delete_blob
 from rapida.session import Session
@@ -159,12 +159,17 @@ def __init__(self, path: str,polygons: str = None,
                     iso3_codes = []
                     for lat, lon in zip(lats, lons):
                         try:
-                            iso3_codes.append(fetch_ccode(lat=lat, lon=lon))
+                            result = rg.search((lat, lon))[0]
+                            iso2_cc = result.get('cc', '')
+                            country = coco.convert(names=iso2_cc, to='ISO3')
+                            iso3_codes.append(country)
                         except Exception as e:
                             logger.warning(f"Failed to fetch ISO3 for point ({lat}, {lon}): {e}")
                             iso3_codes.append(None)
                     gdf["iso3"] = iso3_codes
+
                     self.countries = tuple(sorted(set(filter(lambda x: x in COUNTRY_CODES, gdf["iso3"]))))
+
                     gdf.to_file(
                         filename=self.geopackage_file_path,
                         driver="GPKG",

rapida/stats/raster_zonal_stats.py

@@ -257,4 +257,5 @@ def progress_callback(completed, message, progress=progress, task=task):
             if vname in egdf.columns.tolist():
                 egdf.drop(columns=[vname], inplace=True)
         combined = egdf.merge(combined, on='geometry', how='inner')
+        combined = combined.sort_values(by='geometry', key=lambda geom: geom.to_geo_index().area, ascending=False)
         return combined