# -*- coding: utf-8 -*-
#
# Isoliner - грид и изолинии (QGIS).
# © 2026 ООО «Информ++» (www.informpp.ru).
# SPDX-License-Identifier: GPL-2.0-or-later
#
# Это свободная программа: вы можете распространять её и/или изменять на
# условиях Стандартной общественной лицензии GNU (GNU GPL), опубликованной
# Фондом свободного ПО (FSF), - либо версии 2 Лицензии, либо (на ваше
# усмотрение) любой более поздней версии.
#
# Программа распространяется в надежде на полезность, но БЕЗ КАКИХ-ЛИБО
# ГАРАНТИЙ. Подробнее см. GNU GPL. Полный текст - в файле LICENSE.
"""
Гидрогеология: градиент напора и направление потока (чистый NumPy, без QGIS).

По растру напора (пьезометрической поверхности) строит две производные карты:
гидравлический градиент i = |∇h| и азимут направления потока. Это геометрия
поля напора, без проницаемости: скорость фильтрации по Дарси (v = −K·∇h)
требует коэффициента фильтрации K и здесь не считается.

Соглашение о координатах совпадает с движком kb2d и выводом GeoTIFF: строка 0
растра - север, индекс строки растёт на юг (как в GDAL GeoTransform с
отрицательным шагом по Y). Поэтому производная «на север» берётся с обратным
знаком к производной по строкам.
"""
import numpy as np

EPS = 1e-12


def head_gradient(h, cellx, celly=None, nodata=-9999.0):
    """Градиент напора. Возвращает (модуль, азимут) - два float32-массива.

    Модуль - гидравлический градиент i = |∇h|, безразмерный (м/м): крутизна
    поверхности напора, скорость изменения напора на единицу длины.

    Азимут - компасный, в градусах (0 = север, по часовой), указывает
    направление ПОТОКА, то есть вниз по градиенту (−∇h): подземная вода течёт
    от большего напора к меньшему.

    cellx, celly - размер ячейки по X и Y (если celly не задан - квадратная).
    Ячейки с nodata и края неполного шаблона дифференцирования дают nodata.
    На плоских участках (нулевой градиент) направление не определено, и азимут
    там тоже nodata, а модуль равен нулю."""
    if celly is None:
        celly = cellx
    h = np.asarray(h, float)
    z = np.where(h == nodata, np.nan, h)
    gx = np.gradient(z, axis=1) / float(cellx)      # ∂h/∂x, восток+
    gy = -np.gradient(z, axis=0) / float(celly)     # ∂h/∂y, север+ (строка вниз = юг)
    mag = np.sqrt(gx * gx + gy * gy)
    az = np.degrees(np.arctan2(-gx, -gy)) % 360.0   # вектор потока = −∇h
    valid = np.isfinite(mag)
    mag_out = np.where(valid, mag, nodata).astype(np.float32)
    az_out = np.where(valid & (mag > EPS), az, nodata).astype(np.float32)
    return mag_out, az_out


def flow_samples(mag, az, geotransform, step, nodata=-9999.0):
    """Прореживание для векторного поля стрелок: каждая step-я валидная ячейка.

    Возвращает (xs, ys, az_deg, grad) - массивы координат центров ячеек и
    значений (азимут потока в градусах, гидравлический градиент). geotransform -
    GDAL GeoTransform исходного растра (originX, dx, 0, originY, 0, dy), dy<0.
    Точки с nodata по любому из растров пропускаются."""
    ny, nx = mag.shape
    ox, dx, _, oy, _, dy = geotransform
    step = max(int(step), 1)
    rr = np.arange(0, ny, step)
    cc = np.arange(0, nx, step)
    R, C = np.meshgrid(rr, cc, indexing="ij")
    M = mag[R, C]
    A = az[R, C]
    ok = np.isfinite(M) & (M != nodata) & (A != nodata)
    X = ox + (C + 0.5) * dx
    Y = oy + (R + 0.5) * dy                          # dy<0 -> строка вниз = юг
    return X[ok].ravel(), Y[ok].ravel(), A[ok].ravel(), M[ok].ravel()
