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