Line data    Source code

       1             : /******************************************************************************
       2             :  * Name:     gdalmultidim_gltorthorectification.cpp
       3             :  * Project:  GDAL Core
       4             :  * Purpose:  GLT orthorectification implementation
       5             :  * Author:   Even Rouault <even.rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2023, Even Rouault <even.rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdal_priv.h"
      14             : #include "gdal_pam.h"
      15             : 
      16             : #include <algorithm>
      17             : #include <limits>
      18             : 
      19             : /************************************************************************/
      20             : /*                        GLTOrthoRectifiedArray                        */
      21             : /************************************************************************/
      22             : 
      23             : class GLTOrthoRectifiedArray final : public GDALPamMDArray
      24             : {
      25             :   private:
      26             :     std::shared_ptr<GDALMDArray> m_poParent{};
      27             :     std::vector<std::shared_ptr<GDALDimension>> m_apoDims;
      28             :     std::vector<GUInt64> m_anBlockSize;
      29             :     GDALExtendedDataType m_dt;
      30             :     std::shared_ptr<OGRSpatialReference> m_poSRS{};
      31             :     std::shared_ptr<GDALMDArray> m_poVarX{};
      32             :     std::shared_ptr<GDALMDArray> m_poVarY{};
      33             :     std::shared_ptr<GDALMDArray> m_poGLTX{};
      34             :     std::shared_ptr<GDALMDArray> m_poGLTY{};
      35             :     int m_nGLTIndexOffset = 0;
      36             :     std::vector<GByte> m_abyBandValidity{};
      37             : 
      38             :   protected:
      39           9 :     GLTOrthoRectifiedArray(
      40             :         const std::shared_ptr<GDALMDArray> &poParent,
      41             :         const std::vector<std::shared_ptr<GDALDimension>> &apoDims,
      42             :         const std::vector<GUInt64> &anBlockSize)
      43          18 :         : GDALAbstractMDArray(std::string(), "GLTOrthoRectifiedArray view of " +
      44           9 :                                                  poParent->GetFullName()),
      45          18 :           GDALPamMDArray(std::string(),
      46          18 :                          "GLTOrthoRectifiedArray view of " +
      47           9 :                              poParent->GetFullName(),
      48          18 :                          GDALPamMultiDim::GetPAM(poParent)),
      49           9 :           m_poParent(std::move(poParent)), m_apoDims(apoDims),
      50          45 :           m_anBlockSize(anBlockSize), m_dt(m_poParent->GetDataType())
      51             :     {
      52           9 :         CPLAssert(apoDims.size() == m_poParent->GetDimensionCount());
      53           9 :         CPLAssert(anBlockSize.size() == m_poParent->GetDimensionCount());
      54           9 :     }
      55             : 
      56             :     bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
      57             :                const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
      58             :                const GDALExtendedDataType &bufferDataType,
      59             :                void *pDstBuffer) const override;
      60             : 
      61             :   public:
      62             :     static std::shared_ptr<GDALMDArray>
      63           9 :     Create(const std::shared_ptr<GDALMDArray> &poParent,
      64             :            const std::shared_ptr<GDALGroup> &poRootGroup,
      65             :            const std::shared_ptr<GDALMDArray> &poGLTX,
      66             :            const std::shared_ptr<GDALMDArray> &poGLTY, int nGLTIndexOffset,
      67             :            const std::vector<double> &adfGeoTransform,
      68             :            CSLConstList papszOptions)
      69             :     {
      70          18 :         std::vector<std::shared_ptr<GDALDimension>> apoNewDims;
      71             : 
      72             :         auto poDimY = std::make_shared<GDALDimensionWeakIndexingVar>(
      73           0 :             std::string(), "lat", GDAL_DIM_TYPE_HORIZONTAL_Y, "NORTH",
      74          18 :             poGLTX->GetDimensions()[0]->GetSize());
      75             :         auto varY = GDALMDArrayRegularlySpaced::Create(
      76          27 :             std::string(), poDimY->GetName(), poDimY,
      77          36 :             adfGeoTransform[3] + adfGeoTransform[5] / 2, adfGeoTransform[5], 0);
      78           9 :         poDimY->SetIndexingVariable(varY);
      79           9 :         apoNewDims.emplace_back(poDimY);
      80             : 
      81             :         auto poDimX = std::make_shared<GDALDimensionWeakIndexingVar>(
      82           0 :             std::string(), "lon", GDAL_DIM_TYPE_HORIZONTAL_X, "EAST",
      83          18 :             poGLTX->GetDimensions()[1]->GetSize());
      84             :         auto varX = GDALMDArrayRegularlySpaced::Create(
      85          27 :             std::string(), poDimX->GetName(), poDimX,
      86          36 :             adfGeoTransform[0] + adfGeoTransform[1] / 2, adfGeoTransform[1], 0);
      87           9 :         poDimX->SetIndexingVariable(varX);
      88           9 :         apoNewDims.emplace_back(poDimX);
      89             : 
      90           9 :         if (poParent->GetDimensionCount() == 3)
      91           8 :             apoNewDims.emplace_back(poParent->GetDimensions()[2]);
      92             : 
      93             :         std::vector<GUInt64> anBlockSize = std::vector<GUInt64>{
      94           9 :             std::min<GUInt64>(apoNewDims[0]->GetSize(), 512),
      95          18 :             std::min<GUInt64>(apoNewDims[1]->GetSize(), 512)};
      96           9 :         if (poParent->GetDimensionCount() == 3)
      97             :         {
      98           8 :             anBlockSize.push_back(poParent->GetDimensions()[2]->GetSize());
      99             :         }
     100             : 
     101             :         auto newAr(std::shared_ptr<GLTOrthoRectifiedArray>(
     102          18 :             new GLTOrthoRectifiedArray(poParent, apoNewDims, anBlockSize)));
     103           9 :         newAr->SetSelf(newAr);
     104           9 :         newAr->m_poVarX = varX;
     105           9 :         newAr->m_poVarY = varY;
     106           9 :         newAr->m_poGLTX = poGLTX;
     107           9 :         newAr->m_poGLTY = poGLTY;
     108           9 :         newAr->m_nGLTIndexOffset = nGLTIndexOffset;
     109          18 :         OGRSpatialReference oSRS;
     110           9 :         oSRS.importFromEPSG(4326);
     111           9 :         newAr->m_poSRS.reset(oSRS.Clone());
     112           9 :         newAr->m_poSRS->SetDataAxisToSRSAxisMapping(
     113             :             {/*latIdx = */ 1, /* lonIdx = */ 2});
     114           9 :         if (CPLTestBool(CSLFetchNameValueDef(papszOptions,
     115          17 :                                              "USE_GOOD_WAVELENGTHS", "YES")) &&
     116           8 :             newAr->m_poParent->GetDimensionCount() == 3)
     117             :         {
     118             :             const auto poGoodWaveLengths = poRootGroup->OpenMDArrayFromFullname(
     119          21 :                 "/sensor_band_parameters/good_wavelengths");
     120          10 :             if (poGoodWaveLengths &&
     121          13 :                 poGoodWaveLengths->GetDimensionCount() == 1 &&
     122           3 :                 poGoodWaveLengths->GetDimensions()[0]->GetSize() ==
     123           6 :                     newAr->m_poParent->GetDimensions()[2]->GetSize() &&
     124           3 :                 poGoodWaveLengths->GetDimensions()[0]->GetSize() <
     125          10 :                     1000 * 1000 &&
     126           3 :                 poGoodWaveLengths->GetDataType().GetClass() == GEDTC_NUMERIC)
     127             :             {
     128           3 :                 const GUInt64 arrayStartIdx[] = {0};
     129             :                 const size_t count[] = {static_cast<size_t>(
     130           3 :                     poGoodWaveLengths->GetDimensions()[0]->GetSize())};
     131           3 :                 const GInt64 arrayStep[] = {1};
     132           3 :                 const GPtrDiff_t bufferStride[] = {1};
     133           3 :                 newAr->m_abyBandValidity.resize(count[0], 1);
     134           6 :                 CPL_IGNORE_RET_VAL(poGoodWaveLengths->Read(
     135             :                     arrayStartIdx, count, arrayStep, bufferStride,
     136           6 :                     GDALExtendedDataType::Create(GDT_Byte),
     137           3 :                     newAr->m_abyBandValidity.data()));
     138             :             }
     139             :         }
     140          18 :         return newAr;
     141             :     }
     142             : 
     143           2 :     bool IsWritable() const override
     144             :     {
     145           2 :         return false;
     146             :     }
     147             : 
     148          18 :     const std::string &GetFilename() const override
     149             :     {
     150          18 :         return m_poParent->GetFilename();
     151             :     }
     152             : 
     153             :     const std::vector<std::shared_ptr<GDALDimension>> &
     154         105 :     GetDimensions() const override
     155             :     {
     156         105 :         return m_apoDims;
     157             :     }
     158             : 
     159          77 :     const GDALExtendedDataType &GetDataType() const override
     160             :     {
     161          77 :         return m_dt;
     162             :     }
     163             : 
     164           4 :     std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override
     165             :     {
     166           4 :         return m_poSRS;
     167             :     }
     168             : 
     169           5 :     std::vector<GUInt64> GetBlockSize() const override
     170             :     {
     171           5 :         return m_anBlockSize;
     172             :     }
     173             : 
     174             :     std::shared_ptr<GDALAttribute>
     175           3 :     GetAttribute(const std::string &osName) const override
     176             :     {
     177           3 :         return m_poParent->GetAttribute(osName);
     178             :     }
     179             : 
     180             :     std::vector<std::shared_ptr<GDALAttribute>>
     181           2 :     GetAttributes(CSLConstList papszOptions = nullptr) const override
     182             :     {
     183           2 :         return m_poParent->GetAttributes(papszOptions);
     184             :     }
     185             : 
     186           3 :     const std::string &GetUnit() const override
     187             :     {
     188           3 :         return m_poParent->GetUnit();
     189             :     }
     190             : 
     191          19 :     const void *GetRawNoDataValue() const override
     192             :     {
     193          19 :         return m_poParent->GetRawNoDataValue();
     194             :     }
     195             : 
     196           0 :     double GetOffset(bool *pbHasOffset,
     197             :                      GDALDataType *peStorageType) const override
     198             :     {
     199           0 :         return m_poParent->GetOffset(pbHasOffset, peStorageType);
     200             :     }
     201             : 
     202           0 :     double GetScale(bool *pbHasScale,
     203             :                     GDALDataType *peStorageType) const override
     204             :     {
     205           0 :         return m_poParent->GetScale(pbHasScale, peStorageType);
     206             :     }
     207             : };
     208             : 
     209             : /************************************************************************/
     210             : /*                   GDALMDArrayResampled::IRead()                      */
     211             : /************************************************************************/
     212             : 
     213          16 : bool GLTOrthoRectifiedArray::IRead(const GUInt64 *arrayStartIdx,
     214             :                                    const size_t *count, const GInt64 *arrayStep,
     215             :                                    const GPtrDiff_t *bufferStride,
     216             :                                    const GDALExtendedDataType &bufferDataType,
     217             :                                    void *pDstBuffer) const
     218             : {
     219          16 :     if (bufferDataType.GetClass() != GEDTC_NUMERIC)
     220           0 :         return false;
     221             : 
     222          16 :     const auto eBufferDT = bufferDataType.GetNumericDataType();
     223          16 :     auto pRawNoDataValue = GetRawNoDataValue();
     224          32 :     std::vector<GByte> abyNoData(16);
     225          16 :     if (pRawNoDataValue)
     226          16 :         GDALCopyWords(pRawNoDataValue, GetDataType().GetNumericDataType(), 0,
     227          16 :                       abyNoData.data(), eBufferDT, 0, 1);
     228             : 
     229          16 :     const auto nBufferDTSize = bufferDataType.GetSize();
     230             :     const int nCopyWordsDstStride =
     231          16 :         m_apoDims.size() == 3
     232          16 :             ? static_cast<int>(bufferStride[2] * nBufferDTSize)
     233          16 :             : 0;
     234             :     const int nCopyWordsCount =
     235          16 :         m_apoDims.size() == 3 ? static_cast<int>(count[2]) : 1;
     236             : 
     237             :     const auto FillBufferWithNodata =
     238           5 :         [count, bufferStride, pDstBuffer, nBufferDTSize, &eBufferDT,
     239          49 :          nCopyWordsDstStride, nCopyWordsCount, &abyNoData]()
     240             :     {
     241          12 :         for (size_t iY = 0; iY < count[0]; ++iY)
     242             :         {
     243          22 :             for (size_t iX = 0; iX < count[1]; ++iX)
     244             :             {
     245          15 :                 GByte *pabyDstBuffer =
     246             :                     static_cast<GByte *>(pDstBuffer) +
     247          15 :                     (iY * bufferStride[0] + iX * bufferStride[1]) *
     248          15 :                         static_cast<int>(nBufferDTSize);
     249          15 :                 GDALCopyWords(abyNoData.data(), eBufferDT, 0, pabyDstBuffer,
     250             :                               eBufferDT, nCopyWordsDstStride, nCopyWordsCount);
     251             :             }
     252             :         }
     253          21 :     };
     254             : 
     255          16 :     bool bSomeBandInvalids = false;
     256          16 :     if (!m_abyBandValidity.empty())
     257             :     {
     258           3 :         size_t nCountInvalid = 0;
     259           7 :         for (size_t i = 0; i < count[2]; ++i)
     260             :         {
     261           4 :             const size_t iBand =
     262           4 :                 static_cast<size_t>(arrayStartIdx[2] + i * arrayStep[2]);
     263           4 :             CPLAssert(iBand < m_abyBandValidity.size());
     264           4 :             if (!m_abyBandValidity[iBand])
     265           2 :                 nCountInvalid++;
     266             :         }
     267           3 :         if (nCountInvalid == count[2])
     268             :         {
     269           1 :             FillBufferWithNodata();
     270           1 :             return true;
     271             :         }
     272           2 :         bSomeBandInvalids = true;
     273             :     }
     274             : 
     275          15 :     const size_t nXYValsCount = count[0] * count[1];
     276          30 :     const auto eInt32DT = GDALExtendedDataType::Create(GDT_Int32);
     277          30 :     std::vector<int32_t> anGLTX;
     278          30 :     std::vector<int32_t> anGLTY;
     279             :     try
     280             :     {
     281          15 :         anGLTX.resize(nXYValsCount);
     282          15 :         anGLTY.resize(nXYValsCount);
     283             :     }
     284           0 :     catch (const std::bad_alloc &e)
     285             :     {
     286           0 :         CPLError(CE_Failure, CPLE_OutOfMemory,
     287           0 :                  "GLTOrthoRectifiedArray::IRead(): %s", e.what());
     288           0 :         return false;
     289             :     }
     290          30 :     if (!m_poGLTX->Read(arrayStartIdx, count, arrayStep, nullptr, eInt32DT,
     291          45 :                         anGLTX.data()) ||
     292          30 :         !m_poGLTY->Read(arrayStartIdx, count, arrayStep, nullptr, eInt32DT,
     293          15 :                         anGLTY.data()))
     294             :     {
     295           0 :         return false;
     296             :     }
     297             : 
     298          15 :     int32_t nMinX = std::numeric_limits<int32_t>::max();
     299          15 :     int32_t nMaxX = std::numeric_limits<int32_t>::min();
     300          15 :     const auto nXSize = m_poParent->GetDimensions()[0]->GetSize();
     301         104 :     for (size_t i = 0; i < nXYValsCount; ++i)
     302             :     {
     303             :         const int64_t nX64 =
     304          89 :             static_cast<int64_t>(anGLTX[i]) + m_nGLTIndexOffset;
     305          89 :         if (nX64 >= 0 && static_cast<uint64_t>(nX64) < nXSize)
     306             :         {
     307          38 :             const int32_t nX = static_cast<int32_t>(nX64);
     308          38 :             if (nX < nMinX)
     309          11 :                 nMinX = nX;
     310          38 :             if (nX > nMaxX)
     311          20 :                 nMaxX = nX;
     312             :         }
     313             :     }
     314             : 
     315          15 :     int32_t nMinY = std::numeric_limits<int32_t>::max();
     316          15 :     int32_t nMaxY = std::numeric_limits<int32_t>::min();
     317          15 :     const auto nYSize = m_poParent->GetDimensions()[0]->GetSize();
     318         104 :     for (size_t i = 0; i < nXYValsCount; ++i)
     319             :     {
     320             :         const int64_t nY64 =
     321          89 :             static_cast<int64_t>(anGLTY[i]) + m_nGLTIndexOffset;
     322          89 :         if (nY64 >= 0 && static_cast<uint64_t>(nY64) < nYSize)
     323             :         {
     324          38 :             const int32_t nY = static_cast<int32_t>(nY64);
     325          38 :             if (nY < nMinY)
     326          20 :                 nMinY = nY;
     327          38 :             if (nY > nMaxY)
     328          11 :                 nMaxY = nY;
     329             :         }
     330             :     }
     331             : 
     332          15 :     if (nMinX > nMaxX || nMinY > nMaxY)
     333             :     {
     334           4 :         FillBufferWithNodata();
     335           4 :         return true;
     336             :     }
     337             : 
     338             :     GUInt64 parentArrayIdxStart[3] = {
     339          11 :         static_cast<GUInt64>(nMinY), static_cast<GUInt64>(nMinX),
     340          11 :         m_apoDims.size() == 3 ? arrayStartIdx[2] : 0};
     341          11 :     size_t parentCount[3] = {static_cast<size_t>(nMaxY - nMinY + 1),
     342          11 :                              static_cast<size_t>(nMaxX - nMinX + 1),
     343          11 :                              m_apoDims.size() == 3 ? count[2] : 1};
     344          11 :     GInt64 parentArrayStep[3] = {1, 1,
     345          11 :                                  m_apoDims.size() == 3 ? arrayStep[2] : 0};
     346             : 
     347          11 :     size_t nParentValueSize = nBufferDTSize;
     348          44 :     for (int i = 0; i < 3; ++i)
     349             :     {
     350          66 :         if (parentCount[i] >
     351          33 :             std::numeric_limits<size_t>::max() / nParentValueSize)
     352             :         {
     353           0 :             CPLError(
     354             :                 CE_Failure, CPLE_OutOfMemory,
     355             :                 "GLTOrthoRectifiedArray::IRead(): too big temporary array");
     356           0 :             return false;
     357             :         }
     358          33 :         nParentValueSize *= parentCount[i];
     359             :     }
     360             : 
     361          11 :     GPtrDiff_t parentStride[3] = {
     362          11 :         static_cast<GPtrDiff_t>(parentCount[1] * parentCount[2]),
     363          11 :         static_cast<GPtrDiff_t>(parentCount[2]), 1};
     364          22 :     std::vector<GByte> parentValues;
     365             :     try
     366             :     {
     367          11 :         parentValues.resize(nParentValueSize);
     368             :     }
     369           0 :     catch (const std::bad_alloc &e)
     370             :     {
     371           0 :         CPLError(CE_Failure, CPLE_OutOfMemory,
     372           0 :                  "GLTOrthoRectifiedArray::IRead(): %s", e.what());
     373           0 :         return false;
     374             :     }
     375          22 :     if (!m_poParent->Read(parentArrayIdxStart, parentCount, parentArrayStep,
     376          11 :                           parentStride, bufferDataType, parentValues.data()))
     377             :     {
     378           0 :         return false;
     379             :     }
     380             : 
     381          11 :     size_t iGLTIndex = 0;
     382          11 :     const size_t nXCount = parentCount[1];
     383          11 :     const size_t nBandCount = m_apoDims.size() == 3 ? parentCount[2] : 1;
     384          40 :     for (size_t iY = 0; iY < count[0]; ++iY)
     385             :     {
     386         112 :         for (size_t iX = 0; iX < count[1]; ++iX, ++iGLTIndex)
     387             :         {
     388             :             const int64_t nX64 =
     389          83 :                 static_cast<int64_t>(anGLTX[iGLTIndex]) + m_nGLTIndexOffset;
     390             :             const int64_t nY64 =
     391          83 :                 static_cast<int64_t>(anGLTY[iGLTIndex]) + m_nGLTIndexOffset;
     392          83 :             GByte *pabyDstBuffer =
     393             :                 static_cast<GByte *>(pDstBuffer) +
     394          83 :                 (iY * bufferStride[0] + iX * bufferStride[1]) *
     395          83 :                     static_cast<int>(nBufferDTSize);
     396          83 :             if (nX64 >= nMinX && nX64 <= nMaxX && nY64 >= nMinY &&
     397          38 :                 nY64 <= nMaxY)
     398             :             {
     399          38 :                 const int32_t iSrcX = static_cast<int32_t>(nX64) - nMinX;
     400          38 :                 const int32_t iSrcY = static_cast<int32_t>(nY64) - nMinY;
     401             :                 const GByte *pabySrcBuffer =
     402          38 :                     parentValues.data() +
     403          38 :                     (iSrcY * nXCount + iSrcX) * nBandCount * nBufferDTSize;
     404          38 :                 GDALCopyWords(pabySrcBuffer, eBufferDT,
     405             :                               static_cast<int>(nBufferDTSize), pabyDstBuffer,
     406             :                               eBufferDT, nCopyWordsDstStride, nCopyWordsCount);
     407          38 :                 if (bSomeBandInvalids)
     408             :                 {
     409          20 :                     for (size_t i = 0; i < count[2]; ++i)
     410             :                     {
     411          12 :                         const size_t iBand = static_cast<size_t>(
     412          12 :                             arrayStartIdx[2] + i * arrayStep[2]);
     413          12 :                         if (!m_abyBandValidity[iBand])
     414             :                         {
     415           4 :                             GDALCopyWords(abyNoData.data(), eBufferDT, 0,
     416           4 :                                           pabyDstBuffer +
     417           4 :                                               i * nCopyWordsDstStride,
     418             :                                           eBufferDT, 0, 1);
     419             :                         }
     420             :                     }
     421          38 :                 }
     422             :             }
     423             :             else
     424             :             {
     425          45 :                 GDALCopyWords(abyNoData.data(), eBufferDT, 0, pabyDstBuffer,
     426             :                               eBufferDT, nCopyWordsDstStride, nCopyWordsCount);
     427             :             }
     428             :         }
     429             :     }
     430             : 
     431          11 :     return true;
     432             : }
     433             : 
     434             : /************************************************************************/
     435             : /*                     CreateGLTOrthorectified()                        */
     436             : /************************************************************************/
     437             : 
     438             : //! @cond Doxygen_Suppress
     439             : 
     440           9 : /* static */ std::shared_ptr<GDALMDArray> GDALMDArray::CreateGLTOrthorectified(
     441             :     const std::shared_ptr<GDALMDArray> &poParent,
     442             :     const std::shared_ptr<GDALGroup> &poRootGroup,
     443             :     const std::shared_ptr<GDALMDArray> &poGLTX,
     444             :     const std::shared_ptr<GDALMDArray> &poGLTY, int nGLTIndexOffset,
     445             :     const std::vector<double> &adfGeoTransform, CSLConstList papszOptions)
     446             : {
     447             :     return GLTOrthoRectifiedArray::Create(poParent, poRootGroup, poGLTX, poGLTY,
     448             :                                           nGLTIndexOffset, adfGeoTransform,
     449           9 :                                           papszOptions);
     450             : }
     451             : 
     452             : //! @endcond