diff --git a/deploy.md b/deploy.md index f976ae2..3927738 100644 --- a/deploy.md +++ b/deploy.md @@ -22,9 +22,13 @@ curl -O https://raw.githubusercontent.com/UNDP-Data/rapida/refs/heads/main/deplo curl.exe -O https://raw.githubusercontent.com/UNDP-Data/rapida/refs/heads/main/deploy/pixi.toml ``` -## .env file -Create a .env file and define following environmental variables. -Alternatively ensure the environmental variables defined below exists in your system +## Environmental variables + +There are two types of environmental variables: mandatory and optional. +The **mandatory** ones contain information that is required to operate specific parts of rapida like +downloading imagery from earthdata through a token or authenticating to [space-track.org](www.space-track.org). +Rapida throws and error if these variables are not defined and needed. + ```shell # for uploading to azure TENANT_ID= @@ -34,8 +38,14 @@ EARTHDATA_TOKEN= # for predicting precisely VIIRS orbits SPACETRACK_USER= SPACETRACK_PASSWORD= +``` +The **optional** environmental variables do not generate errors if not defined. Rapida lets the user know whenever these +variables arte detected and detected and override/change default behaviour with a custom one. + +```shell # override road type default seed for connectivity analysis +#Default values #"motorway": 105, #"trunk": 90, @@ -45,10 +55,22 @@ SPACETRACK_PASSWORD= #"unclassified": 40, #"residential": 35, #"service": 25 +# variables that change default values + +#MJOLNIR_MOTORWAY_SPEED=60 +#MJOLNIR_TRUNK_SPEED=50 +#MJOLNIR_PRIMARY_SPEED=35 +#MJOLNIR_SECONDARY_SPEED=25 +#MJOLNIR_UNCLASSIFIED_SPEED=15 +#MJOLNIR_RESIDENTIAL_SPEED=12 +#MJOLNIR_SERVICE_SPEED=8 +#CONNECTIVITY_OSM_SOURCE=MOVISDA # MOVISDA, GEOFABRIK -MJOLNIR_SECONDARY_SPEED=40 ``` + + + ## run NTL ```shell pixi run rapida ntl detect --help diff --git a/deploy/pixi.toml b/deploy/pixi.toml index 84985f4..380e983 100644 --- a/deploy/pixi.toml +++ b/deploy/pixi.toml @@ -8,23 +8,36 @@ platforms = ["linux-64", "win-64", "osx-arm64"] python = ">=3.12" numpy = "<2" gdal = ">=3.12.0" -fiona = "*" -rasterio = "*" -geopandas = "*" -shapely = "*" -exactextract = "*" -rio-cogeo = "*" -pyarrow = "*" -h5py = ">=3.16.0" -satpy = ">=0.59.0" -matplotlib = ">=3.10.9" -netcdf4 = ">=1.7.3" -h5netcdf = ">=1.8.1" +libgdal-netcdf = "*" +osmium-tool = ">=1.19.1,<2" +#playwright="*" +#fiona = "*" +#rasterio = "*" +#geopandas = "*" +#shapely = "*" +#exactextract = "*" +#rio-cogeo = "*" +#pyarrow = "*" +#h5py = ">=3.16.0" +#satpy = ">=0.59.0" +#matplotlib = ">=3.10.9" +#netcdf4 = ">=1.7.3" +#h5netcdf = ">=1.8.1" +#marimo = ">=0.23.14,<0.24" [pypi-dependencies] -# Pixi natively understands Git URLs in the pypi section! +# Pixi natively understands Git URLs in the pypi section in pyproject.toml rapida = { git = "https://github.com/UNDP-Data/rapida.git", branch = "main" } dotenv-cli="*" +#maplibre = { version = ">=0.3.6, <0.4", extras = ["ipywidget"] } +#openlayers = ">=0.1.6, <0.2" + +[activation.env] +# Forces Playwright to install browsers inside the local project folder +PLAYWRIGHT_BROWSERS_PATH = ".pixi/playwright-browsers" + [tasks] -rapida = "dotenv -e .env rapida" \ No newline at end of file +setup = {cmd = "dotenv -e .env playwright install chromium --only-shell"} +rapida = {cmd = "dotenv -e .env rapida", depends-on=["setup"]} +#rapida = "dotenv -e .env rapida" \ No newline at end of file diff --git a/rapida/cli/connectivity.py b/rapida/cli/connectivity.py index 93ab218..383bc8a 100644 --- a/rapida/cli/connectivity.py +++ b/rapida/cli/connectivity.py @@ -98,7 +98,11 @@ def parse_intervals(ctx, param, value): type=str, callback=validate_variables, help=f"One or more RAPIDA population variable to compute zonal stats for withing the connectivity zones" ) - +@click.option('-l','--stats-admin-level', + type=click.IntRange(min=1, max=2, clamp=False), + help='Admin level from where to extract the admin features when computing stats for --popvar/s.'\ + 'The statistics will be computed from the results of intersecting the isochrones with admin data' + ) @click.option( "--dst-dir", "-d", # Short option @@ -134,16 +138,16 @@ def parse_intervals(ctx, param, value): ) -@click.option( - '--smooth', - is_flag=True, - help=( - "By default Valhalla routing engine create isochrones using a grid and the edges of polygons are squarish." - "Use this flag to smoothen them out. in some instances this can alter the isochrones significantly." - "In general it is save to use smoothing on small areas (towns, cities)" - ), - default=False -) +# @click.option( +# '--smooth', +# is_flag=True, +# help=( +# "By default Valhalla routing engine create isochrones using a grid and the edges of polygons are squarish." +# "Use this flag to smoothen them out. in some instances this can alter the isochrones significantly." +# "In general it is save to use smoothing on small areas (towns, cities)" +# ), +# default=True +# ) @click.option( '--disjoint', @@ -161,9 +165,9 @@ def parse_intervals(ctx, param, value): 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, stats_admin_level:int=None, max_snap_distance:float=None, clip_country:str=None, - disjoint:bool=False, smooth:bool=False + disjoint:bool=False, smooth:bool=True ): logger.info(f'Running connectivity analysis') progress = ctx.obj.get('progress') @@ -171,6 +175,6 @@ async def connectivity(ctx, bbox:tuple[float, float, float, float]=None, travel_ return await run_connectivity_analysis( 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, + sites_dataset=sites_dataset, sites_layer=sites_layer, pop_vars=popvar,stats_admin_level=stats_admin_level, progress=progress, radius=max_snap_distance, clip_country=clip_country, disjoint=disjoint, smooth=smooth ) \ No newline at end of file diff --git a/rapida/components/population/__init__.py b/rapida/components/population/__init__.py index e4a5efb..32ba035 100644 --- a/rapida/components/population/__init__.py +++ b/rapida/components/population/__init__.py @@ -511,7 +511,7 @@ def import_raster(self, source=None, **kwargs): **kwargs ) - os.remove(source) + if os.path.exists(source):os.remove(source) os.rename(imported_local_path, source) return source diff --git a/rapida/connectivity/__init__.py b/rapida/connectivity/__init__.py index 535b261..d1bcd5c 100644 --- a/rapida/connectivity/__init__.py +++ b/rapida/connectivity/__init__.py @@ -1,24 +1,37 @@ import datetime -import json + import os.path +from unittest.mock import inplace + from rapida.util.bbox_param_type import get_best_semantic_label import geopandas as gpd +import pandas as pd import logging from rich.progress import Progress -from rapida.connectivity.io import prepare_osm_pbf,extract_health_sites, extract_origins_from_geojson, extract_origins, extract_water_bodies +from rapida.connectivity.io import (prepare_osm_pbf,extract_health_sites, + extract_origins_from_geojson, extract_origins, + extract_water_bodies, extract_roads, filter_polygons) 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 +import gc +from math import nan +import pyogrio 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, + sites_dataset:str=None, sites_layer:str=None,pop_vars:str|tuple[str]=None, stats_admin_level:int=None, progress:Progress=None, year=datetime.datetime.now().year, disjoint:bool=False, radius:float=None, clip_country:str=None, smooth:bool=False ): @@ -45,12 +58,13 @@ async def run_connectivity_analysis( dag_tar_path = await compile_valhalla_graph(pbf_path=bbox_pbf,dst_dir=dest_dir, progress=progress) origins = extract_origins(sites_dataset=sites, src_layer=sites_layer) - + logger.info(f'Computing isochrones for {len(origins)} sites') isochrones_gdf = await connectivity_areas( tar_path=dag_tar_path, origins=origins, travel_mode=travel_mode, intervals_minutes=time_intervals, radius=radius, disjoint=disjoint, smooth=smooth) if clip_country: + logger.info('Clipping 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(): @@ -64,6 +78,7 @@ async def run_connectivity_analysis( 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: + logger.info('Removing water bodies from isochrones') if isochrones_gdf.crs != water_gdf.crs: water_gdf = water_gdf.to_crs(isochrones_gdf.crs) @@ -72,13 +87,50 @@ async def run_connectivity_analysis( # 3. Clean up empty/exploded geometries if any were cut into pieces isochrones_gdf = isochrones_gdf[~isochrones_gdf.is_empty].explode(index_parts=False) + #TD + # ADD THIS: Dissolve by contour to merge overlaps from different sites + isochrones_gdf = isochrones_gdf.dissolve(by='contour', as_index=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 + isochrones_gdf.to_file(isochrones_path, driver="GeoJSON", engine='pyogrio') + if pop_vars: + logger.info(f'Computing zonal stats for regular isochrones') + if stats_admin_level: + url = f"/vsicurl/https://undpngddlsgeohubdev01.blob.core.windows.net/admin/cgaz/geoBoundariesCGAZ_ADM{stats_admin_level}.fgb" + adm_gdf = gpd.read_file(url, bbox=bbox, engine="pyogrio") + if clip_country: + adm_gdf = adm_gdf[adm_gdf['iso3'] == clip_country] + if 'iso3' in isochrones_gdf.columns.tolist(): + adm_gdf.drop(columns=['iso3'], inplace=True) + if isochrones_gdf.crs != adm_gdf.crs: + adm_gdf.to_crs(isochrones_gdf.crs, inplace=True) + results = [] + for i, unit in adm_gdf.iterrows(): + unit_gdf = gpd.GeoDataFrame([unit], crs=adm_gdf.crs, geometry=adm_gdf.geometry.name) + # 2. Fast pre-filter: Skip admin units that don't even touch the isochrones' bounding boxes + if isochrones_gdf.sindex.intersection(unit_gdf.total_bounds).size == 0: + continue + # 3. Run the overlay chunk + chunk_result = gpd.overlay(isochrones_gdf, unit_gdf, how="intersection", keep_geom_type=True) + + if not chunk_result.empty: + results.append(chunk_result) + + # 3. Recombine into the final Spatially-enabled DataFrame + split_isochrones = gpd.GeoDataFrame( + pd.concat(results, ignore_index=True), + crs=isochrones_gdf.crs + ) + split_isochrones.to_file(isochrones_path, driver="GeoJSON") + del split_isochrones + del isochrones_gdf + del results + 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: @@ -112,8 +164,13 @@ async def run_connectivity_analysis( if barriers_dataset is not None and pop_vars: - logger.info(f'Computing isochrones with barriers') + info = pyogrio.read_info(barriers_dataset) + 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, barriers_dataset=barriers_dataset, barriers_layer=barriers_layer, barriers_buffer=barriers_buffer, disjoint=disjoint, radius=radius, @@ -133,14 +190,55 @@ async def run_connectivity_analysis( 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') + 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) + # ADD THIS: Dissolve by contour to merge overlaps from different sites + barrier_isochrones_gdf = barrier_isochrones_gdf.dissolve(by='contour', as_index=False) + + barrier_isochrones_path = os.path.join(dest_dir, 'isochrones_with_barriers.geojson') - barrier_isochrones_gdf.to_file(barrier_isochrones_path, driver="GeoJSON") + barrier_isochrones_gdf.to_file(barrier_isochrones_path, driver="GeoJSON", engine="pyogrio") if pop_vars: logger.info(f'Computing zonal stats for barrier isochrones') + + if stats_admin_level: + url = f"/vsicurl/https://undpngddlsgeohubdev01.blob.core.windows.net/admin/cgaz/geoBoundariesCGAZ_ADM{stats_admin_level}.fgb" + adm_gdf = gpd.read_file(url, bbox=bbox, engine="pyogrio") + if clip_country: + adm_gdf = adm_gdf[adm_gdf['iso3'] == clip_country] + if 'iso3' in barrier_isochrones_gdf.columns.tolist(): + adm_gdf.drop(columns=['iso3'], inplace=True) + if barrier_isochrones_gdf.crs != adm_gdf.crs: + adm_gdf.to_crs(barrier_isochrones_gdf.crs, inplace=True) + results = [] + for i, unit in adm_gdf.iterrows(): + unit_gdf = gpd.GeoDataFrame([unit], crs=adm_gdf.crs, geometry=adm_gdf.geometry.name) + # 2. Fast pre-filter: Skip admin units that don't even touch the isochrones' bounding boxes + if barrier_isochrones_gdf.sindex.intersection(unit_gdf.total_bounds).size == 0: + continue + # 3. Run the overlay chunk + chunk_result = gpd.overlay(barrier_isochrones_gdf, unit_gdf, how="intersection", keep_geom_type=True) + + if not chunk_result.empty: + results.append(chunk_result) + + # 3. Recombine into the final Spatially-enabled DataFrame + split_barrier_isochrones = gpd.GeoDataFrame( + pd.concat(results, ignore_index=True), + crs=barrier_isochrones_gdf.crs + ) + + + split_barrier_isochrones.to_file(barrier_isochrones_path, driver="GeoJSON", engine='pyogrio', promote_to_multi=True, + index=False) + + del barrier_isochrones_gdf + del split_barrier_isochrones + del results + + 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: @@ -156,32 +254,130 @@ async def run_connectivity_analysis( 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.iloc[barrier_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[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, - driver="GeoJSON", - engine="pyogrio", - mode="w", - layer='barrier_isochrones', - promote_to_multi=True, - index=False - ) + + if not stats_admin_level: + 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 + ) + else: + + logger.info('Aggregating zonal stats independently and pivoting to wide format...') + + 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)) + admin_col_name = f'admin{stats_admin_level}_name' + assert admin_col_name in pop_stat_gdf.columns.tolist() + + pop_stat_gdf['contour'] = pop_stat_gdf['contour'].astype(float).astype(int) + barrier_pop_stat_gdf['contour'] = barrier_pop_stat_gdf['contour'].astype(float).astype(int) + + pop_stat_gdf[admin_col_name] = pop_stat_gdf[admin_col_name].astype(str).str.strip() + barrier_pop_stat_gdf[admin_col_name] = barrier_pop_stat_gdf[admin_col_name].astype( + str).str.strip() + + # 1. Group and sum the base stats independently (vectorized aggregation) + base_agg = pop_stat_gdf.groupby([admin_col_name, 'contour'])[pop_col_names].sum().reset_index() + + # 2. Group and sum the barrier stats independently (vectorized aggregation) + barrier_agg = barrier_pop_stat_gdf.groupby([admin_col_name, 'contour'])[ + pop_col_names].sum().reset_index() + barrier_agg.rename(columns=col_name_dict, inplace=True) + + # 3. Merge the tiny aggregated tables + # outer join ensures we don't lose contours if one dataset has contours the other doesn't + grouped_df = base_agg.merge(barrier_agg, on=[admin_col_name, 'contour'], how='outer') + + # 4. Handle NAs and compute differences vector-wise + val_columns = [] + for pvar, bar_pvar in col_name_dict.items(): + grouped_df[pvar] = grouped_df[pvar].fillna(0) + grouped_df[bar_pvar] = grouped_df[bar_pvar].fillna(0) + + # Clean column names for the differences + diff_col = f'{pvar}_diff' + perc_col = f'{pvar}_perc_diff' + + # Calculate differences + grouped_df[diff_col] = (grouped_df[pvar] - grouped_df[bar_pvar]).clip(lower=0) + + # Calculate percentage (safeguarding against division by zero) + safe_div = grouped_df[pvar].replace(0, nan) + grouped_df[perc_col] = (grouped_df[diff_col] / safe_div).fillna(0) * 100 + + val_columns.extend([pvar, bar_pvar, diff_col, perc_col]) + + # 5. Pivot to Wide Format (Option 2: contours become columns) + pivot_df = grouped_df.pivot(index=admin_col_name, columns='contour', values=val_columns) + + # Flatten the MultiIndex columns (e.g., ('male_total_2026', 15.0) -> 'male_total_2026_15min') + pivot_df.columns = [f"{col[0]}_{int(col[1])}min" for col in pivot_df.columns] + wide_df = pivot_df.reset_index() + + # 6. Attach geometries from the original adm_gdf + logger.info('Merging aggregated stats back onto original admin boundaries...') + + # Handle case where CGAZ admin column might natively be 'shapeName' + adm_join_col = admin_col_name if admin_col_name in adm_gdf.columns else 'shapeName' + + final_gdf = adm_gdf.merge( + wide_df, + left_on=adm_join_col, + right_on=admin_col_name, + how='inner' # Use 'inner' to only keep admin units that actually had isochrones + ) + + gc.collect() + + with TemporaryDirectory(dir=dest_dir, delete=True) as admin_project_folder: + logger.info(f'Computing zonal stats for total population ') + adm_ds_path = os.path.join(dest_dir, f'admin_{stats_admin_level}.fgb') + adm_gdf.to_file(adm_ds_path, driver="FlatGeobuf", engine="pyogrio") + admin_project = Project(path=admin_project_folder, polygons=adm_ds_path, + comment='temp project for admin stats') + 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=['total'], + year=year, + project=admin_project.path, + force=False + ) + + admin_stat_gpkg_path = os.path.join(admin_project_folder, 'data', f'{admin_project.name}.gpkg') + admin_pop_stat_gdf = gpd.read_file(admin_stat_gpkg_path, layer='stats.population') + final_gdf = final_gdf.merge(admin_pop_stat_gdf[[admin_col_name, f'total_{year}']],on=admin_col_name) + + if os.path.exists(adm_ds_path):os.remove(adm_ds_path) + + admin_iso_stats = os.path.join(dest_dir, f"admin{stats_admin_level}_iso_stats.geojson") + logger.info(f"Writing final aggregated admin boundaries to {admin_iso_stats}") + + final_gdf.to_file( + filename=admin_iso_stats, # Fixed: this was pointing to barrier_isochrones_path + driver="GeoJSON", + engine="pyogrio", + mode="w", + layer=f"admin{stats_admin_level}_iso_stats", + promote_to_multi=True, + index=False + ) return \ No newline at end of file diff --git a/rapida/connectivity/io.py b/rapida/connectivity/io.py index 0cf2150..31eb6a4 100644 --- a/rapida/connectivity/io.py +++ b/rapida/connectivity/io.py @@ -661,4 +661,125 @@ def process_water_geometries(): if progress: progress.console.print(f"[bold green]✓ Water bodies successfully extracted to: {final_geojson}[/bold green]") - return str(final_geojson) \ No newline at end of file + return str(final_geojson) + + +import asyncio +from pathlib import Path +import geopandas as gpd + + +async def extract_roads(pbf_path: str, dst_dir: str, progress=None) -> gpd.GeoDataFrame: + """ + Extracts road network via Osmium and loads it directly into a GeoDataFrame. + """ + dst_path = Path(dst_dir) + filtered_pbf = dst_path / "roads_only.osm.pbf" + roads_geojsonseq = dst_path / "roads.geojsonseq" + + tags_to_keep = [ + "w/highway=motorway,trunk,primary,secondary,tertiary,unclassified,residential" + ] + + if progress: + progress.console.print("[cyan]Filtering roads from OSM via Osmium...[/cyan]") + + # Step 1: Filter PBF down to specified road types + await asyncio.to_thread( + run_cli, + ["osmium", "tags-filter", str(pbf_path)] + tags_to_keep + ["-o", str(filtered_pbf), "--overwrite"] + ) + + if progress: + progress.console.print("[cyan]Exporting roads to GeoJSONSeq...[/cyan]") + + # Step 2: Export to geojsonseq (LineStrings only) + await asyncio.to_thread( + run_cli, + ["osmium", "export", str(filtered_pbf), "-o", str(roads_geojsonseq), "-f", "geojsonseq", + "--geometry-type=linestring", "--overwrite"] + ) + + # if progress: + # progress.console.print("[cyan]Loading roads into GeoDataFrame...[/cyan]") + return roads_geojsonseq + # # Step 3: Load directly into GeoDataFrame + # def load_gdf(): + # # Using pyogrio to parse the geojsonseq significantly faster than Fiona + # return gpd.read_file(roads_geojsonseq, engine="pyogrio") + # + # roads_gdf = await asyncio.to_thread(load_gdf) + # + # # Step 4: Clean up intermediate files + # for path in [filtered_pbf, roads_geojsonseq]: + # if path.exists(): + # path.unlink() + # + # if progress: + # progress.console.print("[bold green]✓ Roads successfully extracted into GeoDataFrame[/bold green]") + # + # return roads_gdf + + +def filter_polygons(poly_ds_path:str=None, poly_layer:str=None, lines_ds_path:str=None, dst_dir:str=None ): + dst_path = Path(dst_dir) + dst_filtered_poly_ds_path = dst_path / 'filtered_barriers.fgb' + # 1. Load and explode MultiPolygons to single parts + polys = gpd.read_file(poly_ds_path, layer=poly_layer, engine="pyogrio").explode(index_parts=False) + lines = gpd.read_file(lines_ds_path, engine="pyogrio") + + if polys.crs != lines.crs: + polys.to_crs(lines.crs, inplace=True) + + # 2. Intersect with lines and drop the join index + joined = polys.sjoin(lines, how="inner", predicate="intersects") + + # SPEEDUP 1: Drop duplicated polygons to prevent unary_union from choking on overlaps + unique_polys = joined[~joined.index.duplicated(keep='first')] + + # SPEEDUP 2: Isolate just the geometry column so pandas doesn't waste time aggregating attributes + unique_polys[['geometry']].dissolve() \ + .to_file(dst_filtered_poly_ds_path, engine="pyogrio", promote_to_multi=True) + + if dst_filtered_poly_ds_path.exists():return dst_filtered_poly_ds_path + + +async def extract_admin_boundaries(pbf_path: str, dst_dir: str, admin_level: int, progress=None) -> str: + """ + Extracts administrative boundaries via Osmium and exports to a GeoJSONSeq. + Maps standard ADM levels (0, 1, 2) to OSM admin_levels (2, 4, 6). + """ + dst_path = Path(dst_dir) + + # Map ADM 0, 1, 2 to proper OSM admin_levels (Country=2, Province/State=4, District/County=6) + osm_level = {0: 2, 1: 4, 2: 6}.get(admin_level, admin_level) + + filtered_pbf = dst_path / f"admin_{osm_level}_only.osm.pbf" + admin_geojsonseq = dst_path / f"admin_{osm_level}.geojsonseq" + + # Filter for relations and ways matching the specific admin level + tags_to_keep = [f"r/admin_level={osm_level}", f"w/admin_level={osm_level}"] + + if progress: + progress.console.print(f"[cyan]Filtering admin_level={osm_level} from OSM via Osmium...[/cyan]") + + # Step 1: Filter PBF down to specified administrative boundaries + await asyncio.to_thread( + run_cli, + ["osmium", "tags-filter", str(pbf_path)] + tags_to_keep + ["-o", str(filtered_pbf), "--overwrite"] + ) + + if progress: + progress.console.print("[cyan]Exporting admin boundaries to GeoJSONSeq...[/cyan]") + + # Step 2: Export to geojsonseq (Requires polygon geometry type for boundaries) + await asyncio.to_thread( + run_cli, + ["osmium", "export", str(filtered_pbf), "-o", str(admin_geojsonseq), "-f", "geojsonseq", + "--geometry-type=polygon", "--overwrite"] + ) + + if progress: + progress.console.print("[cyan]Loading admin boundaries completed...[/cyan]") + + return str(admin_geojsonseq) \ No newline at end of file