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 :         GDALCopyWords64(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 :                 GDALCopyWords64(abyNoData.data(), eBufferDT, 0, pabyDstBuffer,
     250             :                                 eBufferDT, nCopyWordsDstStride,
     251             :                                 nCopyWordsCount);
     252             :             }
     253             :         }
     254          21 :     };
     255             : 
     256          16 :     bool bSomeBandInvalids = false;
     257          16 :     if (!m_abyBandValidity.empty())
     258             :     {
     259           3 :         size_t nCountInvalid = 0;
     260           7 :         for (size_t i = 0; i < count[2]; ++i)
     261             :         {
     262           4 :             const size_t iBand =
     263           4 :                 static_cast<size_t>(arrayStartIdx[2] + i * arrayStep[2]);
     264           4 :             CPLAssert(iBand < m_abyBandValidity.size());
     265           4 :             if (!m_abyBandValidity[iBand])
     266           2 :                 nCountInvalid++;
     267             :         }
     268           3 :         if (nCountInvalid == count[2])
     269             :         {
     270           1 :             FillBufferWithNodata();
     271           1 :             return true;
     272             :         }
     273           2 :         bSomeBandInvalids = true;
     274             :     }
     275             : 
     276          15 :     const size_t nXYValsCount = count[0] * count[1];
     277          30 :     const auto eInt32DT = GDALExtendedDataType::Create(GDT_Int32);
     278          30 :     std::vector<int32_t> anGLTX;
     279          30 :     std::vector<int32_t> anGLTY;
     280             :     try
     281             :     {
     282          15 :         anGLTX.resize(nXYValsCount);
     283          15 :         anGLTY.resize(nXYValsCount);
     284             :     }
     285           0 :     catch (const std::bad_alloc &e)
     286             :     {
     287           0 :         CPLError(CE_Failure, CPLE_OutOfMemory,
     288           0 :                  "GLTOrthoRectifiedArray::IRead(): %s", e.what());
     289           0 :         return false;
     290             :     }
     291          30 :     if (!m_poGLTX->Read(arrayStartIdx, count, arrayStep, nullptr, eInt32DT,
     292          45 :                         anGLTX.data()) ||
     293          30 :         !m_poGLTY->Read(arrayStartIdx, count, arrayStep, nullptr, eInt32DT,
     294          15 :                         anGLTY.data()))
     295             :     {
     296           0 :         return false;
     297             :     }
     298             : 
     299          15 :     int32_t nMinX = std::numeric_limits<int32_t>::max();
     300          15 :     int32_t nMaxX = std::numeric_limits<int32_t>::min();
     301          15 :     const auto nXSize = m_poParent->GetDimensions()[0]->GetSize();
     302         104 :     for (size_t i = 0; i < nXYValsCount; ++i)
     303             :     {
     304             :         const int64_t nX64 =
     305          89 :             static_cast<int64_t>(anGLTX[i]) + m_nGLTIndexOffset;
     306          89 :         if (nX64 >= 0 && static_cast<uint64_t>(nX64) < nXSize)
     307             :         {
     308          38 :             const int32_t nX = static_cast<int32_t>(nX64);
     309          38 :             if (nX < nMinX)
     310          11 :                 nMinX = nX;
     311          38 :             if (nX > nMaxX)
     312          20 :                 nMaxX = nX;
     313             :         }
     314             :     }
     315             : 
     316          15 :     int32_t nMinY = std::numeric_limits<int32_t>::max();
     317          15 :     int32_t nMaxY = std::numeric_limits<int32_t>::min();
     318          15 :     const auto nYSize = m_poParent->GetDimensions()[0]->GetSize();
     319         104 :     for (size_t i = 0; i < nXYValsCount; ++i)
     320             :     {
     321             :         const int64_t nY64 =
     322          89 :             static_cast<int64_t>(anGLTY[i]) + m_nGLTIndexOffset;
     323          89 :         if (nY64 >= 0 && static_cast<uint64_t>(nY64) < nYSize)
     324             :         {
     325          38 :             const int32_t nY = static_cast<int32_t>(nY64);
     326          38 :             if (nY < nMinY)
     327          20 :                 nMinY = nY;
     328          38 :             if (nY > nMaxY)
     329          11 :                 nMaxY = nY;
     330             :         }
     331             :     }
     332             : 
     333          15 :     if (nMinX > nMaxX || nMinY > nMaxY)
     334             :     {
     335           4 :         FillBufferWithNodata();
     336           4 :         return true;
     337             :     }
     338             : 
     339             :     GUInt64 parentArrayIdxStart[3] = {
     340          11 :         static_cast<GUInt64>(nMinY), static_cast<GUInt64>(nMinX),
     341          11 :         m_apoDims.size() == 3 ? arrayStartIdx[2] : 0};
     342          11 :     size_t parentCount[3] = {static_cast<size_t>(nMaxY - nMinY + 1),
     343          11 :                              static_cast<size_t>(nMaxX - nMinX + 1),
     344          11 :                              m_apoDims.size() == 3 ? count[2] : 1};
     345          11 :     GInt64 parentArrayStep[3] = {1, 1,
     346          11 :                                  m_apoDims.size() == 3 ? arrayStep[2] : 0};
     347             : 
     348          11 :     size_t nParentValueSize = nBufferDTSize;
     349          44 :     for (int i = 0; i < 3; ++i)
     350             :     {
     351          66 :         if (parentCount[i] >
     352          33 :             std::numeric_limits<size_t>::max() / nParentValueSize)
     353             :         {
     354           0 :             CPLError(
     355             :                 CE_Failure, CPLE_OutOfMemory,
     356             :                 "GLTOrthoRectifiedArray::IRead(): too big temporary array");
     357           0 :             return false;
     358             :         }
     359          33 :         nParentValueSize *= parentCount[i];
     360             :     }
     361             : 
     362          11 :     GPtrDiff_t parentStride[3] = {
     363          11 :         static_cast<GPtrDiff_t>(parentCount[1] * parentCount[2]),
     364          11 :         static_cast<GPtrDiff_t>(parentCount[2]), 1};
     365          22 :     std::vector<GByte> parentValues;
     366             :     try
     367             :     {
     368          11 :         parentValues.resize(nParentValueSize);
     369             :     }
     370           0 :     catch (const std::bad_alloc &e)
     371             :     {
     372           0 :         CPLError(CE_Failure, CPLE_OutOfMemory,
     373           0 :                  "GLTOrthoRectifiedArray::IRead(): %s", e.what());
     374           0 :         return false;
     375             :     }
     376          22 :     if (!m_poParent->Read(parentArrayIdxStart, parentCount, parentArrayStep,
     377          11 :                           parentStride, bufferDataType, parentValues.data()))
     378             :     {
     379           0 :         return false;
     380             :     }
     381             : 
     382          11 :     size_t iGLTIndex = 0;
     383          11 :     const size_t nXCount = parentCount[1];
     384          11 :     const size_t nBandCount = m_apoDims.size() == 3 ? parentCount[2] : 1;
     385          40 :     for (size_t iY = 0; iY < count[0]; ++iY)
     386             :     {
     387         112 :         for (size_t iX = 0; iX < count[1]; ++iX, ++iGLTIndex)
     388             :         {
     389             :             const int64_t nX64 =
     390          83 :                 static_cast<int64_t>(anGLTX[iGLTIndex]) + m_nGLTIndexOffset;
     391             :             const int64_t nY64 =
     392          83 :                 static_cast<int64_t>(anGLTY[iGLTIndex]) + m_nGLTIndexOffset;
     393          83 :             GByte *pabyDstBuffer =
     394             :                 static_cast<GByte *>(pDstBuffer) +
     395          83 :                 (iY * bufferStride[0] + iX * bufferStride[1]) *
     396          83 :                     static_cast<int>(nBufferDTSize);
     397          83 :             if (nX64 >= nMinX && nX64 <= nMaxX && nY64 >= nMinY &&
     398          38 :                 nY64 <= nMaxY)
     399             :             {
     400          38 :                 const int32_t iSrcX = static_cast<int32_t>(nX64) - nMinX;
     401          38 :                 const int32_t iSrcY = static_cast<int32_t>(nY64) - nMinY;
     402             :                 const GByte *pabySrcBuffer =
     403          38 :                     parentValues.data() +
     404          38 :                     (iSrcY * nXCount + iSrcX) * nBandCount * nBufferDTSize;
     405          38 :                 GDALCopyWords64(pabySrcBuffer, eBufferDT,
     406             :                                 static_cast<int>(nBufferDTSize), pabyDstBuffer,
     407             :                                 eBufferDT, nCopyWordsDstStride,
     408             :                                 nCopyWordsCount);
     409          38 :                 if (bSomeBandInvalids)
     410             :                 {
     411          20 :                     for (size_t i = 0; i < count[2]; ++i)
     412             :                     {
     413          12 :                         const size_t iBand = static_cast<size_t>(
     414          12 :                             arrayStartIdx[2] + i * arrayStep[2]);
     415          12 :                         if (!m_abyBandValidity[iBand])
     416             :                         {
     417           4 :                             GDALCopyWords64(abyNoData.data(), eBufferDT, 0,
     418           4 :                                             pabyDstBuffer +
     419           4 :                                                 i * nCopyWordsDstStride,
     420             :                                             eBufferDT, 0, 1);
     421             :                         }
     422             :                     }
     423          38 :                 }
     424             :             }
     425             :             else
     426             :             {
     427          45 :                 GDALCopyWords64(abyNoData.data(), eBufferDT, 0, pabyDstBuffer,
     428             :                                 eBufferDT, nCopyWordsDstStride,
     429             :                                 nCopyWordsCount);
     430             :             }
     431             :         }
     432             :     }
     433             : 
     434          11 :     return true;
     435             : }
     436             : 
     437             : /************************************************************************/
     438             : /*                     CreateGLTOrthorectified()                        */
     439             : /************************************************************************/
     440             : 
     441             : //! @cond Doxygen_Suppress
     442             : 
     443           9 : /* static */ std::shared_ptr<GDALMDArray> GDALMDArray::CreateGLTOrthorectified(
     444             :     const std::shared_ptr<GDALMDArray> &poParent,
     445             :     const std::shared_ptr<GDALGroup> &poRootGroup,
     446             :     const std::shared_ptr<GDALMDArray> &poGLTX,
     447             :     const std::shared_ptr<GDALMDArray> &poGLTY, int nGLTIndexOffset,
     448             :     const std::vector<double> &adfGeoTransform, CSLConstList papszOptions)
     449             : {
     450             :     return GLTOrthoRectifiedArray::Create(poParent, poRootGroup, poGLTX, poGLTY,
     451             :                                           nGLTIndexOffset, adfGeoTransform,
     452           9 :                                           papszOptions);
     453             : }
     454             : 
     455             : //! @endcond