Harden priority geospatial compatibility paths

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

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()

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

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,
     )
 

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,13 +30,18 @@
     "get_hdx_dataset",
     "hdx_fb_population",
     "hdx_dataset",
+    "HDXProviderError",
 ]
 
 hdx_config = Configuration.create(
     hdx_site="prod", user_agent="urbanpy", hdx_read_only=True
 )
 
 
+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

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)