gdal_translate: transform entire bounds of -projwin, not just two corners (#12127)

References #12100

Author: dbaston

apps/gdal_translate_lib.cpp

@@ -956,10 +956,11 @@ GDALDatasetH GDALTranslate(const char *pszDest, GDALDatasetH hSrcDataset,
                     OGRCoordinateTransformation *poCT =
                         OGRCreateCoordinateTransformation(&oSRSIn, &oSRSDS);
                     if (!(poCT &&
-                          poCT->Transform(1, &psOptions->dfULX,
-                                          &psOptions->dfULY) &&
-                          poCT->Transform(1, &psOptions->dfLRX,
-                                          &psOptions->dfLRY)))
+                          poCT->TransformBounds(
+                              psOptions->dfULX, psOptions->dfLRY,
+                              psOptions->dfLRX, psOptions->dfULY,
+                              &psOptions->dfULX, &psOptions->dfLRY,
+                              &psOptions->dfLRX, &psOptions->dfULY, 21)))
                     {
                         OGRCoordinateTransformation::DestroyCT(poCT);
 

autotest/gcore/vrt_read.py

@@ -1626,9 +1626,9 @@ def test_vrt_protocol():
         "vrt://data/float32.tif?projwin_srs=OGC:CRS84&projwin=-117.6407,33.90027,-117.6292,33.89181"
     )
 
-    assert ds.GetGeoTransform()[0] == 440840.0
+    assert ds.GetGeoTransform()[0] == 440780.0
     assert ds.GetGeoTransform()[3] == 3751140.0
-    assert ds.GetRasterBand(1).XSize == 18
+    assert ds.GetRasterBand(1).XSize == 19
     assert ds.GetRasterBand(1).YSize == 17
 
     with pytest.raises(Exception):

autotest/utilities/test_gdal_translate.py

@@ -847,7 +847,7 @@ def test_gdal_translate_31(gdal_translate_path, tmp_path):
     dst_tif = str(tmp_path / "test_gdal_translate_30.tif")
 
     gdaltest.runexternal(
-        f"{gdal_translate_path} -projwin_srs EPSG:4267 -projwin -117.641168620797 33.9023526904262 -117.628110837847 33.8915970129613 ../gcore/data/byte.tif {dst_tif}"
+        f"{gdal_translate_path} -projwin_srs EPSG:4267 -projwin -117.6408 33.9023 -117.6282 33.8920 ../gcore/data/byte.tif {dst_tif}"
     )
 
     ds = gdal.Open(dst_tif)
@@ -861,8 +861,6 @@ def test_gdal_translate_31(gdal_translate_path, tmp_path):
         1e-6,
     )
 
-    ds = None
-
 
 ###############################################################################
 # Test subsetting a file with a RPC

autotest/utilities/test_gdal_translate_lib.py

@@ -584,6 +584,49 @@ def test_gdal_translate_lib_projwin_expand():
     assert ds.RasterYSize == 20
 
 
+###############################################################################
+# Test -projwin_srs option
+
+
+def test_gdal_translate_lib_31():
+
+    ds = gdal.Translate(
+        "",
+        "../gcore/data/byte.tif",
+        projWin=(-117.6408, 33.9023, -117.6282, 33.8920),
+        projWinSRS="EPSG:4267",
+        format="MEM",
+    )
+
+    assert ds.GetRasterBand(1).Checksum() == 4672, "Bad checksum"
+
+    gdaltest.check_geotransform(
+        gdal.Open("../gcore/data/byte.tif").GetGeoTransform(),
+        ds.GetGeoTransform(),
+        1e-6,
+    )
+
+
+def test_gdal_translate_lib_projwin_polar(tmp_vsimem):
+
+    src = gdal.GetDriverByName("GTiff").Create(
+        tmp_vsimem / "in.tif", 16620, 30000, options={"SPARSE_OK": True}
+    )
+    src.SetGeoTransform((-640000.0, 90.0, 0.0, -655550.0, 0.0, -90.0))
+    src.SetProjection("EPSG:3413")
+
+    dst = gdal.Translate(
+        "", src, projWin=(-20, 80, -19, 79), projWinSRS="EPSG:4326", format="VRT"
+    )
+
+    dst_gt = dst.GetGeoTransform()
+    assert dst_gt[0] == 458900
+    assert dst_gt[3] == -975950
+
+    assert dst.RasterXSize == 723
+    assert dst.RasterYSize == 1192
+
+
 ###############################################################################
 # Test translate with a MEM source to a anonymous VRT