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.

[1]:
import pyproj
import rioxarray # for the extension to load
import xarray
from shapely.geometry import box

import matplotlib.pyplot as plt

%matplotlib inline
[2]:
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

[3]:
ax = plt.subplot()
xds.plot(ax=ax)
ax.plot(
    *box(*xds.rio.bounds()).exterior.xy,
    color="red",
    linewidth=3,
)
[3]:
[<matplotlib.lines.Line2D at 0x7f324f4456c0>]
../_images/examples_transform_bounds_4_1.png

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.

[4]:
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.

[5]:
reprojected_raster_box = box(*reprojected_raster.rio.bounds())
[6]:
ax = plt.subplot()
reprojected_raster.plot(ax=ax)
ax.plot(
    *reprojected_raster_box.exterior.xy,
    color="red",
    linewidth=3,
)
[6]:
[<matplotlib.lines.Line2D at 0x7f324735e4d0>]
../_images/examples_transform_bounds_9_1.png

Boundary calculated from original corners (incorrect)

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

[7]:
transform_box = box(*transformer.transform(*xds.rio.bounds()))
[8]:
ax = plt.subplot()
reprojected_raster.plot(ax=ax)
ax.plot(
    *transform_box.exterior.xy,
    color="red",
    linewidth=3,
)
[8]:
[<matplotlib.lines.Line2D at 0x7f3245a3bee0>]
../_images/examples_transform_bounds_12_1.png

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.

[9]:
transform_bounds_box = box(*xds.rio.transform_bounds("EPSG:4326"))
[10]:
ax = plt.subplot()
reprojected_raster.plot(ax=ax)
ax.plot(
    *transform_bounds_box.exterior.xy,
    color="red",
    linewidth=3,
)
[10]:
[<matplotlib.lines.Line2D at 0x7f32459316c0>]
../_images/examples_transform_bounds_15_1.png

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

[11]:
pyproj_transform_bounds_box = box(*transformer.transform_bounds(*xds.rio.bounds()))
[12]:
ax = plt.subplot()
reprojected_raster.plot(ax=ax)
ax.plot(
    *transform_bounds_box.exterior.xy,
    color="red",
    linewidth=3,
)
[12]:
[<matplotlib.lines.Line2D at 0x7f324580d120>]
../_images/examples_transform_bounds_18_1.png