Use Lapack solve to more accurately compute affine reverse transforms.#3315
Use Lapack solve to more accurately compute affine reverse transforms.#3315snowman2 merged 14 commits intorasterio:mainfrom
Conversation
On certain transforms, directly inverting the transform matrix leads to subtle numerical inaccuracies. Using np.linalg.solve avoids these issues.
This should work more nicely with affine 3.0 and its new array interface.
|
It looks like there are two other locations where we invert the transform. Should we review that code as well in light of the findings in #3313? |
I think this reimplementation is much easier to understand.
rasterio/rio/warp.py
Outdated
| row2 = ceil(row2) | ||
| rows = np.array([top, top, bottom, bottom]) + eps | ||
| cols = np.array([left, right, right, left]) + eps | ||
| rows, cols = rowcol(src.transform, rows, cols) |
There was a problem hiding this comment.
The default for rowcol is to floor results. The tests pass.
Oddly, when trying to exactly match the old behavior (op=float and doing the floor/ceil calls), the test_warp_no_reproject test fails.
The two transforms do not match.
Affine(9.545306945486766, 0.0, -11858134.818413004,
0.0, -9.545306945486766, 4813698.2926443005) = <open DatasetReader name='test_warp_no_reproject0/test.tif' mode='r'>.transform
Affine(9.5546285343, 0.0, -11858134.818413004,
0.0, -9.5546285343, 4813698.2926443005) = <open DatasetReader name='tests/data/shade.tif' mode='r'>.transform
|
@sgillies I believe this should be ready for a preliminary review.
|
|
The test failures look like they might be related to #3301 |
|
With a moderate risk of an under informed drive-by comment: have you considered using |
|
@geowurster, unfortunately, gdal uses the same numerically unstable method to compute the inverse geotransform as affine. |
|
@groutr ah, thanks. |
|
@groutr mind updating this PR to the latest? |
Updated with the latest changes. |
Co-authored-by: Alan D. Snow <[email protected]>
|
Thanks @groutr 👍 |
3 participants
With certain transforms, directly inverting the transform matrix leads to subtle numerical inaccuracies. Using np.linalg.solve avoids these issues.
Replaces #3313
Fixes #3310