UndefinedRotationError in RotatedCoordinates geotransform roundtrip

[jmilloy]

Description We should be able to recover RotatedCoordinates from their geotransform. But in some cases, the from_geotransform raises an UndefinedRotationError.

Steps to Reproduce

c = RotatedCoordinates(shape=(2, 4), theta=np.pi / 8, origin=[10, 20], step=[-2.0, 1.0], dims=["lat", "lon"])
c2 = RotatedCoordinates.from_geotransform(c.geotransform, c.shape, dims=['lat', 'lon'])
assert c == c2

Expected Behavior The RotatedCoordinates conveltion to geotransform and the back should work and be the same.

Observed Behavior

>>> c = RotatedCoordinates(shape=(2, 4), theta=np.pi / 8, origin=[10, 20], step=[-2.0, 1.0], dims=["lat", "lon"])
>>> c.geotransform
(19.5, 0.9238795325112867, 0.7653668647301796, 11.0, 0.3826834323650898, -1.8477590650225735)
>>> c2 = RotatedCoordinates.from_geotransform(c.geotransform, c.shape, dims=['lat', 'lon'])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/jmilloy/Creare/Pipeline/podpac/podpac/core/coordinates/rotated_coordinates.py", line 121, in from_geotransform
    deg = affine.rotation_angle
  File "/home/jmilloy/Creare/Pipeline/_podpac-38_/lib/python3.8/site-packages/affine/__init__.py", line 385, in rotation_angle
    raise UndefinedRotationError
affine.UndefinedRotationError

[jmilloy]

@mpu-creare This is the RotatedCoordinates used in test_units.py::TestToGeotiff::test_to_geotiff_roundtrip_rotcoords and issue #363. I've fixed a few things for that test and isuse, but this is another reasons that test fails. I'm working on this issue now, but I wonder if you might be able to take a look, too. Any immediate ideas?

[jmilloy]

It looks like the negative step step=[-2.0, 1.0] is part of the problem.

[jmilloy]

Confirmed that using a different parameterization for the rotated coordinates (with a positive step and a different "origin") does work:

>>> c_prime = RotatedCoordinates(shape=(2, 4), theta=np.pi / 8, origin=c[-1, 0].coordinates, step=[2.0, 1.0], dims=["lat", "lon"])
>>> RotatedCoordinates.from_geotransform(c_prime.geotransform, c_prime.shape, dims=['lat', 'lon'])
RotatedCoordinates(('lat', 'lon')): Origin[ 8.15224093 19.23463314], Corner[ 8.8519497 22.7716386], rad[0.3927], shape(2, 4)

So @mpu-creare I can fuss with the geotransform to try to fix the negative step, or we could sort of take a shortcut and only allow RotatedCoordinates with positive steps. What do you think?

[jmilloy]

I guess this could potentially be considered a bug in rasterio? So maybe we can warn about that during RotatedCoordinates initialization and just avoid using negative steps in our unit tests.

[mpu-creare]

Wow you're fast. I was just starting to investigate. Can we do the old try-catch-except negative sign deal? And then just log the error? The negative step size is probably pretty common for i+ == n-->s geotiffs (with relatively small rotation), so being able to handle the case would be nice... and probably important for some datasets.

Is it work asking on the Rasterio thread? Perhaps there's a spec or something we're not aware of.

[jmilloy]

Can we do the old try-catch-except negative sign deal? And then just log the error?

I'm not sure what you mean. But I can say that we already use a try-except to check for UndefinedRotation exceptions in Coordinates.from_geotransform in order to determine if a geotransform is rotated or not. So I don't know how to tell if an UndefinedRotation exception is one of these spurious ones from rotated coordinates with a negative step or if it's actually Uniform. I guess we could explicitly check (mathematically) if a geotransform is rotated or not rather than rely on the exception, which I would generally prefer anyways. Are you thinking there is a way to manually get the rotation and step once we know we have one of these rotated negative step geotransforms?

Is it work asking on the Rasterio thread? Perhaps there's a spec or something we're not aware of.

Yes, I think so. This is something I thought you might be a lot better at than me

[jmilloy]

Let me spend a few minutes putting together a code snippet that only uses rasterio.affine

[jmilloy]

Comes down to:

>>> r = rasterio.Affine.rotation(5)
>>> s = rasterio.Affine.scale(1, -2)
>>> a = r * s
>>> a.rotation_angle
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/jmilloy/Creare/Pipeline/_podpac-38_/lib/python3.8/site-packages/affine/__init__.py", line 385, in rotation_angle
    raise UndefinedRotationError
affine.UndefinedRotationError

[jmilloy]

Here's the code for Affine.scale

https://github.com/rasterio/affine/blob/78c20a0cfbb5ece3dfadc1a5550465936a4fa731/affine/__init__.py#L213

Nothing there in the docstring indicates only positive values are supported, but you can see the simplicity of what's going on. Obviously a negative scale isn't actually a "scale" transformation. I think in our code, we need to make a transformation that doesn't have a negative scale operation in it.

[jmilloy]

That's all I've got for now. Dropping support for negative steps just means you have to create your rotated coordinates with the "origin" at the correct corner. I don't think any geotiffs or other rasterio data sources could ever have a rotated geotransform that uses a different corner. If you are making an array data source, you have to reverse your data.

Related: I'm surprised that the UndefinedRotationError is the necessary check for uniform coordinates. I bet I can make a rotated coordinates with a 0 degree rotation that are essentially uniform grid coordinates, but they'll have a 0 degree rotation_angle. I'm going check tomorrow how the GDAL geotransforms compare for uniform coordinates and 0-degree rotated coordinates that parameterize the same coordinates grid. Curious about that.

[jmilloy]

Closing. This might come up again, but let's see how the AffineCoordinates behave in practice.