smoothing and barrier zonal stats

rapida/cli/__init__.py

@@ -18,10 +18,6 @@
 import click
 import nest_asyncio
 nest_asyncio.apply()
-import uvloop
-import asyncio
-
-asyncio.set_event_loop_policy(uvloop.EventLoopPolicy())
 
 
 

rapida/connectivity/__init__.py

@@ -50,20 +50,33 @@ async def run_connectivity_analysis(
         tar_path=dag_tar_path, origins=origins, travel_mode=travel_mode, intervals_minutes=time_intervals,
         radius=radius, disjoint=disjoint, smooth=smooth)
 
+    if clip_country:
+        url = f"/vsicurl/https://undpngddlsgeohubdev01.blob.core.windows.net/admin/cgaz/geoBoundariesCGAZ_ADM0.fgb"
+        a0_gdf = gpd.read_file(url, bbox=bbox, engine="pyogrio")
+        if not clip_country in a0_gdf['iso3'].tolist():
+            raise Exception(f'--clip-country {clip_country} does not intersect bbox {bbox}')
+        a0_gdf = a0_gdf[a0_gdf['iso3'] == clip_country]
+        if isochrones_gdf.crs != a0_gdf.crs:
+            a0_gdf = a0_gdf.to_crs(isochrones_gdf.crs)
+        isochrones_gdf = isochrones_gdf.clip(a0_gdf)
+        isochrones_gdf['iso3'] = clip_country
+
     water_bodies_path = await extract_water_bodies(pbf_path=bbox_pbf,dst_dir=dest_dir, progress=progress)
     water_gdf = gpd.read_file(water_bodies_path)
     if not water_gdf.empty:
         if isochrones_gdf.crs != water_gdf.crs:
             water_gdf = water_gdf.to_crs(isochrones_gdf.crs)
 
-        single_water_geom = water_gdf.geometry.union_all()
-        isochrones_gdf["geometry"] = isochrones_gdf.geometry.difference(single_water_geom)
+        water_poly = water_gdf.geometry.union_all()
+        isochrones_gdf["geometry"] = isochrones_gdf.geometry.difference(water_poly)
 
         # 3. Clean up empty/exploded geometries if any were cut into pieces
         isochrones_gdf = isochrones_gdf[~isochrones_gdf.is_empty].explode(index_parts=False)
     # 5. Save the final processed isochrones to GeoJSON
     isochrones_path = os.path.join(dest_dir, "isochrones.geojson")
+
     isochrones_gdf.to_file(isochrones_path, driver="GeoJSON")
+    del isochrones_gdf
 
     if pop_vars:
         with TemporaryDirectory(dir=dest_dir, delete=True) as project_folder:
@@ -97,15 +110,35 @@ async def run_connectivity_analysis(
                     index=False
                 )
 
-    if barriers_dataset is not None:
+
+    if barriers_dataset is not None and pop_vars:
         logger.info(f'Computing isochrones with barriers')
-        barrier_results = await connectivity_areas(
+
+        barrier_isochrones_gdf = 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, disjoint=disjoint
+        barriers_dataset=barriers_dataset, barriers_layer=barriers_layer, barriers_buffer=barriers_buffer, disjoint=disjoint, radius=radius,
+            smooth=smooth, progress=progress
                              )
+
+        if clip_country:
+            logger.info('Clipping barriers isochrones with ADM0')
+            url = f"/vsicurl/https://undpngddlsgeohubdev01.blob.core.windows.net/admin/cgaz/geoBoundariesCGAZ_ADM0.fgb"
+            a0_gdf = gpd.read_file(url, bbox=bbox, engine="pyogrio")
+            if not clip_country in a0_gdf['iso3'].tolist():
+                raise Exception(f'--clip-country {clip_country} does not intersect bbox {bbox}')
+            a0_gdf = a0_gdf[a0_gdf['iso3'] == clip_country]
+            if barrier_isochrones_gdf.crs != a0_gdf.crs:
+                a0_gdf = a0_gdf.to_crs(barrier_isochrones_gdf.crs)
+            barrier_isochrones_gdf = barrier_isochrones_gdf.clip(a0_gdf)
+            barrier_isochrones_gdf['iso3'] = clip_country
+        if not water_gdf.empty:
+            logger.info('Removing water bodies from barrier isochrones')
+            #barrier_isochrones_gdf["geometry"] = barrier_isochrones_gdf.geometry.difference(water_poly)
+            barrier_isochrones_gdf = barrier_isochrones_gdf.overlay(water_gdf, how='difference')
+            # 3. Clean up empty/exploded geometries if any were cut into pieces
+            barrier_isochrones_gdf = barrier_isochrones_gdf[~barrier_isochrones_gdf.is_empty].explode(index_parts=False)
         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)
+        barrier_isochrones_gdf.to_file(barrier_isochrones_path, driver="GeoJSON")
         if pop_vars:
             logger.info(f'Computing zonal stats for barrier isochrones')
             with TemporaryDirectory(dir=dest_dir, delete=True) as project_folder:
@@ -128,17 +161,18 @@ async def run_connectivity_analysis(
                     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]
+                    #
+                    # 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[f'{pvar}_{bar_pvar}_perc_difference'] = barrier_pop_stat_gdf[f'{pvar}_{bar_pvar}_difference'] / barrier_pop_stat_gdf[pvar] * 100
 
                     barrier_pop_stat_gdf.to_file(
                         filename=barrier_isochrones_path,

rapida/connectivity/graph.py

@@ -151,9 +151,9 @@ async def compile_valhalla_graph(pbf_path: str, dst_dir: str, progress=None) ->
 
     # Expand the max_locations limit to allow system-wide bulk routing
     valhalla_conf["service_limits"]["isochrone"]["max_locations"] = 50000
-    valhalla_conf["service_limits"]["isochrone"]["max_distance"] = 2000000
+    valhalla_conf["service_limits"]["isochrone"]["max_distance"] = 20000000
     valhalla_conf["service_limits"]["isochrone"]["max_contours"] = 20
-    valhalla_conf["service_limits"]["max_exclude_polygons_length"] = 500000  # Bump to 500km perimeter length
+    valhalla_conf["service_limits"]["max_exclude_polygons_length"] = 500000000  # Bump to 50000km perimeter length
     valhalla_conf["service_limits"]["allow_hard_exclusions"] = True
 
     # ---------------------------------------------------------

rapida/connectivity/io.py

@@ -170,7 +170,8 @@ async def prepare_osm_pbf(bbox: tuple[float, float, float, float], dst_dir: str
         if clip_country:
             if not clip_country in intersecting['iso3'].tolist():
                 logger.warning(f'cli-country={clip_country} was supplied but it was not found to intersect {bbox}')
-            intersecting = intersecting[intersecting['iso3'] == clip_country]
+            else:
+                intersecting = intersecting[intersecting['iso3'] == clip_country]
 
         pbf_urls = []
         if use_geofabrik:
@@ -533,12 +534,27 @@ def read_barriers(src_path: str, src_layer: str = None, barriers_buffer: float =
 
         if lyr is None:
             raise ValueError(f"Layer '{src_layer}' could not be found.")
+        source_srs = lyr.GetSpatialRef()
+        target_srs = osr.SpatialReference()
+        target_srs.ImportFromEPSG(4326)
+        # CRITICAL for GDAL 3+: Force traditional Longitude/Latitude (X/Y) axis order
+        target_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+
+        # 3. Check if they are different to avoid redundant transformations
+        needs_reprojection = False
+        coord_trans = None
 
+        if source_srs is not None and not source_srs.IsSame(target_srs):
+            needs_reprojection = True
+            coord_trans = osr.CoordinateTransformation(source_srs, target_srs)
         lyr.ResetReading()
         for feat in lyr:
             geom_ogr = feat.GetGeometryRef()
             if geom_ogr is None or geom_ogr.IsEmpty():
                 continue
+            # 4. Reproject the geometry in place
+            if needs_reprojection:
+                geom_ogr.Transform(coord_trans)
 
             wkb_output = geom_ogr.ExportToWkb()
             if wkb_output is None or isinstance(wkb_output, int):

rapida/connectivity/isochrone.py

@@ -7,8 +7,9 @@
 from shapely.geometry import shape, mapping, JOIN_STYLE
 from shapely import make_valid
 from pyproj import Transformer
-from shapely.ops import transform
+from shapely.ops import transform, unary_union
 import logging
+import numpy as np
 from rapida.connectivity.io import read_barriers
 logger = logging.getLogger(__name__)
 
@@ -27,6 +28,40 @@
 project_to_degrees = Transformer.from_crs("EPSG:3857", "EPSG:4326", always_xy=True).transform
 
 
+def get_empirical_radius(geom_meters, fallback_radius):
+    """Dynamically calculates the smoothing radius based on polygon segment lengths."""
+    try:
+        # Handle both Polygons and MultiPolygons safely
+        polys = geom_meters.geoms if geom_meters.geom_type == 'MultiPolygon' else [geom_meters]
+        lengths = []
+
+        for poly in polys:
+            coords = poly.exterior.coords
+            # Fast vectorized distance calculation between consecutive vertices
+            x = np.array([c[0] for c in coords])
+            y = np.array([c[1] for c in coords])
+            dist = np.sqrt(np.diff(x) ** 2 + np.diff(y) ** 2)
+
+            # Keep only meaningful segments (ignoring duplicate vertices < 1 meter)
+            lengths.extend(dist[dist > 1.0])
+
+        if lengths:
+            # The 15th percentile reliably targets the smallest common denominator (the grid step)
+            empirical_cell_size = np.percentile(lengths, 50)
+            # Round to the nearest meter to group floating-point variations
+            rounded_lengths = np.round(lengths, decimals=0)
+
+            # Find the most frequent rounded length (the mode)
+            values, counts = np.unique(rounded_lengths, return_counts=True)
+            empirical_cell_size = values[np.argmax(counts)]
+
+            # Apply the 0.85 multiplier to tightly fuse the grid without ballooning
+            return empirical_cell_size * 1.44
+
+    except Exception as e:
+        logger.warning(f"Empirical radius calculation failed, using fallback: {e}")
+
+    return fallback_radius
 
 def make_isochrones_disjoint(gdf, time_col="time", group_col=None):
     """Converts overlapping concentric isochrones into mutually exclusive rings."""
@@ -70,31 +105,6 @@ async def connectivity_areas(
     tar_file = Path(tar_path)
     build_config_file = tar_file.parent / "valhalla.json"
 
-    #runtime_config_file = tar_file.parent / "valhalla_runtime.json"
-
-    # # 1. Bypass Valhalla's hardcoded multi-origin security limits
-    # with open(build_config_file, "r") as f:
-    #     valhalla_config = json.load(f)
-    #
-    # if "service_limits" not in valhalla_config:
-    #     valhalla_config["service_limits"] = {}
-    # if "isochrone" not in valhalla_config["service_limits"]:
-    #     valhalla_config["service_limits"]["isochrone"] = {}
-    #
-    # # Update top-level orchestration limits
-    # valhalla_config["service_limits"]["isochrone"]["max_locations"] = 50000
-    # valhalla_config["service_limits"]["isochrone"]["max_distance"] = 2000000
-    # valhalla_config["service_limits"]["isochrone"]["max_contours"] = 20
-    #
-    #
-    # # THE EXACT MATCHING KEY FROM YOUR CONFIG:
-    # valhalla_config["service_limits"]["max_exclude_polygons_length"] = 500000  # Bump to 500km perimeter length
-    # valhalla_config["service_limits"]["allow_hard_exclusions"] = True
-    #
-    # with open(runtime_config_file, "w") as f:
-    #     json.dump(valhalla_config, f)
-
-
     contours = [{"time": int(mins)} for mins in intervals_minutes]
     barriers_coords = read_barriers(src_path=barriers_dataset, src_layer=barriers_layer, barriers_buffer=barriers_buffer)
     locations = [{"lon": float(lon), "lat": float(lat), "radius": radius} for lon, lat in origins]
@@ -172,27 +182,52 @@ def run_routing():
                     real_cell_size_meters = valhalla_degree_step * 111320
 
                     # 4. Use this true runtime value for your smoothing radius (e.g., 1.5x to 2x the cell size)
-                    smooth_radius_meters = real_cell_size_meters * 1.5
+                    smooth_radius_meters = real_cell_size_meters * .85
 
                     # 5. Apply the Morphological Opening/Closing (Buffer out, in, out)
                     geom_meters = make_valid(transform(project_to_meters, raw_geom_wgs84))
-
-
-                    # 4. Morphological Closing (Now using actual physical meters)
-                    smooth_geom_meters = geom_meters.buffer(
-                        smooth_radius_meters,
-                        join_style=JOIN_STYLE.round
-                    ).buffer(
-                        -(smooth_radius_meters * 2),
-                        join_style=JOIN_STYLE.round
-                    ).buffer(
-                        smooth_radius_meters,
-                        join_style=JOIN_STYLE.round
-                    )
-
-                    # 5. Metric Simplification (Drop vertices closer than 50 meters to the line)
-                    smooth_geom_meters = smooth_geom_meters.simplify(20, preserve_topology=True)
-
+                    # 6. Explode MultiPolygon into individual Polygons
+                    polys = geom_meters.geoms if geom_meters.geom_type == 'MultiPolygon' else [geom_meters]
+                    smoothed_polys = []
+
+                    for poly in polys:
+                        smooth_radius_meters1 = get_empirical_radius(poly, smooth_radius_meters)
+                        # Morphological Closing on each individual polygon
+                        smoothed = poly.buffer(
+                            smooth_radius_meters1,
+                            join_style=JOIN_STYLE.round
+                        ).buffer(
+                            -smooth_radius_meters1,
+                            join_style=JOIN_STYLE.round
+                        )
+
+                        # Simplification on each individual polygon
+                        smoothed = smoothed.simplify(50, preserve_topology=True)
+
+                        # Ignore any polygons that completely disappeared during the negative buffer
+                        if not smoothed.is_empty:
+                            smoothed_polys.append(smoothed)
+
+                    # 7. Recombine back into a valid MultiPolygon
+                    # Using unary_union safely merges any internal overlaps that the buffering might have caused
+                    smooth_geom_meters = make_valid(unary_union(smoothed_polys))
+
+                    # # 5. Calculate empirical radius
+                    # smooth_radius_meters1 = get_empirical_radius(geom_meters, smooth_radius_meters)
+                    #
+                    # # 4. Morphological Closing (Buffer OUT, then IN by the same amount)
+                    # # This fills the jagged grid gaps and rounds corners without severing thin corridors.
+                    # smooth_geom_meters = geom_meters.buffer(
+                    #     smooth_radius_meters1,
+                    #     join_style=JOIN_STYLE.round
+                    # ).buffer(
+                    #     -smooth_radius_meters1,
+                    #     join_style=JOIN_STYLE.round
+                    # )
+                    #
+                    # # 5. Metric Simplification (Drop vertices closer than 50 meters to the line)
+                    # smooth_geom_meters = smooth_geom_meters.simplify(20, preserve_topology=True)
+                    #
                     # 6. Convert back to WGS84 degrees so the GeoJSON renders on a map properly
                     final_geom_wgs84 = transform(project_to_degrees, smooth_geom_meters)