Merge pull request #1101 from gboeing/feature

warning and optimization

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Unreleased
 
 - fix a bug arising from the save_graph_xml function (#1093)
+- warn user if their query area is significantly larger than max query area size (#1101)
 - refactor utils_geo module and deprecate quadrat_width and min_num function arguments (#1100)
 - under-the-hood code clean-up (#1092 #1099)
 

osmnx/_overpass.py

@@ -202,31 +202,25 @@ def _make_overpass_polygon_coord_strs(polygon):
 
     Returns
     -------
-    polygon_coord_strs : list
-        list of exterior coordinate strings for smaller sub-divided polygons
+    coord_strs : list
+        list of strings of exterior coordinates of polygon(s)
     """
-    geometry_proj, crs_proj = projection.project_geometry(polygon)
-    gpcs = utils_geo._consolidate_subdivide_geometry(geometry_proj)
-    geometry, _ = projection.project_geometry(gpcs, crs=crs_proj, to_latlong=True)
-
-    # extract geometry's exterior coords
-    polygons_coords = []
-    for polygon in geometry.geoms:
-        x, y = polygon.exterior.xy
-        polygons_coords.append(list(zip(x, y)))
-
-    # convert exterior coords to the string format the API expects
-    polygon_coord_strs = []
-    for coords in polygons_coords:
-        s = ""
-        separator = " "
-        for coord in list(coords):
-            # round floating point lats and longs to 6 decimals (ie, ~100 mm)
-            # so we can hash and cache strings consistently
-            s = f"{s}{separator}{coord[1]:.6f}{separator}{coord[0]:.6f}"
-        polygon_coord_strs.append(s.strip(separator))
-
-    return polygon_coord_strs
+    # first subdivide the polygon if its area exceeds max size
+    # this results in a multipolygon of 1+ constituent polygons
+    poly_proj, crs_proj = projection.project_geometry(polygon)
+    multi_poly_proj = utils_geo._consolidate_subdivide_geometry(poly_proj)
+    multi_poly, _ = projection.project_geometry(multi_poly_proj, crs=crs_proj, to_latlong=True)
+
+    # then extract each's exterior coords to the string format Overpass
+    # expects, rounding lats and lons to 6 decimals (ie, ~100 mm) so we
+    # can hash and cache URL strings consistently
+    coord_strs = []
+    for geom in multi_poly.geoms:
+        x, y = geom.exterior.xy
+        coord_list = [f'{xy[1]:.6f}{" "}{xy[0]:.6f}' for xy in zip(x, y)]
+        coord_strs.append(" ".join(coord_list))
+
+    return coord_strs
 
 
 def _create_overpass_query(polygon_coord_str, tags):

osmnx/utils_geo.py

@@ -274,6 +274,17 @@ def _consolidate_subdivide_geometry(geometry):
     ):
         geometry = geometry.convex_hull
 
+    # warn user if they passed a geometry with area much larger than max size
+    ratio = int(geometry.area / mqas)
+    warning_threshold = 10
+    if ratio > warning_threshold:
+        msg = (
+            f"This area is {ratio:,} times your configured Overpass max query "
+            "area size. It will automatically be divided up into multiple "
+            "sub-queries accordingly. This may take a long time."
+        )
+        warn(msg, stacklevel=2)
+
     # if geometry area exceeds max size, subdivide it into smaller subpolygons
     # that are no greater than settings.max_query_area_size in size
     if geometry.area > mqas:

tests/test_osmnx.py

@@ -408,7 +408,7 @@ def test_find_nearest():
     ne2 = ox.distance.nearest_edges(G, X[0], Y[0], interpolate=50, return_dist=True)
 
 
-def test_endpoints():
+def test_api_endpoints():
     """Test different API endpoints."""
     default_timeout = ox.settings.timeout
     default_key = ox.settings.nominatim_key
@@ -563,6 +563,11 @@ def test_graph_save_load():
 
 def test_graph_from_functions():
     """Test downloading graphs from Overpass."""
+    # test subdividing a large geometry (raises a UserWarning)
+    bbox = ox.utils_geo.bbox_from_point((0, 0), dist=1e5, project_utm=True)
+    poly = ox.utils_geo.bbox_to_poly(*bbox)
+    _ = ox.utils_geo._consolidate_subdivide_geometry(poly)
+
     # graph from bounding box
     _ = ox.utils_geo.bbox_from_point(location_point, project_utm=True, return_crs=True)
     north, south, east, west = ox.utils_geo.bbox_from_point(location_point, dist=500)