From 3ac0017240e3afeb78253adea5d2012d6a9a2549 Mon Sep 17 00:00:00 2001 From: Claudio Ortega Date: Sun, 9 Aug 2026 21:45:26 -0700 Subject: [PATCH] fix: harden priority geospatial compatibility paths --- tests/compatibility_test.py | 68 ++++++++++++++++++++++++++ tests/download_test.py | 35 +++++++++++-- tests/geom_test.py | 19 +++++++ urbanpy/accessibility/accessibility.py | 51 +++++++++++-------- urbanpy/download/download.py | 44 ++++++++++++----- urbanpy/geom/geom.py | 6 ++- urbanpy/routing/routing.py | 14 ++++-- 7 files changed, 195 insertions(+), 42 deletions(-) create mode 100644 tests/compatibility_test.py diff --git a/tests/compatibility_test.py b/tests/compatibility_test.py new file mode 100644 index 0000000..5417f7a --- /dev/null +++ b/tests/compatibility_test.py @@ -0,0 +1,68 @@ +from importlib.metadata import version + +import geopandas as gpd +import networkx as nx +from shapely.geometry import Point, Polygon + +from urbanpy import accessibility, geom, routing + + +def test_h3_v4_public_api_generates_valid_crs_aware_hexagons(): + assert int(version("h3").split(".", 1)[0]) >= 4 + city = gpd.GeoDataFrame( + geometry=[ + Polygon( + [ + (-77.05, -12.10), + (-77.00, -12.10), + (-77.00, -12.05), + (-77.05, -12.05), + ] + ) + ], + crs="EPSG:4326", + ) + + result = geom.gen_hexagons(8, city) + + assert not result.empty + assert result.crs == city.crs + assert result["hex"].is_unique + + +def test_osmnx_v2_nearest_nodes_uses_longitude_then_latitude(monkeypatch): + graph = nx.MultiDiGraph() + graph.graph["crs"] = "EPSG:4326" + graph.add_node(1, x=-77.04, y=-12.08) + graph.add_node(2, x=-77.03, y=-12.07) + graph.add_edge(1, 2, length=100.0) + nearest = [] + + def fake_nearest(_graph, x, y): + nearest.append((x, y)) + return 1 + + monkeypatch.setattr(routing.routing.ox, "nearest_nodes", fake_nearest) + + result = routing.isochrone_from_graph(graph, [(-77.04, -12.08)], [5], "walking") + + assert nearest == [(-77.04, -12.08)] + assert result.crs.to_string() == "EPSG:4326" + assert result.geometry.notna().all() + + +def test_pressure_map_preserves_crs_for_already_projected_inputs(): + blocks = gpd.GeoDataFrame( + {"demand": [10]}, + geometry=[Polygon([(0, 0), (100, 0), (100, 100), (0, 0)])], + crs="EPSG:3857", + ) + pois = gpd.GeoDataFrame( + geometry=[Point(25, 25)], + crs="EPSG:3857", + ) + + result = accessibility.pressure_map(blocks, pois, "demand", buffer_size=200) + + assert result.crs == blocks.crs + assert result.loc[0, "ds"] == 5 diff --git a/tests/download_test.py b/tests/download_test.py index 2f5b679..0c05078 100644 --- a/tests/download_test.py +++ b/tests/download_test.py @@ -36,9 +36,7 @@ def test_nominatim_parses_captured_geojson(): ], } with responses.RequestsMock() as captured: - captured.get( - "https://nominatim.openstreetmap.org/search.php", json=payload - ) + captured.get("https://nominatim.openstreetmap.org/search.php", json=payload) result = download.nominatim_osm("Lima, Peru", email="dev@example.org") assert list(result["display_name"]) == ["Lima, Peru"] @@ -63,7 +61,9 @@ def test_search_hdx_dataset_uses_captured_provider_records(monkeypatch): }, ] monkeypatch.setattr( - download_module.Dataset, "search_in_hdx", lambda _query: [object()] + download_module.Dataset, + "search_in_hdx", + lambda _query, **_kwargs: [object()], ) monkeypatch.setattr( download_module.Dataset, "get_all_resources", lambda _datasets: records @@ -81,6 +81,33 @@ def test_search_hdx_dataset_uses_captured_provider_records(monkeypatch): } +def test_hdx_single_resource_list_is_a_valid_download_selection(monkeypatch): + resources = download_module.pd.DataFrame( + {"url": ["https://example.org/population.csv"]}, index=[7] + ) + expected = download_module.pd.DataFrame({"population": [10]}) + read_csv = monkeypatch.setattr( + download_module.pd, "read_csv", lambda _url: expected.copy() + ) + + result = download.get_hdx_dataset(resources, [7]) + + assert result.equals(expected) + assert read_csv is None + + +def test_hdx_provider_errors_are_stable_and_hide_provider_details(monkeypatch): + def fail(*_args, **_kwargs): + raise download_module.HDXError("provider response containing internals") + + monkeypatch.setattr(download_module.Dataset, "search_in_hdx", fail) + + with pytest.raises(download.HDXProviderError) as captured: + download.search_hdx_dataset("Peru") + + assert "provider response" not in str(captured.value) + + def test_osmnx_graph_validates_required_arguments(capsys): assert download.osmnx_graph("polygon") is None assert "provide a polygon" in capsys.readouterr().out.lower() diff --git a/tests/geom_test.py b/tests/geom_test.py index d336f1f..3043498 100644 --- a/tests/geom_test.py +++ b/tests/geom_test.py @@ -62,3 +62,22 @@ def test_hexagon_generation_merge_and_downsampling(): assert merged["population"].sum() == 10 assert coarse["population"].sum() == 10 assert coarse.crs == hexagons.crs + + +def test_merge_shape_hex_can_be_rerun_on_its_own_result(): + hexagons = geom.gen_hexagons(8, _city()) + points = gpd.GeoDataFrame( + {"population": [3, 7]}, + geometry=[Point(-77.03, -12.08), Point(-77.01, -12.06)], + crs="EPSG:4326", + ) + + first = geom.merge_shape_hex( + hexagons, points, {"population": "sum"}, predicate="within" + ) + second = geom.merge_shape_hex( + first, points, {"population": "sum"}, predicate="within" + ) + + assert second["population"].equals(first["population"]) + assert second.crs == first.crs diff --git a/urbanpy/accessibility/accessibility.py b/urbanpy/accessibility/accessibility.py index 8175842..950b5e2 100644 --- a/urbanpy/accessibility/accessibility.py +++ b/urbanpy/accessibility/accessibility.py @@ -17,9 +17,7 @@ # Gaussian friction function for distance decay. Accepts scalar or array input. def friction(dm, d0): dm = np.asarray(dm, dtype=float) - return np.where( - dm > d0, 0.0, np.exp(-0.5 * (dm / d0) ** 2) / (1.0 - np.exp(-0.5)) - ) + return np.where(dm > d0, 0.0, np.exp(-0.5 * (dm / d0) ** 2) / (1.0 - np.exp(-0.5))) def hu_access_map(units, pois, population_column, weight=1, d0=1250): @@ -85,12 +83,14 @@ def hu_access_map(units, pois, population_column, weight=1, d0=1250): units["centroid"] = units.geometry.centroid # Create buffer GeoDataFrame - buffers_poi = gpd.GeoDataFrame(pois["idx"], geometry=pois["buffer"]) - buffers_units = gpd.GeoDataFrame(units["idx"], geometry=units["buffer"]) + buffers_poi = gpd.GeoDataFrame(pois["idx"], geometry=pois["buffer"], crs=pois.crs) + buffers_units = gpd.GeoDataFrame( + units["idx"], geometry=units["buffer"], crs=units.crs + ) # Compute catchment area (Rj) for each poi join = gpd.sjoin( - buffers_poi, units[["idx", "geometry"]], op="intersects", how="left" + buffers_poi, units[["idx", "geometry"]], predicate="intersects", how="left" ) join = join.rename(columns={"idx_left": "idx_poi", "idx_right": "idx_unit"}) merge = pd.merge( @@ -126,7 +126,7 @@ def hu_access_map(units, pois, population_column, weight=1, d0=1250): # Compute block accesibility join = gpd.sjoin( - buffers_units, pois[["idx", "geometry"]], op="intersects", how="left" + buffers_units, pois[["idx", "geometry"]], predicate="intersects", how="left" ) join = join.rename(columns={"idx_left": "idx_unit", "idx_right": "idx_poi"}) @@ -143,8 +143,12 @@ def hu_access_map(units, pois, population_column, weight=1, d0=1250): # Compute friction (vectorized over centroid vs POI geometry) ctr_x = np.fromiter((p.x for p in merge["centroid"]), dtype=float, count=len(merge)) ctr_y = np.fromiter((p.y for p in merge["centroid"]), dtype=float, count=len(merge)) - poi_x = np.fromiter((p.x for p in merge["geometry_y"]), dtype=float, count=len(merge)) - poi_y = np.fromiter((p.y for p in merge["geometry_y"]), dtype=float, count=len(merge)) + poi_x = np.fromiter( + (p.x for p in merge["geometry_y"]), dtype=float, count=len(merge) + ) + poi_y = np.fromiter( + (p.y for p in merge["geometry_y"]), dtype=float, count=len(merge) + ) merge["friction"] = friction(np.hypot(ctr_x - poi_x, ctr_y - poi_y), d0) # Compute accesibility Ai @@ -171,7 +175,7 @@ def hu_access_map(units, pois, population_column, weight=1, d0=1250): del df_ai["idx"] - access_map = gpd.GeoDataFrame(df_ai, geometry=df_ai["geometry"]) + access_map = gpd.GeoDataFrame(df_ai, geometry=df_ai["geometry"], crs=units.crs) return access_map @@ -218,28 +222,30 @@ def pressure_map(blocks, pois, demand_column, operation="intersects", buffer_siz ---------- Van Eck, J. R., & de Jong, T. (1999). Accessibility analysis and spatial competition effects in the context of GIS-supported service location planning. Computers, environment and urban systems, 23(2), 75-89. """ - if not pois.crs.is_projected: - pois_proj = project_gdf(pois) - - if not blocks.crs.is_projected: - blocks_proj = project_gdf(blocks) + pois_proj = pois if pois.crs.is_projected else project_gdf(pois) + blocks_proj = blocks if blocks.crs.is_projected else project_gdf(blocks) idx_blocks = [f"block_{i}" for i in blocks.index] blocks_proj["idx"] = idx_blocks buffers = gpd.GeoDataFrame( - idx_blocks, columns=["idx"], geometry=blocks_proj.geometry.buffer(buffer_size) + idx_blocks, + columns=["idx"], + geometry=blocks_proj.geometry.buffer(buffer_size), + crs=blocks_proj.crs, ) - merge = gpd.sjoin(buffers, pois_proj, op=operation) + merge = gpd.sjoin(buffers, pois_proj, predicate=operation) nj = merge.groupby("idx").count()["index_right"] nj.name = "nj" nj = nj.reset_index() - blocks = pd.merge(blocks, nj, how="left") - blocks["ds"] = blocks[demand_column] / (blocks_proj["nj"] + 1) + result = blocks.copy() + result["idx"] = idx_blocks + result = result.merge(nj, how="left", on="idx") + result["ds"] = result[demand_column] / (result["nj"].fillna(0) + 1) - return blocks + return gpd.GeoDataFrame(result, geometry="geometry", crs=blocks.crs) def travel_times(inputs, pois, col_label="poi", nearest_neighbor_dist="haversine"): @@ -277,7 +283,10 @@ def travel_times(inputs, pois, col_label="poi", nearest_neighbor_dist="haversine centroid_list = list(centroids) nearest_geoms = pois.geometry.iloc[ixs.ravel()].reset_index(drop=True).tolist() distance_duration = pd.DataFrame( - [osrm_route(origin=o, destination=d) for o, d in zip(centroid_list, nearest_geoms)], + [ + osrm_route(origin=o, destination=d) + for o, d in zip(centroid_list, nearest_geoms) + ], index=gdf.index, ) diff --git a/urbanpy/download/download.py b/urbanpy/download/download.py index 49015fb..3f5ad95 100644 --- a/urbanpy/download/download.py +++ b/urbanpy/download/download.py @@ -9,6 +9,7 @@ from geopandas import GeoDataFrame, GeoSeries from hdx.api.configuration import Configuration from hdx.data.dataset import Dataset +from hdx.data.hdxobject import HDXError from pandas import DataFrame from shapely.geometry import MultiPolygon, Polygon, Point @@ -18,6 +19,7 @@ overpass_to_gdf, to_overpass_query, ) +from urbanpy.errors import UrbanPyError __all__ = [ "nominatim_osm", @@ -28,6 +30,7 @@ "get_hdx_dataset", "hdx_fb_population", "hdx_dataset", + "HDXProviderError", ] hdx_config = Configuration.create( @@ -35,6 +38,10 @@ ) +class HDXProviderError(UrbanPyError): + """HDX could not complete a dataset search or returned invalid metadata.""" + + def nominatim_osm( query: str, expected_position: "int | None" = 0, email: str = "" ) -> GeoDataFrame: @@ -323,12 +330,21 @@ def search_hdx_dataset( 12 | 2019-06-11 | PER_youth_15_24_2019-06-01_csv.zip | Youth (ages 15-24) | 16.61 | https://data.humdata.org/dataset/4e74db39-87f1... """ # Get dataset list - datasets = Dataset.search_in_hdx(f"title:{country.lower()}-{repository}") + try: + datasets = Dataset.search_in_hdx( + f"title:{country.lower()}-{repository}", rows=100 + ) + resources_records = Dataset.get_all_resources(datasets) + except HDXError as error: + raise HDXProviderError( + "HDX dataset search failed; retry later or run the live provider check." + ) from error - resources_records = Dataset.get_all_resources(datasets) resources_df = pd.DataFrame.from_records(resources_records) if resources_df.shape[0] == 0: - print("No datasets found") + return pd.DataFrame( + columns=["created", "name", "population", "size_mb", "url"] + ).rename_axis("id") else: resources_csv_df = resources_df[ @@ -387,15 +403,16 @@ def get_hdx_dataset( """ - urls = resources_df.loc[ids, "url"] - - print(urls) - if isinstance(ids, list) and len(ids) > 1: - df = pd.concat([pd.read_csv(url) for url in urls]) + selected = resources_df.loc[ids, "url"] + if isinstance(selected, pd.Series): + urls = selected.tolist() + if not urls: + raise ValueError("ids did not select any HDX resources") + df = pd.concat([pd.read_csv(url) for url in urls], ignore_index=True) else: - df = pd.read_csv(urls) + df = pd.read_csv(selected) - if mask: + if mask is not None: if isinstance(mask, GeoDataFrame): mask = mask.unary_union minx, miny, maxx, maxy = mask.bounds @@ -455,6 +472,11 @@ def hdx_fb_population(country, map_type): resources_df["population"] == HDX_POPULATION_TYPES[map_type] ].index.tolist() + if not dataset_ix: + raise ValueError( + f"No {map_type!r} population resource is available for {country!r}." + ) + population = get_hdx_dataset(resources_df, dataset_ix) return population @@ -507,5 +529,5 @@ def hdx_dataset(resource): hdx_url = resource dataset = pd.read_csv(hdx_url) - + return dataset diff --git a/urbanpy/geom/geom.py b/urbanpy/geom/geom.py index a3a6d5a..f70bd3c 100644 --- a/urbanpy/geom/geom.py +++ b/urbanpy/geom/geom.py @@ -233,7 +233,11 @@ def merge_shape_hex( 888e628debfffff | POLYGON ((-76.67982 -12.18998, -76.68413 -12.1... | NaN 888e6299b3fffff | POLYGON ((-76.78876 -11.97286, -76.79307 -11.9... | 3225.658803 """ - joined = gpd.sjoin(shape, hexs, how=how, predicate=predicate) + # A previous result contains the aggregate columns. Exclude them from the + # right-hand join so a second call cannot create `_left`/`_right` suffixes + # and make the requested aggregation columns disappear. + join_hexs = hexs.drop(columns=list(agg), errors="ignore") + joined = gpd.sjoin(shape, join_hexs, how=how, predicate=predicate) # Uses index right based on the order of points and hex. Right takes hex index hex_merge = joined.groupby("index_right").agg(agg) diff --git a/urbanpy/routing/routing.py b/urbanpy/routing/routing.py index 7a61fc2..f177f81 100644 --- a/urbanpy/routing/routing.py +++ b/urbanpy/routing/routing.py @@ -698,7 +698,9 @@ def isochrone_from_graph(graph, locations, time_range, profile): travel_speed = profiles[profile] center_nodes = [ox.nearest_nodes(graph, x, y) for x, y in locations] - G = ox.project_graph(graph) + # Edge lengths are already expressed in meters by OSMnx. Keep the graph in + # its input CRS so returned polygons and CRS metadata remain consistent. + G = graph.copy() meters_per_minute = travel_speed * 1000 / 60 # km per hour to m per minute for u, v, k, edata in G.edges(data=True, keys=True): @@ -716,13 +718,15 @@ def isochrone_from_graph(graph, locations, time_range, profile): for trip_time in sorted_times: reachable_nodes = [n for n, d in dist.items() if d <= trip_time] node_points = [ - Point((node_data[n]["lon"], node_data[n]["lat"])) - for n in reachable_nodes + Point((node_data[n]["x"], node_data[n]["y"])) for n in reachable_nodes ] bounding_poly = MultiPoint(node_points).convex_hull if node_points else None rows.append([ix, trip_time, bounding_poly]) - isochrones = gpd.GeoDataFrame(rows, columns=["group_index", "contour", "geometry"]) - isochrones.crs = "EPSG:4326" + isochrones = gpd.GeoDataFrame( + rows, + columns=["group_index", "contour", "geometry"], + crs=graph.graph.get("crs", "EPSG:4326"), + ) return isochrones