Line data    Source code

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