Границы регионов без упрощения, чтобы Куркино не считалось по Екатеринбургу (#3527)

tradein-mvp/backend/app/services/regions.py

@@ -23,8 +23,10 @@ pathlib — стандартная библиотека, листовость н
 from __future__ import annotations
 
 import json
+import math
 from dataclasses import dataclass
 from pathlib import Path
+from typing import NamedTuple
 
 # bbox = (lat_min, lat_max, lon_min, lon_max) — тот же порядок, что исторический
 # geocoder.EKB_BBOX_TIGHT (см. is_within_bbox ниже).
@@ -105,6 +107,17 @@ def is_within_bbox(lat: float, lon: float, bbox: BBox) -> bool:
 # Парсим один раз на импорте модуля — 27 точек лишний JSON-parse на каждый
 # вызов region_for_point был бы дороже самого ray casting.
 #
+# ТОЧНОСТЬ (правка регрессии первой версии #3052). Полигоны берутся БЕЗ
+# упрощения (`polygon_threshold=0`, 42 755 вершин), а не с `0.005` (~500 м,
+# 2 388 вершин), как в первой версии. Полукилометровое упрощение спрямляет
+# ровно то, ради чего полигоны и заводились — изрезанную границу Москвы с
+# областью. Сверка резолва против сохранённого region_code по проду (277 582
+# объявления + 48 143 дома) на упрощённой геометрии дала: 69 точек НИ в одном
+# полигоне (агломерация → None → вызывающий код `trade_in.py` подставляет
+# DEFAULT_REGION_CODE, и адрес в Куркино считается по коэффициентам
+# Екатеринбурга), Митино/Косино-Ухтомский/Можайский вне Москвы (27 домов),
+# 731 объявление области внутри Москвы. Без упрощения этих спрямлений нет.
+#
 # Ring — одно кольцо GeoJSON (lon, lat) точек: внешний контур ИЛИ дыра.
 # Полигон/мультиполигон региона хранится как ПЛОСКИЙ список всех его колец
 # (для MultiPolygon 77 — кольца всех 10 частей вперемешку, у Polygon 50 —
@@ -148,14 +161,93 @@ def _load_region_boundaries() -> dict[int, list[Ring]]:
 _REGION_BOUNDARIES: dict[int, list[Ring]] = _load_region_boundaries()
 
 
-def _point_in_ring(x: float, y: float, ring: Ring) -> bool:
+class _IndexedRing(NamedTuple):
-    """PNPOLY (W. R. Franklin) ray casting для одного кольца."""
+    """Кольцо + предпосчёт, которым точная геометрия остаётся дешёвой.
+
+    Наивный ray casting линеен по числу вершин, а вершин после отказа от
+    упрощения в 18 раз больше (42 755 против 2 388) — прямым перебором это
+    десятки-сотни микросекунд на вызов. Два предпосчёта снимают почти всё:
+
+    lat_min/lat_max/lon_min/lon_max — bbox САМОГО кольца (не региона): у
+        Москвы мультиполигон из 10 частей, у области внешний контур + 9 дыр,
+        и точка почти всегда вне bbox большинства колец — такое кольцо
+        отвергается четырьмя сравнениями вместо обхода вершин.
+    band_h/bands — кольцо нарезано на горизонтальные полосы по широте; в
+        полосе лежат ТОЛЬКО те рёбра, чей диапазон широт её задевает (ребро
+        попадает во все свои полосы, поэтому пересечения не теряются).
+        Луч ray casting горизонтален, значит пересекать его могут лишь рёбра
+        своей полосы — обход сокращается с «все вершины кольца» до «рёбра
+        одной полосы» (для внешнего контура области — 23 511 → максимум 724,
+        типично десятки).
+    points — те же вершины ОДНИМ объектом с _REGION_BOUNDARIES (не копия):
+        нужны fallback'у по расстоянию до границы в region_for_point.
+    """
+
+    lat_min: float
+    lat_max: float
+    lon_min: float
+    lon_max: float
+    band_h: float
+    bands: tuple[tuple[tuple[float, float, float, float], ...], ...]
+    points: Ring
+
+
+# Целевая плотность полосы: ~4 вершины кольца на полосу. Мельче — полосы
+# начинают стоить памяти и промахов кэша больше, чем экономят на обходе;
+# крупнее — в полосу попадает слишком много рёбер. Потолок числа полос
+# ограничивает память на самых длинных кольцах (внешний контур области).
+_RING_BAND_DENSITY = 4
+_RING_MAX_BANDS = 1024
+
+
+def _index_ring(points: Ring) -> _IndexedRing:
+    """Посчитать bbox и полосы кольца (один раз на импорте модуля)."""
+    lat_min = min(p[1] for p in points)
+    lat_max = max(p[1] for p in points)
+    lat_span = lat_max - lat_min
+    n_bands = max(1, min(_RING_MAX_BANDS, len(points) // _RING_BAND_DENSITY))
+    # Вырожденное по широте кольцо (lat_span == 0) в полосы не режется:
+    # одна полоса, band_h произвольный ненулевой — делить на 0 нельзя.
+    band_h = lat_span / n_bands if lat_span > 0.0 else 1.0
+    buckets: list[list[tuple[float, float, float, float]]] = [[] for _ in range(n_bands)]
+    x_prev, y_prev = points[-1]
+    for x, y in points:
+        edge = (x, y, x_prev, y_prev)
+        lo = int((min(y, y_prev) - lat_min) / band_h)
+        hi = int((max(y, y_prev) - lat_min) / band_h)
+        lo = 0 if lo < 0 else min(lo, n_bands - 1)
+        hi = 0 if hi < 0 else min(hi, n_bands - 1)
+        for band in range(lo, hi + 1):
+            buckets[band].append(edge)
+        x_prev, y_prev = x, y
+    return _IndexedRing(
+        lat_min=lat_min,
+        lat_max=lat_max,
+        lon_min=min(p[0] for p in points),
+        lon_max=max(p[0] for p in points),
+        band_h=band_h,
+        bands=tuple(tuple(b) for b in buckets),
+        points=points,
+    )
+
+
+_REGION_RING_INDEX: dict[int, tuple[_IndexedRing, ...]] = {
+    code: tuple(_index_ring(ring) for ring in rings) for code, rings in _REGION_BOUNDARIES.items()
+}
+
+
+def _point_in_ring(ring: _IndexedRing, lon: float, lat: float) -> bool:
+    """PNPOLY (W. R. Franklin) ray casting по рёбрам ОДНОЙ полосы кольца."""
+    lat_min, lat_max, lon_min, lon_max, band_h, bands, _points = ring
+    if lat < lat_min or lat > lat_max or lon < lon_min or lon > lon_max:
+        return False
+    band = int((lat - lat_min) / band_h)
+    if band >= len(bands):  # lat == lat_max попадает ровно на верхнюю кромку
+        band = len(bands) - 1
     inside = False
-    xj, yj = ring[-1]
+    for xi, yi, xj, yj in bands[band]:
-    for xi, yi in ring:
+        if (yi > lat) != (yj > lat) and lon < (xj - xi) * (lat - yi) / (yj - yi) + xi:
-        if (yi > y) != (yj > y) and x < (xj - xi) * (y - yi) / (yj - yi) + xi:
             inside = not inside
-        xj, yj = xi, yi
     return inside
 
 
@@ -166,12 +258,55 @@ def _point_in_region_polygon(code: int, lat: float, lon: float) -> bool:
     Региона без загруженной границы — граничит только по bbox (см. вызов в
     region_for_point), сюда такой код не попадает."""
     inside = False
-    for ring in _REGION_BOUNDARIES[code]:
+    for ring in _REGION_RING_INDEX[code]:
-        if _point_in_ring(lon, lat, ring):
+        if _point_in_ring(ring, lon, lat):
             inside = not inside
     return inside
 
 
+def _distance_to_boundary(code: int, lat: float, lon: float) -> float:
+    """Расстояние от точки до ближайшего ребра границы региона `code`.
+
+    В градусах широты (1° ≈ 111 км) по ПЛОСКОМУ приближению с поправкой
+    cos(lat) на сжатие долготы: значение сравнивается только с таким же
+    значением другого региона-кандидата на тех же ~десятках метров, где
+    геодезия и плоскость расходятся на доли процента.
+
+    Путь редкий (только щель между полигонами соседних регионов, см.
+    region_for_point), поэтому обход честный — по всем рёбрам, без полос.
+    Стоимость срезается оценкой снизу по bbox кольца: кольцо, чей bbox уже
+    дальше найденного минимума, не может дать ребра ближе.
+    """
+    lat_scale = math.cos(math.radians(lat))
+    px = lon * lat_scale
+    best = math.inf
+    ranked: list[tuple[float, _IndexedRing]] = []
+    for ring in _REGION_RING_INDEX[code]:
+        dx = max(ring.lon_min - lon, 0.0, lon - ring.lon_max) * lat_scale
+        dy = max(ring.lat_min - lat, 0.0, lat - ring.lat_max)
+        ranked.append((dx * dx + dy * dy, ring))
+    for lower_bound, ring in sorted(ranked, key=lambda item: item[0]):
+        if lower_bound >= best:
+            break
+        ax, ay = ring.points[-1]
+        ax *= lat_scale
+        for bx_raw, by in ring.points:
+            bx = bx_raw * lat_scale
+            dx, dy = bx - ax, by - ay
+            if dx or dy:
+                t = ((px - ax) * dx + (lat - ay) * dy) / (dx * dx + dy * dy)
+                t = 0.0 if t < 0.0 else (1.0 if t > 1.0 else t)
+                ex = ax + t * dx - px
+                ey = ay + t * dy - lat
+            else:  # выродившееся в точку ребро (дубль вершины в OSM)
+                ex, ey = ax - px, ay - lat
+            dist = ex * ex + ey * ey
+            if dist < best:
+                best = dist
+            ax, ay = bx, by
+    return math.sqrt(best)
+
+
 # Тиры обогащения (строковые ключи — по label'ам _with_budget в estimator).
 TIER_AVITO_IMV = "avito_imv"
 TIER_YANDEX_VALUATION = "yandex_valuation"
@@ -382,18 +517,38 @@ def region_for_point(lat: float, lon: float) -> Region | None:
          прямоугольнике области, но не в её реальных границах — то есть,
          например, внутри вырезанной дыры Москвы), идём к следующему
          кандидату, а не возвращаем None сразу.
+
+    None — ТОЛЬКО когда точка вне bbox ВСЕХ регионов (Пермь, Таллин). Если
+    хоть один прямоугольник её принял, а все полигоны отвергли, ответом
+    становится кандидат, к ГРАНИЦЕ которого точка ближе
+    (`_distance_to_boundary`). Иначе такая точка досталась бы вызывающему коду
+    как None, а тот (`trade_in.py`: `target_region.code if target_region else
+    DEFAULT_REGION_CODE`) молча подставил бы регион по умолчанию — адрес в
+    агломерации посчитался бы по коэффициентам Екатеринбурга. Щели между
+    полигонами соседних регионов неизбежны (у Москвы и области это независимо
+    оцифрованные в OSM линии, совпадающие с точностью до метров, а не одна
+    общая), и закрывать их возвратом к прямоугольнику нельзя — прямоугольники
+    перекрываются и врут ровно там, ради чего #3052 и делался.
+
     Региона без загруженного полигона в реестре нет (все REGIONS покрыты
-    boundaries.geojson.json) — на practике до fallback-ветки дело не доходит,
+    boundaries.geojson.json) — на практике до fallback-ветки дело не доходит,
     но она есть, чтобы новый регион без границы не «пропадал» молча, а
     работал по старому bbox-поведению до того, как для него добавят полигон.
     """
+    near_miss: list[int] = []
     for code in _POINT_LOOKUP_ORDER:
         region = REGIONS[code]
         if not is_within_bbox(lat, lon, region.bbox_region):
             continue
-        if code not in _REGION_BOUNDARIES or _point_in_region_polygon(code, lat, lon):
+        if code not in _REGION_RING_INDEX or _point_in_region_polygon(code, lat, lon):
             return region
-    return None
+        near_miss.append(code)
+    if not near_miss:
+        return None
+    if len(near_miss) == 1:
+        return REGIONS[near_miss[0]]
+    nearest = min(near_miss, key=lambda code: _distance_to_boundary(code, lat, lon))
+    return REGIONS[nearest]
 
 
 def region_by_city(city: str | None) -> Region | None:

tradein-mvp/backend/tests/test_3052_region_polygon_boundary.py

@@ -14,6 +14,7 @@ OSM/Nominatim из `region_boundaries/boundaries.geojson.json`), не по
 from __future__ import annotations
 
 import os
+import time
 
 os.environ.setdefault("DATABASE_URL", "postgresql+psycopg://test:test@localhost:5432/test")
 
@@ -66,3 +67,105 @@ def test_bbox_prefilter_never_rejects_what_polygon_would_accept() -> None:
                 assert regions.is_within_bbox(lat, lon, bbox), (
                     f"вершина полигона региона {code} ({lat}, {lon}) вне его bbox_region"
                 )
+
+
+# ── Правка регрессии первой версии #3052 (упрощение 0.005 ≈ 500 м) ──────────
+#
+# Сверка резолва против сохранённого region_code по проду (277 582 объявления
+# + 48 143 дома) на упрощённой геометрии: 69 точек агломерации не принял НИ
+# один полигон (→ None → DEFAULT_REGION_CODE=66 у вызывающего кода), куски
+# Москвы (Митино, Косино-Ухтомский, Можайский) уехали в область, 731
+# объявление области — в Москву. Точки ниже — реальные адреса из прода.
+#
+# Кунцево, «ул. Липовой Рощи 1к3» (55.7990, 37.3787) в этом списке стоит как
+# 50, хотя в исходной сверке значилось ожидание 77: улицы с таким именем в
+# Москве нет вообще (Nominatim /search «улица Липовой Рощи» — единственное
+# попадание: квартал «Спутник», Павшино, городской округ Красногорск,
+# Московская область), а сама точка лежит в 119 м ВНУТРИ области от границы
+# Москвы. Ожидание 77 там — след bbox-эпохи (прямоугольник Москвы тянется по
+# долготе до 38.10 и накрывает Красногорск), а не факт. Оно ещё и неисполнимо
+# вместе с соседней точкой той же сверки: Новоивановское (55.7013, 37.3606)
+# лежит в 72 м снаружи Москвы и ожидается как 50 — никакое правило «ближе N
+# метров к Москве ⇒ Москва» не отдаст 77 точке в 119 м и 50 точке в 72 м.
+_BORDER_CONTROL_POINTS: tuple[tuple[str, float, float, int], ...] = (
+    ("Куркино, Соловьиная Роща 8", 55.8917, 37.4008, 77),
+    ("Митино, Муравская 46к1", 55.8617, 37.3769, 77),
+    ("Косино-Ухтомский, Лухмановская 15к4", 55.7226, 37.8922, 77),
+    ("Можайский, Сколковское ш. 40к4", 55.7022, 37.3920, 77),
+    ("Новые Химки, Панфилова 12", 55.8796, 37.4294, 50),
+    ("Химки, Юннатов 11", 55.8972, 37.4613, 50),
+    ("Красногорск, Красногорский б-р 13к1", 55.8156, 37.3707, 50),
+    ("Новоивановское, Одинцовский о.", 55.7013, 37.3606, 50),
+    ("Красногорск, Липовой Рощи 1к3", 55.7990, 37.3787, 50),
+)
+
+
+def test_border_control_points_resolve_on_unsimplified_boundary() -> None:
+    for name, lat, lon, expected_code in _BORDER_CONTROL_POINTS:
+        res = regions.region_for_point(lat, lon)
+        got_code = res.code if res is not None else None
+        assert got_code == expected_code, (
+            f"{name} ({lat}, {lon}): ожидался регион {expected_code}, получен {got_code}"
+        )
+
+
+def test_no_point_inside_a_bbox_resolves_to_none() -> None:
+    """Главный инвариант правки: None законен ТОЛЬКО вне bbox всех регионов.
+
+    Прогон сеткой 0.005° по прямоугольнику Москвы: любая щель между
+    независимо оцифрованными в OSM линиями 50 и 77 всплывёт здесь как None, а
+    None у вызывающего кода (`trade_in.py`: `code if target_region else
+    DEFAULT_REGION_CODE`) означает молчаливую подстановку региона по
+    умолчанию — коэффициенты Екатеринбурга московскому адресу."""
+    checked = 0
+    for lat_step in range(91):
+        lat = 55.55 + lat_step * 0.005
+        for lon_step in range(121):
+            lon = 37.30 + lon_step * 0.005
+            assert any(
+                regions.is_within_bbox(lat, lon, r.bbox_region) for r in regions.REGIONS.values()
+            ), (lat, lon)
+            assert regions.region_for_point(lat, lon) is not None, (
+                f"точка ({lat}, {lon}) внутри bbox региона, но резолв дал None"
+            )
+            checked += 1
+    assert checked == 91 * 121
+
+
+def test_nearest_boundary_fallback_prefers_the_closer_region() -> None:
+    """Щель закрывается по расстоянию до границы, а не «первым в порядке».
+
+    Проверяем саму метрику — в fallback-ветке region_for_point исход решает
+    она. Точка в Красногорске у самой границы Москвы: линия 50 и линия 77
+    здесь общая, поэтому расстояния до них совпадают до метров, и масштаб —
+    сотня метров, а не километры (плоское приближение с поправкой cos(lat)).
+    """
+    lat, lon = 55.7990, 37.3787
+    d_50 = regions._distance_to_boundary(50, lat, lon)
+    d_77 = regions._distance_to_boundary(77, lat, lon)
+    assert abs(d_50 - d_77) * 111_320 < 20.0
+    assert 50.0 < d_50 * 111_320 < 300.0
+    # Центр Москвы от своей границы — километры (метрика не схлопнута в ноль).
+    assert regions._distance_to_boundary(77, 55.7558, 37.6173) * 111_320 > 5_000.0
+
+
+def _measure_resolve(probes: list[tuple[float, float]]) -> float:
+    started = time.perf_counter()
+    for _ in range(200):
+        for lat, lon in probes:
+            regions.region_for_point(lat, lon)
+    return (time.perf_counter() - started) / (200 * len(probes))
+
+
+def test_region_for_point_stays_under_time_budget() -> None:
+    """Точная геометрия не имеет права стоить дороже bbox-эпохи.
+
+    Потолок 150 мкс/вызов при фактических ~4 мкс (bbox каждого кольца +
+    lat-полосы рёбер, см. regions._IndexedRing): запас в 35 раз — тест
+    краснеет от настоящей регрессии алгоритма, а не от шумного раннера."""
+    probes = [(lat, lon) for _n, lat, lon, _c in _BORDER_CONTROL_POINTS]
+    probes += [(lat, lon) for _n, lat, lon, _c in _CONTROL_POINTS]
+    for lat, lon in probes:  # прогрев
+        regions.region_for_point(lat, lon)
+    best = min(_measure_resolve(probes) for _ in range(3))
+    assert best < 150e-6, f"region_for_point: {best * 1e6:.1f} мкс/вызов, потолок 150"