implement pop zonal stats and disjoint isochrones

rapida/cli/__init__.py

@@ -16,12 +16,13 @@
 from rapida.cli.connectivity import connectivity
 from rich.progress import Progress
 import click
-# import nest_asyncio
-# nest_asyncio.apply()
-import uvloop
-import asyncio
-
-asyncio.set_event_loop_policy(uvloop.EventLoopPolicy())
+import nest_asyncio
+nest_asyncio.apply()
+#import uvloop
+# import asyncio
+#
+# asyncio.set_event_loop_policy(uvloop.EventLoopPolicy())
+#
 
 
 

rapida/cli/aclick.py

@@ -13,6 +13,8 @@ def setup_debug_logging(ctx, param, value):
             handler.setLevel(logging.DEBUG)
         logger.debug("Debug logging enabled globally.")
     return value
+
+
 class AsyncCommand(click.Command):
     """
     Async wrapper designed to work alongside nest_asyncio in Jupyter
@@ -29,16 +31,20 @@ def wrapped_callback(*c_args, **c_kwargs):
                 actual_func = inspect.unwrap(orig_callback)
 
                 if inspect.iscoroutinefunction(actual_func):
-                    # Safely get or create the event loop
+                    is_new_loop = False
                     try:
                         loop = asyncio.get_event_loop()
                     except RuntimeError:
                         loop = asyncio.new_event_loop()
                         asyncio.set_event_loop(loop)
+                        is_new_loop = True
 
-                    return loop.run_until_complete(orig_callback(*c_args, **c_kwargs))
-
-                return orig_callback(*c_args, **c_kwargs)
+                    try:
+                        return loop.run_until_complete(orig_callback(*c_args, **c_kwargs))
+                    finally:
+                        if is_new_loop:
+                            loop.close()
+                            asyncio.set_event_loop(None)
 
             self.callback = wrapped_callback
 class RapidaCommandGroup(click.Group):
@@ -51,22 +57,19 @@ def list_commands(self, ctx):
         return self.commands.keys()
 
     def add_command(self, cmd: click.Command, name: str = None) -> None:
-        # Catch-all: forces help display if a subcommand is run empty
-        has_required_params = any(param.required for param in cmd.params)
-        is_grp = isinstance(cmd, click.Group) or issubclass(cmd.__class__, self.__class__)
-        #print(cmd.name, has_required_params, is_grp, cmd.params)
-        # if not is_grp and has_required_params:
-        #     cmd.no_args_is_help = True
-        # 2. Automatically inject --debug into every registered subcommand
-        cmd.params.append(
-            click.Option(
-                ['--debug'],
-                is_flag=True,
-                help='Enable debug logging.',
-                expose_value=False,
-                callback=setup_debug_logging
+        # Check if the --debug option is already injected
+        has_debug = any(param.name == 'debug' for param in cmd.params)
+
+        if not has_debug:
+            cmd.params.append(
+                click.Option(
+                    ['--debug'],
+                    is_flag=True,
+                    help='Enable debug logging.',
+                    expose_value=False,
+                    callback=setup_debug_logging
+                )
             )
-        )
         super().add_command(cmd, name)
 
     def command(self, *args, **kwargs):

rapida/cli/connectivity.py

@@ -36,9 +36,9 @@ def parse_intervals(ctx, param, value):
 @click.command(short_help='run connectivity analysis', cls=AsyncCommand)
 
 @click.option('-b', '--bbox',
-              required=True,
+              required=False,
               type=BboxParamType(),
-              help='Bounding box xmin/west, ymin/south, xmax/east, ymax/north'
+              help='Bounding box xmin/west, ymin/south, xmax/east, ymax/north. If not supplied bbox is derived from the site dataset'
               )
 @click.option(
         '-m', "--mode", "travel_mode",
@@ -113,11 +113,22 @@ def parse_intervals(ctx, param, value):
 )
 
 
+@click.option(
+    '--disjoint',
+    is_flag=True,
+    help=(
+            "By default Valhalla routing engine create isochrones as overlapping polygons."
+            "Use this flag to separate the isochrones areas or make them disjoint"
+    ),
+    default=False
+)
+
 @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,
-                       sites_dataset:str=None, sites_layer:str=None,popvar:str|tuple[str]=None
+                       sites_dataset:str=None, sites_layer:str=None,popvar:str|tuple[str]=None,
+                       disjoint:bool=False
     ):
     logger.info(f'Running connectivity analysis')
     progress = ctx.obj.get('progress')
@@ -126,5 +137,5 @@ async def connectivity(ctx, bbox:tuple[float, float, float, float]=None, travel_
             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,
             sites_dataset=sites_dataset, sites_layer=sites_layer, pop_vars=popvar,
-            progress=progress
+            progress=progress, disjoint=disjoint
         )

rapida/connectivity/__init__.py

@@ -1,23 +1,36 @@
+import datetime
 import json
 import os.path
 from rapida.util.bbox_param_type import get_best_semantic_label
+import geopandas as gpd
+import logging
 from rich.progress import Progress
 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
-
+from rapida.cli.assess import assess
+import click
+from rapida.project.project import Project
+from tempfile import TemporaryDirectory
 
+logger = logging.getLogger(__name__)
 
 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,
         sites_dataset:str=None, sites_layer:str=None,pop_vars:str|tuple[str]=None,
-        progress:Progress=None
+        progress:Progress=None, year=datetime.datetime.now().year, disjoint:bool=False
     ):
+    if bbox is None:
+        assert sites_dataset is not None, f'site_dataset has to be provided when bbox is not'
+        try:
+            slayer = int(sites_layer)
+        except ValueError:
+            slayer = sites_layer
+        gdf = gpd.read_file(sites_dataset, layers=slayer)
+        if not gdf.crs.is_geographic:
+            gdf.to_crs('EPSG:4326', inplace=True)
+        bbox = gdf.total_bounds
     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)
@@ -31,37 +44,95 @@ async def run_connectivity_analysis(
 
 
     results = await connectivity_areas(
-        tar_path=dag_tar_path, origins=origins, travel_mode=travel_mode, intervals_minutes=time_intervals)
+        tar_path=dag_tar_path, origins=origins, travel_mode=travel_mode, intervals_minutes=time_intervals, disjoint=disjoint)
+
+
     isochrones_path = os.path.join(dest_dir,  'isochrones.geojson')
     with open(isochrones_path, "w") as f:
         json.dump(results, f, indent=2)
+    if pop_vars:
+        with TemporaryDirectory(dir=dest_dir, delete=True) 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.
+                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')
+                if not disjoint:
+                    pop_stat_gdf = pop_stat_gdf.iloc[pop_stat_gdf.geometry.area.sort_values(ascending=False).index]
+                pop_stat_gdf = pop_stat_gdf.to_crs('EPSG:4326')
 
-    # 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
-    #         )
-
+                pop_stat_gdf.to_file(
+                    filename=isochrones_path,
+                    driver="GeoJSON",
+                    engine="pyogrio",
+                    mode="w",
+                    layer='isochrones',
+                    promote_to_multi=True,
+                    index=False
+                )
 
     if barriers_dataset is not None:
+        logger.info(f'Computing isochrones with barriers')
         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
+        barriers_dataset=barriers_dataset, barriers_layer=barriers_layer, barriers_buffer=barriers_buffer, disjoint=disjoint
                              )
-        with open(os.path.join(dest_dir, 'isochrones_with_barriers.geojson'), "w") as f:
+        barrier_isochrones_path = os.path.join(dest_dir, 'isochrones_with_barriers.geojson')
+        with open(barrier_isochrones_path, "w") as f:
             json.dump(barrier_results, f, indent=2)
+        if pop_vars:
+            logger.info(f'Computing zonal stats for barrier isochrones')
+            with TemporaryDirectory(dir=dest_dir, delete=True) as project_folder:
+                project = Project(path=project_folder, polygons=barrier_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.
+                    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')
+                    barrier_pop_stat_gdf = gpd.read_file(stat_gpkg_path, layer='stats.population')
+                    if not disjoint:
+                        barrier_pop_stat_gdf = barrier_pop_stat_gdf.iloc[pop_stat_gdf.geometry.area.sort_values(ascending=False).index]
+                    barrier_pop_stat_gdf = barrier_pop_stat_gdf.to_crs('EPSG:4326')
+
+                    pop_col_names = [f'{popv}_{year}' for popv in pop_vars]
+                    new_pop_col_names = [f'{popv}_{year}_barrier' for popv in pop_vars]
+                    col_name_dict = dict(zip(pop_col_names, new_pop_col_names))
+                    barrier_pop_stat_gdf.rename(columns=col_name_dict, inplace=True)
+
+                    data_cols = ['contour']+pop_col_names
+                    data_frame = pop_stat_gdf[data_cols]
+                    barrier_pop_stat_gdf = barrier_pop_stat_gdf.merge(data_frame, on='contour', how='left')
+                    for pvar, bar_pvar in col_name_dict.items():
+                        barrier_pop_stat_gdf[f'{pvar}_{bar_pvar}_difference'] = barrier_pop_stat_gdf[pvar] - barrier_pop_stat_gdf[bar_pvar]
 
+                    barrier_pop_stat_gdf.to_file(
+                        filename=barrier_isochrones_path,
+                        driver="GeoJSON",
+                        engine="pyogrio",
+                        mode="w",
+                        layer='barrier_isochrones',
+                        promote_to_multi=True,
+                        index=False
+                    )
 
     return

rapida/connectivity/isochrone.py

@@ -1,6 +1,8 @@
 import json
 import asyncio
 from pathlib import Path
+
+import geopandas
 from valhalla import Actor
 from shapely.geometry import shape, mapping, JOIN_STYLE
 from pyproj import Transformer
@@ -21,6 +23,30 @@
 
 
 
+def make_isochrones_disjoint(gdf, time_col="time", group_col=None):
+    """Converts overlapping concentric isochrones into mutually exclusive rings."""
+    # Ensure data is sorted by travel time (smallest/inner first)
+    sort_cols = [group_col, time_col] if group_col else [time_col]
+    gdf = gdf.sort_values(sort_cols).reset_index(drop=True)
+
+    # Copy to store the result
+    result_gdf = gdf.copy()
+
+    # Iterate backwards to subtract the smaller inner polygon from the larger outer one
+    for i in range(len(gdf) - 1, 0, -1):
+        # If tracking multiple starting points, ensure they belong to the same group
+        if group_col and gdf.loc[i, group_col] != gdf.loc[i - 1, group_col]:
+            continue
+
+        # Subtract the inner geometry from the outer geometry
+        result_gdf.loc[i, "geometry"] = gdf.loc[i, "geometry"].difference(
+            gdf.loc[i - 1, "geometry"]
+        )
+
+    return result_gdf
+
+
+
 
 async def connectivity_areas(
         tar_path: str,
@@ -30,6 +56,7 @@ async def connectivity_areas(
         barriers_dataset:str=None,
         barriers_layer:str=None,
         barriers_buffer:int=None,
+        disjoint:bool=False,
         progress=None
 ) -> dict:
     tar_file = Path(tar_path)
@@ -95,7 +122,7 @@ def run_routing():
             isochrone_geojson = json.loads(response_str)
 
             # 3. Intercept Valhalla's output and apply Shapely smoothing
-            for feature in isochrone_geojson.get("features", []):
+            for fid, feature in enumerate(isochrone_geojson.get("features", []), start=1):
                 raw_geom_wgs84 = shape(feature["geometry"])
 
                 # 1. Get the bounding box of the raw WGS84 polygon
@@ -143,15 +170,16 @@ def run_routing():
 
                 # 6. Convert back to WGS84 degrees so the GeoJSON renders on a map properly
                 final_geom_wgs84 = transform(project_to_degrees, smooth_geom_meters)
-
+                feature["id"] = fid
                 feature["geometry"] = mapping(final_geom_wgs84)
 
-                #feature["geometry"] = mapping(raw_geom_wgs84)
+
 
                 feature["properties"].update({
                     "mode": travel_mode,
                     "type": "system_catchment",
-                    "facility_count": len(locations)
+                    "facility_count": len(locations),
+                    "id": fid
                 })
                 results["features"].append(feature)
 
@@ -168,4 +196,12 @@ def run_routing():
 
         return results
 
-    return await asyncio.to_thread(run_routing)
+    isos = await asyncio.to_thread(run_routing)
+
+    if disjoint:
+
+        # Example Usage:
+        gdf = geopandas.GeoDataFrame.from_features(isos, crs='EPSG:4326')
+        r_gdf = make_isochrones_disjoint(gdf, time_col='contour')
+        isos = r_gdf.to_geo_dict()
+    return isos

rapida/stats/raster_zonal_stats.py

@@ -257,5 +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)
+        #combined = combined.iloc[combined.geometry.area.sort_values(ascending=False).index]
         return combined