reproject: add bounds_policy parameter (#2187)

docs/source/reference/index.rst

@@ -23,6 +23,7 @@ Reference
    multispectral
    pathfinding
    proximity
+   reproject
    resample
    surface
    terrain_metrics

docs/source/reference/reproject.rst

@@ -0,0 +1,38 @@
+..  _reference.reproject:
+
+*********
+Reproject
+*********
+
+Reproject a raster to a new coordinate reference system, and merge
+multiple rasters into a unified output grid.
+
+Reproject
+=========
+.. autosummary::
+    :toctree: _autosummary
+
+    xrspatial.reproject.reproject
+    xrspatial.reproject.merge
+
+Bounds policy
+=============
+
+The ``bounds_policy`` parameter on :func:`xrspatial.reproject.reproject`
+and :func:`xrspatial.reproject.merge` controls how the output extent is
+derived from the input extent when no explicit ``bounds`` is supplied.
+The available options are:
+
+- ``"auto"`` (default): clamp geographic source bounds inward by
+  0.01 deg from +/-180 longitude and +/-90 latitude, and fall back to
+  the 2nd/98th percentile of a dense interior grid when the projected
+  extent is more than 50x the source extent. Matches the historical
+  behaviour. Emits a ``UserWarning`` when the clamp or percentile
+  fallback actually alters the bounds.
+- ``"raw"``: use the true projected extent of the source corners and
+  edges. No clamp, no percentile. Pass this when downstream tooling
+  needs the exact projection of the input footprint and you are willing
+  to handle large outputs near projection singularities.
+- ``"clamp"``: apply only the geographic-bounds clamp.
+- ``"percentile"``: always use the 2/98 percentile bounds, even when
+  the projected extent looks well-behaved.

examples/user_guide/55_Reproject_Bounds_Policy.ipynb

@@ -0,0 +1,294 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Xarray-Spatial Reproject: bounds policy at projection singularities\n",
+    "\n",
+    "Reprojecting a raster near the antimeridian or a pole can produce huge or undefined output extents. xrspatial.reproject ships heuristics that crop the output to a sane size, but those heuristics fire silently and can discard real data. The `bounds_policy` parameter lets you pick the trade-off explicitly."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### What you'll build\n",
+    "\n",
+    "1. Project a global geographic raster to Web Mercator under each of the four policies.\n",
+    "2. Compare the resulting output bounds side by side.\n",
+    "3. Surface the `UserWarning` emitted under the default `auto` policy.\n",
+    "4. Walk through the geographic clamp on a near-antimeridian raster.\n",
+    "\n",
+    "![preview](images/55_reproject_bounds_policy_preview.png)\n",
+    "\n",
+    "- [Setup](#Setup)\n",
+    "- [The four policies](#The-four-policies)\n",
+    "- [Auto policy emits a warning](#Auto-policy-emits-a-warning)\n",
+    "- [Geographic clamp on a near-antimeridian raster](#Geographic-clamp-on-a-near-antimeridian-raster)\n",
+    "- [References](#References)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Standard imports. The reproject function is the only thing we need from xrspatial."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import warnings\n",
+    "\n",
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "from xrspatial.reproject import reproject"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Setup\n",
+    "\n",
+    "Build a near-global geographic raster (EPSG:4326). The longitude range stops just short of +/-180 so we still have a meaningful projected extent, and the latitudes climb to +/-85 to put us near the Web Mercator polar singularity."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.random.seed(42)\n",
+    "lats = np.linspace(85, -85, 80)\n",
+    "lons = np.linspace(-178, 178, 160)\n",
+    "yy, xx = np.meshgrid(lats, lons, indexing='ij')\n",
+    "data = (np.sin(np.deg2rad(xx) * 2) * np.cos(np.deg2rad(yy)) + 1.0) / 2.0\n",
+    "\n",
+    "raster = xr.DataArray(\n",
+    "    data.astype(np.float32),\n",
+    "    dims=['y', 'x'],\n",
+    "    coords={'y': lats, 'x': lons},\n",
+    "    attrs={'crs': 'EPSG:4326'},\n",
+    ")\n",
+    "raster"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The source values are smooth east-west bands modulated by latitude. The exact pattern is incidental. What matters is the spatial extent."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "raster.plot.imshow(size=4, aspect=2, cmap='viridis', add_colorbar=False)\n",
+    "plt.title('Source raster (EPSG:4326)')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The four policies\n",
+    "\n",
+    "Run the same reprojection four times, once per policy. Suppress the warnings for now and just compare the resulting output bounds and shapes."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "policies = ['auto', 'raw', 'clamp', 'percentile']\n",
+    "results = {}\n",
+    "for policy in policies:\n",
+    "    with warnings.catch_warnings():\n",
+    "        warnings.simplefilter('ignore', UserWarning)\n",
+    "        out = reproject(raster, 'EPSG:3857', resolution=2e5, bounds_policy=policy)\n",
+    "    results[policy] = out\n",
+    "\n",
+    "import pandas as pd\n",
+    "rows = []\n",
+    "for policy, out in results.items():\n",
+    "    x = out.coords['x'].values\n",
+    "    y = out.coords['y'].values\n",
+    "    rows.append({\n",
+    "        'policy': policy,\n",
+    "        'shape': out.shape,\n",
+    "        'x_min': float(x.min()),\n",
+    "        'x_max': float(x.max()),\n",
+    "        'y_min': float(y.min()),\n",
+    "        'y_max': float(y.max()),\n",
+    "    })\n",
+    "pd.DataFrame(rows)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "`auto` and `percentile` produce the smallest extents because the 2/98 percentile fallback throws away the extremes near the poles. `raw` keeps every projected pixel and ends up the largest. `clamp` sits in between because it trims the source extent inward by 0.01 deg before projecting, but it still keeps all the projected output."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(2, 2, figsize=(12, 9))\n",
+    "for ax, policy in zip(axes.flat, policies):\n",
+    "    out = results[policy]\n",
+    "    out.plot.imshow(ax=ax, cmap='viridis', add_colorbar=False)\n",
+    "    ax.set_title(f\"bounds_policy={policy!r}  shape={out.shape}\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<div class=\"alert alert-block alert-warning\">\n",
+    "<b>Picking a policy.</b> Use <code>raw</code> when you need the true projected extent of your input footprint, for example when mosaicking many tiles or computing area statistics. Use <code>auto</code> (or <code>percentile</code>) when you'd rather lose a few pixels near a singularity than allocate a giant output grid.\n",
+    "</div>"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Auto policy emits a warning\n",
+    "\n",
+    "Under `auto`, when the percentile fallback or the geographic clamp actually alters the bounds, reproject emits a `UserWarning`. The warning names the policy and reports the per-side delta vs the raw projected bounds."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "with warnings.catch_warnings(record=True) as caught:\n",
+    "    warnings.simplefilter('always')\n",
+    "    reproject(raster, 'EPSG:3857', resolution=2e5, bounds_policy='auto')\n",
+    "\n",
+    "for w in caught:\n",
+    "    if issubclass(w.category, UserWarning) and 'bounds_policy' in str(w.message):\n",
+    "        print(w.message)\n",
+    "        print('---')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If you intentionally want the silent cropping back, suppress the warning:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "with warnings.catch_warnings():\n",
+    "    warnings.filterwarnings(\n",
+    "        'ignore', category=UserWarning, message='.*bounds_policy.*'\n",
+    "    )\n",
+    "    out = reproject(raster, 'EPSG:3857', resolution=2e5)\n",
+    "    print(f'output shape: {out.shape}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Geographic clamp on a near-antimeridian raster\n",
+    "\n",
+    "When the source CRS is geographic and the source extent touches +/-180 longitude, the default policy trims it inward by 0.01 deg before projecting. This avoids infinities at the antimeridian, but it also trims a sliver of real data. Surface the trim with the `clamp` policy on a synthetic raster whose extent runs to exactly +/-180 longitude."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "edge_raster = xr.DataArray(\n",
+    "    np.random.RandomState(0).rand(40, 80).astype(np.float32),\n",
+    "    dims=['y', 'x'],\n",
+    "    coords={'y': np.linspace(60, -60, 40),\n",
+    "            'x': np.linspace(-180, 180, 80)},\n",
+    "    attrs={'crs': 'EPSG:4326'},\n",
+    ")\n",
+    "\n",
+    "with warnings.catch_warnings(record=True) as caught:\n",
+    "    warnings.simplefilter('always')\n",
+    "    out_clamp = reproject(\n",
+    "        edge_raster, 'EPSG:3857', resolution=2e5, bounds_policy='clamp'\n",
+    "    )\n",
+    "\n",
+    "for w in caught:\n",
+    "    if issubclass(w.category, UserWarning) and 'clamp' in str(w.message):\n",
+    "        print(w.message)\n",
+    "        print('---')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The warning names the original source extent and the trimmed one. If you need the antimeridian-touching data preserved, switch to `bounds_policy='raw'` and handle any projection infinities downstream."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<div class=\"alert alert-block alert-info\">\n",
+    "<b>Antimeridian wrap.</b> Geographic rasters that genuinely span the antimeridian are still tricky. Even with <code>raw</code>, single-pass reprojection cannot reconstruct a discontinuous footprint. Consider splitting the raster on the antimeridian and reprojecting each half separately, then merging.\n",
+    "</div>"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## References\n",
+    "\n",
+    "- xarray-spatial reproject reference: https://xarray-spatial.org/reference/reproject.html\n",
+    "- Web Mercator (EPSG:3857) polar truncation: https://epsg.io/3857\n",
+    "- Issue #2187: bounds heuristics audit"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.x"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}