Line data Source code
1 : /******************************************************************************
2 : *
3 : * Name: gdalmultidim_array_resampled.cpp
4 : * Project: GDAL Core
5 : * Purpose: GDALMDArray::GetResampled() implementation
6 : * Author: Even Rouault <even.rouault at spatialys.com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2019, Even Rouault <even.rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "gdal_multidim.h"
15 : #include "gdal_pam.h"
16 : #include "gdal_pam_multidim.h"
17 : #include "gdal_rasterband.h"
18 : #include "gdal_utils.h"
19 :
20 : #include <algorithm>
21 :
22 : //! @cond Doxygen_Suppress
23 :
24 : /************************************************************************/
25 : /* GDALMDArrayResampled */
26 : /************************************************************************/
27 :
28 : class GDALMDArrayResampledDataset;
29 :
30 : class GDALMDArrayResampledDatasetRasterBand final : public GDALRasterBand
31 : {
32 : protected:
33 : CPLErr IReadBlock(int, int, void *) override;
34 : CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
35 : int nYSize, void *pData, int nBufXSize, int nBufYSize,
36 : GDALDataType eBufType, GSpacing nPixelSpaceBuf,
37 : GSpacing nLineSpaceBuf,
38 : GDALRasterIOExtraArg *psExtraArg) override;
39 :
40 : public:
41 : explicit GDALMDArrayResampledDatasetRasterBand(
42 : GDALMDArrayResampledDataset *poDSIn);
43 :
44 : double GetNoDataValue(int *pbHasNoData) override;
45 : };
46 :
47 : class GDALMDArrayResampledDataset final : public GDALPamDataset
48 : {
49 : friend class GDALMDArrayResampled;
50 : friend class GDALMDArrayResampledDatasetRasterBand;
51 :
52 : std::shared_ptr<GDALMDArray> m_poArray;
53 : const size_t m_iXDim;
54 : const size_t m_iYDim;
55 : GDALGeoTransform m_gt{};
56 : bool m_bHasGT = false;
57 : mutable std::shared_ptr<OGRSpatialReference> m_poSRS{};
58 :
59 : std::vector<GUInt64> m_anOffset{};
60 : std::vector<size_t> m_anCount{};
61 : std::vector<GPtrDiff_t> m_anStride{};
62 :
63 : std::string m_osFilenameLong{};
64 : std::string m_osFilenameLat{};
65 :
66 : public:
67 24 : GDALMDArrayResampledDataset(const std::shared_ptr<GDALMDArray> &array,
68 : size_t iXDim, size_t iYDim)
69 24 : : m_poArray(array), m_iXDim(iXDim), m_iYDim(iYDim),
70 24 : m_anOffset(m_poArray->GetDimensionCount(), 0),
71 24 : m_anCount(m_poArray->GetDimensionCount(), 1),
72 72 : m_anStride(m_poArray->GetDimensionCount(), 0)
73 : {
74 24 : const auto &dims(m_poArray->GetDimensions());
75 :
76 24 : nRasterYSize = static_cast<int>(
77 24 : std::min(static_cast<GUInt64>(INT_MAX), dims[iYDim]->GetSize()));
78 24 : nRasterXSize = static_cast<int>(
79 24 : std::min(static_cast<GUInt64>(INT_MAX), dims[iXDim]->GetSize()));
80 :
81 24 : m_bHasGT = m_poArray->GuessGeoTransform(m_iXDim, m_iYDim, false, m_gt);
82 :
83 24 : SetBand(1, new GDALMDArrayResampledDatasetRasterBand(this));
84 24 : }
85 :
86 : ~GDALMDArrayResampledDataset() override;
87 :
88 43 : CPLErr GetGeoTransform(GDALGeoTransform >) const override
89 : {
90 43 : gt = m_gt;
91 43 : return m_bHasGT ? CE_None : CE_Failure;
92 : }
93 :
94 105 : const OGRSpatialReference *GetSpatialRef() const override
95 : {
96 105 : m_poSRS = m_poArray->GetSpatialRef();
97 105 : if (m_poSRS)
98 : {
99 79 : m_poSRS.reset(m_poSRS->Clone());
100 158 : auto axisMapping = m_poSRS->GetDataAxisToSRSAxisMapping();
101 237 : for (auto &m : axisMapping)
102 : {
103 158 : if (m == static_cast<int>(m_iXDim) + 1)
104 79 : m = 1;
105 79 : else if (m == static_cast<int>(m_iYDim) + 1)
106 79 : m = 2;
107 : }
108 79 : m_poSRS->SetDataAxisToSRSAxisMapping(axisMapping);
109 : }
110 105 : return m_poSRS.get();
111 : }
112 :
113 5 : void SetGeolocationArray(const std::string &osFilenameLong,
114 : const std::string &osFilenameLat)
115 : {
116 5 : m_osFilenameLong = osFilenameLong;
117 5 : m_osFilenameLat = osFilenameLat;
118 10 : CPLStringList aosGeoLoc;
119 5 : aosGeoLoc.SetNameValue("LINE_OFFSET", "0");
120 5 : aosGeoLoc.SetNameValue("LINE_STEP", "1");
121 5 : aosGeoLoc.SetNameValue("PIXEL_OFFSET", "0");
122 5 : aosGeoLoc.SetNameValue("PIXEL_STEP", "1");
123 5 : aosGeoLoc.SetNameValue("SRS", SRS_WKT_WGS84_LAT_LONG); // FIXME?
124 5 : aosGeoLoc.SetNameValue("X_BAND", "1");
125 5 : aosGeoLoc.SetNameValue("X_DATASET", m_osFilenameLong.c_str());
126 5 : aosGeoLoc.SetNameValue("Y_BAND", "1");
127 5 : aosGeoLoc.SetNameValue("Y_DATASET", m_osFilenameLat.c_str());
128 5 : aosGeoLoc.SetNameValue("GEOREFERENCING_CONVENTION", "PIXEL_CENTER");
129 5 : SetMetadata(aosGeoLoc.List(), "GEOLOCATION");
130 5 : }
131 : };
132 :
133 48 : GDALMDArrayResampledDataset::~GDALMDArrayResampledDataset()
134 : {
135 24 : if (!m_osFilenameLong.empty())
136 5 : VSIUnlink(m_osFilenameLong.c_str());
137 24 : if (!m_osFilenameLat.empty())
138 5 : VSIUnlink(m_osFilenameLat.c_str());
139 48 : }
140 :
141 : /************************************************************************/
142 : /* GDALMDArrayResampledDatasetRasterBand() */
143 : /************************************************************************/
144 :
145 24 : GDALMDArrayResampledDatasetRasterBand::GDALMDArrayResampledDatasetRasterBand(
146 24 : GDALMDArrayResampledDataset *poDSIn)
147 : {
148 24 : const auto &poArray(poDSIn->m_poArray);
149 24 : const auto blockSize(poArray->GetBlockSize());
150 24 : nBlockYSize = (blockSize[poDSIn->m_iYDim])
151 24 : ? static_cast<int>(std::min(static_cast<GUInt64>(INT_MAX),
152 13 : blockSize[poDSIn->m_iYDim]))
153 24 : : 1;
154 24 : nBlockXSize = blockSize[poDSIn->m_iXDim]
155 13 : ? static_cast<int>(std::min(static_cast<GUInt64>(INT_MAX),
156 13 : blockSize[poDSIn->m_iXDim]))
157 24 : : poDSIn->GetRasterXSize();
158 24 : eDataType = poArray->GetDataType().GetNumericDataType();
159 24 : eAccess = poDSIn->eAccess;
160 24 : }
161 :
162 : /************************************************************************/
163 : /* GetNoDataValue() */
164 : /************************************************************************/
165 :
166 54 : double GDALMDArrayResampledDatasetRasterBand::GetNoDataValue(int *pbHasNoData)
167 : {
168 54 : auto l_poDS(cpl::down_cast<GDALMDArrayResampledDataset *>(poDS));
169 54 : const auto &poArray(l_poDS->m_poArray);
170 54 : bool bHasNodata = false;
171 54 : double dfRes = poArray->GetNoDataValueAsDouble(&bHasNodata);
172 54 : if (pbHasNoData)
173 48 : *pbHasNoData = bHasNodata;
174 54 : return dfRes;
175 : }
176 :
177 : /************************************************************************/
178 : /* IReadBlock() */
179 : /************************************************************************/
180 :
181 0 : CPLErr GDALMDArrayResampledDatasetRasterBand::IReadBlock(int nBlockXOff,
182 : int nBlockYOff,
183 : void *pImage)
184 : {
185 0 : const int nDTSize(GDALGetDataTypeSizeBytes(eDataType));
186 0 : const int nXOff = nBlockXOff * nBlockXSize;
187 0 : const int nYOff = nBlockYOff * nBlockYSize;
188 0 : const int nReqXSize = std::min(nRasterXSize - nXOff, nBlockXSize);
189 0 : const int nReqYSize = std::min(nRasterYSize - nYOff, nBlockYSize);
190 : GDALRasterIOExtraArg sExtraArg;
191 0 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
192 0 : return IRasterIO(GF_Read, nXOff, nYOff, nReqXSize, nReqYSize, pImage,
193 : nReqXSize, nReqYSize, eDataType, nDTSize,
194 0 : static_cast<GSpacing>(nDTSize) * nBlockXSize, &sExtraArg);
195 : }
196 :
197 : /************************************************************************/
198 : /* IRasterIO() */
199 : /************************************************************************/
200 :
201 32 : CPLErr GDALMDArrayResampledDatasetRasterBand::IRasterIO(
202 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
203 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
204 : GSpacing nPixelSpaceBuf, GSpacing nLineSpaceBuf,
205 : GDALRasterIOExtraArg *psExtraArg)
206 : {
207 32 : auto l_poDS(cpl::down_cast<GDALMDArrayResampledDataset *>(poDS));
208 32 : const auto &poArray(l_poDS->m_poArray);
209 32 : const int nBufferDTSize(GDALGetDataTypeSizeBytes(eBufType));
210 32 : if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&
211 32 : nBufferDTSize > 0 && (nPixelSpaceBuf % nBufferDTSize) == 0 &&
212 32 : (nLineSpaceBuf % nBufferDTSize) == 0)
213 : {
214 32 : l_poDS->m_anOffset[l_poDS->m_iXDim] = static_cast<GUInt64>(nXOff);
215 32 : l_poDS->m_anCount[l_poDS->m_iXDim] = static_cast<size_t>(nXSize);
216 64 : l_poDS->m_anStride[l_poDS->m_iXDim] =
217 32 : static_cast<GPtrDiff_t>(nPixelSpaceBuf / nBufferDTSize);
218 :
219 32 : l_poDS->m_anOffset[l_poDS->m_iYDim] = static_cast<GUInt64>(nYOff);
220 32 : l_poDS->m_anCount[l_poDS->m_iYDim] = static_cast<size_t>(nYSize);
221 64 : l_poDS->m_anStride[l_poDS->m_iYDim] =
222 32 : static_cast<GPtrDiff_t>(nLineSpaceBuf / nBufferDTSize);
223 :
224 64 : return poArray->Read(l_poDS->m_anOffset.data(),
225 32 : l_poDS->m_anCount.data(), nullptr,
226 32 : l_poDS->m_anStride.data(),
227 64 : GDALExtendedDataType::Create(eBufType), pData)
228 32 : ? CE_None
229 32 : : CE_Failure;
230 : }
231 0 : return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
232 : pData, nBufXSize, nBufYSize, eBufType,
233 0 : nPixelSpaceBuf, nLineSpaceBuf, psExtraArg);
234 : }
235 :
236 : class GDALMDArrayResampled final : public GDALPamMDArray
237 : {
238 : private:
239 : std::shared_ptr<GDALMDArray> m_poParent{};
240 : std::vector<std::shared_ptr<GDALDimension>> m_apoDims;
241 : std::vector<GUInt64> m_anBlockSize;
242 : GDALExtendedDataType m_dt;
243 : std::shared_ptr<OGRSpatialReference> m_poSRS{};
244 : std::shared_ptr<GDALMDArray> m_poVarX{};
245 : std::shared_ptr<GDALMDArray> m_poVarY{};
246 : std::unique_ptr<GDALMDArrayResampledDataset> m_poParentDS{};
247 : std::unique_ptr<GDALDataset> m_poReprojectedDS{};
248 :
249 : protected:
250 21 : GDALMDArrayResampled(
251 : const std::shared_ptr<GDALMDArray> &poParent,
252 : const std::vector<std::shared_ptr<GDALDimension>> &apoDims,
253 : const std::vector<GUInt64> &anBlockSize)
254 42 : : GDALAbstractMDArray(std::string(),
255 42 : "Resampled view of " + poParent->GetFullName()),
256 : GDALPamMDArray(
257 42 : std::string(), "Resampled view of " + poParent->GetFullName(),
258 42 : GDALPamMultiDim::GetPAM(poParent), poParent->GetContext()),
259 21 : m_poParent(std::move(poParent)), m_apoDims(apoDims),
260 105 : m_anBlockSize(anBlockSize), m_dt(m_poParent->GetDataType())
261 : {
262 21 : CPLAssert(apoDims.size() == m_poParent->GetDimensionCount());
263 21 : CPLAssert(anBlockSize.size() == m_poParent->GetDimensionCount());
264 21 : }
265 :
266 : bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
267 : const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
268 : const GDALExtendedDataType &bufferDataType,
269 : void *pDstBuffer) const override;
270 :
271 : public:
272 : static std::shared_ptr<GDALMDArray>
273 : Create(const std::shared_ptr<GDALMDArray> &poParent,
274 : const std::vector<std::shared_ptr<GDALDimension>> &apoNewDims,
275 : GDALRIOResampleAlg resampleAlg,
276 : const OGRSpatialReference *poTargetSRS, CSLConstList papszOptions);
277 :
278 42 : ~GDALMDArrayResampled() override
279 21 : {
280 : // First close the warped VRT
281 21 : m_poReprojectedDS.reset();
282 21 : m_poParentDS.reset();
283 42 : }
284 :
285 11 : bool IsWritable() const override
286 : {
287 11 : return false;
288 : }
289 :
290 63 : const std::string &GetFilename() const override
291 : {
292 63 : return m_poParent->GetFilename();
293 : }
294 :
295 : const std::vector<std::shared_ptr<GDALDimension>> &
296 283 : GetDimensions() const override
297 : {
298 283 : return m_apoDims;
299 : }
300 :
301 109 : const GDALExtendedDataType &GetDataType() const override
302 : {
303 109 : return m_dt;
304 : }
305 :
306 21 : std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override
307 : {
308 21 : return m_poSRS;
309 : }
310 :
311 12 : std::vector<GUInt64> GetBlockSize() const override
312 : {
313 12 : return m_anBlockSize;
314 : }
315 :
316 : std::shared_ptr<GDALAttribute>
317 1 : GetAttribute(const std::string &osName) const override
318 : {
319 1 : return m_poParent->GetAttribute(osName);
320 : }
321 :
322 : std::vector<std::shared_ptr<GDALAttribute>>
323 12 : GetAttributes(CSLConstList papszOptions = nullptr) const override
324 : {
325 12 : return m_poParent->GetAttributes(papszOptions);
326 : }
327 :
328 1 : const std::string &GetUnit() const override
329 : {
330 1 : return m_poParent->GetUnit();
331 : }
332 :
333 1 : const void *GetRawNoDataValue() const override
334 : {
335 1 : return m_poParent->GetRawNoDataValue();
336 : }
337 :
338 1 : double GetOffset(bool *pbHasOffset,
339 : GDALDataType *peStorageType) const override
340 : {
341 1 : return m_poParent->GetOffset(pbHasOffset, peStorageType);
342 : }
343 :
344 1 : double GetScale(bool *pbHasScale,
345 : GDALDataType *peStorageType) const override
346 : {
347 1 : return m_poParent->GetScale(pbHasScale, peStorageType);
348 : }
349 : };
350 :
351 : /************************************************************************/
352 : /* GDALMDArrayResampled::Create() */
353 : /************************************************************************/
354 :
355 29 : std::shared_ptr<GDALMDArray> GDALMDArrayResampled::Create(
356 : const std::shared_ptr<GDALMDArray> &poParent,
357 : const std::vector<std::shared_ptr<GDALDimension>> &apoNewDimsIn,
358 : GDALRIOResampleAlg resampleAlg, const OGRSpatialReference *poTargetSRS,
359 : CSLConstList /* papszOptions */)
360 : {
361 29 : const char *pszResampleAlg = "nearest";
362 29 : bool unsupported = false;
363 29 : switch (resampleAlg)
364 : {
365 16 : case GRIORA_NearestNeighbour:
366 16 : pszResampleAlg = "nearest";
367 16 : break;
368 2 : case GRIORA_Bilinear:
369 2 : pszResampleAlg = "bilinear";
370 2 : break;
371 5 : case GRIORA_Cubic:
372 5 : pszResampleAlg = "cubic";
373 5 : break;
374 1 : case GRIORA_CubicSpline:
375 1 : pszResampleAlg = "cubicspline";
376 1 : break;
377 1 : case GRIORA_Lanczos:
378 1 : pszResampleAlg = "lanczos";
379 1 : break;
380 1 : case GRIORA_Average:
381 1 : pszResampleAlg = "average";
382 1 : break;
383 1 : case GRIORA_Mode:
384 1 : pszResampleAlg = "mode";
385 1 : break;
386 1 : case GRIORA_Gauss:
387 1 : unsupported = true;
388 1 : break;
389 0 : case GRIORA_RESERVED_START:
390 0 : unsupported = true;
391 0 : break;
392 0 : case GRIORA_RESERVED_END:
393 0 : unsupported = true;
394 0 : break;
395 1 : case GRIORA_RMS:
396 1 : pszResampleAlg = "rms";
397 1 : break;
398 : }
399 29 : if (unsupported)
400 : {
401 1 : CPLError(CE_Failure, CPLE_NotSupported,
402 : "Unsupported resample method for GetResampled()");
403 1 : return nullptr;
404 : }
405 :
406 28 : if (poParent->GetDimensionCount() < 2)
407 : {
408 1 : CPLError(CE_Failure, CPLE_NotSupported,
409 : "GetResampled() only supports 2 dimensions or more");
410 1 : return nullptr;
411 : }
412 :
413 27 : const auto &aoParentDims = poParent->GetDimensions();
414 27 : if (apoNewDimsIn.size() != aoParentDims.size())
415 : {
416 2 : CPLError(CE_Failure, CPLE_AppDefined,
417 : "GetResampled(): apoNewDims size should be the same as "
418 : "GetDimensionCount()");
419 2 : return nullptr;
420 : }
421 :
422 50 : std::vector<std::shared_ptr<GDALDimension>> apoNewDims;
423 25 : apoNewDims.reserve(apoNewDimsIn.size());
424 :
425 50 : std::vector<GUInt64> anBlockSize;
426 25 : anBlockSize.reserve(apoNewDimsIn.size());
427 50 : const auto &anParentBlockSize = poParent->GetBlockSize();
428 :
429 50 : auto apoParentDims = poParent->GetDimensions();
430 : // Special case for NASA EMIT datasets
431 30 : const bool bYXBandOrder = (apoParentDims.size() == 3 &&
432 7 : apoParentDims[0]->GetName() == "downtrack" &&
433 32 : apoParentDims[1]->GetName() == "crosstrack" &&
434 2 : apoParentDims[2]->GetName() == "bands");
435 :
436 : const size_t iYDimParent =
437 25 : bYXBandOrder ? 0 : poParent->GetDimensionCount() - 2;
438 : const size_t iXDimParent =
439 25 : bYXBandOrder ? 1 : poParent->GetDimensionCount() - 1;
440 :
441 77 : for (unsigned i = 0; i < apoNewDimsIn.size(); ++i)
442 : {
443 53 : if (i == iYDimParent || i == iXDimParent)
444 48 : continue;
445 5 : if (apoNewDimsIn[i] == nullptr)
446 : {
447 3 : apoNewDims.emplace_back(aoParentDims[i]);
448 : }
449 3 : else if (apoNewDimsIn[i]->GetSize() != aoParentDims[i]->GetSize() ||
450 1 : apoNewDimsIn[i]->GetName() != aoParentDims[i]->GetName())
451 : {
452 1 : CPLError(CE_Failure, CPLE_AppDefined,
453 : "GetResampled(): apoNewDims[%u] should be the same "
454 : "as its parent",
455 : i);
456 1 : return nullptr;
457 : }
458 : else
459 : {
460 1 : apoNewDims.emplace_back(aoParentDims[i]);
461 : }
462 4 : anBlockSize.emplace_back(anParentBlockSize[i]);
463 : }
464 :
465 : std::unique_ptr<GDALMDArrayResampledDataset> poParentDS(
466 48 : new GDALMDArrayResampledDataset(poParent, iXDimParent, iYDimParent));
467 :
468 24 : double dfXStart = 0.0;
469 24 : double dfXSpacing = 0.0;
470 24 : bool gotXSpacing = false;
471 48 : auto poNewDimX = apoNewDimsIn[iXDimParent];
472 24 : if (poNewDimX)
473 : {
474 4 : if (poNewDimX->GetSize() > static_cast<GUInt64>(INT_MAX))
475 : {
476 0 : CPLError(CE_Failure, CPLE_NotSupported,
477 : "Too big size for X dimension");
478 0 : return nullptr;
479 : }
480 4 : auto var = poNewDimX->GetIndexingVariable();
481 4 : if (var)
482 : {
483 2 : if (var->GetDimensionCount() != 1 ||
484 2 : var->GetDimensions()[0]->GetSize() != poNewDimX->GetSize() ||
485 1 : !var->IsRegularlySpaced(dfXStart, dfXSpacing))
486 : {
487 0 : CPLError(CE_Failure, CPLE_NotSupported,
488 : "New X dimension should be indexed by a regularly "
489 : "spaced variable");
490 0 : return nullptr;
491 : }
492 1 : gotXSpacing = true;
493 : }
494 : }
495 :
496 24 : double dfYStart = 0.0;
497 24 : double dfYSpacing = 0.0;
498 48 : auto poNewDimY = apoNewDimsIn[iYDimParent];
499 24 : bool gotYSpacing = false;
500 24 : if (poNewDimY)
501 : {
502 4 : if (poNewDimY->GetSize() > static_cast<GUInt64>(INT_MAX))
503 : {
504 0 : CPLError(CE_Failure, CPLE_NotSupported,
505 : "Too big size for Y dimension");
506 0 : return nullptr;
507 : }
508 4 : auto var = poNewDimY->GetIndexingVariable();
509 4 : if (var)
510 : {
511 2 : if (var->GetDimensionCount() != 1 ||
512 2 : var->GetDimensions()[0]->GetSize() != poNewDimY->GetSize() ||
513 1 : !var->IsRegularlySpaced(dfYStart, dfYSpacing))
514 : {
515 0 : CPLError(CE_Failure, CPLE_NotSupported,
516 : "New Y dimension should be indexed by a regularly "
517 : "spaced variable");
518 0 : return nullptr;
519 : }
520 1 : gotYSpacing = true;
521 : }
522 : }
523 :
524 : // This limitation could probably be removed
525 24 : if ((gotXSpacing && !gotYSpacing) || (!gotXSpacing && gotYSpacing))
526 : {
527 0 : CPLError(CE_Failure, CPLE_NotSupported,
528 : "Either none of new X or Y dimension should have an indexing "
529 : "variable, or both should both should have one.");
530 0 : return nullptr;
531 : }
532 :
533 48 : std::string osDstWKT;
534 24 : if (poTargetSRS)
535 : {
536 2 : char *pszDstWKT = nullptr;
537 2 : if (poTargetSRS->exportToWkt(&pszDstWKT) != OGRERR_NONE)
538 : {
539 0 : CPLFree(pszDstWKT);
540 0 : return nullptr;
541 : }
542 2 : osDstWKT = pszDstWKT;
543 2 : CPLFree(pszDstWKT);
544 : }
545 :
546 : // Use coordinate variables for geolocation array
547 48 : const auto apoCoordinateVars = poParent->GetCoordinateVariables();
548 24 : bool useGeolocationArray = false;
549 24 : if (apoCoordinateVars.size() >= 2)
550 : {
551 0 : std::shared_ptr<GDALMDArray> poLongVar;
552 0 : std::shared_ptr<GDALMDArray> poLatVar;
553 15 : for (const auto &poCoordVar : apoCoordinateVars)
554 : {
555 10 : const auto &osName = poCoordVar->GetName();
556 30 : const auto poAttr = poCoordVar->GetAttribute("standard_name");
557 20 : std::string osStandardName;
558 12 : if (poAttr && poAttr->GetDataType().GetClass() == GEDTC_STRING &&
559 2 : poAttr->GetDimensionCount() == 0)
560 : {
561 2 : const char *pszStandardName = poAttr->ReadAsString();
562 2 : if (pszStandardName)
563 2 : osStandardName = pszStandardName;
564 : }
565 21 : if (osName == "lon" || osName == "longitude" ||
566 21 : osName == "Longitude" || osStandardName == "longitude")
567 : {
568 5 : poLongVar = poCoordVar;
569 : }
570 6 : else if (osName == "lat" || osName == "latitude" ||
571 6 : osName == "Latitude" || osStandardName == "latitude")
572 : {
573 5 : poLatVar = poCoordVar;
574 : }
575 : }
576 5 : if (poLatVar != nullptr && poLongVar != nullptr)
577 : {
578 5 : const auto longDimCount = poLongVar->GetDimensionCount();
579 5 : const auto &longDims = poLongVar->GetDimensions();
580 5 : const auto latDimCount = poLatVar->GetDimensionCount();
581 5 : const auto &latDims = poLatVar->GetDimensions();
582 5 : const auto xDimSize = aoParentDims[iXDimParent]->GetSize();
583 5 : const auto yDimSize = aoParentDims[iYDimParent]->GetSize();
584 0 : if (longDimCount == 1 && longDims[0]->GetSize() == xDimSize &&
585 5 : latDimCount == 1 && latDims[0]->GetSize() == yDimSize)
586 : {
587 : // Geolocation arrays are 1D, and of consistent size with
588 : // the variable
589 0 : useGeolocationArray = true;
590 : }
591 1 : else if ((longDimCount == 2 ||
592 6 : (longDimCount == 3 && longDims[0]->GetSize() == 1)) &&
593 10 : longDims[longDimCount - 2]->GetSize() == yDimSize &&
594 10 : longDims[longDimCount - 1]->GetSize() == xDimSize &&
595 1 : (latDimCount == 2 ||
596 6 : (latDimCount == 3 && latDims[0]->GetSize() == 1)) &&
597 15 : latDims[latDimCount - 2]->GetSize() == yDimSize &&
598 5 : latDims[latDimCount - 1]->GetSize() == xDimSize)
599 :
600 : {
601 : // Geolocation arrays are 2D (or 3D with first dimension of
602 : // size 1, as found in Sentinel 5P products), and of consistent
603 : // size with the variable
604 5 : useGeolocationArray = true;
605 : }
606 : else
607 : {
608 0 : CPLDebug(
609 : "GDAL",
610 : "Longitude and latitude coordinate variables found, "
611 : "but their characteristics are not compatible of using "
612 : "them as geolocation arrays");
613 : }
614 5 : if (useGeolocationArray)
615 : {
616 10 : CPLDebug("GDAL",
617 : "Setting geolocation array from variables %s and %s",
618 5 : poLongVar->GetName().c_str(),
619 5 : poLatVar->GetName().c_str());
620 : const std::string osFilenameLong =
621 5 : VSIMemGenerateHiddenFilename("longitude.tif");
622 : const std::string osFilenameLat =
623 5 : VSIMemGenerateHiddenFilename("latitude.tif");
624 : std::unique_ptr<GDALDataset> poTmpLongDS(
625 : longDimCount == 1
626 0 : ? poLongVar->AsClassicDataset(0, 0)
627 20 : : poLongVar->AsClassicDataset(longDimCount - 1,
628 15 : longDimCount - 2));
629 5 : auto hTIFFLongDS = GDALTranslate(
630 : osFilenameLong.c_str(),
631 : GDALDataset::ToHandle(poTmpLongDS.get()), nullptr, nullptr);
632 : std::unique_ptr<GDALDataset> poTmpLatDS(
633 0 : latDimCount == 1 ? poLatVar->AsClassicDataset(0, 0)
634 20 : : poLatVar->AsClassicDataset(
635 15 : latDimCount - 1, latDimCount - 2));
636 5 : auto hTIFFLatDS = GDALTranslate(
637 : osFilenameLat.c_str(),
638 : GDALDataset::ToHandle(poTmpLatDS.get()), nullptr, nullptr);
639 5 : const bool bError =
640 5 : (hTIFFLatDS == nullptr || hTIFFLongDS == nullptr);
641 5 : GDALClose(hTIFFLongDS);
642 5 : GDALClose(hTIFFLatDS);
643 5 : if (bError)
644 : {
645 0 : VSIUnlink(osFilenameLong.c_str());
646 0 : VSIUnlink(osFilenameLat.c_str());
647 0 : return nullptr;
648 : }
649 :
650 5 : poParentDS->SetGeolocationArray(osFilenameLong, osFilenameLat);
651 : }
652 : }
653 : else
654 : {
655 0 : CPLDebug("GDAL",
656 : "Coordinate variables available for %s, but "
657 : "longitude and/or latitude variables were not identified",
658 0 : poParent->GetName().c_str());
659 : }
660 : }
661 :
662 : // Build gdalwarp arguments
663 48 : CPLStringList aosArgv;
664 :
665 24 : aosArgv.AddString("-of");
666 24 : aosArgv.AddString("VRT");
667 :
668 24 : aosArgv.AddString("-r");
669 24 : aosArgv.AddString(pszResampleAlg);
670 :
671 24 : if (!osDstWKT.empty())
672 : {
673 2 : aosArgv.AddString("-t_srs");
674 2 : aosArgv.AddString(osDstWKT.c_str());
675 : }
676 :
677 24 : if (useGeolocationArray)
678 5 : aosArgv.AddString("-geoloc");
679 :
680 24 : if (gotXSpacing && gotYSpacing)
681 : {
682 1 : const double dfXMin = dfXStart - dfXSpacing / 2;
683 : const double dfXMax =
684 1 : dfXMin + dfXSpacing * static_cast<double>(poNewDimX->GetSize());
685 1 : const double dfYMax = dfYStart - dfYSpacing / 2;
686 : const double dfYMin =
687 1 : dfYMax + dfYSpacing * static_cast<double>(poNewDimY->GetSize());
688 1 : aosArgv.AddString("-te");
689 1 : aosArgv.AddString(CPLSPrintf("%.17g", dfXMin));
690 1 : aosArgv.AddString(CPLSPrintf("%.17g", dfYMin));
691 1 : aosArgv.AddString(CPLSPrintf("%.17g", dfXMax));
692 1 : aosArgv.AddString(CPLSPrintf("%.17g", dfYMax));
693 : }
694 :
695 24 : if (poNewDimX && poNewDimY)
696 : {
697 3 : aosArgv.AddString("-ts");
698 : aosArgv.AddString(
699 3 : CPLSPrintf("%d", static_cast<int>(poNewDimX->GetSize())));
700 : aosArgv.AddString(
701 3 : CPLSPrintf("%d", static_cast<int>(poNewDimY->GetSize())));
702 : }
703 21 : else if (poNewDimX)
704 : {
705 1 : aosArgv.AddString("-ts");
706 : aosArgv.AddString(
707 1 : CPLSPrintf("%d", static_cast<int>(poNewDimX->GetSize())));
708 1 : aosArgv.AddString("0");
709 : }
710 20 : else if (poNewDimY)
711 : {
712 1 : aosArgv.AddString("-ts");
713 1 : aosArgv.AddString("0");
714 : aosArgv.AddString(
715 1 : CPLSPrintf("%d", static_cast<int>(poNewDimY->GetSize())));
716 : }
717 :
718 : // Create a warped VRT dataset
719 : GDALWarpAppOptions *psOptions =
720 24 : GDALWarpAppOptionsNew(aosArgv.List(), nullptr);
721 24 : GDALDatasetH hSrcDS = GDALDataset::ToHandle(poParentDS.get());
722 : std::unique_ptr<GDALDataset> poReprojectedDS(GDALDataset::FromHandle(
723 48 : GDALWarp("", nullptr, 1, &hSrcDS, psOptions, nullptr)));
724 24 : GDALWarpAppOptionsFree(psOptions);
725 24 : if (poReprojectedDS == nullptr)
726 3 : return nullptr;
727 :
728 : int nBlockXSize;
729 : int nBlockYSize;
730 21 : poReprojectedDS->GetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
731 21 : anBlockSize.emplace_back(nBlockYSize);
732 21 : anBlockSize.emplace_back(nBlockXSize);
733 :
734 21 : GDALGeoTransform gt;
735 21 : CPLErr eErr = poReprojectedDS->GetGeoTransform(gt);
736 21 : CPLAssert(eErr == CE_None);
737 21 : CPL_IGNORE_RET_VAL(eErr);
738 :
739 : auto poDimY = std::make_shared<GDALDimensionWeakIndexingVar>(
740 0 : std::string(), "dimY", GDAL_DIM_TYPE_HORIZONTAL_Y, "NORTH",
741 42 : poReprojectedDS->GetRasterYSize());
742 : auto varY = GDALMDArrayRegularlySpaced::Create(
743 63 : std::string(), poDimY->GetName(), poDimY, gt.yorig + gt.yscale / 2,
744 84 : gt.yscale, 0);
745 21 : poDimY->SetIndexingVariable(varY);
746 :
747 : auto poDimX = std::make_shared<GDALDimensionWeakIndexingVar>(
748 0 : std::string(), "dimX", GDAL_DIM_TYPE_HORIZONTAL_X, "EAST",
749 42 : poReprojectedDS->GetRasterXSize());
750 : auto varX = GDALMDArrayRegularlySpaced::Create(
751 63 : std::string(), poDimX->GetName(), poDimX, gt.xorig + gt.xscale / 2,
752 84 : gt.xscale, 0);
753 21 : poDimX->SetIndexingVariable(varX);
754 :
755 21 : apoNewDims.emplace_back(poDimY);
756 21 : apoNewDims.emplace_back(poDimX);
757 : auto newAr(std::shared_ptr<GDALMDArrayResampled>(
758 42 : new GDALMDArrayResampled(poParent, apoNewDims, anBlockSize)));
759 21 : newAr->SetSelf(newAr);
760 21 : if (poTargetSRS)
761 : {
762 2 : newAr->m_poSRS.reset(poTargetSRS->Clone());
763 : }
764 : else
765 : {
766 19 : newAr->m_poSRS = poParent->GetSpatialRef();
767 : }
768 21 : newAr->m_poVarX = varX;
769 21 : newAr->m_poVarY = varY;
770 21 : newAr->m_poReprojectedDS = std::move(poReprojectedDS);
771 21 : newAr->m_poParentDS = std::move(poParentDS);
772 :
773 : // If the input array is y,x,band ordered, the above newAr is
774 : // actually band,y,x ordered as it is more convenient for
775 : // GDALMDArrayResampled::IRead() implementation. But transpose that
776 : // array to the order of the input array
777 21 : if (bYXBandOrder)
778 4 : return newAr->Transpose(std::vector<int>{1, 2, 0});
779 :
780 19 : return newAr;
781 : }
782 :
783 : /************************************************************************/
784 : /* GDALMDArrayResampled::IRead() */
785 : /************************************************************************/
786 :
787 29 : bool GDALMDArrayResampled::IRead(const GUInt64 *arrayStartIdx,
788 : const size_t *count, const GInt64 *arrayStep,
789 : const GPtrDiff_t *bufferStride,
790 : const GDALExtendedDataType &bufferDataType,
791 : void *pDstBuffer) const
792 : {
793 29 : if (bufferDataType.GetClass() != GEDTC_NUMERIC)
794 0 : return false;
795 :
796 : struct Stack
797 : {
798 : size_t nIters = 0;
799 : GByte *dst_ptr = nullptr;
800 : GPtrDiff_t dst_inc_offset = 0;
801 : };
802 :
803 29 : const auto nDims = GetDimensionCount();
804 58 : std::vector<Stack> stack(nDims + 1); // +1 to avoid -Wnull-dereference
805 29 : const size_t nBufferDTSize = bufferDataType.GetSize();
806 92 : for (size_t i = 0; i < nDims; i++)
807 : {
808 63 : stack[i].dst_inc_offset =
809 63 : static_cast<GPtrDiff_t>(bufferStride[i] * nBufferDTSize);
810 : }
811 29 : stack[0].dst_ptr = static_cast<GByte *>(pDstBuffer);
812 :
813 29 : size_t dimIdx = 0;
814 29 : const size_t iDimY = nDims - 2;
815 29 : const size_t iDimX = nDims - 1;
816 : // Use an array to avoid a false positive warning from CLang Static
817 : // Analyzer about flushCaches being never read
818 29 : bool flushCaches[] = {false};
819 : const bool bYXBandOrder =
820 29 : m_poParentDS->m_iYDim == 0 && m_poParentDS->m_iXDim == 1;
821 :
822 38 : lbl_next_depth:
823 38 : if (dimIdx == iDimY)
824 : {
825 33 : if (flushCaches[0])
826 : {
827 5 : flushCaches[0] = false;
828 : // When changing of 2D slice, flush GDAL 2D buffers
829 5 : m_poParentDS->FlushCache(false);
830 5 : m_poReprojectedDS->FlushCache(false);
831 : }
832 :
833 33 : if (!GDALMDRasterIOFromBand(m_poReprojectedDS->GetRasterBand(1),
834 : GF_Read, iDimX, iDimY, arrayStartIdx, count,
835 : arrayStep, bufferStride, bufferDataType,
836 33 : stack[dimIdx].dst_ptr))
837 : {
838 0 : return false;
839 : }
840 : }
841 : else
842 : {
843 5 : stack[dimIdx].nIters = count[dimIdx];
844 5 : if (m_poParentDS->m_anOffset[bYXBandOrder ? 2 : dimIdx] !=
845 5 : arrayStartIdx[dimIdx])
846 : {
847 1 : flushCaches[0] = true;
848 : }
849 5 : m_poParentDS->m_anOffset[bYXBandOrder ? 2 : dimIdx] =
850 5 : arrayStartIdx[dimIdx];
851 : while (true)
852 : {
853 9 : dimIdx++;
854 9 : stack[dimIdx].dst_ptr = stack[dimIdx - 1].dst_ptr;
855 9 : goto lbl_next_depth;
856 9 : lbl_return_to_caller:
857 9 : dimIdx--;
858 9 : if ((--stack[dimIdx].nIters) == 0)
859 5 : break;
860 4 : flushCaches[0] = true;
861 4 : ++m_poParentDS->m_anOffset[bYXBandOrder ? 2 : dimIdx];
862 4 : stack[dimIdx].dst_ptr += stack[dimIdx].dst_inc_offset;
863 : }
864 : }
865 38 : if (dimIdx > 0)
866 9 : goto lbl_return_to_caller;
867 :
868 29 : return true;
869 : }
870 :
871 : //! @endcond
872 :
873 : /************************************************************************/
874 : /* GetResampled() */
875 : /************************************************************************/
876 :
877 : /** Return an array that is a resampled / reprojected view of the current array
878 : *
879 : * This is the same as the C function GDALMDArrayGetResampled().
880 : *
881 : * Currently this method can only resample along the last 2 dimensions, unless
882 : * orthorectifying a NASA EMIT dataset.
883 : *
884 : * For NASA EMIT datasets, if apoNewDims[] and poTargetSRS is NULL, the
885 : * geometry lookup table (GLT) is used by default for fast orthorectification.
886 : *
887 : * Options available are:
888 : * <ul>
889 : * <li>EMIT_ORTHORECTIFICATION=YES/NO: defaults to YES for a NASA EMIT dataset.
890 : * Can be set to NO to use generic reprojection method.
891 : * </li>
892 : * <li>USE_GOOD_WAVELENGTHS=YES/NO: defaults to YES. Only used for EMIT
893 : * orthorectification to take into account the value of the
894 : * /sensor_band_parameters/good_wavelengths array to decide if slices of the
895 : * current array along the band dimension are valid.</li>
896 : * </ul>
897 : *
898 : * @param apoNewDims New dimensions. Its size should be GetDimensionCount().
899 : * apoNewDims[i] can be NULL to let the method automatically
900 : * determine it.
901 : * @param resampleAlg Resampling algorithm
902 : * @param poTargetSRS Target SRS, or nullptr
903 : * @param papszOptions NULL-terminated list of options, or NULL.
904 : *
905 : * @return a new array, that holds a reference to the original one, and thus is
906 : * a view of it (not a copy), or nullptr in case of error.
907 : *
908 : * @since 3.4
909 : */
910 38 : std::shared_ptr<GDALMDArray> GDALMDArray::GetResampled(
911 : const std::vector<std::shared_ptr<GDALDimension>> &apoNewDims,
912 : GDALRIOResampleAlg resampleAlg, const OGRSpatialReference *poTargetSRS,
913 : CSLConstList papszOptions) const
914 : {
915 76 : auto self = std::dynamic_pointer_cast<GDALMDArray>(m_pSelf.lock());
916 38 : if (!self)
917 : {
918 0 : CPLError(CE_Failure, CPLE_AppDefined,
919 : "Driver implementation issue: m_pSelf not set !");
920 0 : return nullptr;
921 : }
922 38 : if (GetDataType().GetClass() != GEDTC_NUMERIC)
923 : {
924 0 : CPLError(CE_Failure, CPLE_AppDefined,
925 : "GetResampled() only supports numeric data type");
926 0 : return nullptr;
927 : }
928 :
929 : // Special case for NASA EMIT datasets
930 76 : auto apoDims = GetDimensions();
931 36 : if (poTargetSRS == nullptr &&
932 59 : ((apoDims.size() == 3 && apoDims[0]->GetName() == "downtrack" &&
933 20 : apoDims[1]->GetName() == "crosstrack" &&
934 10 : apoDims[2]->GetName() == "bands" &&
935 48 : (apoNewDims == std::vector<std::shared_ptr<GDALDimension>>(3) ||
936 1 : apoNewDims ==
937 42 : std::vector<std::shared_ptr<GDALDimension>>{nullptr, nullptr,
938 30 : apoDims[2]})) ||
939 51 : (apoDims.size() == 2 && apoDims[0]->GetName() == "downtrack" &&
940 3 : apoDims[1]->GetName() == "crosstrack" &&
941 77 : apoNewDims == std::vector<std::shared_ptr<GDALDimension>>(2))) &&
942 13 : CPLTestBool(CSLFetchNameValueDef(papszOptions,
943 : "EMIT_ORTHORECTIFICATION", "YES")))
944 : {
945 9 : auto poRootGroup = GetRootGroup();
946 9 : if (poRootGroup)
947 : {
948 18 : auto poAttrGeotransform = poRootGroup->GetAttribute("geotransform");
949 18 : auto poLocationGroup = poRootGroup->OpenGroup("location");
950 9 : if (poAttrGeotransform &&
951 9 : poAttrGeotransform->GetDataType().GetClass() == GEDTC_NUMERIC &&
952 9 : poAttrGeotransform->GetDimensionCount() == 1 &&
953 27 : poAttrGeotransform->GetDimensionsSize()[0] == 6 &&
954 9 : poLocationGroup)
955 : {
956 18 : auto poGLT_X = poLocationGroup->OpenMDArray("glt_x");
957 18 : auto poGLT_Y = poLocationGroup->OpenMDArray("glt_y");
958 27 : if (poGLT_X && poGLT_X->GetDimensionCount() == 2 &&
959 18 : poGLT_X->GetDimensions()[0]->GetName() == "ortho_y" &&
960 18 : poGLT_X->GetDimensions()[1]->GetName() == "ortho_x" &&
961 27 : poGLT_Y && poGLT_Y->GetDimensionCount() == 2 &&
962 27 : poGLT_Y->GetDimensions()[0]->GetName() == "ortho_y" &&
963 9 : poGLT_Y->GetDimensions()[1]->GetName() == "ortho_x")
964 : {
965 : return CreateGLTOrthorectified(
966 : self, poRootGroup, poGLT_X, poGLT_Y,
967 : /* nGLTIndexOffset = */ -1,
968 18 : poAttrGeotransform->ReadAsDoubleArray(), papszOptions);
969 : }
970 : }
971 : }
972 : }
973 :
974 29 : if (CPLTestBool(CSLFetchNameValueDef(papszOptions,
975 : "EMIT_ORTHORECTIFICATION", "NO")))
976 : {
977 0 : CPLError(CE_Failure, CPLE_AppDefined,
978 : "EMIT_ORTHORECTIFICATION required, but dataset and/or "
979 : "parameters are not compatible with it");
980 0 : return nullptr;
981 : }
982 :
983 : return GDALMDArrayResampled::Create(self, apoNewDims, resampleAlg,
984 29 : poTargetSRS, papszOptions);
985 : }
|
