"""
Grid Convergence
================

Grid convergence is a spatial coordinate-reference-system (CRS) concept. It
describes how the north direction of a projected coordinate grid is rotated
relative to geographic true north at a location. The calculation does not
depend on any wind-climate workflow.

WindKit exposes this calculation through :py:func:`windkit.spatial.grid_convergence`.
There are two useful forms:

* **Natural grid convergence**: the angle between true north and the grid north
  of a dataset's own CRS.
* **Effective grid convergence**: the rotation between the grid north of a
  source CRS and the grid north of the dataset CRS at the dataset
  locations.
"""

# %%
# True North, Grid North, and Sign Convention
# -------------------------------------------
# True north, or geographic north, follows a geographic meridian toward the
# North Pole. Grid north follows the northing axis of a projected coordinate
# system. These directions are identical in a geographic CRS, but they usually
# diverge away from a projection's central meridian.
#
# WindKit uses the PROJ/pyproj sign convention: positive grid convergence means
# grid north is clockwise, or east, of true north. Negative values mean grid
# north is counter-clockwise, or west, of true north.

# %%
# Compute Natural Grid Convergence
# --------------------------------
# First, create a point dataset in geographic coordinates (EPSG:4326). The
# coordinates are longitude (``west_east``) and latitude (``south_north``).

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pathlib import Path
import geopandas as gpd
import windkit as wk

ds = wk.spatial.create_point(
    west_east=[12.0, 12.1],
    south_north=[55.0, 55.1],
    height=[100.0, 100.0],
    crs=4326,
)

print(ds)

gc_natural = wk.spatial.grid_convergence(ds)
print("Natural grid convergence for geographic coordinates:")
print(gc_natural)

# %%
# In a projected CRS the same locations can have non-zero natural convergence.
# Here the points are reprojected into EPSG:3035, the European Lambert Azimuthal
# Equal Area projection.

ds_projected = wk.spatial.reproject(ds, 3035)
gc_projected = wk.spatial.grid_convergence(ds_projected)
print("Natural grid convergence after reprojection to EPSG:3035:")
print(gc_projected)

# %%
# Spatial Variation Across a Projection
# -------------------------------------
# Grid convergence varies over a projected map. To visualize this, create a grid
# of points over Europe in EPSG:3035, calculate the natural convergence at each
# point, and plot the result.
#
# Positive values (red/warm) mean grid north tilts clockwise (east) relative to
# true north, while negative values (blue/cool) mean grid north tilts
# counter-clockwise (west). At the central meridian (10 degrees east), grid
# convergence is zero. Values grow in magnitude farther east or west because
# projected grid lines no longer follow geographic meridians.

# Define a grid over Europe in EPSG:3035 coordinate space
x = np.linspace(2.5e6, 5.5e6, 100)
y = np.linspace(1.5e6, 4.5e6, 100)
xx, yy = np.meshgrid(x, y)

# Create a WindKit point dataset with the grid locations
ds_grid = wk.spatial.create_point(
    west_east=xx.ravel(),
    south_north=yy.ravel(),
    height=np.ones_like(xx.ravel()) * 100.0,
    crs=3035,
)

# Compute natural grid convergence
gc_grid = wk.spatial.grid_convergence(ds_grid)

# Reshape back to the 2D grid shape for plotting
gc_2d = gc_grid.values.reshape(xx.shape)

# sphinx_gallery_start_ignore
fig, ax = plt.subplots(figsize=(8, 6))

im = ax.imshow(
    gc_2d,
    extent=[x.min() / 1e3, x.max() / 1e3, y.min() / 1e3, y.max() / 1e3],
    origin="lower",
    cmap="RdBu_r",
    vmin=-20,
    vmax=20,
)

geojson_path = None
try:
    geojson_path = Path(__file__).parent / "europe_boundary.geojson"
except NameError:
    for candidate in [
        Path("europe_boundary.geojson"),
        Path("docs/source/examples/europe_boundary.geojson"),
        Path("modules/windkit/docs/source/examples/europe_boundary.geojson"),
    ]:
        if candidate.exists():
            geojson_path = candidate
            break

if geojson_path is not None and geojson_path.exists():
    world = gpd.read_file(geojson_path)
    europe_3035 = world.to_crs(epsg=3035)
    europe_km = europe_3035.scale(xfact=0.001, yfact=0.001, origin=(0, 0))
    europe_km.plot(ax=ax, facecolor="none", edgecolor="black", linewidth=1.0, zorder=2)

plt.colorbar(im, label="Grid Convergence (degrees)")
plt.title("Natural Grid Convergence in EPSG:3035 across Europe")
plt.xlabel("West-East (kilometers)")
plt.ylabel("South-North (kilometers)")
plt.xlim(x.min() / 1e3, x.max() / 1e3)
plt.ylim(y.min() / 1e3, y.max() / 1e3)
plt.grid(True, linestyle="--", alpha=0.5)
# sphinx_gallery_end_ignore

# %%
# Effective Grid Convergence Close to the Pole
# ---------------------------------------------
# The ``source=`` argument computes the convergence between two CRS frames. The
# value is evaluated at the locations of the first argument (``ds_32632`` below)
# and represents the frame rotation from the source CRS grid north into the
# dataset CRS grid north.
#
# A high-latitude example makes the distinction clear. Consider longitude
# 15.0 degrees east, latitude 78.0 degrees north, near Svalbard. Compare:
#
# 1. EPSG:3035 (Lambert Azimuthal Equal Area, central meridian 10 degrees east)
# 2. EPSG:32632 (UTM Zone 32N, central meridian 9 degrees east)
#
# At 15.0 degrees east, both points lie east of their respective central
# meridians. Therefore, both projections have positive natural grid convergence
# angles. The natural angles can be large, while the effective convergence
# between the two projected frames is much smaller because the frames are
# similarly rotated.

# Create a high-latitude point close to the pole in geographic coordinates
ds_high_lat = wk.spatial.create_point(
    west_east=[15.0],
    south_north=[78.0],
    height=[100.0],
    crs=4326,
)

# Project the point into both coordinate systems
ds_3035 = wk.spatial.reproject(ds_high_lat, 3035)
ds_32632 = wk.spatial.reproject(ds_high_lat, 32632)

# Compute natural grid convergence for both
gc_natural_3035 = wk.spatial.grid_convergence(ds_3035)
gc_natural_32632 = wk.spatial.grid_convergence(ds_32632)

# Compute the effective grid convergence between them. This tells us the frame
# rotation from EPSG:3035 grid north to EPSG:32632 grid north, evaluated at the
# EPSG:32632 locations.
gc_effective_from_crs = wk.spatial.grid_convergence(ds_32632, source=3035)

# Note that you can pass a dataset as the source and the projection will be inferred.
gc_effective_from_ds = wk.spatial.grid_convergence(ds_32632, source=ds_3035)

print("Latitude: 78 degrees north, Longitude: 15 degrees east")
print(f"Natural GC in EPSG:3035 (LAEA): {gc_natural_3035.values[0]:.4f} degrees")
print(f"Natural GC in EPSG:32632 (UTM 32N): {gc_natural_32632.values[0]:.4f} degrees")
print(
    "Effective GC (source dataset 3035 to dataset CRS 32632): "
    f"{gc_effective_from_ds.values[0]:.4f} degrees"
)

# %%
# Bin True-North Time Series into a Grid Frame
# -----------------------------------------------
# When a time series is known to be true-north-relative, pass a ``wind_dir_crs``
# to `bwc_from_tswc`. If the TSWC already carries ``attrs["wind_dir_crs"]``,
# WindKit treats the directions as already expressed in that frame and avoids
# rotating them again. A projected CRS rotates directions into that grid
# frame; a geographic CRS retains the true-north frame. WindKit records the
# resolved CRS WKT on the BWC.

# %%
tswc = wk.create_tswc(
    ds_high_lat,
    date_range=pd.date_range("2001-01-01", periods=24, freq="h"),
)

# Mark the input as already expressed in EPSG:32632 grid north.
tswc_marked = tswc.assign_attrs(wind_dir_crs=32632)

bwc_marked = wk.bwc_from_tswc(tswc_marked)
bwc_same_frame = wk.bwc_from_tswc(tswc_marked, wind_dir_crs=32632)
bwc_other_frame = wk.bwc_from_tswc(tswc_marked, wind_dir_crs=3035)

print(bwc_marked.attrs["wind_dir_crs"])
print(bwc_same_frame.attrs["wind_dir_crs"])
print(bwc_other_frame.attrs["wind_dir_crs"])

# %%
# When to Use This
# ----------------
# Use :py:func:`~windkit.spatial.grid_convergence` when reasoning about CRS frames, diagnosing map
# projection behavior, or checking direction reference frames. If you are using
# PyWAsP wind climates, see the PyWAsP meridian-convergence tutorial for how
# these CRS rotations affect BWC, GWC, and WWC workflows.