Example - Transform Bounds

The rio.transform_bounds() method allows you to correctly estimate the bounds of your raster in a different CRS without needing to re-project it. If you simply calculate the bounds by transforming the bounds, there are often situations when this is incorrect due to nonlinear transformations.

import pyproj
import rioxarray # for the extension to load
import xarray
from shapely.geometry import box

import matplotlib.pyplot as plt

%matplotlib inline

xds = xarray.open_dataarray("../../test/test_data/input/MODIS_ARRAY.nc")
transformer = pyproj.Transformer.from_crs(xds.rio.crs, "EPSG:4326", always_xy=True)

Original Raster & Bounds

ax = plt.subplot()
xds.plot(ax=ax)
ax.plot(
    *box(*xds.rio.bounds()).exterior.xy,
    color="red",
    linewidth=3,
)

[<matplotlib.lines.Line2D at 0x7f324f4456c0>]

Determine bounds of re-projected raster

The rio.transform_bounds() method allows you to safely convert a bounding box into another projection taking into account the effects of nonlinear transformations.

reprojected_raster = xds.rio.reproject("EPSG:4326")

Boundary calculated from the re-projected raster (inefficient)

This is the benchmark. However, this method is computationally inefficient. So, if you don’t need to re-project, rio.transform_bounds() is a more efficent method.

reprojected_raster_box = box(*reprojected_raster.rio.bounds())

ax = plt.subplot()
reprojected_raster.plot(ax=ax)
ax.plot(
    *reprojected_raster_box.exterior.xy,
    color="red",
    linewidth=3,
)

[<matplotlib.lines.Line2D at 0x7f324735e4d0>]

Boundary calculated from original corners (incorrect)

Directly transforming the corners is an incorrect method to calculate the new boundary.

transform_box = box(*transformer.transform(*xds.rio.bounds()))

ax = plt.subplot()
reprojected_raster.plot(ax=ax)
ax.plot(
    *transform_box.exterior.xy,
    color="red",
    linewidth=3,
)

[<matplotlib.lines.Line2D at 0x7f3245a3bee0>]

Boundary calculates using transform_bounds

rio.transform_bounds() is both computationally efficient and a correct method for calculating the bounds of your raster in the new projection.

transform_bounds_box = box(*xds.rio.transform_bounds("EPSG:4326"))

ax = plt.subplot()
reprojected_raster.plot(ax=ax)
ax.plot(
    *transform_bounds_box.exterior.xy,
    color="red",
    linewidth=3,
)

[<matplotlib.lines.Line2D at 0x7f32459316c0>]

As seen below, this is equivalent to the Transformer.transform_bounds method in pyproj:

pyproj_transform_bounds_box = box(*transformer.transform_bounds(*xds.rio.bounds()))

ax = plt.subplot()
reprojected_raster.plot(ax=ax)
ax.plot(
    *transform_bounds_box.exterior.xy,
    color="red",
    linewidth=3,
)

[<matplotlib.lines.Line2D at 0x7f324580d120>]