diff --git a/docs/source/reference/index.rst b/docs/source/reference/index.rst
index d851b0523..fa9315866 100644
--- a/docs/source/reference/index.rst
+++ b/docs/source/reference/index.rst
@@ -23,6 +23,7 @@ Reference
multispectral
pathfinding
proximity
+ reproject
resample
surface
terrain_metrics
diff --git a/docs/source/reference/reproject.rst b/docs/source/reference/reproject.rst
new file mode 100644
index 000000000..1d2932431
--- /dev/null
+++ b/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.
diff --git a/examples/user_guide/55_Reproject_Bounds_Policy.ipynb b/examples/user_guide/55_Reproject_Bounds_Policy.ipynb
new file mode 100644
index 000000000..cb4090a6e
--- /dev/null
+++ b/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",
+ "\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": [
+ "\n",
+ "Picking a policy. Use raw when you need the true projected extent of your input footprint, for example when mosaicking many tiles or computing area statistics. Use auto (or percentile) when you'd rather lose a few pixels near a singularity than allocate a giant output grid.\n",
+ "
"
+ ]
+ },
+ {
+ "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": [
+ "\n",
+ "Antimeridian wrap. Geographic rasters that genuinely span the antimeridian are still tricky. Even with raw, single-pass reprojection cannot reconstruct a discontinuous footprint. Consider splitting the raster on the antimeridian and reprojecting each half separately, then merging.\n",
+ "
"
+ ]
+ },
+ {
+ "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
+}
diff --git a/examples/user_guide/images/55_reproject_bounds_policy_preview.png b/examples/user_guide/images/55_reproject_bounds_policy_preview.png
new file mode 100644
index 000000000..2ec410110
Binary files /dev/null and b/examples/user_guide/images/55_reproject_bounds_policy_preview.png differ
diff --git a/xrspatial/reproject/__init__.py b/xrspatial/reproject/__init__.py
index bee09534d..4922e1987 100644
--- a/xrspatial/reproject/__init__.py
+++ b/xrspatial/reproject/__init__.py
@@ -655,6 +655,7 @@ def reproject(
max_memory=None,
src_vertical_crs=None,
tgt_vertical_crs=None,
+ bounds_policy="auto",
):
"""Reproject a raster DataArray to a new coordinate reference system.
@@ -706,6 +707,34 @@ def reproject(
tgt_vertical_crs : str or None
Target vertical datum. Same options as *src_vertical_crs*.
Both must be set to trigger a vertical transformation.
+ bounds_policy : {"auto", "raw", "clamp", "percentile"}, default "auto"
+ How to derive the output extent from the source extent when
+ ``bounds`` is not supplied. Only relevant when projecting near a
+ singularity (antimeridian, pole, projection edge):
+
+ - ``"raw"``: use the true projected extent of the source corners
+ and edges. No clamp, no percentile, no heuristic. The output
+ may be very large if the input straddles a projection
+ singularity. Use this when you want a true projection of the
+ source extent.
+ - ``"clamp"``: trim geographic source bounds inward by 0.01 deg
+ from +/-180 longitude and +/-90 latitude before projecting.
+ Avoids infinities at singularities. No percentile fallback.
+ No-op on projected source CRSes (UTM, Mercator, etc.) since
+ the clamp only applies in degrees.
+ - ``"percentile"``: project a dense interior grid of the source
+ extent and use the 2nd/98th percentiles of the result as the
+ output bounds. Rejects projection outliers at the cost of
+ trimming valid pixels.
+ - ``"auto"`` (default): apply ``"clamp"`` for geographic source
+ CRSes and fall back to ``"percentile"`` when the projected
+ extent is more than 50x the source extent. Matches the
+ historical behaviour.
+
+ When ``"auto"``, ``"clamp"``, or ``"percentile"`` actually alters
+ the bounds, a ``UserWarning`` is emitted naming the policy and
+ reporting the per-side delta versus the raw projected bounds.
+ Filter with ``warnings.filterwarnings`` if the crop is intentional.
Returns
-------
@@ -773,6 +802,9 @@ def reproject(
_validate_resampling(resampling)
+ from ._grid import _validate_bounds_policy
+ _validate_bounds_policy(bounds_policy, func_name='reproject')
+
# Reject unknown vertical-datum tokens at the API boundary so we never
# write None into attrs['vertical_crs'] for typos / unsupported values.
for _name, _val in (('src_vertical_crs', src_vertical_crs),
@@ -830,6 +862,7 @@ def reproject(
src_bounds, src_shape, src_crs, tgt_crs,
resolution=resolution, bounds=bounds,
width=width, height=height,
+ bounds_policy=bounds_policy,
)
out_bounds = grid['bounds']
out_shape = grid['shape']
@@ -1850,6 +1883,7 @@ def merge(
chunk_size=None,
transform_precision=16,
name=None,
+ bounds_policy="auto",
):
"""Merge multiple rasters into a single mosaic.
@@ -1883,6 +1917,10 @@ def merge(
Set to 0 for exact per-pixel transforms matching GDAL/rasterio.
name : str or None
Name for the output DataArray.
+ bounds_policy : {"auto", "raw", "clamp", "percentile"}, default "auto"
+ How to derive the unified output extent from the input rasters
+ when ``bounds`` is not supplied. See :func:`reproject` for the
+ full description of each option. Ignored when ``bounds`` is given.
Returns
-------
@@ -1953,6 +1991,9 @@ def merge(
_validate_resampling(resampling)
_validate_strategy(strategy)
+ from ._grid import _validate_bounds_policy
+ _validate_bounds_policy(bounds_policy, func_name='merge')
+
# Resolve target CRS
tgt_crs = _resolve_crs(target_crs)
if tgt_crs is None:
@@ -2002,15 +2043,43 @@ def merge(
# Compute unified output grid
if bounds is None:
- # Union of all raster bounds in target CRS
+ # Union of all raster bounds in target CRS.
+ # Capture bounds_policy warnings during the per-input gather and
+ # emit a single deduplicated message afterwards. Otherwise a
+ # mosaic of N near-antimeridian rasters yields N identical
+ # warnings, which is noise the caller cannot act on individually.
+ import warnings as _warnings
all_bounds = []
- for info in raster_infos:
- grid = _compute_output_grid(
- info['src_bounds'], info['src_shape'],
- info['src_crs'], tgt_crs,
- resolution=resolution,
+ with _warnings.catch_warnings(record=True) as _caught:
+ _warnings.simplefilter('always', UserWarning)
+ for info in raster_infos:
+ grid = _compute_output_grid(
+ info['src_bounds'], info['src_shape'],
+ info['src_crs'], tgt_crs,
+ resolution=resolution,
+ bounds_policy=bounds_policy,
+ )
+ all_bounds.append(grid['bounds'])
+ _policy_msgs = [
+ str(w.message) for w in _caught
+ if issubclass(w.category, UserWarning)
+ and 'bounds_policy' in str(w.message)
+ ]
+ if _policy_msgs:
+ _unique = list(dict.fromkeys(_policy_msgs))
+ _summary = (
+ f"merge: bounds_policy={bounds_policy!r} altered the "
+ f"projected extent for {len(_policy_msgs)} input "
+ f"raster(s); {len(_unique)} unique trigger(s). "
+ f"First trigger: {_unique[0]}"
)
- all_bounds.append(grid['bounds'])
+ _warnings.warn(_summary, UserWarning, stacklevel=2)
+ # Re-emit non-bounds_policy warnings that we captured.
+ for _w in _caught:
+ if 'bounds_policy' not in str(_w.message):
+ _warnings.warn_explicit(
+ _w.message, _w.category, _w.filename, _w.lineno,
+ )
left = min(b[0] for b in all_bounds)
bottom = min(b[1] for b in all_bounds)
right = max(b[2] for b in all_bounds)
@@ -2025,6 +2094,7 @@ def merge(
info0['src_bounds'], info0['src_shape'],
info0['src_crs'], tgt_crs,
resolution=resolution, bounds=merged_bounds,
+ bounds_policy=bounds_policy,
)
out_bounds = grid['bounds']
out_shape = grid['shape']
diff --git a/xrspatial/reproject/_grid.py b/xrspatial/reproject/_grid.py
index 5eeea1171..312db7126 100644
--- a/xrspatial/reproject/_grid.py
+++ b/xrspatial/reproject/_grid.py
@@ -1,9 +1,33 @@
"""Output grid computation and chunk layout for reprojection."""
from __future__ import annotations
+import warnings
+
import numpy as np
+_VALID_BOUNDS_POLICIES = ("auto", "raw", "clamp", "percentile")
+
+
+def _validate_bounds_policy(bounds_policy, func_name):
+ """Reject unknown bounds_policy tokens at the API boundary."""
+ if bounds_policy not in _VALID_BOUNDS_POLICIES:
+ raise ValueError(
+ f"{func_name}(): bounds_policy must be one of "
+ f"{list(_VALID_BOUNDS_POLICIES)}, got {bounds_policy!r}"
+ )
+
+
+def _bounds_delta(raw, adjusted):
+ """Return per-side deltas (adjusted - raw) for a bounds 4-tuple."""
+ return (
+ adjusted[0] - raw[0],
+ adjusted[1] - raw[1],
+ adjusted[2] - raw[2],
+ adjusted[3] - raw[3],
+ )
+
+
def _validate_grid_params(*, resolution, bounds, width, height,
transform_precision, func_name):
"""Range-check user-supplied grid parameters."""
@@ -121,7 +145,8 @@ def _transform_boundary(source_crs, target_crs, xs, ys):
def _compute_output_grid(source_bounds, source_shape, source_crs, target_crs,
- resolution=None, bounds=None, width=None, height=None):
+ resolution=None, bounds=None, width=None, height=None,
+ bounds_policy="auto"):
"""Compute the output raster grid parameters.
Parameters
@@ -138,28 +163,49 @@ def _compute_output_grid(source_bounds, source_shape, source_crs, target_crs,
Explicit (left, bottom, right, top) in target CRS.
width, height : int or None
Explicit output dimensions.
+ bounds_policy : {"auto", "raw", "clamp", "percentile"}
+ How to handle output bounds near projection singularities.
+ See :func:`reproject` for the full description. Ignored when
+ ``bounds`` is supplied.
Returns
-------
dict with keys: bounds, shape, res_x, res_y
"""
+ _validate_bounds_policy(bounds_policy, func_name='_compute_output_grid')
+
if bounds is None:
# Transform source corners and edges to target CRS
src_left, src_bottom, src_right, src_top = source_bounds
+ src_left_raw, src_bottom_raw = src_left, src_bottom
+ src_right_raw, src_top_raw = src_right, src_top
# Clamp geographic coordinates away from singularities.
# Exactly +/-180 longitude produces inf in many projections;
# latitudes beyond +/-90 are invalid.
- if source_crs.is_geographic:
+ clamp_applied = False
+ if source_crs.is_geographic and bounds_policy in ("auto", "clamp"):
clamp = 0.01 # degrees
- src_left = max(src_left, -180 + clamp)
- src_right = min(src_right, 180 - clamp)
- src_bottom = max(src_bottom, -90 + clamp)
- src_top = min(src_top, 90 - clamp)
+ new_left = max(src_left, -180 + clamp)
+ new_right = min(src_right, 180 - clamp)
+ new_bottom = max(src_bottom, -90 + clamp)
+ new_top = min(src_top, 90 - clamp)
+ clamp_applied = (
+ new_left != src_left or new_right != src_right
+ or new_bottom != src_bottom or new_top != src_top
+ )
+ src_left, src_right = new_left, new_right
+ src_bottom, src_top = new_bottom, new_top
+ # Defensive: if the input was so narrow that the 0.02 deg
+ # clamp inverted the range, restore the originals. Upstream
+ # _validate_grid_params already rejects degenerate bounds,
+ # so this branch should be unreachable in practice.
if src_left >= src_right:
src_left, src_right = source_bounds[0], source_bounds[2]
+ clamp_applied = False
if src_bottom >= src_top:
src_bottom, src_top = source_bounds[1], source_bounds[3]
+ clamp_applied = False
# Sample edges densely plus an interior grid so that
# projections with curvature (e.g. UTM near zone edges)
@@ -185,6 +231,9 @@ def _compute_output_grid(source_bounds, source_shape, source_crs, target_crs,
ixx, iyy = np.meshgrid(ix, iy)
xs = np.concatenate([edge_xs, ixx.ravel()])
ys = np.concatenate([edge_ys, iyy.ravel()])
+ # Under policy='raw' the clamp branch above is skipped, so the
+ # edge samples already start from the original src_*_raw values.
+ # No extra corner injection needed.
tx, ty = _transform_boundary(source_crs, target_crs, xs, ys)
tx = np.asarray(tx)
ty = np.asarray(ty)
@@ -213,7 +262,13 @@ def _compute_output_grid(source_bounds, source_shape, source_crs, target_crs,
# to a dense interior grid to get a tighter bounding box.
ratio_x = out_w_crs / (abs(src_w_crs) + 1e-30)
ratio_y = out_h_crs / (abs(src_h_crs) + 1e-30)
- if ratio_x > 50 or ratio_y > 50:
+ auto_blowup = (ratio_x > 50 or ratio_y > 50)
+ use_percentile = (
+ bounds_policy == "percentile"
+ or (bounds_policy == "auto" and auto_blowup)
+ )
+ percentile_applied = False
+ if use_percentile:
# Sample a dense interior grid instead
n_dense = 51
ix = np.linspace(src_left, src_right, n_dense)
@@ -236,11 +291,44 @@ def _compute_output_grid(source_bounds, source_shape, source_crs, target_crs,
float(np.percentile(itx, q_hi)),
float(np.percentile(ity, q_hi)),
)
+ percentile_applied = True
else:
bounds = raw_bounds
else:
bounds = raw_bounds
+ # Emit a warning when the policy actually altered the bounds vs.
+ # the raw projection. Silent cropping is the bug this knob fixes.
+ if percentile_applied and bounds != raw_bounds:
+ delta = _bounds_delta(raw_bounds, bounds)
+ trigger = "auto blow-up heuristic" if (
+ bounds_policy == "auto" and auto_blowup
+ ) else "percentile policy"
+ warnings.warn(
+ f"reproject bounds_policy={bounds_policy!r}: "
+ f"{trigger} replaced raw projected bounds "
+ f"{raw_bounds} with 2/98 percentile bounds {bounds}; "
+ f"per-side delta (adjusted - raw) = {delta}. Pass "
+ f"bounds_policy='raw' to disable this crop.",
+ UserWarning,
+ stacklevel=2,
+ )
+ elif clamp_applied and not percentile_applied:
+ # The clamp affects which input points get sampled, not the
+ # output bounds directly. Still surface it so callers know
+ # we trimmed source coordinates before projecting.
+ warnings.warn(
+ f"reproject bounds_policy={bounds_policy!r}: "
+ f"geographic clamp trimmed source extent from "
+ f"({src_left_raw}, {src_bottom_raw}, "
+ f"{src_right_raw}, {src_top_raw}) to "
+ f"({src_left}, {src_bottom}, {src_right}, {src_top}) "
+ f"before projecting. Pass bounds_policy='raw' to "
+ f"disable this clamp.",
+ UserWarning,
+ stacklevel=2,
+ )
+
left, bottom, right, top = bounds
# Determine resolution
diff --git a/xrspatial/tests/test_reproject.py b/xrspatial/tests/test_reproject.py
index a2d4ab6ca..ef934372e 100644
--- a/xrspatial/tests/test_reproject.py
+++ b/xrspatial/tests/test_reproject.py
@@ -4655,6 +4655,345 @@ def test_band_first_dask_cupy(self):
assert computed.shape[0] == 3
+# ---------------------------------------------------------------------------
+# Issue #2187: bounds_policy parameter
+# ---------------------------------------------------------------------------
+
+class TestBoundsPolicy:
+ """reproject(): bounds_policy controls the bounds-derivation heuristics.
+
+ Without the policy knob, _compute_output_grid silently clamps
+ geographic bounds and falls back to 2/98 percentile bounds when the
+ projected extent blows up. These tests pin the four policy options:
+ auto (default, current behaviour with warnings), raw (no heuristic),
+ clamp (geographic clamp only), and percentile (force 2/98 fallback).
+ """
+
+ @staticmethod
+ def _global_geographic():
+ """Global lat/lon raster that triggers the polar / antimeridian
+ blow-up when projected to Web Mercator."""
+ data = np.random.RandomState(0).rand(50, 100).astype(np.float32)
+ return xr.DataArray(
+ data, dims=['y', 'x'],
+ coords={'y': np.linspace(90, -90, 50),
+ 'x': np.linspace(-180, 180, 100)},
+ attrs={'crs': 'EPSG:4326'},
+ )
+
+ @staticmethod
+ def _benign_geographic():
+ """Mid-latitude raster well away from any singularity."""
+ data = np.random.RandomState(0).rand(32, 32).astype(np.float32)
+ return xr.DataArray(
+ data, dims=['y', 'x'],
+ coords={'y': np.linspace(55, 45, 32),
+ 'x': np.linspace(-5, 5, 32)},
+ attrs={'crs': 'EPSG:4326'},
+ )
+
+ def test_raw_skips_clamp_and_percentile(self):
+ """bounds_policy='raw' returns un-cropped bounds for a blow-up case.
+
+ A global geographic raster projected to Web Mercator hits the
+ polar singularity. Under 'auto' the percentile fallback fires
+ and crops the output extent. Under 'raw' the caller gets the
+ true projected extent of the corners/edges.
+ """
+ from xrspatial.reproject import reproject
+
+ r = self._global_geographic()
+ import warnings
+ with warnings.catch_warnings():
+ warnings.simplefilter('ignore', UserWarning)
+ auto_result = reproject(r, 'EPSG:3857', bounds_policy='auto')
+ raw_result = reproject(r, 'EPSG:3857', bounds_policy='raw')
+
+ auto_x = auto_result.coords['x'].values
+ raw_x = raw_result.coords['x'].values
+ auto_y = auto_result.coords['y'].values
+ raw_y = raw_result.coords['y'].values
+
+ auto_x_span = auto_x.max() - auto_x.min()
+ raw_x_span = raw_x.max() - raw_x.min()
+ auto_y_span = auto_y.max() - auto_y.min()
+ raw_y_span = raw_y.max() - raw_y.min()
+
+ # Raw should be at least as wide as auto on x, and meaningfully
+ # taller on y (the polar blow-up is the y axis under EPSG:3857).
+ assert raw_x_span >= auto_x_span - 1.0
+ assert raw_y_span > auto_y_span * 1.1, (
+ f"raw y span {raw_y_span} should exceed auto {auto_y_span}"
+ )
+
+ def test_percentile_reproduces_98_2_behaviour(self):
+ """bounds_policy='percentile' matches the previous 2/98 fallback
+ even on inputs that wouldn't trigger the blow-up heuristic."""
+ from xrspatial.reproject import reproject
+ from xrspatial.reproject._grid import _compute_output_grid
+ from xrspatial.reproject._crs_utils import _resolve_crs
+
+ r = self._benign_geographic()
+ src_crs = _resolve_crs('EPSG:4326')
+ tgt_crs = _resolve_crs('EPSG:3857')
+
+ import warnings
+ with warnings.catch_warnings():
+ warnings.simplefilter('ignore', UserWarning)
+ grid_percentile = _compute_output_grid(
+ (-5.0, 45.0, 5.0, 55.0), (32, 32),
+ src_crs, tgt_crs, bounds_policy='percentile',
+ )
+ grid_raw = _compute_output_grid(
+ (-5.0, 45.0, 5.0, 55.0), (32, 32),
+ src_crs, tgt_crs, bounds_policy='raw',
+ )
+
+ # Percentile bounds should be strictly inside raw bounds (or
+ # equal to floating-point precision) for a benign input.
+ pl, pb, pr, pt = grid_percentile['bounds']
+ rl, rb, rr, rt = grid_raw['bounds']
+ assert pl >= rl - 1.0
+ assert pr <= rr + 1.0
+ assert pb >= rb - 1.0
+ assert pt <= rt + 1.0
+
+ def test_warns_when_percentile_fires_under_auto(self):
+ """auto policy emits UserWarning when the 2/98 fallback triggers."""
+ from xrspatial.reproject import reproject
+
+ r = self._global_geographic()
+ import warnings
+ with warnings.catch_warnings(record=True) as w:
+ warnings.simplefilter('always')
+ reproject(r, 'EPSG:3857', bounds_policy='auto')
+
+ matched = [
+ wi for wi in w
+ if issubclass(wi.category, UserWarning)
+ and 'bounds_policy' in str(wi.message)
+ ]
+ assert matched, "expected a bounds_policy warning under auto"
+ assert any('blow-up' in str(wi.message) or 'percentile' in str(wi.message)
+ for wi in matched)
+
+ def test_warns_when_clamp_actually_trims(self):
+ """clamp policy emits UserWarning when source bounds get trimmed."""
+ from xrspatial.reproject._grid import _compute_output_grid
+ from xrspatial.reproject._crs_utils import _resolve_crs
+
+ src_crs = _resolve_crs('EPSG:4326')
+ tgt_crs = _resolve_crs('EPSG:3857')
+
+ import warnings
+ with warnings.catch_warnings(record=True) as w:
+ warnings.simplefilter('always')
+ _compute_output_grid(
+ (-180.0, -90.0, 180.0, 90.0), (50, 100),
+ src_crs, tgt_crs, bounds_policy='clamp',
+ )
+
+ matched = [
+ wi for wi in w
+ if issubclass(wi.category, UserWarning)
+ and 'clamp' in str(wi.message)
+ ]
+ assert matched, "expected a clamp warning on full-globe input"
+
+ def test_no_warning_on_benign_input(self):
+ """auto policy stays silent on inputs that don't trigger heuristics.
+
+ Same-units projections (UTM->UTM) have comparable spans so the
+ blow-up ratio stays below the 50x threshold and the geographic
+ clamp doesn't apply. No warning should fire.
+ """
+ from xrspatial.reproject import reproject
+
+ data = np.random.RandomState(0).rand(32, 32).astype(np.float32)
+ r = xr.DataArray(
+ data, dims=['y', 'x'],
+ coords={'y': np.linspace(5000000, 4000000, 32),
+ 'x': np.linspace(400000, 600000, 32)},
+ attrs={'crs': 'EPSG:32633'},
+ )
+ import warnings
+ with warnings.catch_warnings(record=True) as w:
+ warnings.simplefilter('always')
+ reproject(r, 'EPSG:32632', bounds_policy='auto')
+
+ matched = [
+ wi for wi in w
+ if issubclass(wi.category, UserWarning)
+ and 'bounds_policy' in str(wi.message)
+ ]
+ assert not matched, (
+ f"unexpected bounds_policy warning(s): {[str(m.message) for m in matched]}"
+ )
+
+ def test_invalid_policy_rejected(self):
+ """Unknown bounds_policy tokens raise ValueError at the API boundary."""
+ from xrspatial.reproject import reproject
+
+ r = self._benign_geographic()
+ with pytest.raises(ValueError, match=r"bounds_policy"):
+ reproject(r, 'EPSG:3857', bounds_policy='bogus')
+
+ def test_invalid_policy_rejected_in_merge(self):
+ from xrspatial.reproject import merge
+
+ r = self._benign_geographic()
+ with pytest.raises(ValueError, match=r"bounds_policy"):
+ merge([r], target_crs='EPSG:3857', bounds_policy='bogus')
+
+ def test_explicit_bounds_skips_policy_logic(self):
+ """When the caller passes bounds, the policy heuristics don't run."""
+ from xrspatial.reproject import reproject
+
+ r = self._global_geographic()
+ # Mercator-y bounds chosen well inside the projection envelope.
+ explicit = (-2.0e7, -2.0e7, 2.0e7, 2.0e7)
+ import warnings
+ with warnings.catch_warnings(record=True) as w:
+ warnings.simplefilter('always')
+ out = reproject(
+ r, 'EPSG:3857',
+ bounds=explicit,
+ resolution=2e5,
+ bounds_policy='auto',
+ )
+
+ # No bounds_policy warning fires when bounds are explicit.
+ matched = [
+ wi for wi in w
+ if issubclass(wi.category, UserWarning)
+ and 'bounds_policy' in str(wi.message)
+ ]
+ assert not matched
+ # Output extent reflects the explicit bounds (within one pixel).
+ out_x = out.coords['x'].values
+ out_y = out.coords['y'].values
+ assert abs(out_x.min() - explicit[0]) < 2.5e5
+ assert abs(out_x.max() - explicit[2]) < 2.5e5
+ assert abs(out_y.min() - explicit[1]) < 2.5e5
+ assert abs(out_y.max() - explicit[3]) < 2.5e5
+
+ def test_merge_passes_policy_through(self):
+ """merge() plumbs bounds_policy to _compute_output_grid."""
+ from xrspatial.reproject import merge
+
+ r = self._global_geographic()
+ import warnings
+ with warnings.catch_warnings(record=True) as w:
+ warnings.simplefilter('always')
+ merge([r], target_crs='EPSG:3857', bounds_policy='auto')
+
+ matched = [
+ wi for wi in w
+ if issubclass(wi.category, UserWarning)
+ and 'bounds_policy' in str(wi.message)
+ ]
+ # merge() on a single global geographic raster should also
+ # trigger the percentile fallback warning under auto.
+ assert matched
+
+ @pytest.mark.skipif(not HAS_DASK, reason="dask required")
+ def test_raw_policy_dask_backend(self):
+ """bounds_policy='raw' works with a dask-backed input."""
+ from xrspatial.reproject import reproject
+
+ r = self._benign_geographic()
+ r = r.chunk({'y': 16, 'x': 16})
+ out = reproject(r, 'EPSG:3857', bounds_policy='raw')
+ # Result is also dask-backed (lazy).
+ assert hasattr(out.data, 'dask')
+ # Compute and confirm we got finite output.
+ arr = out.compute()
+ assert np.isfinite(arr.data).any()
+
+ def test_clamp_policy_noop_on_benign_geographic(self):
+ """bounds_policy='clamp' is silent on a mid-latitude geographic
+ input whose extent does not touch +/-180 or +/-90.
+
+ The clamp branch runs but trims nothing, so no warning should
+ fire. This pins the behaviour so a future change that always
+ emits a clamp warning shows up here.
+ """
+ from xrspatial.reproject import reproject
+
+ r = self._benign_geographic()
+ import warnings
+ with warnings.catch_warnings(record=True) as w:
+ warnings.simplefilter('always')
+ reproject(r, 'EPSG:3857', bounds_policy='clamp')
+
+ matched = [
+ wi for wi in w
+ if issubclass(wi.category, UserWarning)
+ and 'bounds_policy' in str(wi.message)
+ ]
+ assert not matched, (
+ f"clamp on benign geographic input should not warn; got "
+ f"{[str(m.message) for m in matched]}"
+ )
+
+ def test_clamp_policy_noop_on_projected_source(self):
+ """bounds_policy='clamp' is a no-op when source CRS is projected.
+
+ The clamp condition is gated on `source_crs.is_geographic`, so
+ a UTM input under 'clamp' should run without trimming or
+ warning regardless of how close to a singularity the extent is.
+ """
+ from xrspatial.reproject import reproject
+
+ data = np.random.RandomState(0).rand(32, 32).astype(np.float32)
+ # UTM-style coords, mid-latitudes.
+ r = xr.DataArray(
+ data, dims=['y', 'x'],
+ coords={'y': np.linspace(5000000, 4000000, 32),
+ 'x': np.linspace(400000, 600000, 32)},
+ attrs={'crs': 'EPSG:32633'},
+ )
+ import warnings
+ with warnings.catch_warnings(record=True) as w:
+ warnings.simplefilter('always')
+ reproject(r, 'EPSG:4326', bounds_policy='clamp')
+
+ matched = [
+ wi for wi in w
+ if issubclass(wi.category, UserWarning)
+ and 'bounds_policy' in str(wi.message)
+ ]
+ assert not matched
+
+ def test_merge_dedupes_per_input_warnings(self):
+ """merge() collapses per-input bounds_policy warnings into one.
+
+ When several inputs all trigger the percentile fallback, the
+ caller should see a single summary warning rather than N
+ near-identical messages.
+ """
+ from xrspatial.reproject import merge
+
+ r = self._global_geographic()
+ import warnings
+ with warnings.catch_warnings(record=True) as w:
+ warnings.simplefilter('always')
+ merge([r, r, r], target_crs='EPSG:3857', bounds_policy='auto')
+
+ matched = [
+ wi for wi in w
+ if issubclass(wi.category, UserWarning)
+ and 'bounds_policy' in str(wi.message)
+ ]
+ # Three identical inputs should yield exactly one summary
+ # warning from merge(), not three.
+ summary = [m for m in matched if 'merge:' in str(m.message)]
+ assert len(summary) == 1, (
+ f"expected one merge summary warning, got "
+ f"{[str(m.message) for m in matched]}"
+ )
+
+
# ---------------------------------------------------------------------------
# Integer dtype nodata handling (#2185)
# ---------------------------------------------------------------------------