Repository: rasterio/rasterio (language: python)

Issue:
`reproject` with `src_geoloc_array` introduces unexpected shift for group of pixels
<!--

WELCOME ABOARD!

Hi and welcome to the Rasterio project. We appreciate bug reports, questions
about documentation, and suggestions for new features. This issue template
isn't intended to ward you off; only to intercept and redirect some particular
categories of reports, and to collect a few important facts that issue reporters
often omit.

The primary forum for questions about installation and usage of Rasterio is 
https://rasterio.groups.io/g/main. The authors and other users will answer 
questions when they have expertise to share and time to explain. Please take the
time to craft a clear question and be patient about responses. Please do not
bring these questions to Rasterio's issue tracker, which we want to reserve for
bug reports and other actionable issues.

Questions about development of Rasterio, brainstorming, requests for comment,
and not-yet-actionable proposals are welcome in the project's developers 
discussion group https://rasterio.groups.io/g/dev. Issues opened in Rasterio's
GitHub repo which haven't been socialized there may be perfunctorily closed.

Please note: Rasterio contains extension modules and is thus susceptible to
C library compatibility issues. If you are reporting an installation or module
import issue, please note that this project only accepts reports about problems
with packages downloaded from the Python Package Index. Conda users should take
issues to one of the following trackers:

- https://github.com/ContinuumIO/anaconda-issues/issues
- https://github.com/conda-forge/rasterio-feedstock

You think you've found a bug? We believe you!
-->

## Expected behavior and actual behavior.

When using the new `src_geoloc_array` functionality in 1.4a2, I occasionally see a noticeably shifted block of pixels within the image. As my coordinate array is smooth (I think!), I do not expect to see any stepping effects in the output image.

My test pattern looks like this:

![input](https://github.com/rasterio/rasterio/assets/28623784/07662a4f-e00c-4275-acca-cdc06d682eb3)

And the output like this, with a distinct step around y=512:

![output](https://github.com/rasterio/rasterio/assets/28623784/b6b54e2d-7178-48d5-9fd1-42d7f327bb1f)

The `src_geoloc_array` is as follows, with no apparent step:
![x_coords](https://github.com/rasterio/rasterio/assets/28623784/57228418-3a7d-4122-8dde-23dc8cd708a0)
![y_coords](https://github.com/rasterio/rasterio/assets/28623784/cb185e76-fef4-4731-98ec-0286f15b5718)
And e.g. a slice through the `y` coordinate array at `x=256` showing a smooth curve:
![image](https://github.com/rasterio/rasterio/assets/28623784/279bf27b-24d8-4d98-821c-c703059788d2)

I have wondered if this was an artifact related to shoving in a dummy CRS to run the warp, but given that it appears in the matplotlib plots of the underlying data arrays, and doesn't seem to be affected by different choices of CRS (4326 vs 32639), I'm tempted to discount that.

Another possibility is that I've somehow introduced a data type that causes a rounding error, although I can't quite explain why that would cause the observed pattern.

## Steps to reproduce the problem.

Script below. It generates an array containing a grid pattern, and a dummy transformation grid with some radial distortion. The data is warped and then displayed.

```
import numpy as np
import rasterio as rio
from rasterio.warp import calculate_default_transform, reproject
from rasterio.enums import Resampling
import matplotlib.pyplot as plt


WIDTH = 1280
HEIGHT = 1024
X_CENTER = WIDTH // 2
Y_CENTER = HEIGHT // 2
K1 = -3.8e-8
CRS = 4326


def generate_warp_grid(xys):
    out_x = np.zeros(len(xys))
    out_y = np.zeros(len(xys))
    for i, (x, y) in enumerate(xys):
        x_bar = x - X_CENTER
        y_bar = y - Y_CENTER
        r = np.sqrt(x_bar**2 + y_bar**2)

        x_prime = x + x_bar * (K1 * r**2)
        y_prime = y + y_bar * (K1 * r**2)

        out_x[i] = x_prime
        out_y[i] = y_prime

    return np.array([out_x, out_y]).T


def main():
    data = np.zeros((HEIGHT, WIDTH), dtype=np.uint16)
    for i in np.arange(0, HEIGHT, 64, dtype=np.int32):
        data[i:i+4, :] = 1
    for j in np.arange(0, WIDTH, 64, dtype=np.int32):
        data[:, j:j+4] = 1

    xx, yy = np.meshgrid(
        np.arange(0, WIDTH),
        np.arange(0, HEIGHT),
    )

    stacked = np.vstack([xx.reshape((-1)), yy.reshape((-1))]).T
    coordinate_grid = generate_warp_grid(stacked).T.reshape((2, HEIGHT, WIDTH))

    crs = rio.CRS.from_epsg(CRS)
    out_transform, out_width, out_height = calculate_default_transform(
        crs,
        crs,
        WIDTH,
        HEIGHT,
        src_geoloc_array=coordinate_grid,
    )

    out_arr = np.zeros((out_height, out_width), dtype=data.dtype)

    reproject(
        data,
        destination=out_arr,
        src_geoloc_array=coordinate_grid,
        src_crs=crs,
        dst_crs=crs,
        dst_transform=out_transform,
        resampling=Resampling.cubic,
    )

    plt.imsave("x_coords.png", coordinate_grid[0])
    plt.imsave("y_coords.png", coordinate_grid[1])
    plt.imsave("input.png", data)
    plt.imsave("output.png", out_arr)

if __name__ == "__main__":
    main()
```

#### Environment Information
```
rasterio info:
  rasterio: 1.4a2
      GDAL: 3.8.4
      PROJ: 9.3.1
      GEOS: 3.11.1
 PROJ DATA: /home/danielevans/.cache/pypoetry/virtualenvs/image-processing-LWpwzOMN-py3.10/lib/python3.10/site-packages/rasterio/proj_data
 GDAL DATA: /home/danielevans/.cache/pypoetry/virtualenvs/image-processing-LWpwzOMN-py3.10/lib/python3.10/site-packages/rasterio/gdal_data

System:
    python: 3.10.12 (main, Nov 20 2023, 15:14:05) [GCC 11.4.0]
executable: /home/danielevans/.cache/pypoetry/virtualenvs/image-processing-LWpwzOMN-py3.10/bin/python
   machine: Linux-5.15.146.1-microsoft-standard-WSL2-x86_64-with-glibc2.35

Python deps:
    affine: 2.4.0
     attrs: 23.1.0
   certifi: 2023.07.22
     click: 8.1.7
     cligj: 0.7.2
    cython: None
     numpy: 1.26.4
click-plugins: None
setuptools: 68.2.2
```

## Installation Method

1.4a2 wheel from PyPI, using Poetry/Pip



Discussion / hints:
Perhaps useful for context, this is the effect on a real-world raster, where the roads become significantly misaligned:
![image](https://github.com/rasterio/rasterio/assets/28623784/6e95d4a3-9b79-4371-b82a-1d13f489014a)

This bug seems to be caused by rasterio choosing to use an Approximate Transformer - the `bUseApproxTransformer` flag is always `True` here unless using RPCs:

https://github.com/rasterio/rasterio/blob/842743055124b1b4ef0f43907467a3a95b8fdb82/rasterio/_warp.pyx#L503-L512

With the current main branch, with `cdef bint bUseApproxTransformer = True`:
![output_main_approx_transformer](https://github.com/rasterio/rasterio/assets/28623784/e9ff3472-edae-473e-8ead-06543cbe4733)

And changing it to `False`
![output_main_full_transformer](https://github.com/rasterio/rasterio/assets/28623784/67bd507e-0541-43a0-9071-3389814e9d3a)

My naive suggestion for a fix would be adding `bUseApproxTransformer = False` inside the `elif src_geoloc_array is not None` branch in the _reproject function, although I'm not fully sure of what other side effects using a non-approximate transformer might give (slightly worse performance?):
https://github.com/rasterio/rasterio/blob/842743055124b1b4ef0f43907467a3a95b8fdb82/rasterio/_warp.pyx#L470-L481

Fix patch:
diff --git a/rasterio/_warp.pyx b/rasterio/_warp.pyx
index 84ff2907d..fc627bcb8 100644
--- a/rasterio/_warp.pyx
+++ b/rasterio/_warp.pyx
@@ -479,6 +479,7 @@ def _reproject(
         imgProjOptions = CSLSetNameValue(
             imgProjOptions, "SRC_SRS", src_crs.to_string().encode("utf-8")
         )
+        bUseApproxTransformer = False
 
     # See https://gdal.org/doxygen/gdal__alg_8h.html#a94cd172f78dbc41d6f407d662914f2e3
     # for a list of supported options. I (Sean) don't see harm in
@@ -727,6 +728,7 @@ def _calculate_default_transform(
             imgProjOptions = CSLSetNameValue(
                 imgProjOptions, "SRC_SRS", src_crs.to_string().encode("utf-8")
             )
+            bUseApproxTransformer = False
 
         hTransformArg = exc_wrap_pointer(
             GDALCreateGenImgProjTransformer2(hds, NULL, imgProjOptions)
@@ -1319,14 +1321,14 @@ def _suggested_proxy_vrt_doc(width, height, transform=None, crs=None, gcps=None,
             "line_off", "samp_off", "height_off", "lat_off", "long_off",
             "line_scale", "samp_scale", "height_scale", "lat_scale", "long_scale",
             "err_bias", "err_rand",
-        ): 
+        ):
             item = ET.SubElement(metadata, "MDI")
             item.attrib["key"] = key.upper()
             item.text = str(getattr(rpcs, key))
         for key in ("line_num_coeff", "line_den_coeff", "samp_num_coeff", "samp_den_coeff"):
             item = ET.SubElement(metadata, "MDI")
             item.attrib["key"] = key.upper()
-            item.text = " ".join(["{:+e}".format(val) for val in getattr(rpcs, key)]) 
+            item.text = " ".join(["{:+e}".format(val) for val in getattr(rpcs, key)])
 
     return ET.tostring(vrtdataset)
 


Test patch:
diff --git a/tests/test_warp.py b/tests/test_warp.py
index 49659c92e..b6b7199e3 100644
--- a/tests/test_warp.py
+++ b/tests/test_warp.py
@@ -2317,5 +2317,5 @@ def test_geoloc_warp_array_subsampled(path_rgb_byte_tif, tmp_path):
     )
 
     # This value is specific to DST_TRANSFORM and an 800 x 880 array.
-    assert np.count_nonzero(output[0]) in [471492]
+    assert np.count_nonzero(output[0]) in [471533]