From 1a4eca8b72b88a84e1ddb7c7570bd4a3a6d4e2ef Mon Sep 17 00:00:00 2001 From: Ioan Ferencik Date: Thu, 30 Jul 2026 17:17:05 +0200 Subject: [PATCH] smoothing and barrier zonal stats --- rapida/cli/__init__.py | 4 - rapida/connectivity/__init__.py | 70 +++++++++++++----- rapida/connectivity/graph.py | 4 +- rapida/connectivity/io.py | 18 ++++- rapida/connectivity/isochrone.py | 123 ++++++++++++++++++++----------- rapida/project/project.py | 94 ++++++++++++++++------- 6 files changed, 219 insertions(+), 94 deletions(-) diff --git a/rapida/cli/__init__.py b/rapida/cli/__init__.py index 4d5181a..e36b2eb 100644 --- a/rapida/cli/__init__.py +++ b/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()) diff --git a/rapida/connectivity/__init__.py b/rapida/connectivity/__init__.py index e9482c1..535b261 100644 --- a/rapida/connectivity/__init__.py +++ b/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, diff --git a/rapida/connectivity/graph.py b/rapida/connectivity/graph.py index 37cb5d0..01c6da9 100644 --- a/rapida/connectivity/graph.py +++ b/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 # --------------------------------------------------------- diff --git a/rapida/connectivity/io.py b/rapida/connectivity/io.py index 1a5f349..0cf2150 100644 --- a/rapida/connectivity/io.py +++ b/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): diff --git a/rapida/connectivity/isochrone.py b/rapida/connectivity/isochrone.py index ebbdb5e..bf9d141 100644 --- a/rapida/connectivity/isochrone.py +++ b/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) diff --git a/rapida/project/project.py b/rapida/project/project.py index 721e07a..acff8d8 100644 --- a/rapida/project/project.py +++ b/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 +import geopandas as gpd from rapida import constants from rapida.az.blobstorage import check_blob_exists, delete_blob from rapida.session import Session @@ -107,8 +107,10 @@ def __init__(self, path: str,polygons: str = None, if polygons is not None: l = geopandas.list_layers(polygons) + lnames = l.name.tolist() lcount = len(lnames) + if lcount > 1: click.echo(f'{polygons} contains {lcount} layers: {",".join(lnames)}') layer_name = click.prompt( @@ -142,34 +144,38 @@ def __init__(self, path: str,polygons: str = None, if "iso3" in cols and gdf['iso3'].isna().any(): gdf.drop(columns=["iso3"], inplace=True) cols = gdf.columns.tolist() - if 'iso3' not in cols: + + if not 'iso3' in cols: logger.info( 'ISO3 column missing or some values are empty in ISO3 column. ' - 'Going to add country codes into "iso3" column' + 'Going to add country codes into "iso3" column via spatial join.' ) - geo_srs = osr.SpatialReference() - geo_srs.ImportFromEPSG(4326) - geo_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) - coord_trans = osr.CoordinateTransformation(self.target_srs, geo_srs) - centroids = gdf.geometry.centroid - transformed = [coord_trans.TransformPoint(x, y) for x, y in zip(centroids.x, centroids.y)] - lats = [t[1] for t in transformed] - lons = [t[0] for t in transformed] - - iso3_codes = [] - for lat, lon in zip(lats, lons): - try: - 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"])))) + url = f"/vsicurl/https://undpngddlsgeohubdev01.blob.core.windows.net/admin/cgaz/geoBoundariesCGAZ_ADM0.fgb" + # 1. Extract centroids first + centroids = gdf.copy() + if 'iso3' in centroids.columns: + centroids = centroids.drop(columns=['iso3']) + centroids.geometry = centroids.geometry.centroid.representative_point() + + # 2. Project centroids to EPSG:4326 for the remote mask query + centroids_4326 = centroids.to_crs(epsg=4326) + centroids_4326.to_file('/tmp/isocen.fgb', engine='pyogrio') + + # 3. Use 'mask' instead of 'bbox'. This asks GDAL to only download + # features that intersect your specific points, ignoring empty space. + a0_gdf = gpd.read_file(url, mask=centroids_4326, engine="pyogrio") + # 4. Perform vectorized spatial join (both are now in EPSG:4326) + joined = gpd.sjoin(centroids_4326, a0_gdf, how='left', predicate='within') + + # 5. Drop potential duplicates from borders and assign back to original gdf + joined = joined[~joined.index.duplicated(keep='first')] + gdf['iso3'] = joined['iso3'] + + # # 6. Generate the tuple of valid country codes + # self.countries = tuple(sorted(set(filter(lambda x: x in COUNTRY_CODES, gdf["iso3"])))) + + # 7. Export to GeoPackage gdf.to_file( filename=self.geopackage_file_path, driver="GPKG", @@ -180,6 +186,44 @@ def __init__(self, path: str,polygons: str = None, index=False ) + # if 'iso3' not in cols: + # logger.info( + # 'ISO3 column missing or some values are empty in ISO3 column. ' + # 'Going to add country codes into "iso3" column' + # ) + # geo_srs = osr.SpatialReference() + # geo_srs.ImportFromEPSG(4326) + # geo_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) + # coord_trans = osr.CoordinateTransformation(self.target_srs, geo_srs) + # centroids = gdf.geometry.centroid + # transformed = [coord_trans.TransformPoint(x, y) for x, y in zip(centroids.x, centroids.y)] + # lats = [t[1] for t in transformed] + # lons = [t[0] for t in transformed] + # + # iso3_codes = [] + # for lat, lon in zip(lats, lons): + # try: + # 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", + # engine="pyogrio", + # mode="w", + # layer=self.polygons_layer_name, + # promote_to_multi=True, + # index=False + # ) + if 'h3id' in cols: h3ids = gdf['h3id'].tolist() # check duplicated