Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Name:     gdalmultidim_gridded.cpp
       4             :  * Project:  GDAL Core
       5             :  * Purpose:  GDALMDArray::GetGridded() implementation
       6             :  * Author:   Even Rouault <even.rouault at spatialys.com>
       7             :  *
       8             :  ******************************************************************************
       9             :  * Copyright (c) 2023, Even Rouault <even.rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "gdal_alg.h"
      15             : #include "gdalgrid.h"
      16             : #include "gdal_priv.h"
      17             : #include "gdal_pam.h"
      18             : #include "ogrsf_frmts.h"
      19             : 
      20             : #include <algorithm>
      21             : #include <cassert>
      22             : #include <limits>
      23             : #include <new>
      24             : 
      25             : /************************************************************************/
      26             : /*                         GDALMDArrayGridded                           */
      27             : /************************************************************************/
      28             : 
      29             : class GDALMDArrayGridded final : public GDALPamMDArray
      30             : {
      31             :   private:
      32             :     std::shared_ptr<GDALMDArray> m_poParent{};
      33             :     std::vector<std::shared_ptr<GDALDimension>> m_apoDims{};
      34             :     std::shared_ptr<GDALMDArray> m_poVarX{};
      35             :     std::shared_ptr<GDALMDArray> m_poVarY{};
      36             :     std::unique_ptr<GDALDataset> m_poVectorDS{};
      37             :     GDALGridAlgorithm m_eAlg;
      38             :     std::unique_ptr<void, VSIFreeReleaser> m_poGridOptions;
      39             :     const GDALExtendedDataType m_dt;
      40             :     std::vector<GUInt64> m_anBlockSize{};
      41             :     const double m_dfNoDataValue;
      42             :     const double m_dfMinX;
      43             :     const double m_dfResX;
      44             :     const double m_dfMinY;
      45             :     const double m_dfResY;
      46             :     const double m_dfRadius;
      47             :     mutable std::vector<GUInt64> m_anLastStartIdx{};
      48             :     mutable std::vector<double> m_adfZ{};
      49             : 
      50             :   protected:
      51           4 :     explicit GDALMDArrayGridded(
      52             :         const std::shared_ptr<GDALMDArray> &poParent,
      53             :         const std::vector<std::shared_ptr<GDALDimension>> &apoDims,
      54             :         const std::shared_ptr<GDALMDArray> &poVarX,
      55             :         const std::shared_ptr<GDALMDArray> &poVarY,
      56             :         std::unique_ptr<GDALDataset> &&poVectorDS, GDALGridAlgorithm eAlg,
      57             :         std::unique_ptr<void, VSIFreeReleaser> &&poGridOptions,
      58             :         double dfNoDataValue, double dfMinX, double dfResX, double dfMinY,
      59             :         double dfResY, double dfRadius)
      60           8 :         : GDALAbstractMDArray(std::string(),
      61           8 :                               "Gridded view of " + poParent->GetFullName()),
      62             :           GDALPamMDArray(
      63           8 :               std::string(), "Gridded view of " + poParent->GetFullName(),
      64           8 :               GDALPamMultiDim::GetPAM(poParent), poParent->GetContext()),
      65           4 :           m_poParent(std::move(poParent)), m_apoDims(apoDims), m_poVarX(poVarX),
      66           4 :           m_poVarY(poVarY), m_poVectorDS(std::move(poVectorDS)), m_eAlg(eAlg),
      67           4 :           m_poGridOptions(std::move(poGridOptions)),
      68             :           m_dt(GDALExtendedDataType::Create(GDT_Float64)),
      69             :           m_dfNoDataValue(dfNoDataValue), m_dfMinX(dfMinX), m_dfResX(dfResX),
      70          20 :           m_dfMinY(dfMinY), m_dfResY(dfResY), m_dfRadius(dfRadius)
      71             :     {
      72           4 :         const auto anParentBlockSize = m_poParent->GetBlockSize();
      73           4 :         m_anBlockSize.resize(m_apoDims.size());
      74          12 :         for (size_t i = 0; i + 1 < m_apoDims.size(); ++i)
      75           8 :             m_anBlockSize[i] = anParentBlockSize[i];
      76           4 :         m_anBlockSize[m_apoDims.size() - 2] = 256;
      77           4 :         m_anBlockSize[m_apoDims.size() - 1] = 256;
      78           4 :     }
      79             : 
      80             :     bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
      81             :                const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
      82             :                const GDALExtendedDataType &bufferDataType,
      83             :                void *pDstBuffer) const override;
      84             : 
      85             :   public:
      86             :     static std::shared_ptr<GDALMDArrayGridded>
      87           4 :     Create(const std::shared_ptr<GDALMDArray> &poParent,
      88             :            const std::vector<std::shared_ptr<GDALDimension>> &apoDims,
      89             :            const std::shared_ptr<GDALMDArray> &poVarX,
      90             :            const std::shared_ptr<GDALMDArray> &poVarY,
      91             :            std::unique_ptr<GDALDataset> &&poVectorDS, GDALGridAlgorithm eAlg,
      92             :            std::unique_ptr<void, VSIFreeReleaser> &&poGridOptions,
      93             :            double dfNoDataValue, double dfMinX, double dfResX, double dfMinY,
      94             :            double dfResY, double dfRadius)
      95             :     {
      96             :         auto newAr(std::shared_ptr<GDALMDArrayGridded>(new GDALMDArrayGridded(
      97           4 :             poParent, apoDims, poVarX, poVarY, std::move(poVectorDS), eAlg,
      98           4 :             std::move(poGridOptions), dfNoDataValue, dfMinX, dfResX, dfMinY,
      99           4 :             dfResY, dfRadius)));
     100           4 :         newAr->SetSelf(newAr);
     101           4 :         return newAr;
     102             :     }
     103             : 
     104           1 :     bool IsWritable() const override
     105             :     {
     106           1 :         return false;
     107             :     }
     108             : 
     109           9 :     const std::string &GetFilename() const override
     110             :     {
     111           9 :         return m_poParent->GetFilename();
     112             :     }
     113             : 
     114             :     const std::vector<std::shared_ptr<GDALDimension>> &
     115          63 :     GetDimensions() const override
     116             :     {
     117          63 :         return m_apoDims;
     118             :     }
     119             : 
     120           1 :     const void *GetRawNoDataValue() const override
     121             :     {
     122           1 :         return &m_dfNoDataValue;
     123             :     }
     124             : 
     125           2 :     std::vector<GUInt64> GetBlockSize() const override
     126             :     {
     127           2 :         return m_anBlockSize;
     128             :     }
     129             : 
     130          25 :     const GDALExtendedDataType &GetDataType() const override
     131             :     {
     132          25 :         return m_dt;
     133             :     }
     134             : 
     135           1 :     const std::string &GetUnit() const override
     136             :     {
     137           1 :         return m_poParent->GetUnit();
     138             :     }
     139             : 
     140           1 :     std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override
     141             :     {
     142           1 :         return m_poParent->GetSpatialRef();
     143             :     }
     144             : 
     145             :     std::shared_ptr<GDALAttribute>
     146           1 :     GetAttribute(const std::string &osName) const override
     147             :     {
     148           1 :         return m_poParent->GetAttribute(osName);
     149             :     }
     150             : 
     151             :     std::vector<std::shared_ptr<GDALAttribute>>
     152           2 :     GetAttributes(CSLConstList papszOptions = nullptr) const override
     153             :     {
     154           2 :         return m_poParent->GetAttributes(papszOptions);
     155             :     }
     156             : };
     157             : 
     158             : /************************************************************************/
     159             : /*                             IRead()                                  */
     160             : /************************************************************************/
     161             : 
     162           7 : bool GDALMDArrayGridded::IRead(const GUInt64 *arrayStartIdx,
     163             :                                const size_t *count, const GInt64 *arrayStep,
     164             :                                const GPtrDiff_t *bufferStride,
     165             :                                const GDALExtendedDataType &bufferDataType,
     166             :                                void *pDstBuffer) const
     167             : {
     168           7 :     if (bufferDataType.GetClass() != GEDTC_NUMERIC)
     169             :     {
     170           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     171             :                  "GDALMDArrayGridded::IRead() only support numeric "
     172             :                  "bufferDataType");
     173           0 :         return false;
     174             :     }
     175           7 :     const auto nDims = GetDimensionCount();
     176             : 
     177          14 :     std::vector<GUInt64> anStartIdx;
     178          13 :     for (size_t i = 0; i + 2 < nDims; ++i)
     179             :     {
     180           7 :         anStartIdx.push_back(arrayStartIdx[i]);
     181           7 :         if (count[i] != 1)
     182             :         {
     183           1 :             CPLError(CE_Failure, CPLE_NotSupported,
     184             :                      "GDALMDArrayGridded::IRead() only support count = 1 in "
     185             :                      "the first dimensions, except the last 2 Y,X ones");
     186           1 :             return false;
     187             :         }
     188             :     }
     189             : 
     190           6 :     const auto iDimX = nDims - 1;
     191           6 :     const auto iDimY = nDims - 2;
     192             : 
     193           6 :     if (arrayStep[iDimX] < 0)
     194             :     {
     195           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     196             :                  "GDALMDArrayGridded::IRead(): "
     197             :                  "arrayStep[iDimX] < 0 not supported");
     198           1 :         return false;
     199             :     }
     200             : 
     201           5 :     if (arrayStep[iDimY] < 0)
     202             :     {
     203           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     204             :                  "GDALMDArrayGridded::IRead(): "
     205             :                  "arrayStep[iDimY] < 0 not supported");
     206           0 :         return false;
     207             :     }
     208             : 
     209             :     // Load the values taken by the variable at the considered slice
     210             :     // (if not already done)
     211           5 :     if (m_adfZ.empty() || m_anLastStartIdx != anStartIdx)
     212             :     {
     213           4 :         std::vector<GUInt64> anTempStartIdx(anStartIdx);
     214           4 :         anTempStartIdx.push_back(0);
     215             :         const std::vector<GInt64> anTempArrayStep(
     216           4 :             m_poParent->GetDimensionCount(), 1);
     217             :         std::vector<GPtrDiff_t> anTempBufferStride(
     218           4 :             m_poParent->GetDimensionCount() - 1, 0);
     219           4 :         anTempBufferStride.push_back(1);
     220           4 :         std::vector<size_t> anTempCount(m_poParent->GetDimensionCount() - 1, 1);
     221           4 :         anTempCount.push_back(
     222           4 :             static_cast<size_t>(m_poParent->GetDimensions().back()->GetSize()));
     223           4 :         CPLAssert(anTempStartIdx.size() == m_poParent->GetDimensionCount());
     224           4 :         CPLAssert(anTempCount.size() == m_poParent->GetDimensionCount());
     225           4 :         CPLAssert(anTempArrayStep.size() == m_poParent->GetDimensionCount());
     226           4 :         CPLAssert(anTempBufferStride.size() == m_poParent->GetDimensionCount());
     227             : 
     228             :         try
     229             :         {
     230           4 :             m_adfZ.resize(anTempCount.back());
     231             :         }
     232           0 :         catch (const std::bad_alloc &e)
     233             :         {
     234           0 :             CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
     235           0 :             return false;
     236             :         }
     237             : 
     238           8 :         if (!m_poParent->Read(anTempStartIdx.data(), anTempCount.data(),
     239           4 :                               anTempArrayStep.data(), anTempBufferStride.data(),
     240           4 :                               m_dt, m_adfZ.data()))
     241             :         {
     242           0 :             return false;
     243             :         }
     244             :         // GCC 13.1 warns here. Definitely a false positive.
     245             : #if defined(__GNUC__) && __GNUC__ >= 13
     246             : #pragma GCC diagnostic push
     247             : #pragma GCC diagnostic ignored "-Wnull-dereference"
     248             : #endif
     249           4 :         m_anLastStartIdx = std::move(anStartIdx);
     250             : #if defined(__GNUC__) && __GNUC__ >= 13
     251             : #pragma GCC diagnostic pop
     252             : #endif
     253             :     }
     254             : 
     255             :     // Determine the X,Y spatial extent of the request
     256           5 :     const double dfX1 = m_dfMinX + arrayStartIdx[iDimX] * m_dfResX;
     257           5 :     const double dfX2 = m_dfMinX + (arrayStartIdx[iDimX] +
     258           5 :                                     (count[iDimX] - 1) * arrayStep[iDimX]) *
     259           5 :                                        m_dfResX;
     260           5 :     const double dfMinX = std::min(dfX1, dfX2) - m_dfResX / 2;
     261           5 :     const double dfMaxX = std::max(dfX1, dfX2) + m_dfResX / 2;
     262             : 
     263           5 :     const double dfY1 = m_dfMinY + arrayStartIdx[iDimY] * m_dfResY;
     264           5 :     const double dfY2 = m_dfMinY + (arrayStartIdx[iDimY] +
     265           5 :                                     (count[iDimY] - 1) * arrayStep[iDimY]) *
     266           5 :                                        m_dfResY;
     267           5 :     const double dfMinY = std::min(dfY1, dfY2) - m_dfResY / 2;
     268           5 :     const double dfMaxY = std::max(dfY1, dfY2) + m_dfResY / 2;
     269             : 
     270             :     // Extract relevant variable values
     271           5 :     auto poLyr = m_poVectorDS->GetLayer(0);
     272           9 :     if (!(arrayStartIdx[iDimX] == 0 &&
     273           4 :           arrayStartIdx[iDimX] + (count[iDimX] - 1) * arrayStep[iDimX] ==
     274           4 :               m_apoDims[iDimX]->GetSize() - 1 &&
     275           4 :           arrayStartIdx[iDimY] == 0 &&
     276           4 :           arrayStartIdx[iDimY] + (count[iDimY] - 1) * arrayStep[iDimY] ==
     277           4 :               m_apoDims[iDimY]->GetSize() - 1))
     278             :     {
     279           1 :         poLyr->SetSpatialFilterRect(dfMinX - m_dfRadius, dfMinY - m_dfRadius,
     280           1 :                                     dfMaxX + m_dfRadius, dfMaxY + m_dfRadius);
     281             :     }
     282             :     else
     283             :     {
     284           4 :         poLyr->SetSpatialFilter(nullptr);
     285             :     }
     286          10 :     std::vector<double> adfX;
     287          10 :     std::vector<double> adfY;
     288          10 :     std::vector<double> adfZ;
     289             :     try
     290             :     {
     291          35 :         for (auto &&poFeat : poLyr)
     292             :         {
     293          30 :             const auto poGeom = poFeat->GetGeometryRef();
     294          30 :             CPLAssert(poGeom);
     295          30 :             CPLAssert(poGeom->getGeometryType() == wkbPoint);
     296          30 :             adfX.push_back(poGeom->toPoint()->getX());
     297          30 :             adfY.push_back(poGeom->toPoint()->getY());
     298          30 :             const auto nIdxInDataset = poFeat->GetFieldAsInteger64(0);
     299          30 :             assert(nIdxInDataset >= 0 &&
     300             :                    static_cast<size_t>(nIdxInDataset) < m_adfZ.size());
     301          30 :             adfZ.push_back(m_adfZ[static_cast<size_t>(nIdxInDataset)]);
     302             :         }
     303             :     }
     304           0 :     catch (const std::bad_alloc &e)
     305             :     {
     306           0 :         CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
     307           0 :         return false;
     308             :     }
     309             : 
     310           5 :     const size_t nXSize =
     311           5 :         static_cast<size_t>((count[iDimX] - 1) * arrayStep[iDimX] + 1);
     312           5 :     const size_t nYSize =
     313           5 :         static_cast<size_t>((count[iDimY] - 1) * arrayStep[iDimY] + 1);
     314           5 :     if (nXSize > std::numeric_limits<GUInt32>::max() / nYSize)
     315             :     {
     316           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     317             :                  "Too many points queried at once");
     318           0 :         return false;
     319             :     }
     320          10 :     std::vector<double> adfRes;
     321             :     try
     322             :     {
     323           5 :         adfRes.resize(nXSize * nYSize, m_dfNoDataValue);
     324             :     }
     325           0 :     catch (const std::bad_alloc &e)
     326             :     {
     327           0 :         CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
     328           0 :         return false;
     329             :     }
     330             : #if 0
     331             :     CPLDebug("GDAL",
     332             :              "dfMinX=%f, dfMaxX=%f, dfMinY=%f, dfMaxY=%f",
     333             :              dfMinX, dfMaxX, dfMinY, dfMaxY);
     334             : #endif
     335             : 
     336             :     // Finally do the gridded interpolation
     337          10 :     if (!adfX.empty() &&
     338           5 :         GDALGridCreate(
     339           5 :             m_eAlg, m_poGridOptions.get(), static_cast<GUInt32>(adfX.size()),
     340           5 :             adfX.data(), adfY.data(), adfZ.data(), dfMinX, dfMaxX, dfMinY,
     341             :             dfMaxY, static_cast<GUInt32>(nXSize), static_cast<GUInt32>(nYSize),
     342           5 :             GDT_Float64, adfRes.data(), nullptr, nullptr) != CE_None)
     343             :     {
     344           0 :         return false;
     345             :     }
     346             : 
     347             :     // Copy interpolated data into destination buffer
     348           5 :     GByte *pabyDestBuffer = static_cast<GByte *>(pDstBuffer);
     349           5 :     const auto eBufferDT = bufferDataType.GetNumericDataType();
     350           5 :     const auto nBufferDTSize = GDALGetDataTypeSizeBytes(eBufferDT);
     351          15 :     for (size_t iY = 0; iY < count[iDimY]; ++iY)
     352             :     {
     353          20 :         GDALCopyWords64(
     354          10 :             adfRes.data() + iY * arrayStep[iDimY] * nXSize, GDT_Float64,
     355          10 :             static_cast<int>(sizeof(double) * arrayStep[iDimX]),
     356          10 :             pabyDestBuffer + iY * bufferStride[iDimY] * nBufferDTSize,
     357          10 :             eBufferDT, static_cast<int>(bufferStride[iDimX] * nBufferDTSize),
     358          10 :             static_cast<GPtrDiff_t>(count[iDimX]));
     359             :     }
     360             : 
     361           5 :     return true;
     362             : }
     363             : 
     364             : /************************************************************************/
     365             : /*                            GetGridded()                              */
     366             : /************************************************************************/
     367             : 
     368             : /** Return a gridded array from scattered point data, that is from an array
     369             :  * whose last dimension is the indexing variable of X and Y arrays.
     370             :  *
     371             :  * The gridding is done in 2D, using GDALGridCreate(), on-the-fly at Read()
     372             :  * time, taking into account the spatial extent of the request to limit the
     373             :  * gridding. The results got on the whole extent or a subset of it might not be
     374             :  * strictly identical depending on the gridding algorithm and its radius.
     375             :  * Setting a radius in osGridOptions is recommended to improve performance.
     376             :  * For arrays which have more dimensions than the dimension of the indexing
     377             :  * variable of the X and Y arrays, Read() must be called on slices of the extra
     378             :  * dimensions (ie count[i] must be set to 1, except for the X and Y dimensions
     379             :  * of the array returned by this method).
     380             :  *
     381             :  * This is the same as the C function GDALMDArrayGetGridded().
     382             :  *
     383             :  * @param osGridOptions Gridding algorithm and options.
     384             :  * e.g. "invdist:nodata=nan:radius1=1:radius2=1:max_points=5".
     385             :  * See documentation of the <a href="/programs/gdal_grid.html">gdal_grid</a>
     386             :  * utility for all options.
     387             :  * @param poXArrayIn Single-dimension array containing X values, and whose
     388             :  * dimension is the last one of this array. If set to nullptr, the "coordinates"
     389             :  * attribute must exist on this array, and the X variable will be the (N-1)th one
     390             :  * mentioned in it, unless there is a "x" or "lon" variable in "coordinates".
     391             :  * @param poYArrayIn Single-dimension array containing Y values, and whose
     392             :  * dimension is the last one of this array. If set to nullptr, the "coordinates"
     393             :  * attribute must exist on this array, and the Y variable will be the (N-2)th one
     394             :  * mentioned in it,  unless there is a "y" or "lat" variable in "coordinates".
     395             :  * @param papszOptions NULL terminated list of options, or nullptr. Supported
     396             :  * options are:
     397             :  * <ul>
     398             :  * <li>RESOLUTION=val: Spatial resolution of the returned array. If not set,
     399             :  * will be guessed from the typical spacing of (X,Y) points.</li>
     400             :  * </ul>
     401             :  *
     402             :  * @return gridded array, or nullptr in case of error.
     403             :  *
     404             :  * @since GDAL 3.7
     405             :  */
     406             : std::shared_ptr<GDALMDArray>
     407          23 : GDALMDArray::GetGridded(const std::string &osGridOptions,
     408             :                         const std::shared_ptr<GDALMDArray> &poXArrayIn,
     409             :                         const std::shared_ptr<GDALMDArray> &poYArrayIn,
     410             :                         CSLConstList papszOptions) const
     411             : {
     412          46 :     auto self = std::dynamic_pointer_cast<GDALMDArray>(m_pSelf.lock());
     413          23 :     if (!self)
     414             :     {
     415           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     416             :                  "Driver implementation issue: m_pSelf not set !");
     417           0 :         return nullptr;
     418             :     }
     419             : 
     420             :     GDALGridAlgorithm eAlg;
     421          23 :     void *pOptions = nullptr;
     422          23 :     if (GDALGridParseAlgorithmAndOptions(osGridOptions.c_str(), &eAlg,
     423          23 :                                          &pOptions) != CE_None)
     424             :     {
     425           1 :         return nullptr;
     426             :     }
     427             : 
     428          44 :     std::unique_ptr<void, VSIFreeReleaser> poGridOptions(pOptions);
     429             : 
     430          22 :     if (GetDataType().GetClass() != GEDTC_NUMERIC)
     431             :     {
     432           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     433             :                  "GetDataType().GetClass() != GEDTC_NUMERIC");
     434           1 :         return nullptr;
     435             :     }
     436             : 
     437          21 :     if (GetDimensionCount() == 0)
     438             :     {
     439           1 :         CPLError(CE_Failure, CPLE_NotSupported, "GetDimensionCount() == 0");
     440           1 :         return nullptr;
     441             :     }
     442             : 
     443          20 :     if (poXArrayIn && !poYArrayIn)
     444             :     {
     445           1 :         CPLError(
     446             :             CE_Failure, CPLE_AppDefined,
     447             :             "As poXArrayIn is specified, poYArrayIn must also be specified");
     448           1 :         return nullptr;
     449             :     }
     450          19 :     else if (!poXArrayIn && poYArrayIn)
     451             :     {
     452           1 :         CPLError(
     453             :             CE_Failure, CPLE_AppDefined,
     454             :             "As poYArrayIn is specified, poXArrayIn must also be specified");
     455           1 :         return nullptr;
     456             :     }
     457          36 :     std::shared_ptr<GDALMDArray> poXArray = poXArrayIn;
     458          36 :     std::shared_ptr<GDALMDArray> poYArray = poYArrayIn;
     459             : 
     460          18 :     if (!poXArray)
     461             :     {
     462           8 :         const auto aoCoordVariables = GetCoordinateVariables();
     463           8 :         if (aoCoordVariables.size() < 2)
     464             :         {
     465           1 :             CPLError(CE_Failure, CPLE_NotSupported,
     466             :                      "aoCoordVariables.size() < 2");
     467           1 :             return nullptr;
     468             :         }
     469             : 
     470           7 :         if (aoCoordVariables.size() != GetDimensionCount() + 1)
     471             :         {
     472           2 :             CPLError(CE_Failure, CPLE_NotSupported,
     473             :                      "aoCoordVariables.size() != GetDimensionCount() + 1");
     474           2 :             return nullptr;
     475             :         }
     476             : 
     477             :         // Default choice for X and Y arrays
     478           5 :         poYArray = aoCoordVariables[aoCoordVariables.size() - 2];
     479           5 :         poXArray = aoCoordVariables[aoCoordVariables.size() - 1];
     480             : 
     481             :         // Detect X and Y array from coordinate variables
     482          20 :         for (const auto &poVar : aoCoordVariables)
     483             :         {
     484          15 :             const auto &osVarName = poVar->GetName();
     485             : #ifdef disabled
     486             :             if (poVar->GetDimensionCount() != 1)
     487             :             {
     488             :                 CPLError(CE_Failure, CPLE_NotSupported,
     489             :                          "Coordinate variable %s is not 1-dimensional",
     490             :                          osVarName.c_str());
     491             :                 return nullptr;
     492             :             }
     493             :             const auto &osVarDimName = poVar->GetDimensions()[0]->GetFullName();
     494             :             bool bFound = false;
     495             :             for (const auto &poDim : GetDimensions())
     496             :             {
     497             :                 if (osVarDimName == poDim->GetFullName())
     498             :                 {
     499             :                     bFound = true;
     500             :                     break;
     501             :                 }
     502             :             }
     503             :             if (!bFound)
     504             :             {
     505             :                 CPLError(CE_Failure, CPLE_NotSupported,
     506             :                          "Dimension %s of coordinate variable %s is not a "
     507             :                          "dimension of this array",
     508             :                          osVarDimName.c_str(), osVarName.c_str());
     509             :                 return nullptr;
     510             :             }
     511             : #endif
     512          15 :             if (osVarName == "x" || osVarName == "lon")
     513             :             {
     514           0 :                 poXArray = poVar;
     515             :             }
     516          15 :             else if (osVarName == "y" || osVarName == "lat")
     517             :             {
     518           0 :                 poYArray = poVar;
     519             :             }
     520             :         }
     521             :     }
     522             : 
     523          15 :     if (poYArray->GetDimensionCount() != 1)
     524             :     {
     525           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     526             :                  "aoCoordVariables[aoCoordVariables.size() - "
     527             :                  "2]->GetDimensionCount() != 1");
     528           1 :         return nullptr;
     529             :     }
     530          14 :     if (poYArray->GetDataType().GetClass() != GEDTC_NUMERIC)
     531             :     {
     532           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     533             :                  "poYArray->GetDataType().GetClass() != GEDTC_NUMERIC");
     534           1 :         return nullptr;
     535             :     }
     536             : 
     537          13 :     if (poXArray->GetDimensionCount() != 1)
     538             :     {
     539           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     540             :                  "aoCoordVariables[aoCoordVariables.size() - "
     541             :                  "1]->GetDimensionCount() != 1");
     542           1 :         return nullptr;
     543             :     }
     544          12 :     if (poXArray->GetDataType().GetClass() != GEDTC_NUMERIC)
     545             :     {
     546           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     547             :                  "poXArray->GetDataType().GetClass() != GEDTC_NUMERIC");
     548           1 :         return nullptr;
     549             :     }
     550             : 
     551          11 :     if (poYArray->GetDimensions()[0]->GetFullName() !=
     552          11 :         poXArray->GetDimensions()[0]->GetFullName())
     553             :     {
     554           3 :         CPLError(CE_Failure, CPLE_NotSupported,
     555             :                  "poYArray->GetDimensions()[0]->GetFullName() != "
     556             :                  "poXArray->GetDimensions()[0]->GetFullName()");
     557           3 :         return nullptr;
     558             :     }
     559             : 
     560           8 :     if (poXArray->GetDimensions()[0]->GetFullName() !=
     561           8 :         GetDimensions().back()->GetFullName())
     562             :     {
     563           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     564             :                  "poYArray->GetDimensions()[0]->GetFullName() != "
     565             :                  "GetDimensions().back()->GetFullName()");
     566           0 :         return nullptr;
     567             :     }
     568             : 
     569           8 :     if (poXArray->GetTotalElementsCount() <= 2)
     570             :     {
     571           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     572             :                  "poXArray->GetTotalElementsCount() <= 2");
     573           1 :         return nullptr;
     574             :     }
     575             : 
     576           7 :     if (poXArray->GetTotalElementsCount() >
     577           7 :         std::numeric_limits<size_t>::max() / 2)
     578             :     {
     579           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     580             :                  "poXArray->GetTotalElementsCount() > "
     581             :                  "std::numeric_limits<size_t>::max() / 2");
     582           0 :         return nullptr;
     583             :     }
     584             : 
     585           8 :     if (poXArray->GetTotalElementsCount() > 10 * 1024 * 1024 &&
     586           1 :         !CPLTestBool(CSLFetchNameValueDef(
     587             :             papszOptions, "ACCEPT_BIG_SPATIAL_INDEXING_VARIABLE", "NO")))
     588             :     {
     589           1 :         CPLError(CE_Failure, CPLE_AppDefined,
     590             :                  "The spatial indexing variable has " CPL_FRMT_GIB " elements. "
     591             :                  "Set the ACCEPT_BIG_SPATIAL_INDEXING_VARIABLE=YES option of "
     592             :                  "GetGridded() to mean you want to continue and are aware of "
     593             :                  "big RAM and CPU time requirements",
     594           1 :                  static_cast<GIntBig>(poXArray->GetTotalElementsCount()));
     595           1 :         return nullptr;
     596             :     }
     597             : 
     598          12 :     std::vector<double> adfXVals;
     599          12 :     std::vector<double> adfYVals;
     600             :     try
     601             :     {
     602           6 :         adfXVals.resize(static_cast<size_t>(poXArray->GetTotalElementsCount()));
     603           6 :         adfYVals.resize(static_cast<size_t>(poXArray->GetTotalElementsCount()));
     604             :     }
     605           0 :     catch (const std::bad_alloc &e)
     606             :     {
     607           0 :         CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
     608           0 :         return nullptr;
     609             :     }
     610             : 
     611             :     // Ingest X and Y arrays
     612           6 :     const GUInt64 arrayStartIdx[] = {0};
     613           6 :     const size_t count[] = {adfXVals.size()};
     614           6 :     const GInt64 arrayStep[] = {1};
     615           6 :     const GPtrDiff_t bufferStride[] = {1};
     616          12 :     if (!poXArray->Read(arrayStartIdx, count, arrayStep, bufferStride,
     617          12 :                         GDALExtendedDataType::Create(GDT_Float64),
     618           6 :                         adfXVals.data()))
     619             :     {
     620           0 :         CPLError(CE_Failure, CPLE_AppDefined, "poXArray->Read() failed");
     621           0 :         return nullptr;
     622             :     }
     623          12 :     if (!poYArray->Read(arrayStartIdx, count, arrayStep, bufferStride,
     624          12 :                         GDALExtendedDataType::Create(GDT_Float64),
     625           6 :                         adfYVals.data()))
     626             :     {
     627           0 :         CPLError(CE_Failure, CPLE_AppDefined, "poYArray->Read() failed");
     628           0 :         return nullptr;
     629             :     }
     630             : 
     631           6 :     const char *pszExt = "fgb";
     632           6 :     GDALDriver *poDrv = GetGDALDriverManager()->GetDriverByName("FlatGeoBuf");
     633           6 :     if (!poDrv)
     634             :     {
     635           0 :         pszExt = "gpkg";
     636           0 :         poDrv = GetGDALDriverManager()->GetDriverByName("GPKG");
     637           0 :         if (!poDrv)
     638             :         {
     639           0 :             pszExt = "mem";
     640           0 :             poDrv = GetGDALDriverManager()->GetDriverByName("Memory");
     641           0 :             if (!poDrv)
     642             :             {
     643           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     644             :                          "Cannot get driver FlatGeoBuf, GPKG or Memory");
     645           0 :                 return nullptr;
     646             :             }
     647             :         }
     648             :     }
     649             : 
     650             :     // Create a in-memory vector layer with (X,Y) points
     651             :     const std::string osTmpFilename(VSIMemGenerateHiddenFilename(
     652          18 :         std::string("tmp.").append(pszExt).c_str()));
     653             :     auto poDS = std::unique_ptr<GDALDataset>(
     654          12 :         poDrv->Create(osTmpFilename.c_str(), 0, 0, 0, GDT_Unknown, nullptr));
     655           6 :     if (!poDS)
     656           0 :         return nullptr;
     657           6 :     auto poLyr = poDS->CreateLayer("layer", nullptr, wkbPoint);
     658           6 :     if (!poLyr)
     659           0 :         return nullptr;
     660          12 :     OGRFieldDefn oFieldDefn("IDX", OFTInteger64);
     661           6 :     poLyr->CreateField(&oFieldDefn);
     662           6 :     if (poLyr->StartTransaction() != OGRERR_NONE)
     663           0 :         return nullptr;
     664          12 :     OGRFeature oFeat(poLyr->GetLayerDefn());
     665          42 :     for (size_t i = 0; i < adfXVals.size(); ++i)
     666             :     {
     667          36 :         auto poPoint = new OGRPoint(adfXVals[i], adfYVals[i]);
     668          36 :         oFeat.SetFID(OGRNullFID);
     669          36 :         oFeat.SetGeometryDirectly(poPoint);
     670          36 :         oFeat.SetField(0, static_cast<GIntBig>(i));
     671          36 :         if (poLyr->CreateFeature(&oFeat) != OGRERR_NONE)
     672           0 :             return nullptr;
     673             :     }
     674           6 :     if (poLyr->CommitTransaction() != OGRERR_NONE)
     675           0 :         return nullptr;
     676           6 :     OGREnvelope sEnvelope;
     677           6 :     CPL_IGNORE_RET_VAL(poLyr->GetExtent(&sEnvelope));
     678           6 :     if (!EQUAL(pszExt, "mem"))
     679             :     {
     680           6 :         if (poDS->Close() != OGRERR_NONE)
     681           0 :             return nullptr;
     682           6 :         poDS.reset(GDALDataset::Open(osTmpFilename.c_str(), GDAL_OF_VECTOR));
     683           6 :         if (!poDS)
     684           0 :             return nullptr;
     685           6 :         poDS->MarkSuppressOnClose();
     686             :     }
     687             : 
     688             :     /* Set of constraints:
     689             :     nX * nY = nCount
     690             :     nX * res = MaxX - MinX
     691             :     nY * res = MaxY - MinY
     692             :     */
     693             : 
     694             :     double dfRes;
     695           6 :     const char *pszRes = CSLFetchNameValue(papszOptions, "RESOLUTION");
     696           6 :     if (pszRes)
     697             :     {
     698           3 :         dfRes = CPLAtofM(pszRes);
     699             :     }
     700             :     else
     701             :     {
     702           3 :         const double dfTotalArea = (sEnvelope.MaxY - sEnvelope.MinY) *
     703           3 :                                    (sEnvelope.MaxX - sEnvelope.MinX);
     704           3 :         dfRes = sqrt(dfTotalArea / static_cast<double>(adfXVals.size()));
     705             :         // CPLDebug("GDAL", "dfRes = %f", dfRes);
     706             : 
     707             :         // Take 10 "random" points in the set, and find the minimum distance from
     708             :         // each to the closest one. And take the geometric average of those minimum
     709             :         // distances as the resolution.
     710           3 :         const size_t nNumSamplePoints = std::min<size_t>(10, adfXVals.size());
     711             : 
     712           3 :         poLyr = poDS->GetLayer(0);
     713           3 :         double dfSumDist2Min = 0;
     714           3 :         int nCountDistMin = 0;
     715          21 :         for (size_t i = 0; i < nNumSamplePoints; ++i)
     716             :         {
     717          18 :             const auto nIdx = i * adfXVals.size() / nNumSamplePoints;
     718          90 :             poLyr->SetSpatialFilterRect(
     719          18 :                 adfXVals[nIdx] - 2 * dfRes, adfYVals[nIdx] - 2 * dfRes,
     720          18 :                 adfXVals[nIdx] + 2 * dfRes, adfYVals[nIdx] + 2 * dfRes);
     721          18 :             double dfDist2Min = std::numeric_limits<double>::max();
     722         102 :             for (auto &&poFeat : poLyr)
     723             :             {
     724          84 :                 const auto poGeom = poFeat->GetGeometryRef();
     725          84 :                 CPLAssert(poGeom);
     726          84 :                 CPLAssert(poGeom->getGeometryType() == wkbPoint);
     727          84 :                 double dfDX = poGeom->toPoint()->getX() - adfXVals[nIdx];
     728          84 :                 double dfDY = poGeom->toPoint()->getY() - adfYVals[nIdx];
     729          84 :                 double dfDist2 = dfDX * dfDX + dfDY * dfDY;
     730          84 :                 if (dfDist2 > 0 && dfDist2 < dfDist2Min)
     731          24 :                     dfDist2Min = dfDist2;
     732             :             }
     733          18 :             if (dfDist2Min < std::numeric_limits<double>::max())
     734             :             {
     735          18 :                 dfSumDist2Min += dfDist2Min;
     736          18 :                 nCountDistMin++;
     737             :             }
     738             :         }
     739           3 :         poLyr->SetSpatialFilter(nullptr);
     740           3 :         if (nCountDistMin > 0)
     741             :         {
     742           3 :             const double dfNewRes = sqrt(dfSumDist2Min / nCountDistMin);
     743             :             // CPLDebug("GDAL", "dfRes = %f, dfNewRes = %f", dfRes, dfNewRes);
     744           3 :             dfRes = dfNewRes;
     745             :         }
     746             :     }
     747             : 
     748           6 :     if (!(dfRes > 0))
     749             :     {
     750           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Invalid RESOLUTION");
     751           1 :         return nullptr;
     752             :     }
     753             : 
     754           5 :     constexpr double EPS = 1e-8;
     755           5 :     const double dfXSize =
     756           5 :         1 + std::floor((sEnvelope.MaxX - sEnvelope.MinX) / dfRes + EPS);
     757           5 :     if (dfXSize > std::numeric_limits<int>::max())
     758             :     {
     759           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Integer overflow with dfXSize");
     760           1 :         return nullptr;
     761             :     }
     762           4 :     const int nXSize = std::max(2, static_cast<int>(dfXSize));
     763             : 
     764           4 :     const double dfYSize =
     765           4 :         1 + std::floor((sEnvelope.MaxY - sEnvelope.MinY) / dfRes + EPS);
     766           4 :     if (dfYSize > std::numeric_limits<int>::max())
     767             :     {
     768           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Integer overflow with dfYSize");
     769           0 :         return nullptr;
     770             :     }
     771           4 :     const int nYSize = std::max(2, static_cast<int>(dfYSize));
     772             : 
     773           4 :     const double dfResX = (sEnvelope.MaxX - sEnvelope.MinX) / (nXSize - 1);
     774           4 :     const double dfResY = (sEnvelope.MaxY - sEnvelope.MinY) / (nYSize - 1);
     775             :     // CPLDebug("GDAL", "nXSize = %d, nYSize = %d", nXSize, nYSize);
     776             : 
     777           8 :     std::vector<std::shared_ptr<GDALDimension>> apoNewDims;
     778           4 :     const auto &apoSelfDims = GetDimensions();
     779           8 :     for (size_t i = 0; i < GetDimensionCount() - 1; ++i)
     780           4 :         apoNewDims.emplace_back(apoSelfDims[i]);
     781             : 
     782             :     auto poDimY = std::make_shared<GDALDimensionWeakIndexingVar>(
     783           8 :         std::string(), "dimY", GDAL_DIM_TYPE_HORIZONTAL_Y, "NORTH", nYSize);
     784             :     auto varY = GDALMDArrayRegularlySpaced::Create(
     785          12 :         std::string(), poDimY->GetName(), poDimY, sEnvelope.MinY, dfResY, 0);
     786           4 :     poDimY->SetIndexingVariable(varY);
     787             : 
     788             :     auto poDimX = std::make_shared<GDALDimensionWeakIndexingVar>(
     789           8 :         std::string(), "dimX", GDAL_DIM_TYPE_HORIZONTAL_X, "EAST", nXSize);
     790             :     auto varX = GDALMDArrayRegularlySpaced::Create(
     791          12 :         std::string(), poDimX->GetName(), poDimX, sEnvelope.MinX, dfResX, 0);
     792           4 :     poDimX->SetIndexingVariable(varX);
     793             : 
     794           4 :     apoNewDims.emplace_back(poDimY);
     795           4 :     apoNewDims.emplace_back(poDimX);
     796             : 
     797             :     const CPLStringList aosTokens(
     798           4 :         CSLTokenizeString2(osGridOptions.c_str(), ":", FALSE));
     799             : 
     800             :     // Extract nodata value from gridding options
     801           4 :     const char *pszNoDataValue = aosTokens.FetchNameValue("nodata");
     802           4 :     double dfNoDataValue = 0;
     803           4 :     if (pszNoDataValue != nullptr)
     804             :     {
     805           1 :         dfNoDataValue = CPLAtofM(pszNoDataValue);
     806             :     }
     807             : 
     808             :     // Extract radius from gridding options
     809           4 :     double dfRadius = 5 * std::max(dfResX, dfResY);
     810           4 :     const char *pszRadius = aosTokens.FetchNameValue("radius");
     811           4 :     if (pszRadius)
     812             :     {
     813           0 :         dfRadius = CPLAtofM(pszRadius);
     814             :     }
     815             :     else
     816             :     {
     817           4 :         const char *pszRadius1 = aosTokens.FetchNameValue("radius1");
     818           4 :         if (pszRadius1)
     819             :         {
     820           2 :             dfRadius = CPLAtofM(pszRadius1);
     821           2 :             const char *pszRadius2 = aosTokens.FetchNameValue("radius2");
     822           2 :             if (pszRadius2)
     823             :             {
     824           2 :                 dfRadius = std::max(dfRadius, CPLAtofM(pszRadius2));
     825             :             }
     826             :         }
     827             :     }
     828             : 
     829           8 :     return GDALMDArrayGridded::Create(
     830           4 :         self, apoNewDims, varX, varY, std::move(poDS), eAlg,
     831           4 :         std::move(poGridOptions), dfNoDataValue, sEnvelope.MinX, dfResX,
     832           4 :         sEnvelope.MinY, dfResY, dfRadius);
     833             : }