From cb4cfa9fc426e26bea72ae73d9850d6acd06b37b Mon Sep 17 00:00:00 2001
From: "Alan D. Snow" <alansnow21@gmail.com>
Date: Fri, 18 Jun 2021 10:18:53 -0500
Subject: [PATCH] ENH: Add transformer_options to WarpedVRT (#2195)

* ENH: Add transformer_options to WarpedVRT

* REF: Merge transformer_options into warp_extras
---
 rasterio/_warp.pyx      | 49 ++++++++++++++++++++++++++++++++++++++---
 rasterio/gdal.pxi       | 10 +++++++++
 tests/test_warpedvrt.py | 24 +++++++++++++++++++-
 3 files changed, 79 insertions(+), 4 deletions(-)

diff --git a/rasterio/_warp.pyx b/rasterio/_warp.pyx
index d77c27dd..1ec3b8ba 100644
--- a/rasterio/_warp.pyx
+++ b/rasterio/_warp.pyx
@@ -715,6 +715,38 @@ def _calculate_default_transform(src_crs, dst_crs, width, height,
 DEFAULT_NODATA_FLAG = object()
 
 
+cdef GDALDatasetH auto_create_warped_vrt(
+    GDALDatasetH hSrcDS,
+    const char *pszSrcWKT,
+    const char *pszDstWKT,
+    GDALResampleAlg eResampleAlg,
+    double dfMaxError,
+    const GDALWarpOptions *psOptions,
+    const char **transformer_options,
+) nogil:
+
+    IF (CTE_GDAL_MAJOR_VERSION, CTE_GDAL_MINOR_VERSION) >= (3, 2):
+        return GDALAutoCreateWarpedVRTEx(
+            hSrcDS,
+            pszSrcWKT,
+            pszDstWKT,
+            eResampleAlg,
+            dfMaxError,
+            psOptions,
+            transformer_options,
+        )
+
+    ELSE:
+        return GDALAutoCreateWarpedVRT(
+            hSrcDS,
+            pszSrcWKT,
+            pszDstWKT,
+            eResampleAlg,
+            dfMaxError,
+            psOptions,
+        )
+
+
 cdef class WarpedVRTReaderBase(DatasetReaderBase):
 
     def __init__(self, src_dataset, src_crs=None, crs=None,
@@ -773,9 +805,12 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
             means 64 MB with GDAL 2.2.
         dtype : str, optional
             The working data type for warp operation and output.
-        warp_extras : dict
-            GDAL extra warp options. See
+        warp_extras : dict, optional
+            GDAL extra warp options. See:
             https://gdal.org/doxygen/structGDALWarpOptions.html.
+            Also, GDALCreateGenImgProjTransformer2() options.
+            Requires rasterio 1.3+, GDAL 3.2+. See:
+            https://gdal.org/doxygen/gdal__alg_8h.html#a94cd172f78dbc41d6f407d662914f2e3
 
         Returns
         -------
@@ -950,7 +985,15 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
         if not self.dst_transform and (not self.dst_width or not self.dst_height):
 
             with nogil:
-                hds_warped = GDALAutoCreateWarpedVRT(hds, src_crs_wkt, dst_crs_wkt, c_resampling, c_tolerance, psWOptions)
+                hds_warped = auto_create_warped_vrt(
+                    hds,
+                    src_crs_wkt,
+                    dst_crs_wkt,
+                    c_resampling,
+                    c_tolerance,
+                    psWOptions,
+                    <const char **>c_warp_extras,
+                )
 
         else:
 
diff --git a/rasterio/gdal.pxi b/rasterio/gdal.pxi
index 0170226d..958ddf0c 100644
--- a/rasterio/gdal.pxi
+++ b/rasterio/gdal.pxi
@@ -564,6 +564,16 @@ cdef extern from "gdalwarper.h" nogil:
          double *padfGeoTransform, const GDALWarpOptions *psOptionsIn)
 
 
+IF (CTE_GDAL_MAJOR_VERSION, CTE_GDAL_MINOR_VERSION) >= (3, 2):
+    cdef extern from "gdalwarper.h" nogil:
+
+        GDALDatasetH GDALAutoCreateWarpedVRTEx(
+            GDALDatasetH hSrcDS, const char *pszSrcWKT, const char *pszDstWKT,
+            GDALResampleAlg eResampleAlg, double dfMaxError,
+            const GDALWarpOptions *psOptions, char** papszTransformerOptions)
+
+
+
 cdef extern from "gdal_alg.h" nogil:
     void *GDALCreateRPCTransformer( GDALRPCInfo *psRPC, int bReversed,
                           double dfPixErrThreshold,
diff --git a/tests/test_warpedvrt.py b/tests/test_warpedvrt.py
index 77f431ef..24e17f3a 100644
--- a/tests/test_warpedvrt.py
+++ b/tests/test_warpedvrt.py
@@ -19,7 +19,7 @@ from rasterio.vrt import WarpedVRT
 from rasterio.warp import transform_bounds
 from rasterio.windows import Window
 
-from .conftest import requires_gdal21, requires_gdal2
+from .conftest import gdal_version, requires_gdal21, requires_gdal2
 
 log = logging.getLogger(__name__)
 
@@ -205,6 +205,28 @@ def test_warp_extras(path_rgb_byte_tif):
             assert (rgb[:, 0, 0] == 255).all()
 
 
+def test_transformer_options(path_rgb_byte_tif):
+    transform = affine.Affine(
+        1.0003577499128138, 0.0, -8848646.496183893,
+        0.0, -1.0003577499128138, 720.9022505759253,
+    )
+    transformer_options = {
+        "SRC_METHOD": "NO_GEOTRANSFORM",
+        "DST_METHOD": "NO_GEOTRANSFORM",
+    }
+    with rasterio.open(path_rgb_byte_tif) as src, WarpedVRT(
+        src,
+        crs=DST_CRS,
+        **transformer_options,
+    ) as vrt:
+        for key, value in transformer_options.items():
+            assert vrt.warp_extras[key] == value
+        if gdal_version.at_least("3.2"):
+            assert vrt.transform.almost_equals(transform)
+        else:
+            assert not vrt.transform.almost_equals(transform)
+
+
 @requires_gdal21(reason="S3 raster access requires GDAL 2.1+")
 @credentials
 @pytest.mark.network