# -*- coding: utf-8 -*-
"""Tests for the pyflwdir.gis_utils module."""

import pytest
import numpy as np
from pyflwdir import gis_utils as gis
from affine import Affine

# glob total area
glob_area = 4 * np.pi * gis._R**2
glob_circ = 2 * np.pi * gis._R


## TRANSFORM
# Adapted from https://github.com/mapbox/rasterio/blob/master/tests/test_transform.py
prof = {
    "width": 6120,
    "height": 4920,
    "res": 1 / 1200.0,
    "bounds": [-10.5, 51.4, -5.4, 55.5],
    "transform": Affine(1 / 1200.0, 0.0, -10.5, 0.0, -1 / 1200.0, 55.5),
}


def test_from_origin():
    w, _, _, n = prof["bounds"]
    tr = gis.transform_from_origin(w, n, prof["res"], prof["res"])
    assert [round(v, 7) for v in tr] == [round(v, 7) for v in prof["transform"]]


def test_from_bounds():
    w, s, e, n = prof["bounds"]
    tr = gis.transform_from_bounds(w, s, e, n, prof["width"], prof["height"])
    assert [round(v, 7) for v in tr] == [round(v, 7) for v in prof["transform"]]


def test_array_bounds():
    bounds0 = np.asarray(prof["bounds"])
    bounds = gis.array_bounds(prof["height"], prof["width"], prof["transform"])
    assert np.all(bounds0 == np.asarray(bounds).round(7))


def test_xy():
    aff = prof["transform"]
    ul_x, ul_y = aff * (0, 0)
    xoff = aff.a
    yoff = aff.e
    assert gis.xy(aff, 0, 0, offset="ul") == (ul_x, ul_y)
    assert gis.xy(aff, 0, 0, offset="ur") == (ul_x + xoff, ul_y)
    assert gis.xy(aff, 0, 0, offset="ll") == (ul_x, ul_y + yoff)
    expected = (ul_x + xoff, ul_y + yoff)
    assert gis.xy(aff, 0, 0, offset="lr") == expected
    expected = (ul_x + xoff / 2, ul_y + yoff / 2)
    assert gis.xy(aff, 0, 0, offset="center") == expected
    assert (
        gis.xy(aff, 0, 0, offset="lr")
        == gis.xy(aff, 0, 1, offset="ll")
        == gis.xy(aff, 1, 1, offset="ul")
        == gis.xy(aff, 1, 0, offset="ur")
    )


def test_rowcol():
    aff = gis.IDENTITY  # N->S changed in version 0.5
    left, bottom, right, top = (0, -200, 100, 0)
    assert gis.rowcol(aff, left, top) == (top, left)
    assert gis.rowcol(aff, right, top) == (top, right)
    assert gis.rowcol(aff, right, bottom) == (-bottom, right)
    assert gis.rowcol(aff, left, bottom) == (-bottom, left)


def test_idxs_to_coords():
    shape = (10, 8)
    idxs = np.arange(shape[0] * shape[1]).reshape(shape)
    transform = gis.Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0)
    xs, ys = gis.idxs_to_coords(idxs, transform, shape)
    assert np.all(ys == (np.arange(shape[0]) + 0.5)[:, None])
    assert np.all(xs == np.arange(shape[1]) + 0.5)
    with pytest.raises(IndexError):
        gis.idxs_to_coords(np.array([-1]), transform, shape)


def test_coords_to_idxs():
    shape = (10, 8)
    idxs0 = np.arange(shape[0] * shape[1])
    transform = gis.Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0)
    xs, ys = np.meshgrid(np.arange(shape[1]) + 0.5, np.arange(shape[0]) + 0.5)
    idxs = gis.coords_to_idxs(xs, ys, transform, shape)
    assert np.all(idxs.ravel() == idxs0)
    with pytest.raises(IndexError):
        gis.coords_to_idxs(ys, xs, transform, shape)


def test_affine_to_coords():
    shape = (10, 8)
    transform = gis.Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0)
    xs, ys = gis.affine_to_coords(transform, shape)
    assert np.all(ys == np.arange(shape[0]) + 0.5)
    assert np.all(xs == np.arange(shape[1]) + 0.5)


def test_reggrid_dx():
    # also tests degree_metres_x
    # area of glob in 1 degree cells
    lats = np.array([0.0])
    lons = np.arange(-179.5, 180)
    dx = gis.reggrid_dx(lats, lons)
    assert dx.shape == (lats.size, lons.size)
    assert dx.sum().round(3) == 40075004.88


def test_reggrid_dy():
    # also tests degree_metres_y
    # area of glob in 1 degree cells
    lats = np.arange(-89.5, 90)
    lons = np.array([0.0])
    dy = gis.reggrid_dy(lats, lons)
    assert dy.shape == (lats.size, lons.size)
    assert dy.sum().round(3) == 20003925.600


def test_cellarea():
    # area of whole sphere
    assert gis.cellarea(0, 360, 180) == glob_area
    # area of 1 degree cell
    assert gis.cellarea(0, 1, 1) == 12364154779.389229


def test_reggrid_area():
    # area of glob in 1 degree cells
    lats = np.arange(-89.5, 90)
    lons = np.arange(-179.5, 180)
    assert gis.reggrid_area(lats, lons).sum().round() == np.round(glob_area)


def test_distance():
    # transform=gis.IDENTITY
    assert gis.distance(0, 1, 3) == 1  # horizontal
    assert gis.distance(0, 3, 3) == 1  # vertical
    assert gis.distance(4, 0, 3) == np.hypot(1, 1)  # diagonal
    assert gis.distance(0, 4, 3, True) == gis.distance(4, 0, 3, True)
    assert gis.distance(0, 1, 3, False) == gis.distance(7, 8, 3, False)
    assert gis.distance(0, 1, 3, True) != gis.distance(7, 8, 3, True)


def test_edge():
    a = np.ones((5, 5), dtype=bool)
    b = a.copy()
    b[1:-1, 1:-1] = False
    assert np.all(gis.get_edge(a) == b)
    a[np.diag_indices(5)] = False
    assert np.all(gis.get_edge(a) == a)
    b = a.copy()
    b[1, 3], b[3, 1] = False, False
    d4 = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
    assert np.all(gis.get_edge(a, structure=d4) == b)


def test_spread():
    a = np.zeros((5, 5))
    a[2, 2] = 1
    out, src, dst = gis.spread2d(a, nodata=0)
    assert np.all(out == 1)
    assert np.all(src == 12)
    assert np.isclose(np.max(dst), 2 * np.hypot(1, 1))
    a[-1, -1] = 2
    out, src, dst = gis.spread2d(a, nodata=0, msk=a != 2, latlon=True)
    assert np.all(out[a != 2] == 1)
    assert np.all(out.flat[src] == out)
    assert dst[-1, -1] == 0

# -*- coding: utf-8 -*-
"""Tests for the pyflwdir.gis_utils module."""

import pytest
import numpy as np
from pyflwdir import gis_utils as gis
from affine import Affine

# glob total area
glob_area = 4 * np.pi * gis._R**2
glob_circ = 2 * np.pi * gis._R


## TRANSFORM
# Adapted from https://github.com/mapbox/rasterio/blob/master/tests/test_transform.py
prof = {
    "width": 6120,
    "height": 4920,
    "res": 1 / 1200.0,
    "bounds": [-10.5, 51.4, -5.4, 55.5],
    "transform": Affine(1 / 1200.0, 0.0, -10.5, 0.0, -1 / 1200.0, 55.5),
}


def test_from_origin():
    w, _, _, n = prof["bounds"]
    tr = gis.transform_from_origin(w, n, prof["res"], prof["res"])
    assert [round(v, 7) for v in tr] == [round(v, 7) for v in prof["transform"]]


def test_from_bounds():
    w, s, e, n = prof["bounds"]
    tr = gis.transform_from_bounds(w, s, e, n, prof["width"], prof["height"])
    assert [round(v, 7) for v in tr] == [round(v, 7) for v in prof["transform"]]


def test_array_bounds():
    bounds0 = np.asarray(prof["bounds"])
    bounds = gis.array_bounds(prof["height"], prof["width"], prof["transform"])
    assert np.all(bounds0 == np.asarray(bounds).round(7))


def test_xy():
    aff = prof["transform"]
    ul_x, ul_y = aff * (0, 0)
    xoff = aff.a
    yoff = aff.e
    assert gis.xy(aff, 0, 0, offset="ul") == (ul_x, ul_y)
    assert gis.xy(aff, 0, 0, offset="ur") == (ul_x + xoff, ul_y)
    assert gis.xy(aff, 0, 0, offset="ll") == (ul_x, ul_y + yoff)
    expected = (ul_x + xoff, ul_y + yoff)
    assert gis.xy(aff, 0, 0, offset="lr") == expected
    expected = (ul_x + xoff / 2, ul_y + yoff / 2)
    assert gis.xy(aff, 0, 0, offset="center") == expected
    assert (
        gis.xy(aff, 0, 0, offset="lr")
        == gis.xy(aff, 0, 1, offset="ll")
        == gis.xy(aff, 1, 1, offset="ul")
        == gis.xy(aff, 1, 0, offset="ur")
    )


def test_rowcol():
    aff = gis.IDENTITY  # N->S changed in version 0.5
    left, bottom, right, top = (0, -200, 100, 0)
    assert gis.rowcol(aff, left, top) == (top, left)
    assert gis.rowcol(aff, right, top) == (top, right)
    assert gis.rowcol(aff, right, bottom) == (-bottom, right)
    assert gis.rowcol(aff, left, bottom) == (-bottom, left)


def test_idxs_to_coords():
    shape = (10, 8)
    idxs = np.arange(shape[0] * shape[1]).reshape(shape)
    transform = gis.Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0)
    xs, ys = gis.idxs_to_coords(idxs, transform, shape)
    assert np.all(ys == (np.arange(shape[0]) + 0.5)[:, None])
    assert np.all(xs == np.arange(shape[1]) + 0.5)
    with pytest.raises(IndexError):
        gis.idxs_to_coords(np.array([-1]), transform, shape)


def test_coords_to_idxs():
    shape = (10, 8)
    idxs0 = np.arange(shape[0] * shape[1])
    transform = gis.Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0)
    xs, ys = np.meshgrid(np.arange(shape[1]) + 0.5, np.arange(shape[0]) + 0.5)
    idxs = gis.coords_to_idxs(xs, ys, transform, shape)
    assert np.all(idxs.ravel() == idxs0)
    with pytest.raises(IndexError):
        gis.coords_to_idxs(ys, xs, transform, shape)


def test_affine_to_coords():
    shape = (10, 8)
    transform = gis.Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0)
    xs, ys = gis.affine_to_coords(transform, shape)
    assert np.all(ys == np.arange(shape[0]) + 0.5)
    assert np.all(xs == np.arange(shape[1]) + 0.5)


def test_reggrid_dx():
    # also tests degree_metres_x
    # area of glob in 1 degree cells
    lats = np.array([0.0])
    lons = np.arange(-179.5, 180)
    dx = gis.reggrid_dx(lats, lons)
    assert dx.shape == (lats.size, lons.size)
    assert dx.sum().round(3) == 40075004.88


def test_reggrid_dy():
    # also tests degree_metres_y
    # area of glob in 1 degree cells
    lats = np.arange(-89.5, 90)
    lons = np.array([0.0])
    dy = gis.reggrid_dy(lats, lons)
    assert dy.shape == (lats.size, lons.size)
    assert dy.sum().round(3) == 20003925.600


def test_cellarea():
    # area of whole sphere
    assert gis.cellarea(0, 360, 180) == glob_area
    # area of 1 degree cell
    assert gis.cellarea(0, 1, 1) == 12364154779.389229


def test_reggrid_area():
    # area of glob in 1 degree cells
    lats = np.arange(-89.5, 90)
    lons = np.arange(-179.5, 180)
    assert gis.reggrid_area(lats, lons).sum().round() == np.round(glob_area)


def test_distance():
    # transform=gis.IDENTITY
    assert gis.distance(0, 1, 3) == 1  # horizontal
    assert gis.distance(0, 3, 3) == 1  # vertical
    assert gis.distance(4, 0, 3) == np.hypot(1, 1)  # diagonal
    assert gis.distance(0, 4, 3, True) == gis.distance(4, 0, 3, True)
    assert gis.distance(0, 1, 3, False) == gis.distance(7, 8, 3, False)
    assert gis.distance(0, 1, 3, True) != gis.distance(7, 8, 3, True)


def test_edge():
    a = np.ones((5, 5), dtype=bool)
    b = a.copy()
    b[1:-1, 1:-1] = False
    assert np.all(gis.get_edge(a) == b)
    a[np.diag_indices(5)] = False
    assert np.all(gis.get_edge(a) == a)
    b = a.copy()
    b[1, 3], b[3, 1] = False, False
    d4 = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
    assert np.all(gis.get_edge(a, structure=d4) == b)


def test_spread():
    a = np.zeros((5, 5))
    a[2, 2] = 1
    out, src, dst = gis.spread2d(a, nodata=0)
    assert np.all(out == 1)
    assert np.all(src == 12)
    assert np.isclose(np.max(dst), 2 * np.hypot(1, 1))
    a[-1, -1] = 2
    out, src, dst = gis.spread2d(a, nodata=0, msk=a != 2, latlon=True)
    assert np.all(out[a != 2] == 1)
    assert np.all(out.flat[src] == out)
    assert dst[-1, -1] == 0

# -*- coding: utf-8 -*- Tests for the pyflwdir.gis_utils module. # glob total area ## TRANSFORM # Adapted from https://github.com/mapbox/rasterio/blob/master/tests/test_transform.py # N->S changed in version 0.5 # also tests degree_metres_x # area of glob in 1 degree cells # also tests degree_metres_y # area of glob in 1 degree cells # area of whole sphere # area of 1 degree cell # area of glob in 1 degree cells # transform=gis.IDENTITY # horizontal # vertical # diagonal