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