vector/services/osm_local.py

"""Локальный OSM-слой (B16): пространственные запросы к PostGIS вместо Overpass.

Таблицы planet_osm_* заливаются osm2pgsql из дампа Беларуси
(scripts/b16-import.sh; SRID 3857 — дефолт osm2pgsql, проекция Web Mercator).
Обновление — разово/по мере надобности. Никакого интернета.

Доступ к БД — свой sync-движок из DATABASE_URL (asyncpg → psycopg2).
"""
from __future__ import annotations

import logging
import math
import os
from typing import Any, Optional

from sqlalchemy import create_engine, text
from sqlalchemy.engine import Engine

logger = logging.getLogger('vector.geo.osm')

_PG_URL = os.getenv('DATABASE_URL', 'postgresql://postgres:postgres@postgres:5432/vector_mchs')
if _PG_URL.startswith('postgresql+asyncpg://'):
    _PG_URL = _PG_URL.replace('postgresql+asyncpg://', 'postgresql://')

_engine: Optional[Engine] = None
_available: Optional[bool] = None


def _get_engine() -> Engine:
    global _engine
    if _engine is None:
        _engine = create_engine(_PG_URL, pool_pre_ping=True)
    return _engine


def osm_available() -> bool:
    """True, если локальные OSM-таблицы залиты (иначе geo_service падает на Overpass)."""
    global _available
    if _available is not None:
        return _available
    try:
        with _get_engine().connect() as conn:
            row = conn.execute(text(
                "SELECT count(*) FROM information_schema.tables "
                "WHERE table_name = 'planet_osm_line'"
            )).scalar()
            _available = bool(row)
    except Exception as e:
        logger.warning('OSM PostGIS недоступен: %s', e)
        _available = False
    return _available


def _sector_wkt_3857(lat: float, lon: float, radius_m: int, direction: str) -> str:
    """Сектор-«пирог» 45° в Web Mercator (проекция данных planet_osm_*).

    Строим в WGS84 (локальные метрики), потом трансформируем — mercator-метры
    на широте Беларуси раздуты, строить дугу в них нельзя.
    """
    angles = {'N': 0, 'NE': 45, 'E': 90, 'SE': 135, 'S': 180, 'SW': 225, 'W': 270, 'NW': 315}
    az = angles.get(direction, 0)
    half = 22.5
    seg = 24
    pts: list[tuple[float, float]] = []
    for i in range(seg + 1):
        bearing = math.radians(az - half + (2 * half) * i / seg)
        dlat = (radius_m / 111320.0) * math.cos(bearing)
        dlon = (radius_m / (111320.0 * max(0.2, math.cos(math.radians(lat))))) * math.sin(bearing)
        pts.append((lon + dlon, lat + dlat))
    # Замыкаем пирог через центр
    pts_wkt = f'{lon:.7f} {lat:.7f}, ' + ','.join(f'{p[0]:.7f} {p[1]:.7f}' for p in pts)
    return f'POLYGON(({pts_wkt}, {lon:.7f} {lat:.7f}))'


def get_zone_features_postgis(lat: float, lon: float, direction: str, radius_m: int) -> dict[str, Any]:
    """Характеристики сектора — 1 запрос к локальному PostGIS.

    Контракт тот же, что get_zone_features (Overpass): roads_km, road_density,
    water_distance_km, settlement_distance_km, forest_pct.
    """
    sector_wkt = _sector_wkt_3857(lat, lon, radius_m, direction)
    center_pt = f'POINT({lon} {lat})'

    q = text("""
        WITH s AS (
            SELECT ST_Transform(ST_SetSRID(ST_GeomFromText(:sector_wkt), 4326), 3857) AS sector,
                   ST_Transform(ST_SetSRID(ST_GeomFromText(:center_pt), 4326), 3857) AS center_m
        ),
        road_len AS (
            SELECT COALESCE(SUM(ST_Length(l.way)), 0) AS meters
            FROM planet_osm_line l, s
            WHERE l.highway IS NOT NULL
              AND l.way && s.sector
              AND ST_Intersects(l.way, s.sector)
        ),
        forest_area AS (
            SELECT COALESCE(SUM(ST_Area(ST_Intersection(p.way, s.sector))), 0) AS m2
            FROM planet_osm_polygon p, s
            WHERE (p.landuse = 'forest' OR p."natural" = 'wood')
              AND p.way && s.sector
        ),
        water_near AS (
            SELECT MIN(ST_Distance(p.way, s.center_m)) AS meters
            FROM planet_osm_polygon p, s
            WHERE (p."natural" = 'water' OR p.landuse = 'reservoir')
              AND p.way && s.sector
              AND ST_Intersects(p.way, s.sector)
        ),
        waterway_near AS (
            SELECT MIN(ST_Distance(l.way, s.center_m)) AS meters
            FROM planet_osm_line l, s
            WHERE l.waterway IN ('river','stream','canal','ditch','drain')
              AND l.way && s.sector
              AND ST_Intersects(l.way, s.sector)
        ),
        settle_near AS (
            SELECT MIN(ST_Distance(p.way, s.center_m)) AS meters
            FROM planet_osm_point p, s
            WHERE p.place IN ('village','town','city')
              AND p.way && s.sector
              AND ST_Intersects(p.way, s.sector)
        )
        SELECT
            (SELECT meters FROM road_len),
            (SELECT m2 FROM forest_area),
            (SELECT meters FROM water_near),
            (SELECT meters FROM waterway_near),
            (SELECT meters FROM settle_near)
    """)

    # Масштаб 3857 на широте Минска (~53.9°): cos(lat)
    scale = max(0.2, math.cos(math.radians(lat)))

    with _get_engine().connect() as conn:
        row = conn.execute(q, {
            'sector_wkt': sector_wkt, 'center_pt': center_pt, 'scale': scale,
        }).one()

    # 3857-метры на широте lat раздуты в 1/cos(lat) раз:
    # real = mercator * cos(lat) (проверено: merc/geo = 1/cos(lat) точно).
    roads_km = float(row[0] or 0) / 1000.0 * scale
    forest_m2 = float(row[1] or 0) * (scale * scale)   # mercator-м² → реальные м²
    # Дистанции: mercator-метры → реальные метры (та же поправка, что у дорог:
    # раньше деление на cos здесь ДВОИЛО раздувание — 850 «км» до озера в Браславе).
    water_poly_m = float(row[2]) * scale if row[2] is not None else None
    waterway_m = float(row[3]) * scale if row[3] is not None else None
    settle_m = float(row[4]) * scale if row[4] is not None else None

    search_area_km2 = math.pi * (radius_m / 1000.0) ** 2
    road_density = roads_km / search_area_km2 if search_area_km2 > 0 else 0.0

    # Контракт Zone.forest_pct — ДОЛЯ леса 0..1 (так его читает
    # scoring_service.score_zone: weight × forest_pct; тесты скоринга
    # используют 0.3–0.8 как доли). Раньше отдались проценты 0..100 —
    # лес в 100 раз перевешивал остальные факторы (зоны «вверх-вправо»).
    forest_pct = 0.0
    if forest_m2 > 0 and search_area_km2 > 0:
        forest_pct = round(min(1.0, (forest_m2 / 1e6) / search_area_km2), 3)

    water_candidates = [float(m) / 1000.0 for m in (water_poly_m, waterway_m) if m is not None]
    water_distance_km = round(min(water_candidates), 3) if water_candidates else None
    settlement_distance_km = round(float(settle_m) / 1000.0, 3) if settle_m is not None else None

    return {
        'roads_km': round(roads_km, 3),
        'road_density': round(road_density, 2),
        'water_distance_km': water_distance_km,
        'settlement_distance_km': settlement_distance_km,
        'forest_pct': forest_pct,
    }