Fix/barriers: remove isochrones inside large barriers areas

rapida/connectivity/__init__.py

@@ -1,8 +1,8 @@
 import datetime
 
 import os.path
-from unittest.mock import inplace
 
+from shapely.geometry import Point
 from rapida.util.bbox_param_type import get_best_semantic_label
 import geopandas as gpd
 import pandas as pd
@@ -243,16 +243,34 @@ async def run_connectivity_analysis(
 
 
 
-    if barriers_dataset is not None and pop_vars:
+    if barriers_dataset is not None:
         info = pyogrio.read_info(barriers_dataset)
+        # 1. Load barriers and align CRS with the EPSG:4326 origins
+        barriers_gdf = gpd.read_file(barriers_dataset, engine="pyogrio")
+        if barriers_gdf.crs and barriers_gdf.crs != "EPSG:4326":
+            barriers_gdf = barriers_gdf.to_crs("EPSG:4326")
+
+        # 2. Merge all barrier polygons into a single geometry for fast intersection checking
+        unified_barriers = barriers_gdf.geometry.union_all()
+
+        logger.info(f'Filtering out origins located inside barriers...')
+
+        # 3. Filter the list of (lon, lat) tuples using Shapely
+        valid_origins = [
+            (lon, lat) for lon, lat in origins
+            if not unified_barriers.contains(Point(lon, lat))
+        ]
+
+        logger.info(f'Proceeding with {len(valid_origins)} valid origins outside barrier zones.')
+
         if 'polygon' in info['geometry_type'].lower(): # keep only polys that actually intersect the roads
             logger.info(f'Removing barrier polygons that do not intersect roads...')
             roads_dataset = await extract_roads(pbf_path=bbox_pbf, dst_dir=dest_dir, progress=progress)
             barriers_dataset = filter_polygons(poly_ds_path=barriers_dataset, lines_ds_path=roads_dataset, dst_dir=dest_dir)
 
         logger.info(f'Computing isochrones with barriers')
         barrier_isochrones_gdf = await connectivity_areas(
-        tar_path=dag_tar_path, origins=origins, travel_mode=travel_mode, intervals_minutes=time_intervals,
+        tar_path=dag_tar_path, origins=valid_origins, travel_mode=travel_mode, intervals_minutes=time_intervals,
         barriers_dataset=barriers_dataset, barriers_layer=barriers_layer, barriers_buffer=barriers_buffer, disjoint=disjoint, radius=radius,
             smooth=smooth, progress=progress
                              )
@@ -268,6 +286,14 @@ async def run_connectivity_analysis(
                 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
+
+        logger.info('Removing barrier areas from barrier isochrones')
+        if barrier_isochrones_gdf.crs != barriers_gdf.crs:
+            barriers_gdf = barriers_gdf.to_crs(barrier_isochrones_gdf.crs)
+
+        barriers_poly = barriers_gdf.geometry.union_all()
+        barrier_isochrones_gdf["geometry"] = barrier_isochrones_gdf.geometry.difference(barriers_poly)
+
         if not water_gdf.empty:
             logger.info('Removing water bodies from barrier isochrones')
             barrier_isochrones_gdf["geometry"] = barrier_isochrones_gdf.geometry.difference(water_poly)

rapida/connectivity/isochrone.py

@@ -28,7 +28,7 @@
 project_to_degrees = Transformer.from_crs("EPSG:3857", "EPSG:4326", always_xy=True).transform
 
 
-def get_empirical_radius(geom_meters, fallback_radius):
+def get_empirical_radius_old(geom_meters, fallback_radius):
     """Dynamically calculates the smoothing radius based on polygon segment lengths."""
     try:
         # Handle both Polygons and MultiPolygons safely
@@ -62,6 +62,48 @@ def get_empirical_radius(geom_meters, fallback_radius):
         logger.warning(f"Empirical radius calculation failed, using fallback: {e}")
 
     return fallback_radius
+def get_empirical_radius(geom_meters, fallback_radius):
+    """Dynamically calculates the smoothing radius based on polygon segment lengths."""
+    try:
+        # Safely extract only Polygon geometries, handling GeometryCollections and MultiPolygons
+        def extract_polygons(geom):
+            if geom.geom_type == 'Polygon':
+                return [geom]
+            elif hasattr(geom, 'geoms'):
+                polys = []
+                for part in geom.geoms:
+                    polys.extend(extract_polygons(part))
+                return polys
+            return []
+
+        polys = extract_polygons(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:
+            # 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."""