xrspatial.reproject.reproject#

xrspatial.reproject.reproject(raster, target_crs, *, source_crs=None, resolution=None, bounds=None, width=None, height=None, resampling='bilinear', nodata=None, transform_precision=16, chunk_size=None, name=None, max_memory=None, source_vertical_crs=None, target_vertical_crs=None, bounds_policy='auto', src_vertical_crs=<object object>, tgt_vertical_crs=<object object>)[source]#

Reproject a raster DataArray to a new coordinate reference system.

Supports numpy, cupy, dask+numpy, and dask+cupy backends. For dask inputs, the computation is fully lazy: each output chunk independently reads only the source pixels it needs.

Parameters:

  • raster (xr.DataArray) – Input raster with y/x coordinates.

  • target_crs – Target CRS in any format accepted by pyproj.CRS().

  • source_crs (optional) – Source CRS. Auto-detected from raster if None.

  • resolution (float or (float, float) or None) – Output pixel size in target CRS units.

  • bounds ((left, bottom, right, top) or None) – Explicit output extent in target CRS.

  • width (int or None) – Explicit output grid dimensions.

  • height (int or None) – Explicit output grid dimensions.

  • resampling (str) – One of ‘nearest’, ‘bilinear’, ‘cubic’.

  • nodata (float or None) – Nodata value. Auto-detected if None. For integer input dtypes, an explicit value that does not fit the dtype range raises ValueError (e.g. nodata=-9999 with a uint8 raster). Attrs/rioxarray-derived out-of-range values emit a UserWarning and fall back to dtype.min for signed or dtype.max for unsigned so legacy files still load (#2572).

  • transform_precision (int) – Control-grid subdivisions for the coordinate transform (default 16). Higher values increase accuracy at the cost of more pyproj calls. Set to 0 for exact per-pixel transforms matching GDAL/rasterio.

  • chunk_size (int or (int, int) or None) – Output chunk size for dask. If None, defaults to 512 for the standard dask path and 2048 for the in-memory streaming and dask+cupy paths (chosen to amortize kernel launch overhead).

  • name (str or None) – Name for the output DataArray.

  • max_memory (int or str or None) – Maximum memory budget for the reprojection working set. Accepts bytes (int) or human-readable strings like '4GB', '512MB'. Controls how many output tiles are processed in parallel for large-dataset streaming mode. Default None uses 1GB. Has no effect for small datasets that fit in memory.

  • source_vertical_crs (str or None) –

    Source vertical datum for height values. One of:

    • 'EGM96' – orthometric heights relative to EGM96 geoid (MSL)

    • 'EGM2008' – orthometric heights relative to EGM2008 geoid

    • 'ellipsoidal' – heights relative to the WGS84 ellipsoid

    • None – no vertical transformation (default)

  • target_vertical_crs (str or None) – Target vertical datum. Same options as source_vertical_crs. Both must be set to trigger a vertical transformation.

  • src_vertical_crs (str or None) – Deprecated alias for source_vertical_crs. Passing it emits a DeprecationWarning.

  • tgt_vertical_crs (str or None) – Deprecated alias for target_vertical_crs. Passing it emits a DeprecationWarning.

  • 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:

The output attrs['crs'] is in WKT format. Whenever target_vertical_crs is set, attrs['vertical_crs'] records the target vertical datum’s EPSG code (5773 for EGM96, 3855 for EGM2008, 4979 for ellipsoidal WGS84) to match the convention used by xrspatial.geotiff. The friendly string token ('EGM96' etc.) is preserved under attrs['vertical_datum']. Both attrs are written even when no shift is applied (e.g. when source_vertical_crs equals target_vertical_crs, or when only the target is given), so the output’s vertical reference is always explicit.

The output y coordinate is always emitted in descending order (top-down, north-up) and the output x coordinate is always emitted in ascending order (left-to-right) regardless of the input directions. This matches the standard raster convention and the output of common GIS libraries. Inputs with descending x are detected from the x coordinate values and handled the same way as descending y: the pixel-index mapping is mirrored so the output values stay correct.

Non-spatial coords from the input (such as a scalar time coord or a non-dimension coord that is not aligned to the spatial dims) are carried through to the output. Coords that are aligned to the input y or x dims are dropped because their values do not apply to the rebuilt grid.

Return type:

xr.DataArray

Examples

>>> import xarray as xr
>>> import numpy as np
>>> from xrspatial.reproject import reproject
>>> raster = xr.DataArray(
...     np.random.rand(64, 64),
...     dims=['y', 'x'],
...     coords={'y': np.linspace(50, 40, 64),
...             'x': np.linspace(-5, 5, 64)},
...     attrs={'crs': 'EPSG:4326'},
... )
>>> result = reproject(raster, 'EPSG:3857')
>>> result.attrs['crs'].startswith(('PROJCRS', 'PROJCS'))
True