help with rotating a data array

[adamancer]

I'm trying to use rioxarray to rotate some AVIRIS data and have run into an issue I can't figure out. The function I've come up with to do this does okay with the rotation, but the resulting shape is smaller than and offset from the footprint from the flight:

I don't really understand affine math, but my best guess is that this has something to do with rio.transform(). The original dataset has a rotation component, but the rotation is missing from the transform generated by the rio.transform() method:

• GeoTransform: [14.26, 9.26, 305827.93, 9.26, -14.26, 5223236.60]
• xda.rio.transform(): [14.26, 0.00, 305827.93, 0.00, -14.26, 5223236.60]

The missing rotation component might explain why the resolution/scale after an additional rotation is wrong, but I'm at a loss for how to fix it. Does anyone know how to get the rotated object to scale properly?

Here's my code:

import earthpy.plot as ep
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
from rasterio.plot import plotting_extent
import rioxarray as rxr
from shapely.geometry import Polygon


def rotation(xobj, angle):
    """Rotates an xarray object created using rioxarray around its center

    Parameters
    ----------
    xobj: xarray.DataArray or xarray.Dataset
        xarry object to be rotated
    angle: float
        angle in degrees to rotate the object

    Returns
    -------
    xarray.DataArray or xarray.Dataset
        the rotated xarray object
    """

    # Set nodata to np.nan if it isn't set
    if xobj.rio.nodata is None:
        xobj = xobj.rio.set_nodata(np.nan, inplace=False)

    # Pad the object based on its diagonal to get a shape that will accomodate
    # any rotation around the center. If the orientation of the data is not
    # the same as the array (for example, if the array has already been
    # rotated), this will overestimate the size of the canvas required. Rows
    # and columns consisting of all NaNs are trimmed later to account for this.
    x1, y1, x2, y2 = (0, 0, xobj.sizes["x"], xobj.sizes["y"])
    width = x2 - x1
    height = y2 - y1
    maxdim = int((width ** 2 + height ** 2) ** 0.5)

    # Calculate pivot based on array size
    pivot = [d // 2 for d in (xobj.sizes["x"], xobj.sizes["y"])]

    # Rotate object and translate it to the center of the new, larger canvas
    transform = xobj.rio.transform()
    transform *= transform.rotation(angle, pivot)
    transform *= transform.translation((width - maxdim) / 2,
                                       (height - maxdim) / 2)

    rotated = xobj.rio.reproject(
        xobj.rio.crs, shape=(maxdim, maxdim), transform=transform
    )

    # Trim all-NaN rows and columns before returning
    return rotated.dropna("x", "all").dropna("y", "all")


# Read raster data for one band extracted from an AVIRIS dataset
fn = "f190808t01p00r07_corr_v111_1293nm.tif"
band = rxr.open_rasterio(fn, masked=True)

# Create a GeoDataFrame of the polygon for the same flight
geom = Polygon((
    (-118.53101, 46.10141),
    (-118.39009, 46.16763),
    (-119.41823, 47.20150),
    (-119.56051, 47.13400),
    (-118.53101, 46.10141)
))
flight_poly = gpd.GeoDataFrame(
    {"name": ["f190808t01p00r07"], "geometry": [geom]},
    crs="EPSG:4326"
).to_crs(band.rio.crs)

# Rotate the band using the custom function above
rtd = rotation(band, 33)

# The earthpy plot functions require a masked np.array and an extent
arr = rtd.values
masked = np.ma.masked_array(arr, np.isnan(arr))
extent = plotting_extent(rtd[0], rtd[0].rio.transform())

# Plot raster and vector data
fig, ax = plt.subplots(figsize=(8, 4))
ep.plot_bands(masked,
              ax=ax,
              extent=extent,
              cbar=False,
              title="Rotated AVIRIS raster\nwith expected footprint")
flight_poly.plot(ax=ax, facecolor="none", edgecolor="lightgray", linewidth=2)
plt.tight_layout()
plt.show()

Raster file: data.zip

[snowman2]

I don't see the rotation component in your dataset using rasterio:

with rasterio.open("f190808t01p00r07_corr_v111_1293nm.tif") as rds:
    print(rds.transform.to_gdal())

(305827.93, 14.25739965507221, 0.0, 5223236.6, 0.0, -14.25739965507221)

Or with gdal:

ds = gdal.Open("f190808t01p00r07_corr_v111_1293nm.tif")
print(ds.GetGeoTransform())

(305827.93, 14.25739965507221, 0.0, 5223236.6, 0.0, -14.25739965507221)

Version information:

• rasterio 1.2
• GDAL 3.1.4

What are you using to get the GeoTransform with the rotation component?

[adamancer]

The problem is with the file. Here's a new one with the correct transform: corrected.zip.

Both files are a single band extracted from a larger array. I saved the first one using rio.to_raster(), but that seems to have overwritten the original transform. For this one I just used rasterio.

Version info since I didn't include it in the original post:

• GDAL 3.1.4
• rasterio 1.2.0
• rioxarray 0.3.1

[snowman2]

Thanks for the corrected file. I was able to verify that rioxarray does not return the correct transform. Non-rectilinear support is not great at the moment in rioxarray, but rioxarray should definitely return the correct transform at a minimum. Will need to address this bug.

[snowman2]

Forgot to mention, added fix in #285 to return the correct transform.

[adamancer]

That did it. The rotation above now works as expected. Thanks!

[snowman2]

Great 👍