From c180f4be6363f0b127689afc8cd41589a7a7779f Mon Sep 17 00:00:00 2001 From: Daniel Baston Date: Thu, 9 Jan 2025 13:28:19 -0500 Subject: [PATCH 1/4] VRTProcessedDataset: unscale source raster values --- autotest/gdrivers/vrtprocesseddataset.py | 112 +++++++++++++- .../drivers/raster/vrt_processed_dataset.rst | 2 +- frmts/vrt/vrtdataset.h | 3 + frmts/vrt/vrtprocesseddataset.cpp | 137 +++++++++++++++--- 4 files changed, 233 insertions(+), 21 deletions(-) diff --git a/autotest/gdrivers/vrtprocesseddataset.py b/autotest/gdrivers/vrtprocesseddataset.py index 633f2158a654..0d1f15d03837 100755 --- a/autotest/gdrivers/vrtprocesseddataset.py +++ b/autotest/gdrivers/vrtprocesseddataset.py @@ -21,7 +21,7 @@ ) np = pytest.importorskip("numpy") -pytest.importorskip("osgeo.gdal_array") +gdal_array = pytest.importorskip("osgeo.gdal_array") ############################################################################### # Test error cases in general VRTProcessedDataset XML structure @@ -73,6 +73,16 @@ def test_vrtprocesseddataset_errors(tmp_vsimem): src_ds.GetRasterBand(3).Fill(3) src_ds.Close() + with pytest.raises(Exception, match="Invalid value of 'unscale'"): + gdal.Open( + f""" + + {src_filename} + + + """ + ) + with pytest.raises(Exception, match="ProcessingSteps element missing"): gdal.Open( f""" @@ -1216,7 +1226,7 @@ def test_vrtprocesseddataset_serialize(tmp_vsimem): vrt_filename = str(tmp_vsimem / "the.vrt") content = f""" - + {src_filename} @@ -1517,3 +1527,101 @@ def test_vrtprocesseddataset_RasterIO(tmp_vsimem): assert ds.GetRasterBand(1).GetBlockSize() == [1, 1] with pytest.raises(Exception): ds.ReadAsArray() + + +############################################################################### +# Test reading input datasets with scale and offset + + +@pytest.mark.parametrize("input_scaled", (True, False)) +@pytest.mark.parametrize("unscale", (True, False, "auto"), ids=lambda x: f"unscale={x}") +@pytest.mark.parametrize( + "dtype", (gdal.GDT_Int16, gdal.GDT_Float32), ids=gdal.GetDataTypeName +) +def test_vrtprocesseddataset_scaled_inputs(tmp_vsimem, input_scaled, dtype, unscale): + + src_filename = tmp_vsimem / "src.tif" + + nx = 2 + ny = 3 + nz = 2 + + if dtype == gdal.GDT_Float32: + nodata = float("nan") + else: + nodata = 99 + + np_type = gdal_array.GDALTypeCodeToNumericTypeCode(dtype) + + data = np.arange(nx * ny * nz, dtype=np_type).reshape(nz, ny, nx) + data[:, 2, 1] = nodata + + if input_scaled: + offsets = [i + 2 for i in range(nz)] + scales = [(i + 1) / 4 for i in range(nz)] + else: + offsets = [0 for i in range(nz)] + scales = [1 for i in range(nz)] + + with gdal.GetDriverByName("GTiff").Create( + src_filename, nx, ny, nz, eType=dtype + ) as src_ds: + src_ds.WriteArray(data) + for i in range(src_ds.RasterCount): + bnd = src_ds.GetRasterBand(i + 1) + bnd.SetOffset(offsets[i]) + bnd.SetScale(scales[i]) + bnd.SetNoDataValue(nodata) + + ds = gdal.Open( + f""" + + + {src_filename} + + + + BandAffineCombination + 0,1,0 + 0,0,1 + + + """ + ) + + assert ds.RasterCount == nz + + if unscale is True or (unscale == "auto" and input_scaled): + for i in range(ds.RasterCount): + bnd = ds.GetRasterBand(i + 1) + assert bnd.DataType == gdal.GDT_Float64 + assert bnd.GetScale() in (None, 1) + assert bnd.GetOffset() in (None, 0) + else: + for i in range(ds.RasterCount): + bnd = ds.GetRasterBand(i + 1) + assert bnd.DataType == dtype + assert bnd.GetScale() == scales[i] + assert bnd.GetOffset() == offsets[i] + assert ( + np.isnan(bnd.GetNoDataValue()) + if np.isnan(nodata) + else bnd.GetNoDataValue() == nodata + ) + + result = np.ma.stack( + [ds.GetRasterBand(i + 1).ReadAsMaskedArray() for i in range(ds.RasterCount)] + ) + + if unscale: + expected = np.ma.masked_array( + np.stack([data[i, :, :] * scales[i] + offsets[i] for i in range(nz)]), + np.isnan(data) if np.isnan(nodata) else data == nodata, + ) + else: + expected = np.ma.masked_array( + data, np.isnan(data) if np.isnan(nodata) else data == nodata + ) + + np.testing.assert_array_equal(result.mask, expected.mask) + np.testing.assert_array_equal(result[~result.mask], expected[~expected.mask]) diff --git a/doc/source/drivers/raster/vrt_processed_dataset.rst b/doc/source/drivers/raster/vrt_processed_dataset.rst index 28f001904af7..b3383126763b 100644 --- a/doc/source/drivers/raster/vrt_processed_dataset.rst +++ b/doc/source/drivers/raster/vrt_processed_dataset.rst @@ -121,7 +121,7 @@ The following child elements of ``VRTDataset`` may be defined: ``SRS``, ``GeoTra The ``VRTDataset`` root element must also have the 2 following child elements: -- ``Input``, which must have one and only one of the following ``SourceFilename`` or ``VRTDataset`` as child elements, to define the input dataset to which to apply the processing steps. +- ``Input``, which must have one and only one of the following ``SourceFilename`` or ``VRTDataset`` as child elements, to define the input dataset to which to apply the processing steps. Starting with GDAL 3.11, values from the input dataset will be automatically unscaled; this can be disabled by setting the ``unscale`` attribute of ``Input`` to ``false``. - ``ProcessingSteps``, with at least one child ``Step`` element. diff --git a/frmts/vrt/vrtdataset.h b/frmts/vrt/vrtdataset.h index 50fbac8e0a46..6354f23878d3 100644 --- a/frmts/vrt/vrtdataset.h +++ b/frmts/vrt/vrtdataset.h @@ -770,6 +770,9 @@ class VRTProcessedDataset final : public VRTDataset //! Value of CPLGetUsablePhysicalRAM() / 10 * 4 GIntBig m_nAllowedRAMUsage = 0; + //! Whether to apply a scale and offset if defined by the input dataset + bool m_bUnscale = true; + CPLErr Init(const CPLXMLNode *, const char *, const VRTProcessedDataset *poParentDS, GDALDataset *poParentSrcDS, int iOvrLevel); diff --git a/frmts/vrt/vrtprocesseddataset.cpp b/frmts/vrt/vrtprocesseddataset.cpp index 318da09023aa..f91971a8bc51 100644 --- a/frmts/vrt/vrtprocesseddataset.cpp +++ b/frmts/vrt/vrtprocesseddataset.cpp @@ -206,6 +206,27 @@ CPLErr VRTProcessedDataset::XMLInit(const CPLXMLNode *psTree, return CE_None; } +static bool HasScaleOffset(GDALDataset &oSrcDS) +{ + for (int i = 1; i <= oSrcDS.GetRasterCount(); i++) + { + int pbSuccess; + GDALRasterBand &oBand = *oSrcDS.GetRasterBand(i); + double scale = oBand.GetScale(&pbSuccess); + if (pbSuccess && scale != 1) + { + return true; + } + double offset = oBand.GetOffset(&pbSuccess); + if (pbSuccess && offset != 0) + { + return true; + } + } + + return false; +} + /** Instantiate object from XML tree */ CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree, const char *pszVRTPathIn, @@ -260,6 +281,27 @@ CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree, if (!m_poSrcDS) return CE_Failure; + const char *pszUnscale = CPLGetXMLValue(psInput, "unscale", "auto"); + if (EQUAL(pszUnscale, "AUTO")) + { + m_bUnscale = HasScaleOffset(*m_poSrcDS); + } + else if (EQUAL(pszUnscale, "NO") || EQUAL(pszUnscale, "OFF") || + EQUAL(pszUnscale, "FALSE") || EQUAL(pszUnscale, "0")) + { + m_bUnscale = false; + } + else if (EQUAL(pszUnscale, "YES") || EQUAL(pszUnscale, "ON") || + EQUAL(pszUnscale, "TRUE") || EQUAL(pszUnscale, "1")) + { + m_bUnscale = true; + } + else + { + CPLError(CE_Failure, CPLE_AppDefined, "Invalid value of 'unscale'"); + return CE_Failure; + } + if (nRasterXSize == 0 && nRasterYSize == 0) { nRasterXSize = m_poSrcDS->GetRasterXSize(); @@ -437,17 +479,23 @@ CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree, return CE_Failure; } - const auto eInDT = poSrcFirstBand->GetRasterDataType(); - for (int i = 1; i < m_poSrcDS->GetRasterCount(); ++i) + const auto eInDT = + m_bUnscale ? GDT_Float64 : poSrcFirstBand->GetRasterDataType(); + if (!m_bUnscale) { - const auto eDT = m_poSrcDS->GetRasterBand(i + 1)->GetRasterDataType(); - if (eDT != eInDT) + for (int i = 1; i < m_poSrcDS->GetRasterCount(); ++i) { - CPLError(CE_Warning, CPLE_AppDefined, - "Not all bands of the input dataset have the same data " - "type. The data type of the first band will be used as " - "the reference one."); - break; + const auto eDT = + m_poSrcDS->GetRasterBand(i + 1)->GetRasterDataType(); + if (eDT != eInDT) + { + CPLError( + CE_Warning, CPLE_AppDefined, + "Not all bands of the input dataset have the same data " + "type. The data type of the first band will be used as " + "the reference one."); + break; + } } } @@ -455,13 +503,21 @@ CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree, int nCurrentBandCount = m_poSrcDS->GetRasterCount(); std::vector adfNoData; - for (int i = 1; i <= nCurrentBandCount; ++i) + if (m_bUnscale) + { + adfNoData.resize(nCurrentBandCount, + std::numeric_limits::quiet_NaN()); + } + else { - int bHasVal = FALSE; - const double dfVal = - m_poSrcDS->GetRasterBand(i)->GetNoDataValue(&bHasVal); - adfNoData.emplace_back( - bHasVal ? dfVal : std::numeric_limits::quiet_NaN()); + for (int i = 1; i <= nCurrentBandCount; ++i) + { + int bHasVal = FALSE; + const double dfVal = + m_poSrcDS->GetRasterBand(i)->GetNoDataValue(&bHasVal); + adfNoData.emplace_back( + bHasVal ? dfVal : std::numeric_limits::quiet_NaN()); + } } int nStepCount = 0; @@ -654,10 +710,18 @@ CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree, { const auto poSrcBand = m_poSrcDS->GetRasterBand(i + 1); const GDALDataType eOutputBandType = - GetOutputBandType(poSrcBand->GetRasterDataType()); + m_bUnscale ? GDT_Float64 + : GetOutputBandType(poSrcBand->GetRasterDataType()); auto poBand = new VRTProcessedRasterBand(this, i + 1, eOutputBandType); poBand->CopyCommonInfoFrom(poSrcBand); + if (m_bUnscale) + { + poBand->SetScale(1); + poBand->SetOffset(0); + poBand->SetNoDataValue( + std::numeric_limits::quiet_NaN()); + } SetBand(i + 1, poBand); } } @@ -1137,7 +1201,8 @@ bool VRTProcessedDataset::ProcessRegion(int nXOff, int nYOff, int nBufXSize, const char *pszInterleave = m_poSrcDS->GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE"); - if (nFirstBandCount > 1 && (!pszInterleave || EQUAL(pszInterleave, "BAND"))) + if (m_bUnscale || (nFirstBandCount > 1 && + (!pszInterleave || EQUAL(pszInterleave, "BAND")))) { // If there are several bands and the source dataset organization // is apparently band interleaved, then first acquire data in @@ -1145,7 +1210,9 @@ bool VRTProcessedDataset::ProcessRegion(int nXOff, int nYOff, int nBufXSize, // data type. // And then transpose it and convert it to the expected data type // of the first step. - const auto eSrcDT = m_poSrcDS->GetRasterBand(1)->GetRasterDataType(); + const auto eSrcDT = + m_bUnscale ? eFirstDT + : m_poSrcDS->GetRasterBand(1)->GetRasterDataType(); try { abyInput.resize(nPixelCount * nFirstBandCount * @@ -1188,6 +1255,40 @@ bool VRTProcessedDataset::ProcessRegion(int nXOff, int nYOff, int nBufXSize, if (!bOK) return false; + if (m_bUnscale) + { + double *pdfValue = reinterpret_cast(abyInput.data()); + for (int nBand = 0; nBand < nFirstBandCount; nBand++) + { + auto poBand = m_poSrcDS->GetRasterBand(nBand + 1); + const double dfScale = poBand->GetScale(nullptr); + const double dfOffset = poBand->GetOffset(nullptr); + int bHasNoData; + const double dfNoData = poBand->GetNoDataValue(&bHasNoData); + const double *pdfEnd = pdfValue + nPixelCount; + + if (bHasNoData) + { + while (pdfValue != pdfEnd) + { + *pdfValue = + (*pdfValue != dfNoData) + ? (*pdfValue * dfScale + dfOffset) + : std::numeric_limits::quiet_NaN(); + pdfValue++; + } + } + else + { + while (pdfValue != pdfEnd) + { + *pdfValue = *pdfValue * dfScale + dfOffset; + pdfValue++; + } + } + } + } + CPLDebugOnly("VRT", "ProcessRegion(): start GDALTranspose2D()"); GDALTranspose2D(abyInput.data(), eSrcDT, abyOutput.data(), eFirstDT, static_cast(nBufXSize) * nBufYSize, From c038111e644b7b4bb3167b42152d52a7740a69c5 Mon Sep 17 00:00:00 2001 From: Daniel Baston Date: Thu, 9 Jan 2025 16:17:09 -0500 Subject: [PATCH 2/4] VRTProcessedDataset: Use GDALTranslate to handle unscaling --- autotest/gdrivers/vrtprocesseddataset.py | 4 +- frmts/vrt/vrtdataset.h | 6 +- frmts/vrt/vrtprocesseddataset.cpp | 138 +++++++++-------------- 3 files changed, 58 insertions(+), 90 deletions(-) diff --git a/autotest/gdrivers/vrtprocesseddataset.py b/autotest/gdrivers/vrtprocesseddataset.py index 0d1f15d03837..0a31d8f3cc8f 100755 --- a/autotest/gdrivers/vrtprocesseddataset.py +++ b/autotest/gdrivers/vrtprocesseddataset.py @@ -1533,7 +1533,9 @@ def test_vrtprocesseddataset_RasterIO(tmp_vsimem): # Test reading input datasets with scale and offset -@pytest.mark.parametrize("input_scaled", (True, False)) +@pytest.mark.parametrize( + "input_scaled", (True, False), ids=lambda x: f"input scaled={x}" +) @pytest.mark.parametrize("unscale", (True, False, "auto"), ids=lambda x: f"unscale={x}") @pytest.mark.parametrize( "dtype", (gdal.GDT_Int16, gdal.GDT_Float32), ids=gdal.GetDataTypeName diff --git a/frmts/vrt/vrtdataset.h b/frmts/vrt/vrtdataset.h index 6354f23878d3..aa6852b6b80c 100644 --- a/frmts/vrt/vrtdataset.h +++ b/frmts/vrt/vrtdataset.h @@ -724,6 +724,9 @@ class VRTProcessedDataset final : public VRTDataset //! Directory of the VRT std::string m_osVRTPath{}; + //! Source of source dataset generated with GDALTranslate + std::unique_ptr m_poVRTSrcDS{}; + //! Source dataset std::unique_ptr m_poSrcDS{}; @@ -770,9 +773,6 @@ class VRTProcessedDataset final : public VRTDataset //! Value of CPLGetUsablePhysicalRAM() / 10 * 4 GIntBig m_nAllowedRAMUsage = 0; - //! Whether to apply a scale and offset if defined by the input dataset - bool m_bUnscale = true; - CPLErr Init(const CPLXMLNode *, const char *, const VRTProcessedDataset *poParentDS, GDALDataset *poParentSrcDS, int iOvrLevel); diff --git a/frmts/vrt/vrtprocesseddataset.cpp b/frmts/vrt/vrtprocesseddataset.cpp index f91971a8bc51..5591039fe430 100644 --- a/frmts/vrt/vrtprocesseddataset.cpp +++ b/frmts/vrt/vrtprocesseddataset.cpp @@ -12,6 +12,7 @@ #include "cpl_minixml.h" #include "cpl_string.h" +#include "gdal_utils.h" #include "vrtdataset.h" #include @@ -281,27 +282,51 @@ CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree, if (!m_poSrcDS) return CE_Failure; - const char *pszUnscale = CPLGetXMLValue(psInput, "unscale", "auto"); + const char *pszUnscale = CPLGetXMLValue(psInput, "unscale", "AUTO"); + bool bUnscale = false; if (EQUAL(pszUnscale, "AUTO")) { - m_bUnscale = HasScaleOffset(*m_poSrcDS); - } - else if (EQUAL(pszUnscale, "NO") || EQUAL(pszUnscale, "OFF") || - EQUAL(pszUnscale, "FALSE") || EQUAL(pszUnscale, "0")) - { - m_bUnscale = false; + if (HasScaleOffset(*m_poSrcDS)) + { + bUnscale = true; + } } else if (EQUAL(pszUnscale, "YES") || EQUAL(pszUnscale, "ON") || EQUAL(pszUnscale, "TRUE") || EQUAL(pszUnscale, "1")) { - m_bUnscale = true; + bUnscale = true; } - else + else if (!(EQUAL(pszUnscale, "NO") || EQUAL(pszUnscale, "OFF") || + EQUAL(pszUnscale, "FALSE") || EQUAL(pszUnscale, "0"))) { CPLError(CE_Failure, CPLE_AppDefined, "Invalid value of 'unscale'"); return CE_Failure; } + if (bUnscale) + { + CPLStringList oArgs; + oArgs.AddString("-unscale"); + oArgs.AddString("-ot"); + oArgs.AddString("Float64"); + oArgs.AddString("-of"); + oArgs.AddString("VRT"); + oArgs.AddString("-a_nodata"); + oArgs.AddString("nan"); + auto *poArgs = GDALTranslateOptionsNew(oArgs.List(), nullptr); + int pbUsageError; + CPLAssert(poArgs); + m_poVRTSrcDS = std::move(m_poSrcDS); + m_poSrcDS.reset(GDALDataset::FromHandle( + GDALTranslate("", m_poVRTSrcDS.get(), poArgs, &pbUsageError))); + GDALTranslateOptionsFree(poArgs); + + if (pbUsageError || !m_poSrcDS) + { + return CE_Failure; + } + } + if (nRasterXSize == 0 && nRasterYSize == 0) { nRasterXSize = m_poSrcDS->GetRasterXSize(); @@ -479,23 +504,17 @@ CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree, return CE_Failure; } - const auto eInDT = - m_bUnscale ? GDT_Float64 : poSrcFirstBand->GetRasterDataType(); - if (!m_bUnscale) + const auto eInDT = poSrcFirstBand->GetRasterDataType(); + for (int i = 1; i < m_poSrcDS->GetRasterCount(); ++i) { - for (int i = 1; i < m_poSrcDS->GetRasterCount(); ++i) + const auto eDT = m_poSrcDS->GetRasterBand(i + 1)->GetRasterDataType(); + if (eDT != eInDT) { - const auto eDT = - m_poSrcDS->GetRasterBand(i + 1)->GetRasterDataType(); - if (eDT != eInDT) - { - CPLError( - CE_Warning, CPLE_AppDefined, - "Not all bands of the input dataset have the same data " - "type. The data type of the first band will be used as " - "the reference one."); - break; - } + CPLError(CE_Warning, CPLE_AppDefined, + "Not all bands of the input dataset have the same data " + "type. The data type of the first band will be used as " + "the reference one."); + break; } } @@ -503,21 +522,13 @@ CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree, int nCurrentBandCount = m_poSrcDS->GetRasterCount(); std::vector adfNoData; - if (m_bUnscale) - { - adfNoData.resize(nCurrentBandCount, - std::numeric_limits::quiet_NaN()); - } - else + for (int i = 1; i <= nCurrentBandCount; ++i) { - for (int i = 1; i <= nCurrentBandCount; ++i) - { - int bHasVal = FALSE; - const double dfVal = - m_poSrcDS->GetRasterBand(i)->GetNoDataValue(&bHasVal); - adfNoData.emplace_back( - bHasVal ? dfVal : std::numeric_limits::quiet_NaN()); - } + int bHasVal = FALSE; + const double dfVal = + m_poSrcDS->GetRasterBand(i)->GetNoDataValue(&bHasVal); + adfNoData.emplace_back( + bHasVal ? dfVal : std::numeric_limits::quiet_NaN()); } int nStepCount = 0; @@ -710,18 +721,10 @@ CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree, { const auto poSrcBand = m_poSrcDS->GetRasterBand(i + 1); const GDALDataType eOutputBandType = - m_bUnscale ? GDT_Float64 - : GetOutputBandType(poSrcBand->GetRasterDataType()); + GetOutputBandType(poSrcBand->GetRasterDataType()); auto poBand = new VRTProcessedRasterBand(this, i + 1, eOutputBandType); poBand->CopyCommonInfoFrom(poSrcBand); - if (m_bUnscale) - { - poBand->SetScale(1); - poBand->SetOffset(0); - poBand->SetNoDataValue( - std::numeric_limits::quiet_NaN()); - } SetBand(i + 1, poBand); } } @@ -1201,8 +1204,7 @@ bool VRTProcessedDataset::ProcessRegion(int nXOff, int nYOff, int nBufXSize, const char *pszInterleave = m_poSrcDS->GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE"); - if (m_bUnscale || (nFirstBandCount > 1 && - (!pszInterleave || EQUAL(pszInterleave, "BAND")))) + if (nFirstBandCount > 1 && (!pszInterleave || EQUAL(pszInterleave, "BAND"))) { // If there are several bands and the source dataset organization // is apparently band interleaved, then first acquire data in @@ -1210,9 +1212,7 @@ bool VRTProcessedDataset::ProcessRegion(int nXOff, int nYOff, int nBufXSize, // data type. // And then transpose it and convert it to the expected data type // of the first step. - const auto eSrcDT = - m_bUnscale ? eFirstDT - : m_poSrcDS->GetRasterBand(1)->GetRasterDataType(); + const auto eSrcDT = m_poSrcDS->GetRasterBand(1)->GetRasterDataType(); try { abyInput.resize(nPixelCount * nFirstBandCount * @@ -1255,40 +1255,6 @@ bool VRTProcessedDataset::ProcessRegion(int nXOff, int nYOff, int nBufXSize, if (!bOK) return false; - if (m_bUnscale) - { - double *pdfValue = reinterpret_cast(abyInput.data()); - for (int nBand = 0; nBand < nFirstBandCount; nBand++) - { - auto poBand = m_poSrcDS->GetRasterBand(nBand + 1); - const double dfScale = poBand->GetScale(nullptr); - const double dfOffset = poBand->GetOffset(nullptr); - int bHasNoData; - const double dfNoData = poBand->GetNoDataValue(&bHasNoData); - const double *pdfEnd = pdfValue + nPixelCount; - - if (bHasNoData) - { - while (pdfValue != pdfEnd) - { - *pdfValue = - (*pdfValue != dfNoData) - ? (*pdfValue * dfScale + dfOffset) - : std::numeric_limits::quiet_NaN(); - pdfValue++; - } - } - else - { - while (pdfValue != pdfEnd) - { - *pdfValue = *pdfValue * dfScale + dfOffset; - pdfValue++; - } - } - } - } - CPLDebugOnly("VRT", "ProcessRegion(): start GDALTranspose2D()"); GDALTranspose2D(abyInput.data(), eSrcDT, abyOutput.data(), eFirstDT, static_cast(nBufXSize) * nBufYSize, From d514839baeee342427d51719de15ae0c0b131fa3 Mon Sep 17 00:00:00 2001 From: Daniel Baston Date: Thu, 9 Jan 2025 19:59:03 -0500 Subject: [PATCH 3/4] vrtprocesseddataset.cpp: Add cppcheck suppression --- frmts/vrt/vrtprocesseddataset.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/frmts/vrt/vrtprocesseddataset.cpp b/frmts/vrt/vrtprocesseddataset.cpp index 5591039fe430..ce84a6ac5c57 100644 --- a/frmts/vrt/vrtprocesseddataset.cpp +++ b/frmts/vrt/vrtprocesseddataset.cpp @@ -316,7 +316,9 @@ CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree, auto *poArgs = GDALTranslateOptionsNew(oArgs.List(), nullptr); int pbUsageError; CPLAssert(poArgs); - m_poVRTSrcDS = std::move(m_poSrcDS); + m_poVRTSrcDS.reset(m_poSrcDS.release()); + // https://trac.cppcheck.net/ticket/11325 + // cppcheck-suppress accessMoved m_poSrcDS.reset(GDALDataset::FromHandle( GDALTranslate("", m_poVRTSrcDS.get(), poArgs, &pbUsageError))); GDALTranslateOptionsFree(poArgs); From 5c60422f6384b39f4713f25394b49b69220d2ce6 Mon Sep 17 00:00:00 2001 From: Daniel Baston Date: Fri, 10 Jan 2025 09:56:31 -0500 Subject: [PATCH 4/4] VRT: Update xsd to allow Input unscale attribute --- .../processed_OutputBands_FROM_LAST_STEP.vrt | 2 +- autotest/gdrivers/vrtprocesseddataset.py | 21 +++++++++++++++++++ frmts/vrt/data/gdalvrt.xsd | 12 ++++++++++- 3 files changed, 33 insertions(+), 2 deletions(-) diff --git a/autotest/gdrivers/data/vrt/processed_OutputBands_FROM_LAST_STEP.vrt b/autotest/gdrivers/data/vrt/processed_OutputBands_FROM_LAST_STEP.vrt index 0f68c6c8e697..c8284f55c73b 100644 --- a/autotest/gdrivers/data/vrt/processed_OutputBands_FROM_LAST_STEP.vrt +++ b/autotest/gdrivers/data/vrt/processed_OutputBands_FROM_LAST_STEP.vrt @@ -1,5 +1,5 @@ - + ../byte.tif diff --git a/autotest/gdrivers/vrtprocesseddataset.py b/autotest/gdrivers/vrtprocesseddataset.py index 0a31d8f3cc8f..cf6c5f038351 100755 --- a/autotest/gdrivers/vrtprocesseddataset.py +++ b/autotest/gdrivers/vrtprocesseddataset.py @@ -10,11 +10,15 @@ # SPDX-License-Identifier: MIT ############################################################################### +import os + import gdaltest import pytest from osgeo import gdal +from .vrtderived import _validate + pytestmark = pytest.mark.skipif( not gdaltest.vrt_has_open_support(), reason="VRT driver open missing", @@ -1529,6 +1533,23 @@ def test_vrtprocesseddataset_RasterIO(tmp_vsimem): ds.ReadAsArray() +############################################################################### +# Validate processed datasets according to xsd + + +@pytest.mark.parametrize( + "fname", + [ + f + for f in os.listdir(os.path.join(os.path.dirname(__file__), "data/vrt")) + if f.startswith("processed") + ], +) +def test_vrt_processeddataset_validate(fname): + with open(os.path.join("data/vrt", fname)) as f: + _validate(f.read()) + + ############################################################################### # Test reading input datasets with scale and offset diff --git a/frmts/vrt/data/gdalvrt.xsd b/frmts/vrt/data/gdalvrt.xsd index ad9d0aa63902..6cc250615091 100644 --- a/frmts/vrt/data/gdalvrt.xsd +++ b/frmts/vrt/data/gdalvrt.xsd @@ -198,6 +198,16 @@ + + + + YES, NO, or AUTO. + If not specified, AUTO is the default and will result in + unscaling all input bands to Float64 if any input band has + a defined scale/offset. + + + @@ -249,7 +259,7 @@ - Allowed names are specific of each processing function + Allowed names are specific to each processing function