Line data    Source code

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