From 4b8ec04163309757e6b4396251390675ed0d4269 Mon Sep 17 00:00:00 2001 From: Ioan Ferencik Date: Wed, 8 Jul 2026 16:29:19 +0200 Subject: [PATCH] implement pop zonal stats and disjoint isochrones --- rapida/cli/__init__.py | 13 +-- rapida/cli/aclick.py | 41 +++++----- rapida/cli/connectivity.py | 19 ++++- rapida/connectivity/__init__.py | 125 ++++++++++++++++++++++------- rapida/connectivity/isochrone.py | 46 +++++++++-- rapida/stats/raster_zonal_stats.py | 2 +- 6 files changed, 184 insertions(+), 62 deletions(-) diff --git a/rapida/cli/__init__.py b/rapida/cli/__init__.py index 481dc4b..278d409 100644 --- a/rapida/cli/__init__.py +++ b/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()) +# diff --git a/rapida/cli/aclick.py b/rapida/cli/aclick.py index bf1ee26..d61ea39 100644 --- a/rapida/cli/aclick.py +++ b/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): diff --git a/rapida/cli/connectivity.py b/rapida/cli/connectivity.py index 36b4996..d0d0067 100644 --- a/rapida/cli/connectivity.py +++ b/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 ) \ No newline at end of file diff --git a/rapida/connectivity/__init__.py b/rapida/connectivity/__init__.py index 1ba6886..c75e40b 100644 --- a/rapida/connectivity/__init__.py +++ b/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 \ No newline at end of file diff --git a/rapida/connectivity/isochrone.py b/rapida/connectivity/isochrone.py index a7b10c4..a555958 100644 --- a/rapida/connectivity/isochrone.py +++ b/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) \ No newline at end of file + 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 \ No newline at end of file diff --git a/rapida/stats/raster_zonal_stats.py b/rapida/stats/raster_zonal_stats.py index 51ac901..3da8d20 100644 --- a/rapida/stats/raster_zonal_stats.py +++ b/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 \ No newline at end of file