Merge pull request #11248 from rouault/fix_rasterio_3070

VRT: fix reading from a CFloat32/CFloat64 ComplexSource in a non-complex-type buffer
OSGeo · Nov 13, 2024 · f2c7ec4 · f2c7ec4
2 parents c17c048 + 934312e
commit f2c7ec4
Show file tree

Hide file tree

Showing 5 changed files with 51 additions and 2 deletions.
diff --git a/autotest/gcore/data/gtiff/complex_non_zero_real_zero_imag.tif b/autotest/gcore/data/gtiff/complex_non_zero_real_zero_imag.tif
diff --git a/autotest/gcore/data/vrt/complex_non_zero_real_zero_imag.vrt b/autotest/gcore/data/vrt/complex_non_zero_real_zero_imag.vrt
@@ -0,0 +1,12 @@
+<VRTDataset rasterXSize="2" rasterYSize="2">
+  <SRS dataAxisToSRSAxisMapping="1,2">PROJCS["WGS 84 / UTM zone 11N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32611"]]</SRS>
+  <GeoTransform>  2.1520000000000000e+05,  3.0000000000000000e+01,  0.0000000000000000e+00,  4.3975000000000000e+06,  0.0000000000000000e+00, -3.0000000000000000e+01</GeoTransform>
+  <VRTRasterBand dataType="CFloat32" band="1">
+    <ColorInterp>Gray</ColorInterp>
+    <ComplexSource>
+      <NODATA>0</NODATA>
+      <SourceFilename relativeToVRT="1">../gtiff/complex_non_zero_real_zero_imag.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+    </ComplexSource>
+  </VRTRasterBand>
+</VRTDataset>
diff --git a/autotest/gcore/data/vrt/complex_non_zero_real_zero_imag_as_float32.vrt b/autotest/gcore/data/vrt/complex_non_zero_real_zero_imag_as_float32.vrt
@@ -0,0 +1,12 @@
+<VRTDataset rasterXSize="2" rasterYSize="2">
+  <SRS dataAxisToSRSAxisMapping="1,2">PROJCS["WGS 84 / UTM zone 11N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32611"]]</SRS>
+  <GeoTransform>  2.1520000000000000e+05,  3.0000000000000000e+01,  0.0000000000000000e+00,  4.3975000000000000e+06,  0.0000000000000000e+00, -3.0000000000000000e+01</GeoTransform>
+  <VRTRasterBand dataType="Float32" band="1">
+    <ColorInterp>Gray</ColorInterp>
+    <ComplexSource>
+      <NODATA>0</NODATA>
+      <SourceFilename relativeToVRT="1">../gtiff/complex_non_zero_real_zero_imag.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+    </ComplexSource>
+  </VRTRasterBand>
+</VRTDataset>
diff --git a/autotest/gcore/vrt_read.py b/autotest/gcore/vrt_read.py
@@ -2837,3 +2837,29 @@ def test_vrt_update_nested_VRTDataset(tmp_vsimem):
     with gdal.Open(vrt_filename) as ds:
         assert ds.GetRasterBand(1).Checksum() == 4672
         assert ds.GetRasterBand(1).GetMinimum() == 74
+
+
+###############################################################################
+# Test reading a VRT with a ComplexSource of type CFloat32
+
+
+def test_vrt_read_cfloat32_complex_source_as_float32():
+
+    ds = gdal.Open("data/vrt/complex_non_zero_real_zero_imag.vrt")
+    assert struct.unpack("f" * 8, ds.ReadRaster()) == (1, 0, 1, 0, 1, 0, 1, 0)
+    assert struct.unpack("f" * 4, ds.ReadRaster(buf_type=gdal.GDT_Float32)) == (
+        1,
+        1,
+        1,
+        1,
+    )
+
+
+###############################################################################
+# Test reading a VRT with a ComplexSource of type Float32 (from a underlying CFloat32 source)
+
+
+def test_vrt_read_float32_complex_source_from_cfloat32():
+
+    ds = gdal.Open("data/vrt/complex_non_zero_real_zero_imag_as_float32.vrt")
+    assert struct.unpack("f" * 4, ds.ReadRaster()) == (1, 1, 1, 1)
diff --git a/frmts/vrt/vrtsources.cpp b/frmts/vrt/vrtsources.cpp
@@ -3068,8 +3068,7 @@ CPLErr VRTComplexSource::RasterIO(GDALDataType eVRTBandDataType, int nXOff,
         }
     }
 
-    const bool bIsComplex =
-        CPL_TO_BOOL(GDALDataTypeIsComplex(eVRTBandDataType));
+    const bool bIsComplex = CPL_TO_BOOL(GDALDataTypeIsComplex(eBufType));
     CPLErr eErr;
     // For Int32, float32 isn't sufficiently precise as working data type
     if (eVRTBandDataType == GDT_CInt32 || eVRTBandDataType == GDT_CFloat64 ||