{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Grid Convergence\n\nGrid convergence is a spatial coordinate-reference-system (CRS) concept. It\ndescribes how the north direction of a projected coordinate grid is rotated\nrelative to geographic true north at a location. The calculation does not\ndepend on any wind-climate workflow.\n\nWindKit exposes this calculation through :py:func:`windkit.spatial.grid_convergence`.\nThere are two useful forms:\n\n* **Natural grid convergence**: the angle between true north and the grid north\n  of a dataset's own CRS.\n* **Effective grid convergence**: the rotation between the grid north of a\n  source CRS and the grid north of the dataset CRS at the dataset\n  locations.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## True North, Grid North, and Sign Convention\nTrue north, or geographic north, follows a geographic meridian toward the\nNorth Pole. Grid north follows the northing axis of a projected coordinate\nsystem. These directions are identical in a geographic CRS, but they usually\ndiverge away from a projection's central meridian.\n\nWindKit uses the PROJ/pyproj sign convention: positive grid convergence means\ngrid north is clockwise, or east, of true north. Negative values mean grid\nnorth is counter-clockwise, or west, of true north.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute Natural Grid Convergence\nFirst, create a point dataset in geographic coordinates (EPSG:4326). The\ncoordinates are longitude (``west_east``) and latitude (``south_north``).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom pathlib import Path\nimport geopandas as gpd\nimport windkit as wk\n\nds = wk.spatial.create_point(\n    west_east=[12.0, 12.1],\n    south_north=[55.0, 55.1],\n    height=[100.0, 100.0],\n    crs=4326,\n)\n\nprint(ds)\n\ngc_natural = wk.spatial.grid_convergence(ds)\nprint(\"Natural grid convergence for geographic coordinates:\")\nprint(gc_natural)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In a projected CRS the same locations can have non-zero natural convergence.\nHere the points are reprojected into EPSG:3035, the European Lambert Azimuthal\nEqual Area projection.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ds_projected = wk.spatial.reproject(ds, 3035)\ngc_projected = wk.spatial.grid_convergence(ds_projected)\nprint(\"Natural grid convergence after reprojection to EPSG:3035:\")\nprint(gc_projected)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Spatial Variation Across a Projection\nGrid convergence varies over a projected map. To visualize this, create a grid\nof points over Europe in EPSG:3035, calculate the natural convergence at each\npoint, and plot the result.\n\nPositive values (red/warm) mean grid north tilts clockwise (east) relative to\ntrue north, while negative values (blue/cool) mean grid north tilts\ncounter-clockwise (west). At the central meridian (10 degrees east), grid\nconvergence is zero. Values grow in magnitude farther east or west because\nprojected grid lines no longer follow geographic meridians.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Define a grid over Europe in EPSG:3035 coordinate space\nx = np.linspace(2.5e6, 5.5e6, 100)\ny = np.linspace(1.5e6, 4.5e6, 100)\nxx, yy = np.meshgrid(x, y)\n\n# Create a WindKit point dataset with the grid locations\nds_grid = wk.spatial.create_point(\n    west_east=xx.ravel(),\n    south_north=yy.ravel(),\n    height=np.ones_like(xx.ravel()) * 100.0,\n    crs=3035,\n)\n\n# Compute natural grid convergence\ngc_grid = wk.spatial.grid_convergence(ds_grid)\n\n# Reshape back to the 2D grid shape for plotting\ngc_2d = gc_grid.values.reshape(xx.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Effective Grid Convergence Close to the Pole\nThe ``source=`` argument computes the convergence between two CRS frames. The\nvalue is evaluated at the locations of the first argument (``ds_32632`` below)\nand represents the frame rotation from the source CRS grid north into the\ndataset CRS grid north.\n\nA high-latitude example makes the distinction clear. Consider longitude\n15.0 degrees east, latitude 78.0 degrees north, near Svalbard. Compare:\n\n1. EPSG:3035 (Lambert Azimuthal Equal Area, central meridian 10 degrees east)\n2. EPSG:32632 (UTM Zone 32N, central meridian 9 degrees east)\n\nAt 15.0 degrees east, both points lie east of their respective central\nmeridians. Therefore, both projections have positive natural grid convergence\nangles. The natural angles can be large, while the effective convergence\nbetween the two projected frames is much smaller because the frames are\nsimilarly rotated.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Create a high-latitude point close to the pole in geographic coordinates\nds_high_lat = wk.spatial.create_point(\n    west_east=[15.0],\n    south_north=[78.0],\n    height=[100.0],\n    crs=4326,\n)\n\n# Project the point into both coordinate systems\nds_3035 = wk.spatial.reproject(ds_high_lat, 3035)\nds_32632 = wk.spatial.reproject(ds_high_lat, 32632)\n\n# Compute natural grid convergence for both\ngc_natural_3035 = wk.spatial.grid_convergence(ds_3035)\ngc_natural_32632 = wk.spatial.grid_convergence(ds_32632)\n\n# Compute the effective grid convergence between them. This tells us the frame\n# rotation from EPSG:3035 grid north to EPSG:32632 grid north, evaluated at the\n# EPSG:32632 locations.\ngc_effective_from_crs = wk.spatial.grid_convergence(ds_32632, source=3035)\n\n# Note that you can pass a dataset as the source and the projection will be inferred.\ngc_effective_from_ds = wk.spatial.grid_convergence(ds_32632, source=ds_3035)\n\nprint(\"Latitude: 78 degrees north, Longitude: 15 degrees east\")\nprint(f\"Natural GC in EPSG:3035 (LAEA): {gc_natural_3035.values[0]:.4f} degrees\")\nprint(f\"Natural GC in EPSG:32632 (UTM 32N): {gc_natural_32632.values[0]:.4f} degrees\")\nprint(\n    \"Effective GC (source dataset 3035 to dataset CRS 32632): \"\n    f\"{gc_effective_from_ds.values[0]:.4f} degrees\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Bin True-North Time Series into a Grid Frame\nWhen a time series is known to be true-north-relative, pass a ``wind_dir_crs``\nto `bwc_from_tswc`. If the TSWC already carries ``attrs[\"wind_dir_crs\"]``,\nWindKit treats the directions as already expressed in that frame and avoids\nrotating them again. A projected CRS rotates directions into that grid\nframe; a geographic CRS retains the true-north frame. WindKit records the\nresolved CRS WKT on the BWC.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "tswc = wk.create_tswc(\n    ds_high_lat,\n    date_range=pd.date_range(\"2001-01-01\", periods=24, freq=\"h\"),\n)\n\n# Mark the input as already expressed in EPSG:32632 grid north.\ntswc_marked = tswc.assign_attrs(wind_dir_crs=32632)\n\nbwc_marked = wk.bwc_from_tswc(tswc_marked)\nbwc_same_frame = wk.bwc_from_tswc(tswc_marked, wind_dir_crs=32632)\nbwc_other_frame = wk.bwc_from_tswc(tswc_marked, wind_dir_crs=3035)\n\nprint(bwc_marked.attrs[\"wind_dir_crs\"])\nprint(bwc_same_frame.attrs[\"wind_dir_crs\"])\nprint(bwc_other_frame.attrs[\"wind_dir_crs\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## When to Use This\nUse :py:func:`~windkit.spatial.grid_convergence` when reasoning about CRS frames, diagnosing map\nprojection behavior, or checking direction reference frames. If you are using\nPyWAsP wind climates, see the PyWAsP meridian-convergence tutorial for how\nthese CRS rotations affect BWC, GWC, and WWC workflows.\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.14.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}