diff --git a/rapida/cli/assess.py b/rapida/cli/assess.py index 112a299..5216883 100644 --- a/rapida/cli/assess.py +++ b/rapida/cli/assess.py @@ -147,7 +147,7 @@ def build_variable_help(): if vars_list: vars_str = ", ".join(vars_list) parts.append(f"{comp} ({vars_str})") - if comp == 'population': print(len(vars_list)) + return "The variable/s to be assessed. Will be filtered by selected components. Available variables per component:\n\n" + "\n\n".join(parts) @@ -211,7 +211,7 @@ def assess(ctx, all=False, components=None, variables=None, year=None, datetime rapida assess -c landuse -dt 2025-02-01/2025-05-31 -cc 10: Search Sentinel 2 item which is less than 10% of cloud cover from February to May 2025. """ - + progress = ctx.obj.get('progress') if not is_rapida_initialized(): return @@ -231,7 +231,7 @@ def assess(ctx, all=False, components=None, variables=None, year=None, datetime sys.exit(0) logger.info(f'Current project/folder: {prj.path}') - with Progress(disable=False, console=None) as progress: + with progress: with Session() as session: all_components = session.get_components() target_components = components diff --git a/rapida/cli/connectivity.py b/rapida/cli/connectivity.py index e207702..36b4996 100644 --- a/rapida/cli/connectivity.py +++ b/rapida/cli/connectivity.py @@ -1,4 +1,4 @@ -from typing import Union +from typing import Union, Iterable import click import logging @@ -7,10 +7,21 @@ from rapida.connectivity import run_connectivity_analysis from rapida.cli.aclick import AsyncCommand from rapida.connectivity.isochrone import MODE_MAP - +from rapida.cli.assess import get_variables_by_components logger = logging.getLogger(__name__) + +def validate_variables(ctx, param, value): + """ + click callback function to validate polulation + """ + valid_vars = get_variables_by_components(['population'])['population'] + invalid = [v for v in value if v not in valid_vars] + if invalid: + raise click.BadParameter(f"Invalid variable{'s' if len(invalid) > 1 else ''}: {', '.join(invalid)} for population . Valid options: {', '.join(valid_vars)}") + return value + def parse_intervals(ctx, param, value): """Parses a comma-separated string of numbers into a list of integers.""" if not value: @@ -46,6 +57,20 @@ def parse_intervals(ctx, param, value): help="Comma-separated time intervals in minutes for the catchment areas." ) + +@click.option( + '-sd', '--sites-dataset', + type=click.Path(exists=True, file_okay=True, dir_okay=False, readable=True), + default=None, + help="Path to an OGR-supported vector data source (e.g., GPKG, Shapefile) containing sites." +) +@click.option( + '-sl', '--sites-layer', + type=str, + default="0", + help="Name or index of the layer to use from sites dataset. Defaults to the first layer (layer 0)." +) + @click.option( '-bd', '--barriers-dataset', type=click.Path(exists=True, file_okay=True, dir_okay=False, readable=True), @@ -59,12 +84,19 @@ def parse_intervals(ctx, param, value): help="Name or index of the layer to use from barriers dataset. Defaults to the first layer (layer 0)." ) + @click.option('-bb', "--barriers-buffer", type=int, default=5, required=False, help="The value in meters to used to buffer the geometries in barriers/dataset/layer in case the barriers are lines" ) + +@click.option('--popvar', required=False, multiple=True, + type=str, callback=validate_variables, + help=f"Open or more RAPIDA population variable to compute zonal stats for withing the connectivity zones" + ) + @click.option( "--dst-dir", "-d", # Short option @@ -84,12 +116,15 @@ def parse_intervals(ctx, param, value): @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 + 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 ): - logger.info(f'Running connectivity analysis ') + logger.info(f'Running connectivity analysis') progress = ctx.obj.get('progress') with progress: 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, progress=progress + 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 ) \ No newline at end of file diff --git a/rapida/connectivity/__init__.py b/rapida/connectivity/__init__.py index 48bb3e0..1ba6886 100644 --- a/rapida/connectivity/__init__.py +++ b/rapida/connectivity/__init__.py @@ -2,27 +2,66 @@ import os.path from rapida.util.bbox_param_type import get_best_semantic_label from rich.progress import Progress -from rapida.connectivity.io import prepare_osm_pbf,extract_health_sites, extract_origins_from_geojson +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 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, progress:Progress=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 ): 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) - health_sites = await extract_health_sites(pbf_path=bbox_pbf, dst_dir=dest_dir, progress=progress) + if sites_dataset is None: + sites = await extract_health_sites(pbf_path=bbox_pbf, dst_dir=dest_dir, progress=progress) + else: + sites = sites_dataset + dag_tar_path = await compile_valhalla_graph(pbf_path=bbox_pbf,dst_dir=dest_dir, progress=progress) - origins = extract_origins_from_geojson(geojson_path=health_sites) + origins = extract_origins(sites_dataset=sites, src_layer=sites_layer) + + results = await connectivity_areas( + tar_path=dag_tar_path, origins=origins, travel_mode=travel_mode, intervals_minutes=time_intervals) + isochrones_path = os.path.join(dest_dir, 'isochrones.geojson') + with open(isochrones_path, "w") as f: + json.dump(results, f, indent=2) + + # 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 + # ) + + + if barriers_dataset is not None: + 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 - ) - with open(os.path.join(dest_dir, 'isochrones.geojson'), "w") as f: - json.dump(results, f, indent=2) + ) + with open(os.path.join(dest_dir, 'isochrones_with_barriers.geojson'), "w") as f: + json.dump(barrier_results, f, indent=2) + return \ No newline at end of file diff --git a/rapida/connectivity/io.py b/rapida/connectivity/io.py index e391689..e27504a 100644 --- a/rapida/connectivity/io.py +++ b/rapida/connectivity/io.py @@ -11,7 +11,7 @@ from rich.progress import Progress import json from shapely.geometry import shape, mapping -from osgeo import gdal +from osgeo import gdal, ogr, osr from shapely.wkb import loads as load_wkb from shapely.ops import orient, transform import numpy as np @@ -58,6 +58,7 @@ async def prepare_osm_pbf(bbox: tuple[float, float, float, float], dst_dir: str # Perform quick spatial overlay check gdf = gpd.GeoDataFrame.from_features(features, crs="EPSG:4326") + intersecting = gdf[gdf.intersects(bbox_geom)] if intersecting.empty: @@ -214,6 +215,67 @@ def extract_origins_from_geojson(geojson_path: str) -> list[tuple[float, float]] return origins +def extract_origins(sites_dataset: str=None, src_layer: str = None) -> list[tuple[float, float]]: + """ + Extracts a list of (longitude, latitude) tuples from a spatial file using OGR. + Handles reprojection to WGS84 (EPSG:4326) if the source is not in lat/lon. + """ + # Open the dataset + with gdal.OpenEx(sites_dataset, gdal.OF_VECTOR) as src_ds: + + try: + layer = src_ds.GetLayer(int(src_layer)) + except ValueError: + layer = src_ds.GetLayerByName(str(src_layer)) + + if layer is None: + raise ValueError(f"Layer '{src_layer}' could not be found in the dataset {sites_dataset}.") + + + + # Set up coordinate transformation to WGS84 (EPSG:4326) + source_srs = layer.GetSpatialRef() + target_srs = osr.SpatialReference() + target_srs.ImportFromEPSG(4326) + + # OGR 3+ strict axis mapping strategy (ensures Longitude/Latitude order) + target_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) + + transform = None + if source_srs and not source_srs.IsSame(target_srs): + transform = osr.CoordinateTransformation(source_srs, target_srs) + + origins = [] + layer.ResetReading() + # Iterate through features + for feature in layer: + # If you still need the filter: if feature.GetField("osm_id") != 80: continue + + geom = feature.GetGeometryRef() + if geom is not None: + # Clone geometry to avoid modifying original layer data during transform + geom_clone = geom.Clone() + + # Reproject if necessary + if transform: + geom_clone.Transform(transform) + + # Helper function to recursively extract points from nested collections + def extract_points(g): + name = g.GetGeometryName() + if name == "POINT": + origins.append((g.GetX(), g.GetY())) + elif name in ("MULTIPOINT", "GEOMETRYCOLLECTION"): + for i in range(g.GetGeometryCount()): + sub_geom = g.GetGeometryRef(i) + if sub_geom is not None: + extract_points(sub_geom) + + extract_points(geom_clone) + + return origins + + def read_barriers_grid(src_path: str, src_layer: str = None, barriers_buffer:float=None) -> list: """Reads a vector source and cuts features into micro-tiles to stay under Valhalla's limit.""" if not src_path: diff --git a/rapida/project/project.py b/rapida/project/project.py index 97e378b..721e07a 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 from rapida import constants from rapida.az.blobstorage import check_blob_exists, delete_blob from rapida.session import Session @@ -159,12 +159,17 @@ def __init__(self, path: str,polygons: str = None, iso3_codes = [] for lat, lon in zip(lats, lons): try: - iso3_codes.append(fetch_ccode(lat=lat, lon=lon)) + 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", diff --git a/rapida/stats/raster_zonal_stats.py b/rapida/stats/raster_zonal_stats.py index b60fc7f..51ac901 100644 --- a/rapida/stats/raster_zonal_stats.py +++ b/rapida/stats/raster_zonal_stats.py @@ -257,4 +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) return combined \ No newline at end of file