gendesign/backend/app/services/cadastre/grid_geometry.py

"""Геометрические утилиты для grid-walk по кадастровому кварталу.

Используется в bulk_harvest.py для вычисления bbox квартала в EPSG:3857
и генерации сетки точек для WMS GetFeatureInfo.
"""

from __future__ import annotations

import logging
import math

from sqlalchemy import text
from sqlalchemy.orm import Session

logger = logging.getLogger(__name__)

# Радиус сферы Web Mercator (EPSG:3857), м. Нужен для инверсии Y → широта.
_WEB_MERCATOR_R = 6378137.0


def quarter_bbox_3857(db: Session, quarter: str) -> tuple[float, float, float, float] | None:
    """Получить bbox квартала из cad_quarters_geom в EPSG:3857.

    Args:
        db: SQLAlchemy session.
        quarter: кадастровый номер квартала, например '66:41:0303161'.

    Returns:
        (xmin, ymin, xmax, ymax) в метрах EPSG:3857, или None если квартал не найден.
    """
    row = (
        db.execute(
            text("""
                SELECT
                    ST_XMin(b) AS xmin,
                    ST_YMin(b) AS ymin,
                    ST_XMax(b) AS xmax,
                    ST_YMax(b) AS ymax
                FROM (
                    SELECT ST_Extent(ST_Transform(geom, 3857)) AS b
                    FROM cad_quarters_geom
                    WHERE cad_number = :q
                ) sub
            """),
            {"q": quarter},
        )
        .mappings()
        .first()
    )
    if not row or row["xmin"] is None:
        logger.warning("quarter_bbox_3857: квартал %s не найден в cad_quarters_geom", quarter)
        return None
    bbox = (
        float(row["xmin"]),
        float(row["ymin"]),
        float(row["xmax"]),
        float(row["ymax"]),
    )
    # Sanity check: cad_quarters_geom имеет ~72% rows с broken micro-precision
    # geometry (width <0.01m). Это вызывает useless 1px tiles в WMS и 500 errors
    # от NSPD. Реальный квартал ЕКБ имеет width 200-4000м. Skip grid-walk для
    # broken — snapshot phase даст 20 per cat, что норм для MVP.
    #
    # Bug #1629: EPSG:3857 (Web Mercator) — конформная, НЕ равнопротяжённая
    # проекция. Линейный масштаб = sec(lat); на широте ЕКБ (~56-60°N) фактор
    # ≈1.8-2.0, т.е. 1 наземный метр ≈ 1.8-2.0 единиц 3857. Пороги ниже заданы
    # в НАЗЕМНЫХ метрах, поэтому переводим протяжённость bbox из единиц 3857 в
    # наземные метры умножением на cos(lat) (lat берём из центра bbox через
    # инверсию Y). Иначе верхняя граница 10000 отсекала бы реальные крупные/
    # вытянутые кварталы (10000 единиц 3857 ≈ всего ~5050 м земли на широте ЕКБ).
    y_mid = (bbox[1] + bbox[3]) / 2.0
    lat_mid = 2.0 * math.atan(math.exp(y_mid / _WEB_MERCATOR_R)) - math.pi / 2.0
    ground_scale = math.cos(lat_mid)  # ground_m = mercator_units * cos(lat)
    width = (bbox[2] - bbox[0]) * ground_scale
    height = (bbox[3] - bbox[1]) * ground_scale
    if width < 100 or height < 100 or width > 10000 or height > 10000:
        logger.warning(
            "quarter_bbox_3857: квартал %s broken geom — width=%.2fm height=%.2fm "
            "(expected 100-10000m). Grid-walk skipped, only snapshot phase.",
            quarter,
            width,
            height,
        )
        return None
    return bbox


def generate_grid_click_points(
    bbox: tuple[float, float, float, float],
    grid_size: int = 15,
    tile_size: int = 512,
) -> list[tuple[tuple[float, float, float, float], tuple[int, int]]]:
    """Генерировать grid_size × grid_size ячеек для bbox квартала.

    Каждая ячейка — свой sub-bbox + центральный пиксель (click_xy) для
    WMS GetFeatureInfo. 15×15 = 225 запросов обеспечивают хорошее покрытие
    при шаге ~50м для типичного квартала ЕКБ (~750м × ~750м).

    Args:
        bbox: (xmin, ymin, xmax, ymax) в EPSG:3857.
        grid_size: количество ячеек по каждой оси (15 → 225 запросов).
        tile_size: размер виртуального WMS тайла в пикселях (квадратный).

    Returns:
        Список (sub_bbox, click_xy) — sub_bbox в EPSG:3857, click_xy в пикселях.
    """
    xmin, ymin, xmax, ymax = bbox
    x_step = (xmax - xmin) / grid_size
    y_step = (ymax - ymin) / grid_size

    result: list[tuple[tuple[float, float, float, float], tuple[int, int]]] = []

    for i in range(grid_size):
        for j in range(grid_size):
            cell_xmin = xmin + i * x_step
            cell_ymin = ymin + j * y_step
            cell_xmax = cell_xmin + x_step
            cell_ymax = cell_ymin + y_step
            # Клик в центр ячейки
            click_x = tile_size // 2
            click_y = tile_size // 2
            result.append(
                (
                    (cell_xmin, cell_ymin, cell_xmax, cell_ymax),
                    (click_x, click_y),
                )
            )

    return result