Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Virtual GDAL Datasets
4 : * Purpose: Tile index based VRT
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2023, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : /*! @cond Doxygen_Suppress */
14 :
15 : #include <array>
16 : #include <algorithm>
17 : #include <atomic>
18 : #include <cmath>
19 : #include <limits>
20 : #include <mutex>
21 : #include <set>
22 : #include <tuple>
23 : #include <utility>
24 : #include <vector>
25 :
26 : #include "cpl_port.h"
27 : #include "cpl_error_internal.h"
28 : #include "cpl_json.h"
29 : #include "cpl_mem_cache.h"
30 : #include "cpl_minixml.h"
31 : #include "cpl_quad_tree.h"
32 : #include "cpl_worker_thread_pool.h"
33 : #include "vrtdataset.h"
34 : #include "vrt_priv.h"
35 : #include "ogrsf_frmts.h"
36 : #include "ogrwarpedlayer.h"
37 : #include "gdal_frmts.h"
38 : #include "gdal_proxy.h"
39 : #include "gdalsubdatasetinfo.h"
40 : #include "gdal_thread_pool.h"
41 : #include "gdal_utils.h"
42 :
43 : #include "gdalalg_raster_index.h"
44 :
45 : #ifdef USE_NEON_OPTIMIZATIONS
46 : #define USE_SSE2_OPTIM
47 : #define USE_SSE41_OPTIM
48 : #include "include_sse2neon.h"
49 : #elif defined(__SSE2__) || defined(_M_X64)
50 : #define USE_SSE2_OPTIM
51 : #include <emmintrin.h>
52 : // MSVC doesn't define __SSE4_1__, but if -arch:AVX2 is enabled, we do have SSE4.1
53 : #if defined(__SSE4_1__) || defined(__AVX2__)
54 : #define USE_SSE41_OPTIM
55 : #include <smmintrin.h>
56 : #endif
57 : #endif
58 :
59 : #ifndef _
60 : #define _(x) (x)
61 : #endif
62 :
63 : constexpr const char *GTI_PREFIX = "GTI:";
64 :
65 : constexpr const char *MD_DS_TILE_INDEX_LAYER = "TILE_INDEX_LAYER";
66 : constexpr const char *MD_DS_TILE_INDEX_SQL = "TILE_INDEX_SQL";
67 : constexpr const char *MD_DS_TILE_INDEX_SPATIAL_SQL = "TILE_INDEX_SPATIAL_SQL";
68 :
69 : constexpr const char *MD_RESX = "RESX";
70 : constexpr const char *MD_RESY = "RESY";
71 : constexpr const char *MD_BAND_COUNT = "BAND_COUNT";
72 : constexpr const char *MD_DATA_TYPE = "DATA_TYPE";
73 : constexpr const char *MD_NODATA = "NODATA";
74 : constexpr const char *MD_MINX = "MINX";
75 : constexpr const char *MD_MINY = "MINY";
76 : constexpr const char *MD_MAXX = "MAXX";
77 : constexpr const char *MD_MAXY = "MAXY";
78 : constexpr const char *MD_GEOTRANSFORM = "GEOTRANSFORM";
79 : constexpr const char *MD_XSIZE = "XSIZE";
80 : constexpr const char *MD_YSIZE = "YSIZE";
81 : constexpr const char *MD_COLOR_INTERPRETATION = "COLOR_INTERPRETATION";
82 : constexpr const char *MD_SRS = "SRS";
83 : constexpr const char *MD_LOCATION_FIELD = "LOCATION_FIELD";
84 : constexpr const char *MD_SORT_FIELD = "SORT_FIELD";
85 : constexpr const char *MD_SORT_FIELD_ASC = "SORT_FIELD_ASC";
86 : constexpr const char *MD_BLOCK_X_SIZE = "BLOCKXSIZE";
87 : constexpr const char *MD_BLOCK_Y_SIZE = "BLOCKYSIZE";
88 : constexpr const char *MD_MASK_BAND = "MASK_BAND";
89 : constexpr const char *MD_RESAMPLING = "RESAMPLING";
90 : constexpr const char *MD_INTERLEAVE = "INTERLEAVE";
91 :
92 : constexpr const char *const apszTIOptions[] = {MD_RESX,
93 : MD_RESY,
94 : MD_BAND_COUNT,
95 : MD_DATA_TYPE,
96 : MD_NODATA,
97 : MD_MINX,
98 : MD_MINY,
99 : MD_MAXX,
100 : MD_MAXY,
101 : MD_GEOTRANSFORM,
102 : MD_XSIZE,
103 : MD_YSIZE,
104 : MD_COLOR_INTERPRETATION,
105 : MD_SRS,
106 : MD_LOCATION_FIELD,
107 : MD_SORT_FIELD,
108 : MD_SORT_FIELD_ASC,
109 : MD_BLOCK_X_SIZE,
110 : MD_BLOCK_Y_SIZE,
111 : MD_MASK_BAND,
112 : MD_RESAMPLING,
113 : MD_INTERLEAVE};
114 :
115 : constexpr const char *const MD_BAND_OFFSET = "OFFSET";
116 : constexpr const char *const MD_BAND_SCALE = "SCALE";
117 : constexpr const char *const MD_BAND_UNITTYPE = "UNITTYPE";
118 : constexpr const char *const apszReservedBandItems[] = {
119 : MD_BAND_OFFSET, MD_BAND_SCALE, MD_BAND_UNITTYPE};
120 :
121 : constexpr const char *GTI_XML_BANDCOUNT = "BandCount";
122 : constexpr const char *GTI_XML_DATATYPE = "DataType";
123 : constexpr const char *GTI_XML_NODATAVALUE = "NoDataValue";
124 : constexpr const char *GTI_XML_COLORINTERP = "ColorInterp";
125 : constexpr const char *GTI_XML_LOCATIONFIELD = "LocationField";
126 : constexpr const char *GTI_XML_SORTFIELD = "SortField";
127 : constexpr const char *GTI_XML_SORTFIELDASC = "SortFieldAsc";
128 : constexpr const char *GTI_XML_MASKBAND = "MaskBand";
129 : constexpr const char *GTI_XML_OVERVIEW_ELEMENT = "Overview";
130 : constexpr const char *GTI_XML_OVERVIEW_DATASET = "Dataset";
131 : constexpr const char *GTI_XML_OVERVIEW_LAYER = "Layer";
132 : constexpr const char *GTI_XML_OVERVIEW_FACTOR = "Factor";
133 :
134 : constexpr const char *GTI_XML_BAND_ELEMENT = "Band";
135 : constexpr const char *GTI_XML_BAND_NUMBER = "band";
136 : constexpr const char *GTI_XML_BAND_DATATYPE = "dataType";
137 : constexpr const char *GTI_XML_BAND_DESCRIPTION = "Description";
138 : constexpr const char *GTI_XML_BAND_OFFSET = "Offset";
139 : constexpr const char *GTI_XML_BAND_SCALE = "Scale";
140 : constexpr const char *GTI_XML_BAND_NODATAVALUE = "NoDataValue";
141 : constexpr const char *GTI_XML_BAND_UNITTYPE = "UnitType";
142 : constexpr const char *GTI_XML_BAND_COLORINTERP = "ColorInterp";
143 : constexpr const char *GTI_XML_CATEGORYNAMES = "CategoryNames";
144 : constexpr const char *GTI_XML_COLORTABLE = "ColorTable";
145 : constexpr const char *GTI_XML_RAT = "GDALRasterAttributeTable";
146 :
147 : /************************************************************************/
148 : /* ENDS_WITH_CI() */
149 : /************************************************************************/
150 :
151 67916 : static inline bool ENDS_WITH_CI(const char *a, const char *b)
152 : {
153 67916 : return strlen(a) >= strlen(b) && EQUAL(a + strlen(a) - strlen(b), b);
154 : }
155 :
156 : /************************************************************************/
157 : /* GTISharedSourceKey */
158 : /************************************************************************/
159 :
160 : struct GTISharedSourceKey
161 : {
162 : std::string osTileName{};
163 : std::vector<int> anBands{};
164 :
165 409 : bool operator==(const GTISharedSourceKey &other) const
166 : {
167 409 : return osTileName == other.osTileName && anBands == other.anBands;
168 : }
169 :
170 : CPL_NOSANITIZE_UNSIGNED_INT_OVERFLOW
171 1420 : size_t getHash() const noexcept
172 : {
173 1420 : size_t h = std::hash<std::string>{}(osTileName);
174 1467 : for (int b : anBands)
175 : {
176 : // Cf https://www.boost.org/doc/libs/1_36_0/doc/html/hash/reference.html#boost.hash_combine
177 47 : h ^= std::hash<int>{}(b) + 0x9e3779b9 + (h << 6) + (h >> 2);
178 : }
179 1420 : return h;
180 : }
181 : };
182 :
183 : namespace std
184 : {
185 : template <> struct hash<GTISharedSourceKey>
186 : {
187 1420 : size_t operator()(const GTISharedSourceKey &val) const noexcept
188 : {
189 1420 : return val.getHash();
190 : }
191 : };
192 : } // namespace std
193 :
194 : /************************************************************************/
195 : /* GDALTileIndexDataset */
196 : /************************************************************************/
197 :
198 : class GDALTileIndexBand;
199 :
200 : class GDALTileIndexDataset final : public GDALPamDataset
201 : {
202 : public:
203 : GDALTileIndexDataset();
204 : ~GDALTileIndexDataset() override;
205 :
206 : bool Open(GDALOpenInfo *poOpenInfo);
207 :
208 : CPLErr FlushCache(bool bAtClosing) override;
209 :
210 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
211 : const OGRSpatialReference *GetSpatialRef() const override;
212 :
213 : CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
214 : int nYSize, void *pData, int nBufXSize, int nBufYSize,
215 : GDALDataType eBufType, int nBandCount,
216 : BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
217 : GSpacing nLineSpace, GSpacing nBandSpace,
218 : GDALRasterIOExtraArg *psExtraArg) override;
219 :
220 : const char *GetMetadataItem(const char *pszName,
221 : const char *pszDomain) override;
222 : CPLErr SetMetadataItem(const char *pszName, const char *pszValue,
223 : const char *pszDomain) override;
224 : CPLErr SetMetadata(CSLConstList papszMD, const char *pszDomain) override;
225 :
226 : void LoadOverviews();
227 :
228 : std::vector<GTISourceDesc> GetSourcesMoreRecentThan(int64_t mTime);
229 :
230 : private:
231 : friend class GDALTileIndexBand;
232 :
233 : //! Optional GTI XML
234 : CPLXMLTreeCloser m_psXMLTree{nullptr};
235 :
236 : //! Whether the GTI XML might be modified (by SetMetadata/SetMetadataItem)
237 : bool m_bXMLUpdatable = false;
238 :
239 : //! Whether the GTI XML has been modified (by SetMetadata/SetMetadataItem)
240 : bool m_bXMLModified = false;
241 :
242 : //! Unique string (without the process) for this tile index. Passed to
243 : //! GDALProxyPoolDataset to ensure that sources are unique for a given owner
244 : const std::string m_osUniqueHandle;
245 :
246 : //! Vector dataset with the sources
247 : std::unique_ptr<GDALDataset> m_poVectorDS{};
248 :
249 : //! Generic SQL request to return features. May be empty.
250 : std::string m_osSQL{};
251 :
252 : //! SQL request to return features with placeholders for spatial filtering. May be empty
253 : std::string m_osSpatialSQL{};
254 :
255 : //! Vector layer with the sources
256 : OGRLayer *m_poLayer = nullptr;
257 :
258 : //! Whether m_poLayer should be freed with m_poVectorDS->ReleaseResultSet()
259 : bool m_bIsSQLResultLayer = false;
260 :
261 : //! When the SRS of m_poLayer is not the one we expose
262 : std::unique_ptr<OGRLayer> m_poWarpedLayerKeeper{};
263 :
264 : //! Geotransform matrix of the tile index
265 : GDALGeoTransform m_gt{};
266 :
267 : //! Index of the "location" (or alternate name given by user) field
268 : //! (within m_poLayer->GetLayerDefn()), that contain source dataset names.
269 : int m_nLocationFieldIndex = -1;
270 :
271 : //! SRS of the tile index.
272 : OGRSpatialReference m_oSRS{};
273 :
274 : struct SharedDataset
275 : {
276 : //! Source dataset (possibly warped).
277 : std::shared_ptr<GDALDataset> poDS{};
278 :
279 : //! Source dataset, raw/unwarped
280 : GDALDataset *poUnreprojectedDS = nullptr;
281 :
282 : //! Whether (nBandCount, panBandMap) is taken into account by poDS
283 : bool bBandMapTakenIntoAccount = false;
284 : };
285 :
286 : //! Cache from dataset name to dataset handle.
287 : //! Note that the dataset objects are ultimately GDALProxyPoolDataset,
288 : //! and that the GDALProxyPoolDataset limits the number of simultaneously
289 : //! opened real datasets (controlled by GDAL_MAX_DATASET_POOL_SIZE). Hence 500 is not too big.
290 : lru11::Cache<GTISharedSourceKey, std::shared_ptr<SharedDataset>>
291 : m_oMapSharedSources{500};
292 :
293 : //! Mask band (e.g. for JPEG compressed + mask band)
294 : std::unique_ptr<GDALTileIndexBand> m_poMaskBand{};
295 :
296 : //! Whether all bands of the tile index have the same data type.
297 : bool m_bSameDataType = true;
298 :
299 : //! Whether all bands of the tile index have the same nodata value.
300 : bool m_bSameNoData = true;
301 :
302 : //! Whether a band interleave must be exposed.
303 : bool m_bBandInterleave = false;
304 :
305 : //! Minimum X of the current pixel request, in georeferenced units.
306 : double m_dfLastMinXFilter = std::numeric_limits<double>::quiet_NaN();
307 :
308 : //! Minimum Y of the current pixel request, in georeferenced units.
309 : double m_dfLastMinYFilter = std::numeric_limits<double>::quiet_NaN();
310 :
311 : //! Maximum X of the current pixel request, in georeferenced units.
312 : double m_dfLastMaxXFilter = std::numeric_limits<double>::quiet_NaN();
313 :
314 : //! Maximum Y of the current pixel request, in georeferenced units.
315 : double m_dfLastMaxYFilter = std::numeric_limits<double>::quiet_NaN();
316 :
317 : //! Bands for which m_aoSourceDesc is valid (only if m_bBandInterleave)
318 : std::vector<int> m_anLastBands{};
319 :
320 : //! Index of the field (within m_poLayer->GetLayerDefn()) used to sort, or -1 if none.
321 : int m_nSortFieldIndex = -1;
322 :
323 : //! Whether sorting must be ascending (true) or descending (false).
324 : bool m_bSortFieldAsc = true;
325 :
326 : //! Resampling method by default for warping or when a source has not
327 : //! the same resolution as the tile index.
328 : std::string m_osResampling = "near";
329 : GDALRIOResampleAlg m_eResampling = GRIORA_NearestNeighbour;
330 :
331 : //! WKT2 representation of the tile index SRS (if needed, typically for on-the-fly warping).
332 : std::string m_osWKT{};
333 :
334 : //! Whether we had to open of the sources at tile index opening.
335 : bool m_bScannedOneFeatureAtOpening = false;
336 :
337 : //! Array of overview descriptors.
338 : //! Each descriptor is a tuple (dataset_name, concatenated_open_options, layer_name, overview_factor).
339 : std::vector<std::tuple<std::string, CPLStringList, std::string, double>>
340 : m_aoOverviewDescriptor{};
341 :
342 : //! Array of overview datasets.
343 : std::vector<std::unique_ptr<GDALDataset>> m_apoOverviews{};
344 :
345 : //! Cache of buffers used by VRTComplexSource to avoid memory reallocation.
346 : VRTSource::WorkingState m_oWorkingState{};
347 :
348 : //! Used by IRasterIO() when using multi-threading
349 : struct QueueWorkingStates
350 : {
351 : std::mutex oMutex{};
352 : std::vector<std::unique_ptr<VRTSource::WorkingState>> oStates{};
353 : };
354 :
355 : //! Used by IRasterIO() when using multi-threading
356 : QueueWorkingStates m_oQueueWorkingStates{};
357 :
358 : //! Structure describing one of the source raster in the tile index.
359 : struct SourceDesc
360 : {
361 : //! Source dataset name.
362 : std::string osName{};
363 :
364 : //! Source dataset (possibly warped).
365 : std::shared_ptr<GDALDataset> poDS{};
366 :
367 : //! Source dataset, raw/unwarped
368 : GDALDataset *poUnreprojectedDS = nullptr;
369 :
370 : //! Whether (nBandCount, panBandMap) is taken into account by poDS
371 : bool bBandMapTakenIntoAccount = false;
372 :
373 : //! VRTSimpleSource or VRTComplexSource for the source.
374 : std::unique_ptr<VRTSimpleSource> poSource{};
375 :
376 : //! OGRFeature corresponding to the source in the tile index.
377 : std::unique_ptr<OGRFeature> poFeature{};
378 :
379 : //! Work buffer containing the value of the mask band for the current pixel query.
380 : mutable std::vector<GByte> abyMask{};
381 :
382 : //! Whether the source covers the whole area of interest of the current pixel query.
383 : bool bCoversWholeAOI = false;
384 :
385 : //! Whether the source has a nodata value at least in one of its band.
386 : bool bHasNoData = false;
387 :
388 : //! Whether all bands of the source have the same nodata value.
389 : bool bSameNoData = false;
390 :
391 : //! Nodata value of all bands (when bSameNoData == true).
392 : double dfSameNoData = 0;
393 :
394 : //! Mask band of the source.
395 : GDALRasterBand *poMaskBand = nullptr;
396 : };
397 :
398 : //! Array of sources participating to the current pixel query.
399 : std::vector<SourceDesc> m_aoSourceDesc{};
400 :
401 : //! Maximum number of threads. Updated by CollectSources().
402 : int m_nNumThreads = -1;
403 :
404 : //! Whereas the multi-threading rendering code path must be used. Updated by CollectSources().
405 : bool m_bLastMustUseMultiThreading = false;
406 :
407 : //! Whether the GTI file is a STAC collection
408 : bool m_bSTACCollection = false;
409 :
410 : std::string m_osWarpMemory{};
411 :
412 : //! From a source dataset name, return its SourceDesc description structure.
413 : bool GetSourceDesc(const std::string &osTileName, SourceDesc &oSourceDesc,
414 : std::mutex *pMutex, int nBandCount,
415 : const int *panBandMap);
416 :
417 : //! Collect sources corresponding to the georeferenced window of interest,
418 : //! and store them in m_aoSourceDesc[].
419 : bool CollectSources(double dfXOff, double dfYOff, double dfXSize,
420 : double dfYSize, int nBandCount, const int *panBandMap,
421 : bool bMultiThreadAllowed);
422 :
423 : //! Sort sources according to m_nSortFieldIndex.
424 : void SortSourceDesc();
425 :
426 : //! Whether the output buffer needs to be nodata initialized, or if
427 : //! sources are fully covering it.
428 : bool NeedInitBuffer(int nBandCount, const int *panBandMap) const;
429 :
430 : //! Nodata initialize the output buffer.
431 : void InitBuffer(void *pData, int nBufXSize, int nBufYSize,
432 : GDALDataType eBufType, int nBandCount,
433 : const int *panBandMap, GSpacing nPixelSpace,
434 : GSpacing nLineSpace, GSpacing nBandSpace) const;
435 :
436 : //! Render one source. Used by IRasterIO()
437 : CPLErr RenderSource(const SourceDesc &oSourceDesc, bool bNeedInitBuffer,
438 : int nBandNrMax, int nXOff, int nYOff, int nXSize,
439 : int nYSize, double dfXOff, double dfYOff,
440 : double dfXSize, double dfYSize, int nBufXSize,
441 : int nBufYSize, void *pData, GDALDataType eBufType,
442 : int nBandCount, BANDMAP_TYPE panBandMap,
443 : GSpacing nPixelSpace, GSpacing nLineSpace,
444 : GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg,
445 : VRTSource::WorkingState &oWorkingState) const;
446 :
447 : //! Whether m_poVectorDS supports SetMetadata()/SetMetadataItem()
448 : bool TileIndexSupportsEditingLayerMetadata() const;
449 :
450 : //! Return number of threads that can be used
451 : int GetNumThreads() const;
452 :
453 : /** Structure used to declare a threaded job to satisfy IRasterIO()
454 : * on a given source.
455 : */
456 : struct RasterIOJob
457 : {
458 : std::atomic<int> *pnCompletedJobs = nullptr;
459 : std::atomic<bool> *pbSuccess = nullptr;
460 : CPLErrorAccumulator *poErrorAccumulator = nullptr;
461 : GDALTileIndexDataset *poDS = nullptr;
462 : GDALTileIndexDataset::QueueWorkingStates *poQueueWorkingStates =
463 : nullptr;
464 : int nBandNrMax = 0;
465 : bool bSTACCollection = false;
466 :
467 : int nXOff = 0;
468 : int nYOff = 0;
469 : int nXSize = 0;
470 : int nYSize = 0;
471 : void *pData = nullptr;
472 : int nBufXSize = 0;
473 : int nBufYSize = 0;
474 : int nBandCount = 0;
475 : BANDMAP_TYPE panBandMap = nullptr;
476 : GDALDataType eBufType = GDT_Unknown;
477 : GSpacing nPixelSpace = 0;
478 : GSpacing nLineSpace = 0;
479 : GSpacing nBandSpace = 0;
480 : GDALRasterIOExtraArg *psExtraArg = nullptr;
481 :
482 : std::string osTileName{};
483 :
484 : static void Func(void *pData);
485 : };
486 :
487 : CPL_DISALLOW_COPY_ASSIGN(GDALTileIndexDataset)
488 : };
489 :
490 : /************************************************************************/
491 : /* GDALTileIndexBand */
492 : /************************************************************************/
493 :
494 : class GDALTileIndexBand final : public GDALPamRasterBand
495 : {
496 : public:
497 : GDALTileIndexBand(GDALTileIndexDataset *poDSIn, int nBandIn,
498 : GDALDataType eDT, int nBlockXSizeIn, int nBlockYSizeIn);
499 :
500 110 : double GetNoDataValue(int *pbHasNoData) override
501 : {
502 110 : if (pbHasNoData)
503 107 : *pbHasNoData = m_bNoDataValueSet;
504 110 : return m_dfNoDataValue;
505 : }
506 :
507 68 : GDALColorInterp GetColorInterpretation() override
508 : {
509 68 : return m_eColorInterp;
510 : }
511 :
512 : CPLErr IReadBlock(int nBlockXOff, int nBlockYOff, void *pData) override;
513 :
514 : CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
515 : int nYSize, void *pData, int nBufXSize, int nBufYSize,
516 : GDALDataType eBufType, GSpacing nPixelSpace,
517 : GSpacing nLineSpace,
518 : GDALRasterIOExtraArg *psExtraArg) override;
519 :
520 : int IGetDataCoverageStatus(int nXOff, int nYOff, int nXSize, int nYSize,
521 : int nMaskFlagStop, double *pdfDataPct) override;
522 :
523 32 : int GetMaskFlags() override
524 : {
525 32 : if (m_poDS->m_poMaskBand && m_poDS->m_poMaskBand.get() != this)
526 4 : return GMF_PER_DATASET;
527 28 : return GDALPamRasterBand::GetMaskFlags();
528 : }
529 :
530 34 : GDALRasterBand *GetMaskBand() override
531 : {
532 34 : if (m_poDS->m_poMaskBand && m_poDS->m_poMaskBand.get() != this)
533 7 : return m_poDS->m_poMaskBand.get();
534 27 : return GDALPamRasterBand::GetMaskBand();
535 : }
536 :
537 13 : double GetOffset(int *pbHasValue) override
538 : {
539 13 : int bHasValue = FALSE;
540 13 : double dfVal = GDALPamRasterBand::GetOffset(&bHasValue);
541 13 : if (bHasValue)
542 : {
543 0 : if (pbHasValue)
544 0 : *pbHasValue = true;
545 0 : return dfVal;
546 : }
547 13 : if (pbHasValue)
548 10 : *pbHasValue = !std::isnan(m_dfOffset);
549 13 : return std::isnan(m_dfOffset) ? 0.0 : m_dfOffset;
550 : }
551 :
552 13 : double GetScale(int *pbHasValue) override
553 : {
554 13 : int bHasValue = FALSE;
555 13 : double dfVal = GDALPamRasterBand::GetScale(&bHasValue);
556 13 : if (bHasValue)
557 : {
558 0 : if (pbHasValue)
559 0 : *pbHasValue = true;
560 0 : return dfVal;
561 : }
562 13 : if (pbHasValue)
563 10 : *pbHasValue = !std::isnan(m_dfScale);
564 13 : return std::isnan(m_dfScale) ? 1.0 : m_dfScale;
565 : }
566 :
567 9 : const char *GetUnitType() override
568 : {
569 9 : const char *pszVal = GDALPamRasterBand::GetUnitType();
570 9 : if (pszVal && *pszVal)
571 0 : return pszVal;
572 9 : return m_osUnit.c_str();
573 : }
574 :
575 5 : char **GetCategoryNames() override
576 : {
577 5 : return m_aosCategoryNames.List();
578 : }
579 :
580 11 : GDALColorTable *GetColorTable() override
581 : {
582 11 : return m_poColorTable.get();
583 : }
584 :
585 5 : GDALRasterAttributeTable *GetDefaultRAT() override
586 : {
587 5 : return m_poRAT.get();
588 : }
589 :
590 : int GetOverviewCount() override;
591 : GDALRasterBand *GetOverview(int iOvr) override;
592 :
593 : char **GetMetadataDomainList() override;
594 : const char *GetMetadataItem(const char *pszName,
595 : const char *pszDomain) override;
596 : CPLErr SetMetadataItem(const char *pszName, const char *pszValue,
597 : const char *pszDomain) override;
598 : CPLErr SetMetadata(CSLConstList papszMD, const char *pszDomain) override;
599 :
600 : private:
601 : friend class GDALTileIndexDataset;
602 :
603 : //! Dataset that owns this band.
604 : GDALTileIndexDataset *m_poDS = nullptr;
605 :
606 : //! Whether a nodata value is set to this band.
607 : bool m_bNoDataValueSet = false;
608 :
609 : //! Nodata value.
610 : double m_dfNoDataValue = 0;
611 :
612 : //! Color interpretation.
613 : GDALColorInterp m_eColorInterp = GCI_Undefined;
614 :
615 : //! Cached value for GetMetadataItem("Pixel_X_Y", "LocationInfo").
616 : std::string m_osLastLocationInfo{};
617 :
618 : //! Scale value (returned by GetScale())
619 : double m_dfScale = std::numeric_limits<double>::quiet_NaN();
620 :
621 : //! Offset value (returned by GetOffset())
622 : double m_dfOffset = std::numeric_limits<double>::quiet_NaN();
623 :
624 : //! Unit type (returned by GetUnitType()).
625 : std::string m_osUnit{};
626 :
627 : //! Category names (returned by GetCategoryNames()).
628 : CPLStringList m_aosCategoryNames{};
629 :
630 : //! Color table (returned by GetColorTable()).
631 : std::unique_ptr<GDALColorTable> m_poColorTable{};
632 :
633 : //! Raster attribute table (returned by GetDefaultRAT()).
634 : std::unique_ptr<GDALRasterAttributeTable> m_poRAT{};
635 :
636 : CPL_DISALLOW_COPY_ASSIGN(GDALTileIndexBand)
637 : };
638 :
639 : /************************************************************************/
640 : /* IsSameNaNAware() */
641 : /************************************************************************/
642 :
643 313 : static inline bool IsSameNaNAware(double a, double b)
644 : {
645 313 : return a == b || (std::isnan(a) && std::isnan(b));
646 : }
647 :
648 : /************************************************************************/
649 : /* GDALTileIndexDataset() */
650 : /************************************************************************/
651 :
652 286 : GDALTileIndexDataset::GDALTileIndexDataset()
653 286 : : m_osUniqueHandle(CPLSPrintf("%p", this))
654 : {
655 286 : }
656 :
657 : /************************************************************************/
658 : /* GetAbsoluteFileName() */
659 : /************************************************************************/
660 :
661 613 : static std::string GetAbsoluteFileName(const char *pszTileName,
662 : const char *pszVRTName,
663 : bool bIsStacCollection)
664 : {
665 : std::string osRet =
666 1839 : bIsStacCollection ? VSIURIToVSIPath(pszTileName) : pszTileName;
667 613 : if (osRet != pszTileName)
668 6 : return osRet;
669 :
670 607 : if (CPLIsFilenameRelative(pszTileName) &&
671 625 : !STARTS_WITH(pszTileName, "<VRTDataset") &&
672 18 : !STARTS_WITH(pszVRTName, "<GDALTileIndexDataset"))
673 : {
674 18 : const auto oSubDSInfo(GDALGetSubdatasetInfo(pszTileName));
675 18 : if (oSubDSInfo && !oSubDSInfo->GetPathComponent().empty())
676 : {
677 4 : const std::string osPath(oSubDSInfo->GetPathComponent());
678 2 : osRet = CPLIsFilenameRelative(osPath.c_str())
679 5 : ? oSubDSInfo->ModifyPathComponent(
680 4 : CPLProjectRelativeFilenameSafe(
681 3 : CPLGetPathSafe(pszVRTName).c_str(),
682 : osPath.c_str()))
683 2 : : std::string(pszTileName);
684 2 : GDALDestroySubdatasetInfo(oSubDSInfo);
685 2 : return osRet;
686 : }
687 :
688 : std::string osRelativeMadeAbsolute = CPLProjectRelativeFilenameSafe(
689 16 : CPLGetPathSafe(pszVRTName).c_str(), pszTileName);
690 : VSIStatBufL sStat;
691 16 : if (VSIStatL(osRelativeMadeAbsolute.c_str(), &sStat) == 0)
692 9 : return osRelativeMadeAbsolute;
693 : }
694 596 : return pszTileName;
695 : }
696 :
697 : /************************************************************************/
698 : /* GTIDoPaletteExpansionIfNeeded() */
699 : /************************************************************************/
700 :
701 : //! Do palette -> RGB(A) expansion
702 : static bool
703 479 : GTIDoPaletteExpansionIfNeeded(std::shared_ptr<GDALDataset> &poTileDS,
704 : int nBandCount)
705 : {
706 479 : bool bRet = true;
707 756 : if (poTileDS->GetRasterCount() == 1 &&
708 758 : (nBandCount == 3 || nBandCount == 4) &&
709 4 : poTileDS->GetRasterBand(1)->GetColorTable() != nullptr)
710 : {
711 :
712 8 : CPLStringList aosOptions;
713 4 : aosOptions.AddString("-of");
714 4 : aosOptions.AddString("VRT");
715 :
716 4 : aosOptions.AddString("-expand");
717 4 : aosOptions.AddString(nBandCount == 3 ? "rgb" : "rgba");
718 :
719 : GDALTranslateOptions *psOptions =
720 4 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
721 4 : int bUsageError = false;
722 : auto poRGBDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(
723 : GDALTranslate("", GDALDataset::ToHandle(poTileDS.get()), psOptions,
724 8 : &bUsageError)));
725 4 : GDALTranslateOptionsFree(psOptions);
726 4 : bRet = poRGBDS != nullptr;
727 4 : poTileDS = std::move(poRGBDS);
728 : }
729 479 : return bRet;
730 : }
731 :
732 : /************************************************************************/
733 : /* Open() */
734 : /************************************************************************/
735 :
736 286 : bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
737 : {
738 286 : eAccess = poOpenInfo->eAccess;
739 :
740 286 : CPLXMLNode *psRoot = nullptr;
741 286 : const char *pszIndexDataset = poOpenInfo->pszFilename;
742 :
743 286 : if (STARTS_WITH(poOpenInfo->pszFilename, GTI_PREFIX))
744 : {
745 14 : pszIndexDataset = poOpenInfo->pszFilename + strlen(GTI_PREFIX);
746 : }
747 272 : else if (STARTS_WITH(poOpenInfo->pszFilename, "<GDALTileIndexDataset"))
748 : {
749 : // CPLParseXMLString() emits an error in case of failure
750 25 : m_psXMLTree.reset(CPLParseXMLString(poOpenInfo->pszFilename));
751 25 : if (m_psXMLTree == nullptr)
752 1 : return false;
753 : }
754 247 : else if (poOpenInfo->nHeaderBytes > 0 &&
755 247 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
756 : "<GDALTileIndexDataset"))
757 : {
758 : // CPLParseXMLFile() emits an error in case of failure
759 6 : m_psXMLTree.reset(CPLParseXMLFile(poOpenInfo->pszFilename));
760 6 : if (m_psXMLTree == nullptr)
761 1 : return false;
762 5 : m_bXMLUpdatable = (poOpenInfo->eAccess == GA_Update);
763 : }
764 :
765 284 : if (m_psXMLTree)
766 : {
767 29 : psRoot = CPLGetXMLNode(m_psXMLTree.get(), "=GDALTileIndexDataset");
768 29 : if (psRoot == nullptr)
769 : {
770 1 : CPLError(CE_Failure, CPLE_AppDefined,
771 : "Missing GDALTileIndexDataset root element.");
772 1 : return false;
773 : }
774 :
775 28 : pszIndexDataset = CPLGetXMLValue(psRoot, "IndexDataset", nullptr);
776 28 : if (!pszIndexDataset)
777 : {
778 1 : CPLError(CE_Failure, CPLE_AppDefined,
779 : "Missing IndexDataset element.");
780 1 : return false;
781 : }
782 : }
783 :
784 282 : if (ENDS_WITH_CI(pszIndexDataset, ".gti.gpkg") &&
785 526 : poOpenInfo->nHeaderBytes >= 100 &&
786 244 : STARTS_WITH(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
787 : "SQLite format 3"))
788 : {
789 240 : const char *const apszAllowedDrivers[] = {"GPKG", nullptr};
790 240 : m_poVectorDS.reset(GDALDataset::Open(
791 480 : std::string("GPKG:\"").append(pszIndexDataset).append("\"").c_str(),
792 240 : GDAL_OF_VECTOR | GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR |
793 240 : ((poOpenInfo->nOpenFlags & GDAL_OF_UPDATE) ? GDAL_OF_UPDATE
794 : : GDAL_OF_READONLY),
795 : apszAllowedDrivers));
796 240 : if (!m_poVectorDS)
797 : {
798 1 : return false;
799 : }
800 242 : if (m_poVectorDS->GetLayerCount() == 0 &&
801 2 : (m_poVectorDS->GetRasterCount() != 0 ||
802 1 : m_poVectorDS->GetMetadata("SUBDATASETS") != nullptr))
803 : {
804 1 : return false;
805 : }
806 : }
807 : else
808 : {
809 42 : m_poVectorDS.reset(GDALDataset::Open(
810 42 : pszIndexDataset, GDAL_OF_VECTOR | GDAL_OF_VERBOSE_ERROR |
811 42 : ((poOpenInfo->nOpenFlags & GDAL_OF_UPDATE)
812 42 : ? GDAL_OF_UPDATE
813 : : GDAL_OF_READONLY)));
814 42 : if (!m_poVectorDS)
815 : {
816 1 : return false;
817 : }
818 : }
819 :
820 280 : if (m_poVectorDS->GetLayerCount() == 0)
821 : {
822 1 : CPLError(CE_Failure, CPLE_AppDefined, "%s has no vector layer",
823 : poOpenInfo->pszFilename);
824 1 : return false;
825 : }
826 :
827 279 : double dfOvrFactor = 1.0;
828 279 : if (const char *pszFactor =
829 279 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "FACTOR"))
830 : {
831 5 : dfOvrFactor = CPLAtof(pszFactor);
832 5 : if (!(dfOvrFactor > 1.0))
833 : {
834 1 : CPLError(CE_Failure, CPLE_AppDefined, "Wrong overview factor");
835 1 : return false;
836 : }
837 : }
838 :
839 278 : m_osSQL = CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "SQL", "");
840 278 : if (m_osSQL.empty())
841 : {
842 275 : if (!psRoot)
843 : {
844 249 : if (const char *pszVal =
845 249 : m_poVectorDS->GetMetadataItem(MD_DS_TILE_INDEX_SQL))
846 0 : m_osSQL = pszVal;
847 : }
848 : else
849 26 : m_osSQL = CPLGetXMLValue(psRoot, "SQL", "");
850 : }
851 :
852 278 : if (!m_osSQL.empty())
853 : {
854 5 : m_osSpatialSQL = CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
855 5 : "SPATIAL_SQL", "");
856 5 : if (m_osSpatialSQL.empty())
857 : {
858 3 : if (!psRoot)
859 : {
860 2 : if (const char *pszVal = m_poVectorDS->GetMetadataItem(
861 1 : MD_DS_TILE_INDEX_SPATIAL_SQL))
862 0 : m_osSpatialSQL = pszVal;
863 : }
864 : else
865 2 : m_osSpatialSQL = CPLGetXMLValue(psRoot, "SpatialSQL", "");
866 : }
867 : }
868 :
869 : const char *pszLayerName;
870 :
871 278 : if ((pszLayerName = CSLFetchNameValue(poOpenInfo->papszOpenOptions,
872 278 : "LAYER")) != nullptr)
873 : {
874 6 : m_poLayer = m_poVectorDS->GetLayerByName(pszLayerName);
875 6 : if (!m_poLayer)
876 : {
877 2 : CPLError(CE_Failure, CPLE_AppDefined, "Layer %s does not exist",
878 : pszLayerName);
879 2 : return false;
880 : }
881 : }
882 272 : else if (psRoot && (pszLayerName = CPLGetXMLValue(psRoot, "IndexLayer",
883 : nullptr)) != nullptr)
884 : {
885 8 : m_poLayer = m_poVectorDS->GetLayerByName(pszLayerName);
886 8 : if (!m_poLayer)
887 : {
888 1 : CPLError(CE_Failure, CPLE_AppDefined, "Layer %s does not exist",
889 : pszLayerName);
890 1 : return false;
891 : }
892 : }
893 514 : else if (!psRoot && (pszLayerName = m_poVectorDS->GetMetadataItem(
894 250 : MD_DS_TILE_INDEX_LAYER)) != nullptr)
895 : {
896 2 : m_poLayer = m_poVectorDS->GetLayerByName(pszLayerName);
897 2 : if (!m_poLayer)
898 : {
899 1 : CPLError(CE_Failure, CPLE_AppDefined, "Layer %s does not exist",
900 : pszLayerName);
901 1 : return false;
902 : }
903 : }
904 262 : else if (!m_osSQL.empty())
905 : {
906 5 : m_poLayer = m_poVectorDS->ExecuteSQL(m_osSQL.c_str(), nullptr, nullptr);
907 5 : if (!m_poLayer)
908 : {
909 1 : CPLError(CE_Failure, CPLE_AppDefined, "SQL request %s failed",
910 : m_osSQL.c_str());
911 1 : return false;
912 : }
913 4 : m_bIsSQLResultLayer = true;
914 : }
915 257 : else if (m_poVectorDS->GetLayerCount() == 1)
916 : {
917 255 : m_poLayer = m_poVectorDS->GetLayer(0);
918 255 : if (!m_poLayer)
919 : {
920 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot open layer 0");
921 0 : return false;
922 : }
923 : }
924 : else
925 : {
926 2 : if (STARTS_WITH(poOpenInfo->pszFilename, GTI_PREFIX))
927 : {
928 1 : CPLError(CE_Failure, CPLE_AppDefined,
929 : "%s has more than one layer. LAYER open option "
930 : "must be defined to specify which one to "
931 : "use as the tile index",
932 : pszIndexDataset);
933 : }
934 1 : else if (psRoot)
935 : {
936 0 : CPLError(CE_Failure, CPLE_AppDefined,
937 : "%s has more than one layer. IndexLayer element must be "
938 : "defined to specify which one to "
939 : "use as the tile index",
940 : pszIndexDataset);
941 : }
942 : else
943 : {
944 1 : CPLError(CE_Failure, CPLE_AppDefined,
945 : "%s has more than one layer. %s "
946 : "metadata item must be defined to specify which one to "
947 : "use as the tile index",
948 : pszIndexDataset, MD_DS_TILE_INDEX_LAYER);
949 : }
950 2 : return false;
951 : }
952 :
953 : // Try to get the metadata from an embedded xml:GTI domain
954 271 : if (!m_psXMLTree)
955 : {
956 247 : CSLConstList papszMD = m_poLayer->GetMetadata("xml:GTI");
957 247 : if (papszMD && papszMD[0])
958 : {
959 1 : m_psXMLTree.reset(CPLParseXMLString(papszMD[0]));
960 1 : if (m_psXMLTree == nullptr)
961 0 : return false;
962 :
963 1 : psRoot = CPLGetXMLNode(m_psXMLTree.get(), "=GDALTileIndexDataset");
964 1 : if (psRoot == nullptr)
965 : {
966 0 : CPLError(CE_Failure, CPLE_AppDefined,
967 : "Missing GDALTileIndexDataset root element.");
968 0 : return false;
969 : }
970 : }
971 : }
972 :
973 : // Get the value of an option.
974 : // The order of lookup is the following one (first to last):
975 : // - open options
976 : // - XML file
977 : // - Layer metadata items.
978 25424 : const auto GetOption = [poOpenInfo, psRoot, this](const char *pszItem)
979 : {
980 : const char *pszVal =
981 8105 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, pszItem);
982 8105 : if (pszVal)
983 10 : return pszVal;
984 :
985 8095 : if (psRoot)
986 : {
987 595 : pszVal = CPLGetXMLValue(psRoot, pszItem, nullptr);
988 595 : if (pszVal)
989 27 : return pszVal;
990 :
991 568 : if (EQUAL(pszItem, MD_BAND_COUNT))
992 22 : pszItem = GTI_XML_BANDCOUNT;
993 546 : else if (EQUAL(pszItem, MD_DATA_TYPE))
994 25 : pszItem = GTI_XML_DATATYPE;
995 521 : else if (EQUAL(pszItem, MD_NODATA))
996 20 : pszItem = GTI_XML_NODATAVALUE;
997 501 : else if (EQUAL(pszItem, MD_COLOR_INTERPRETATION))
998 25 : pszItem = GTI_XML_COLORINTERP;
999 476 : else if (EQUAL(pszItem, MD_LOCATION_FIELD))
1000 25 : pszItem = GTI_XML_LOCATIONFIELD;
1001 451 : else if (EQUAL(pszItem, MD_SORT_FIELD))
1002 25 : pszItem = GTI_XML_SORTFIELD;
1003 426 : else if (EQUAL(pszItem, MD_SORT_FIELD_ASC))
1004 2 : pszItem = GTI_XML_SORTFIELDASC;
1005 424 : else if (EQUAL(pszItem, MD_MASK_BAND))
1006 20 : pszItem = GTI_XML_MASKBAND;
1007 568 : pszVal = CPLGetXMLValue(psRoot, pszItem, nullptr);
1008 568 : if (pszVal)
1009 7 : return pszVal;
1010 : }
1011 :
1012 8061 : return m_poLayer->GetMetadataItem(pszItem);
1013 271 : };
1014 :
1015 271 : const char *pszFilter = GetOption("Filter");
1016 271 : if (pszFilter)
1017 : {
1018 1 : if (m_poLayer->SetAttributeFilter(pszFilter) != OGRERR_NONE)
1019 0 : return false;
1020 : }
1021 :
1022 271 : const OGRFeatureDefn *poLayerDefn = m_poLayer->GetLayerDefn();
1023 :
1024 542 : std::string osLocationFieldName;
1025 : {
1026 271 : const char *pszLocationFieldName = GetOption(MD_LOCATION_FIELD);
1027 271 : if (pszLocationFieldName)
1028 : {
1029 10 : osLocationFieldName = pszLocationFieldName;
1030 : }
1031 : else
1032 : {
1033 : // Is this a https://stac-utils.github.io/stac-geoparquet/latest/spec/stac-geoparquet-spec ?
1034 261 : if (poLayerDefn->GetFieldIndex("assets.data.href") >= 0)
1035 : {
1036 0 : m_bSTACCollection = true;
1037 0 : osLocationFieldName = "assets.data.href";
1038 0 : CPLDebug("GTI", "Using %s as location field",
1039 : osLocationFieldName.c_str());
1040 : }
1041 261 : else if (poLayerDefn->GetFieldIndex("assets.image.href") >= 0)
1042 : {
1043 3 : m_bSTACCollection = true;
1044 3 : osLocationFieldName = "assets.image.href";
1045 3 : CPLDebug("GTI", "Using %s as location field",
1046 : osLocationFieldName.c_str());
1047 : }
1048 512 : else if (poLayerDefn->GetFieldIndex("stac_version") >= 0 ||
1049 254 : poLayerDefn->GetFieldIndex("stac_extensions") >= 0)
1050 : {
1051 4 : m_bSTACCollection = true;
1052 4 : const int nFieldCount = poLayerDefn->GetFieldCount();
1053 : // Look for "assets.xxxxx.href" fields
1054 4 : int nAssetCount = 0;
1055 60 : for (int i = 0; i < nFieldCount; ++i)
1056 : {
1057 56 : const auto poFDefn = poLayerDefn->GetFieldDefn(i);
1058 56 : const char *pszFieldName = poFDefn->GetNameRef();
1059 56 : if (STARTS_WITH(pszFieldName, "assets.") &&
1060 44 : EQUAL(pszFieldName + strlen(pszFieldName) -
1061 : strlen(".href"),
1062 4 : ".href") &&
1063 : // Assets with "metadata" in them are very much likely
1064 : // not rasters... We could potentially confirm that by
1065 : // inspecting the value of the assets.XXX.type or
1066 : // assets.XXX.roles fields of one feature
1067 4 : !strstr(pszFieldName, "metadata"))
1068 : {
1069 4 : ++nAssetCount;
1070 4 : if (!osLocationFieldName.empty())
1071 : {
1072 0 : osLocationFieldName += ", ";
1073 : }
1074 4 : osLocationFieldName += pszFieldName;
1075 : }
1076 : }
1077 4 : if (nAssetCount > 1)
1078 : {
1079 0 : CPLError(CE_Failure, CPLE_AppDefined,
1080 : "Several potential STAC assets. Please select one "
1081 : "among %s with the LOCATION_FIELD open option",
1082 : osLocationFieldName.c_str());
1083 0 : return false;
1084 : }
1085 4 : else if (nAssetCount == 0)
1086 : {
1087 0 : CPLError(CE_Failure, CPLE_AppDefined,
1088 : "File has stac_version or stac_extensions "
1089 : "property but lacks assets");
1090 0 : return false;
1091 : }
1092 : }
1093 : else
1094 : {
1095 254 : constexpr const char *DEFAULT_LOCATION_FIELD_NAME = "location";
1096 254 : osLocationFieldName = DEFAULT_LOCATION_FIELD_NAME;
1097 : }
1098 : }
1099 : }
1100 :
1101 271 : m_nLocationFieldIndex =
1102 271 : poLayerDefn->GetFieldIndex(osLocationFieldName.c_str());
1103 271 : if (m_nLocationFieldIndex < 0)
1104 : {
1105 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find field %s",
1106 : osLocationFieldName.c_str());
1107 1 : return false;
1108 : }
1109 270 : if (poLayerDefn->GetFieldDefn(m_nLocationFieldIndex)->GetType() !=
1110 : OFTString)
1111 : {
1112 1 : CPLError(CE_Failure, CPLE_AppDefined, "Field %s is not of type string",
1113 : osLocationFieldName.c_str());
1114 1 : return false;
1115 : }
1116 :
1117 269 : const char *pszSortFieldName = GetOption(MD_SORT_FIELD);
1118 269 : if (pszSortFieldName)
1119 : {
1120 96 : m_nSortFieldIndex = poLayerDefn->GetFieldIndex(pszSortFieldName);
1121 96 : if (m_nSortFieldIndex < 0)
1122 : {
1123 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find field %s",
1124 : pszSortFieldName);
1125 1 : return false;
1126 : }
1127 :
1128 : const auto eFieldType =
1129 95 : poLayerDefn->GetFieldDefn(m_nSortFieldIndex)->GetType();
1130 95 : if (eFieldType != OFTString && eFieldType != OFTInteger &&
1131 61 : eFieldType != OFTInteger64 && eFieldType != OFTReal &&
1132 19 : eFieldType != OFTDate && eFieldType != OFTDateTime)
1133 : {
1134 1 : CPLError(CE_Failure, CPLE_AppDefined,
1135 : "Unsupported type for field %s", pszSortFieldName);
1136 1 : return false;
1137 : }
1138 :
1139 94 : const char *pszSortFieldAsc = GetOption(MD_SORT_FIELD_ASC);
1140 94 : if (pszSortFieldAsc)
1141 : {
1142 3 : m_bSortFieldAsc = CPLTestBool(pszSortFieldAsc);
1143 : }
1144 : }
1145 :
1146 267 : const char *pszResX = GetOption(MD_RESX);
1147 267 : const char *pszResY = GetOption(MD_RESY);
1148 267 : if (pszResX && !pszResY)
1149 : {
1150 1 : CPLError(CE_Failure, CPLE_AppDefined,
1151 : "%s metadata item defined, but not %s", MD_RESX, MD_RESY);
1152 1 : return false;
1153 : }
1154 266 : if (!pszResX && pszResY)
1155 : {
1156 1 : CPLError(CE_Failure, CPLE_AppDefined,
1157 : "%s metadata item defined, but not %s", MD_RESY, MD_RESX);
1158 1 : return false;
1159 : }
1160 :
1161 265 : const char *pszResampling = GetOption(MD_RESAMPLING);
1162 265 : if (pszResampling)
1163 : {
1164 8 : const auto nErrorCountBefore = CPLGetErrorCounter();
1165 8 : m_eResampling = GDALRasterIOGetResampleAlg(pszResampling);
1166 8 : if (nErrorCountBefore != CPLGetErrorCounter())
1167 : {
1168 0 : return false;
1169 : }
1170 8 : m_osResampling = pszResampling;
1171 : }
1172 :
1173 265 : const char *pszMinX = GetOption(MD_MINX);
1174 265 : const char *pszMinY = GetOption(MD_MINY);
1175 265 : const char *pszMaxX = GetOption(MD_MAXX);
1176 265 : const char *pszMaxY = GetOption(MD_MAXY);
1177 265 : int nCountMinMaxXY = (pszMinX ? 1 : 0) + (pszMinY ? 1 : 0) +
1178 265 : (pszMaxX ? 1 : 0) + (pszMaxY ? 1 : 0);
1179 265 : if (nCountMinMaxXY != 0 && nCountMinMaxXY != 4)
1180 : {
1181 4 : CPLError(CE_Failure, CPLE_AppDefined,
1182 : "None or all of %s, %s, %s and %s must be specified", MD_MINX,
1183 : MD_MINY, MD_MAXX, MD_MAXY);
1184 4 : return false;
1185 : }
1186 :
1187 261 : const char *pszXSize = GetOption(MD_XSIZE);
1188 261 : const char *pszYSize = GetOption(MD_YSIZE);
1189 261 : const char *pszGeoTransform = GetOption(MD_GEOTRANSFORM);
1190 261 : const int nCountXSizeYSizeGT =
1191 261 : (pszXSize ? 1 : 0) + (pszYSize ? 1 : 0) + (pszGeoTransform ? 1 : 0);
1192 261 : if (nCountXSizeYSizeGT != 0 && nCountXSizeYSizeGT != 3)
1193 : {
1194 3 : CPLError(CE_Failure, CPLE_AppDefined,
1195 : "None or all of %s, %s, %s must be specified", MD_XSIZE,
1196 : MD_YSIZE, MD_GEOTRANSFORM);
1197 3 : return false;
1198 : }
1199 :
1200 258 : const char *pszDataType = GetOption(MD_DATA_TYPE);
1201 258 : const char *pszColorInterp = GetOption(MD_COLOR_INTERPRETATION);
1202 258 : int nBandCount = 0;
1203 516 : std::vector<GDALDataType> aeDataTypes;
1204 516 : std::vector<std::pair<bool, double>> aNoData;
1205 516 : std::vector<GDALColorInterp> aeColorInterp;
1206 :
1207 258 : const char *pszSRS = GetOption(MD_SRS);
1208 258 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1209 258 : if (pszSRS)
1210 : {
1211 2 : if (m_oSRS.SetFromUserInput(
1212 : pszSRS,
1213 2 : OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get()) !=
1214 : OGRERR_NONE)
1215 : {
1216 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid %s", MD_SRS);
1217 1 : return false;
1218 : }
1219 : }
1220 256 : else if (const auto poSRS = m_poLayer->GetSpatialRef())
1221 : {
1222 : // Ignore GPKG "Undefined geographic SRS" and "Undefined Cartesian SRS"
1223 132 : if (!STARTS_WITH(poSRS->GetName(), "Undefined "))
1224 132 : m_oSRS = *poSRS;
1225 : }
1226 :
1227 514 : std::vector<const CPLXMLNode *> apoXMLNodeBands;
1228 257 : if (psRoot)
1229 : {
1230 25 : int nExpectedBandNumber = 1;
1231 113 : for (const CPLXMLNode *psIter = psRoot->psChild; psIter;
1232 88 : psIter = psIter->psNext)
1233 : {
1234 91 : if (psIter->eType == CXT_Element &&
1235 91 : strcmp(psIter->pszValue, GTI_XML_BAND_ELEMENT) == 0)
1236 : {
1237 : const char *pszBand =
1238 8 : CPLGetXMLValue(psIter, GTI_XML_BAND_NUMBER, nullptr);
1239 8 : if (!pszBand)
1240 : {
1241 1 : CPLError(CE_Failure, CPLE_AppDefined,
1242 : "%s attribute missing on %s element",
1243 : GTI_XML_BAND_NUMBER, GTI_XML_BAND_ELEMENT);
1244 3 : return false;
1245 : }
1246 7 : const int nBand = atoi(pszBand);
1247 7 : if (nBand <= 0)
1248 : {
1249 1 : CPLError(CE_Failure, CPLE_AppDefined,
1250 : "Invalid band number");
1251 1 : return false;
1252 : }
1253 6 : if (nBand != nExpectedBandNumber)
1254 : {
1255 1 : CPLError(CE_Failure, CPLE_AppDefined,
1256 : "Invalid band number: found %d, expected %d",
1257 : nBand, nExpectedBandNumber);
1258 1 : return false;
1259 : }
1260 5 : apoXMLNodeBands.push_back(psIter);
1261 5 : ++nExpectedBandNumber;
1262 : }
1263 : }
1264 : }
1265 :
1266 254 : const char *pszBandCount = GetOption(MD_BAND_COUNT);
1267 254 : if (pszBandCount)
1268 24 : nBandCount = atoi(pszBandCount);
1269 :
1270 254 : if (!apoXMLNodeBands.empty())
1271 : {
1272 5 : if (!pszBandCount)
1273 4 : nBandCount = static_cast<int>(apoXMLNodeBands.size());
1274 1 : else if (nBandCount != static_cast<int>(apoXMLNodeBands.size()))
1275 : {
1276 1 : CPLError(CE_Failure, CPLE_AppDefined,
1277 : "Inconsistent %s with actual number of %s elements",
1278 : GTI_XML_BANDCOUNT, GTI_XML_BAND_ELEMENT);
1279 1 : return false;
1280 : }
1281 : }
1282 :
1283 : // Take into STAC GeoParquet proj:code / proj:epsg / proj:wkt2 / proj:projjson
1284 : // and proj:transform fields
1285 253 : std::unique_ptr<OGRFeature> poFeature;
1286 506 : std::string osResX, osResY, osMinX, osMinY, osMaxX, osMaxY;
1287 253 : int iProjCode = -1;
1288 253 : int iProjEPSG = -1;
1289 253 : int iProjWKT2 = -1;
1290 253 : int iProjPROJSON = -1;
1291 253 : int iProjTransform = -1;
1292 :
1293 : const bool bIsStacGeoParquet =
1294 260 : STARTS_WITH(osLocationFieldName.c_str(), "assets.") &&
1295 7 : EQUAL(osLocationFieldName.c_str() + osLocationFieldName.size() -
1296 : strlen(".href"),
1297 : ".href");
1298 506 : std::string osAssetName;
1299 253 : if (bIsStacGeoParquet)
1300 : {
1301 7 : m_bSTACCollection = true;
1302 14 : osAssetName = osLocationFieldName.substr(
1303 : strlen("assets."),
1304 14 : osLocationFieldName.size() - strlen("assets.") - strlen(".href"));
1305 : }
1306 :
1307 : const auto GetAssetFieldIndex =
1308 38 : [poLayerDefn, &osAssetName](const char *pszFieldName)
1309 : {
1310 21 : const int idx = poLayerDefn->GetFieldIndex(
1311 21 : CPLSPrintf("assets.%s.%s", osAssetName.c_str(), pszFieldName));
1312 21 : if (idx >= 0)
1313 4 : return idx;
1314 17 : return poLayerDefn->GetFieldIndex(pszFieldName);
1315 253 : };
1316 :
1317 7 : if (bIsStacGeoParquet && !pszSRS && !pszResX && !pszResY && !pszMinX &&
1318 7 : !pszMinY && !pszMaxX && !pszMaxY &&
1319 7 : ((iProjCode = GetAssetFieldIndex("proj:code")) >= 0 ||
1320 4 : (iProjEPSG = GetAssetFieldIndex("proj:epsg")) >= 0 ||
1321 2 : (iProjWKT2 = GetAssetFieldIndex("proj:wkt2")) >= 0 ||
1322 261 : (iProjPROJSON = GetAssetFieldIndex("proj:projjson")) >= 0) &&
1323 7 : ((iProjTransform = GetAssetFieldIndex("proj:transform")) >= 0))
1324 : {
1325 7 : poFeature.reset(m_poLayer->GetNextFeature());
1326 : const auto poProjTransformField =
1327 7 : poLayerDefn->GetFieldDefn(iProjTransform);
1328 7 : if (poFeature &&
1329 7 : ((iProjCode >= 0 && poFeature->IsFieldSet(iProjCode)) ||
1330 4 : (iProjEPSG >= 0 && poFeature->IsFieldSet(iProjEPSG)) ||
1331 2 : (iProjWKT2 >= 0 && poFeature->IsFieldSet(iProjWKT2)) ||
1332 8 : (iProjPROJSON >= 0 && poFeature->IsFieldSet(iProjPROJSON))) &&
1333 25 : poFeature->IsFieldSet(iProjTransform) &&
1334 11 : (poProjTransformField->GetType() == OFTRealList ||
1335 4 : poProjTransformField->GetType() == OFTIntegerList ||
1336 0 : poProjTransformField->GetType() == OFTInteger64List))
1337 : {
1338 14 : OGRSpatialReference oSTACSRS;
1339 7 : oSTACSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1340 :
1341 7 : if (iProjCode >= 0 && poFeature->IsFieldSet(iProjCode))
1342 3 : oSTACSRS.SetFromUserInput(
1343 : poFeature->GetFieldAsString(iProjCode),
1344 : OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get());
1345 :
1346 4 : else if (iProjEPSG >= 0 && poFeature->IsFieldSet(iProjEPSG))
1347 2 : oSTACSRS.importFromEPSG(
1348 : poFeature->GetFieldAsInteger(iProjEPSG));
1349 :
1350 2 : else if (iProjWKT2 >= 0 && poFeature->IsFieldSet(iProjWKT2))
1351 1 : oSTACSRS.SetFromUserInput(
1352 : poFeature->GetFieldAsString(iProjWKT2),
1353 : OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get());
1354 :
1355 1 : else if (iProjPROJSON >= 0 && poFeature->IsFieldSet(iProjPROJSON))
1356 1 : oSTACSRS.SetFromUserInput(
1357 : poFeature->GetFieldAsString(iProjPROJSON),
1358 : OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get());
1359 :
1360 7 : if (!oSTACSRS.IsEmpty())
1361 : {
1362 6 : int nTransformCount = 0;
1363 : // Note: different coefficient ordering than GDAL geotransform
1364 6 : double adfProjTransform[6] = {0, 0, 0, 0, 0, 0};
1365 6 : if (poProjTransformField->GetType() == OFTRealList)
1366 : {
1367 : const auto padfFeatureTransform =
1368 2 : poFeature->GetFieldAsDoubleList(iProjTransform,
1369 : &nTransformCount);
1370 2 : if (nTransformCount >= 6)
1371 2 : memcpy(adfProjTransform, padfFeatureTransform,
1372 : 6 * sizeof(double));
1373 : }
1374 4 : else if (poProjTransformField->GetType() == OFTInteger64List)
1375 : {
1376 : const auto paFeatureTransform =
1377 0 : poFeature->GetFieldAsInteger64List(iProjTransform,
1378 : &nTransformCount);
1379 0 : if (nTransformCount >= 6)
1380 : {
1381 0 : for (int i = 0; i < 6; ++i)
1382 0 : adfProjTransform[i] =
1383 0 : static_cast<double>(paFeatureTransform[i]);
1384 : }
1385 : }
1386 4 : else if (poProjTransformField->GetType() == OFTIntegerList)
1387 : {
1388 : const auto paFeatureTransform =
1389 4 : poFeature->GetFieldAsIntegerList(iProjTransform,
1390 : &nTransformCount);
1391 4 : if (nTransformCount >= 6)
1392 : {
1393 28 : for (int i = 0; i < 6; ++i)
1394 24 : adfProjTransform[i] = paFeatureTransform[i];
1395 : }
1396 : }
1397 6 : OGREnvelope sEnvelope;
1398 12 : if (nTransformCount >= 6 && m_poLayer->GetSpatialRef() &&
1399 6 : m_poLayer->GetExtent(&sEnvelope, /* bForce = */ true) ==
1400 : OGRERR_NONE)
1401 : {
1402 6 : const double dfResX = adfProjTransform[0];
1403 6 : osResX = CPLSPrintf("%.17g", dfResX);
1404 6 : const double dfResY = std::fabs(adfProjTransform[4]);
1405 6 : osResY = CPLSPrintf("%.17g", dfResY);
1406 :
1407 : auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
1408 : OGRCreateCoordinateTransformation(
1409 12 : m_poLayer->GetSpatialRef(), &oSTACSRS));
1410 : auto poInvCT = std::unique_ptr<OGRCoordinateTransformation>(
1411 12 : poCT ? poCT->GetInverse() : nullptr);
1412 6 : double dfOutMinX = 0;
1413 6 : double dfOutMinY = 0;
1414 6 : double dfOutMaxX = 0;
1415 6 : double dfOutMaxY = 0;
1416 12 : if (dfResX > 0 && dfResY > 0 && poCT && poInvCT &&
1417 12 : poCT->TransformBounds(sEnvelope.MinX, sEnvelope.MinY,
1418 : sEnvelope.MaxX, sEnvelope.MaxY,
1419 : &dfOutMinX, &dfOutMinY,
1420 6 : &dfOutMaxX, &dfOutMaxY, 21))
1421 : {
1422 6 : constexpr double EPSILON = 1e-3;
1423 : const bool bTileAlignedOnRes =
1424 6 : (fmod(std::fabs(adfProjTransform[3]), dfResX) <=
1425 12 : EPSILON * dfResX &&
1426 6 : fmod(std::fabs(adfProjTransform[5]), dfResY) <=
1427 6 : EPSILON * dfResY);
1428 :
1429 : osMinX = CPLSPrintf(
1430 : "%.17g",
1431 : !bTileAlignedOnRes
1432 : ? dfOutMinX
1433 6 : : std::floor(dfOutMinX / dfResX) * dfResX);
1434 : osMinY = CPLSPrintf(
1435 : "%.17g",
1436 : !bTileAlignedOnRes
1437 : ? dfOutMinY
1438 6 : : std::floor(dfOutMinY / dfResY) * dfResY);
1439 : osMaxX = CPLSPrintf(
1440 : "%.17g",
1441 : !bTileAlignedOnRes
1442 : ? dfOutMaxX
1443 6 : : std::ceil(dfOutMaxX / dfResX) * dfResX);
1444 : osMaxY = CPLSPrintf(
1445 : "%.17g",
1446 : !bTileAlignedOnRes
1447 : ? dfOutMaxY
1448 6 : : std::ceil(dfOutMaxY / dfResY) * dfResY);
1449 :
1450 6 : m_oSRS = std::move(oSTACSRS);
1451 6 : pszResX = osResX.c_str();
1452 6 : pszResY = osResY.c_str();
1453 6 : pszMinX = osMinX.c_str();
1454 6 : pszMinY = osMinY.c_str();
1455 6 : pszMaxX = osMaxX.c_str();
1456 6 : pszMaxY = osMaxY.c_str();
1457 6 : nCountMinMaxXY = 4;
1458 :
1459 6 : poFeature.reset();
1460 6 : m_poLayer->ResetReading();
1461 :
1462 : m_poWarpedLayerKeeper =
1463 6 : std::make_unique<OGRWarpedLayer>(
1464 0 : m_poLayer, /* iGeomField = */ 0,
1465 6 : /* bTakeOwnership = */ false, std::move(poCT),
1466 12 : std::move(poInvCT));
1467 6 : m_poLayer = m_poWarpedLayerKeeper.get();
1468 6 : poLayerDefn = m_poLayer->GetLayerDefn();
1469 : }
1470 : }
1471 : }
1472 : }
1473 : }
1474 :
1475 253 : OGREnvelope sEnvelope;
1476 253 : if (nCountMinMaxXY == 4)
1477 : {
1478 17 : sEnvelope.MinX = CPLAtof(pszMinX);
1479 17 : sEnvelope.MinY = CPLAtof(pszMinY);
1480 17 : sEnvelope.MaxX = CPLAtof(pszMaxX);
1481 17 : sEnvelope.MaxY = CPLAtof(pszMaxY);
1482 17 : if (!(sEnvelope.MaxX > sEnvelope.MinX))
1483 : {
1484 1 : CPLError(CE_Failure, CPLE_AppDefined,
1485 : "%s metadata item must be > %s", MD_MAXX, MD_MINX);
1486 1 : return false;
1487 : }
1488 16 : if (!(sEnvelope.MaxY > sEnvelope.MinY))
1489 : {
1490 1 : CPLError(CE_Failure, CPLE_AppDefined,
1491 : "%s metadata item must be > %s", MD_MAXY, MD_MINY);
1492 1 : return false;
1493 : }
1494 : }
1495 :
1496 251 : bool bHasMaskBand = false;
1497 251 : std::unique_ptr<GDALColorTable> poSingleColorTable;
1498 278 : if ((!pszBandCount && apoXMLNodeBands.empty()) ||
1499 27 : (!(pszResX && pszResY) && nCountXSizeYSizeGT == 0))
1500 : {
1501 238 : CPLDebug("GTI", "Inspecting one feature due to missing metadata items");
1502 238 : m_bScannedOneFeatureAtOpening = true;
1503 :
1504 238 : if (!poFeature)
1505 237 : poFeature.reset(m_poLayer->GetNextFeature());
1506 474 : if (!poFeature ||
1507 236 : !poFeature->IsFieldSetAndNotNull(m_nLocationFieldIndex))
1508 : {
1509 2 : CPLError(
1510 : CE_Failure, CPLE_AppDefined,
1511 : "BAND_COUNT(+DATA_TYPE+COLOR_INTERPRETATION)+ (RESX+RESY or "
1512 : "XSIZE+YSIZE+GEOTRANSFORM) metadata items "
1513 : "missing");
1514 10 : return false;
1515 : }
1516 :
1517 : const char *pszTileName =
1518 236 : poFeature->GetFieldAsString(m_nLocationFieldIndex);
1519 : const std::string osTileName(GetAbsoluteFileName(
1520 236 : pszTileName, poOpenInfo->pszFilename, m_bSTACCollection));
1521 236 : pszTileName = osTileName.c_str();
1522 :
1523 : auto poTileDS = std::shared_ptr<GDALDataset>(
1524 : GDALDataset::Open(pszTileName,
1525 : GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR),
1526 236 : GDALDatasetUniquePtrReleaser());
1527 236 : if (!poTileDS)
1528 : {
1529 1 : return false;
1530 : }
1531 :
1532 : // do palette -> RGB(A) expansion if needed
1533 235 : if (!GTIDoPaletteExpansionIfNeeded(poTileDS, nBandCount))
1534 0 : return false;
1535 :
1536 235 : const int nTileBandCount = poTileDS->GetRasterCount();
1537 544 : for (int i = 0; i < nTileBandCount; ++i)
1538 : {
1539 309 : auto poTileBand = poTileDS->GetRasterBand(i + 1);
1540 309 : aeDataTypes.push_back(poTileBand->GetRasterDataType());
1541 309 : int bHasNoData = FALSE;
1542 309 : const double dfNoData = poTileBand->GetNoDataValue(&bHasNoData);
1543 309 : aNoData.emplace_back(CPL_TO_BOOL(bHasNoData), dfNoData);
1544 309 : aeColorInterp.push_back(poTileBand->GetColorInterpretation());
1545 309 : if (nTileBandCount == 1)
1546 : {
1547 200 : if (auto poCT = poTileBand->GetColorTable())
1548 : {
1549 : // We assume that this will apply to all tiles...
1550 : // TODO: detect if that it is really the case, and warn
1551 : // if not, or do approximate palette matching like
1552 : // done in GDALRasterBand::GetIndexColorTranslationTo()
1553 0 : poSingleColorTable.reset(poCT->Clone());
1554 : }
1555 : }
1556 :
1557 309 : if (poTileBand->GetMaskFlags() == GMF_PER_DATASET)
1558 1 : bHasMaskBand = true;
1559 : }
1560 235 : if (!pszBandCount && nBandCount == 0)
1561 221 : nBandCount = nTileBandCount;
1562 :
1563 235 : auto poTileSRS = poTileDS->GetSpatialRef();
1564 235 : if (!m_oSRS.IsEmpty() && poTileSRS && !m_oSRS.IsSame(poTileSRS))
1565 : {
1566 11 : CPLStringList aosOptions;
1567 11 : aosOptions.AddString("-of");
1568 11 : aosOptions.AddString("VRT");
1569 :
1570 11 : char *pszWKT = nullptr;
1571 11 : const char *const apszWKTOptions[] = {"FORMAT=WKT2_2019", nullptr};
1572 11 : m_oSRS.exportToWkt(&pszWKT, apszWKTOptions);
1573 11 : if (pszWKT)
1574 11 : m_osWKT = pszWKT;
1575 11 : CPLFree(pszWKT);
1576 :
1577 11 : if (m_osWKT.empty())
1578 : {
1579 0 : CPLError(CE_Failure, CPLE_AppDefined,
1580 : "Cannot export VRT SRS to WKT2");
1581 0 : return false;
1582 : }
1583 :
1584 11 : aosOptions.AddString("-t_srs");
1585 11 : aosOptions.AddString(m_osWKT.c_str());
1586 :
1587 : GDALWarpAppOptions *psWarpOptions =
1588 11 : GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
1589 11 : GDALDatasetH ahSrcDS[] = {GDALDataset::ToHandle(poTileDS.get())};
1590 11 : int bUsageError = false;
1591 : auto poWarpDS =
1592 : std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(GDALWarp(
1593 11 : "", nullptr, 1, ahSrcDS, psWarpOptions, &bUsageError)));
1594 11 : GDALWarpAppOptionsFree(psWarpOptions);
1595 11 : if (!poWarpDS)
1596 : {
1597 0 : return false;
1598 : }
1599 :
1600 11 : poTileDS = std::move(poWarpDS);
1601 11 : poTileSRS = poTileDS->GetSpatialRef();
1602 11 : CPL_IGNORE_RET_VAL(poTileSRS);
1603 : }
1604 :
1605 235 : GDALGeoTransform gtTile;
1606 235 : if (poTileDS->GetGeoTransform(gtTile) != CE_None)
1607 : {
1608 1 : CPLError(CE_Failure, CPLE_AppDefined,
1609 : "Cannot find geotransform on %s", pszTileName);
1610 1 : return false;
1611 : }
1612 234 : if (!(gtTile.xrot == 0))
1613 : {
1614 1 : CPLError(CE_Failure, CPLE_AppDefined,
1615 : "3rd value of GeoTransform of %s must be 0", pszTileName);
1616 1 : return false;
1617 : }
1618 233 : if (!(gtTile.yrot == 0))
1619 : {
1620 1 : CPLError(CE_Failure, CPLE_AppDefined,
1621 : "5th value of GeoTransform of %s must be 0", pszTileName);
1622 1 : return false;
1623 : }
1624 :
1625 232 : const double dfResX = gtTile.xscale;
1626 232 : const double dfResY = gtTile.yscale;
1627 232 : if (!(dfResX > 0))
1628 : {
1629 1 : CPLError(CE_Failure, CPLE_AppDefined,
1630 : "2nd value of GeoTransform of %s must be > 0",
1631 : pszTileName);
1632 1 : return false;
1633 : }
1634 231 : if (!(dfResY != 0))
1635 : {
1636 0 : CPLError(CE_Failure, CPLE_AppDefined,
1637 : "6th value of GeoTransform of %s must be != 0",
1638 : pszTileName);
1639 0 : return false;
1640 : }
1641 :
1642 450 : if (!sEnvelope.IsInit() &&
1643 219 : m_poLayer->GetExtent(&sEnvelope, /* bForce = */ false) ==
1644 : OGRERR_FAILURE)
1645 : {
1646 2 : if (m_poLayer->GetExtent(&sEnvelope, /* bForce = */ true) ==
1647 : OGRERR_FAILURE)
1648 : {
1649 1 : CPLError(CE_Failure, CPLE_AppDefined,
1650 : "Cannot get layer extent");
1651 1 : return false;
1652 : }
1653 1 : CPLError(CE_Warning, CPLE_AppDefined,
1654 : "Could get layer extent, but using a slower method");
1655 : }
1656 :
1657 230 : const double dfXSize = (sEnvelope.MaxX - sEnvelope.MinX) / dfResX;
1658 230 : if (!(dfXSize >= 0 && dfXSize < INT_MAX))
1659 : {
1660 1 : CPLError(CE_Failure, CPLE_AppDefined,
1661 : "Too small %s, or wrong layer extent", MD_RESX);
1662 1 : return false;
1663 : }
1664 :
1665 229 : const double dfYSize =
1666 229 : (sEnvelope.MaxY - sEnvelope.MinY) / std::fabs(dfResY);
1667 229 : if (!(dfYSize >= 0 && dfYSize < INT_MAX))
1668 : {
1669 1 : CPLError(CE_Failure, CPLE_AppDefined,
1670 : "Too small %s, or wrong layer extent", MD_RESY);
1671 1 : return false;
1672 : }
1673 :
1674 228 : m_gt.xorig = sEnvelope.MinX;
1675 228 : m_gt.xscale = dfResX;
1676 228 : m_gt.xrot = 0;
1677 228 : m_gt.yorig = sEnvelope.MaxY;
1678 228 : m_gt.yrot = 0;
1679 228 : m_gt.yscale = -std::fabs(dfResY);
1680 :
1681 228 : nRasterXSize = static_cast<int>(std::ceil(dfXSize));
1682 228 : nRasterYSize = static_cast<int>(std::ceil(dfYSize));
1683 : }
1684 :
1685 241 : if (pszXSize && pszYSize && pszGeoTransform)
1686 : {
1687 12 : const int nXSize = atoi(pszXSize);
1688 12 : if (nXSize <= 0)
1689 : {
1690 1 : CPLError(CE_Failure, CPLE_AppDefined,
1691 : "%s metadata item must be > 0", MD_XSIZE);
1692 6 : return false;
1693 : }
1694 :
1695 11 : const int nYSize = atoi(pszYSize);
1696 11 : if (nYSize <= 0)
1697 : {
1698 1 : CPLError(CE_Failure, CPLE_AppDefined,
1699 : "%s metadata item must be > 0", MD_YSIZE);
1700 1 : return false;
1701 : }
1702 :
1703 : const CPLStringList aosTokens(
1704 10 : CSLTokenizeString2(pszGeoTransform, ",", 0));
1705 10 : if (aosTokens.size() != 6)
1706 : {
1707 1 : CPLError(CE_Failure, CPLE_AppDefined,
1708 : "%s metadata item must be 6 numeric values "
1709 : "separated with comma",
1710 : MD_GEOTRANSFORM);
1711 1 : return false;
1712 : }
1713 63 : for (int i = 0; i < 6; ++i)
1714 : {
1715 54 : m_gt[i] = CPLAtof(aosTokens[i]);
1716 : }
1717 9 : if (!(m_gt.xscale > 0))
1718 : {
1719 0 : CPLError(CE_Failure, CPLE_AppDefined, "2nd value of %s must be > 0",
1720 : MD_GEOTRANSFORM);
1721 0 : return false;
1722 : }
1723 9 : if (!(m_gt.xrot == 0))
1724 : {
1725 1 : CPLError(CE_Failure, CPLE_AppDefined, "3rd value of %s must be 0",
1726 : MD_GEOTRANSFORM);
1727 1 : return false;
1728 : }
1729 8 : if (!(m_gt.yrot == 0))
1730 : {
1731 1 : CPLError(CE_Failure, CPLE_AppDefined, "5th value of %s must be 0",
1732 : MD_GEOTRANSFORM);
1733 1 : return false;
1734 : }
1735 7 : if (!(m_gt.yscale < 0))
1736 : {
1737 1 : CPLError(CE_Failure, CPLE_AppDefined, "6th value of %s must be < 0",
1738 : MD_GEOTRANSFORM);
1739 1 : return false;
1740 : }
1741 6 : nRasterXSize = nXSize;
1742 12 : nRasterYSize = nYSize;
1743 : }
1744 229 : else if (pszResX && pszResY)
1745 : {
1746 22 : const double dfResX = CPLAtof(pszResX);
1747 22 : if (!(dfResX > 0))
1748 : {
1749 1 : CPLError(CE_Failure, CPLE_AppDefined,
1750 : "RESX metadata item must be > 0");
1751 1 : return false;
1752 : }
1753 21 : const double dfResY = CPLAtof(pszResY);
1754 21 : if (!(dfResY > 0))
1755 : {
1756 1 : CPLError(CE_Failure, CPLE_AppDefined,
1757 : "RESY metadata item must be > 0");
1758 1 : return false;
1759 : }
1760 :
1761 20 : if (nCountMinMaxXY == 4)
1762 : {
1763 11 : if (pszXSize || pszYSize || pszGeoTransform)
1764 : {
1765 0 : CPLError(CE_Warning, CPLE_AppDefined,
1766 : "Ignoring %s, %s and %s when %s, "
1767 : "%s, %s and %s are specified",
1768 : MD_XSIZE, MD_YSIZE, MD_GEOTRANSFORM, MD_MINX, MD_MINY,
1769 : MD_MAXX, MD_MAXY);
1770 : }
1771 : }
1772 13 : else if (!sEnvelope.IsInit() &&
1773 4 : m_poLayer->GetExtent(&sEnvelope, /* bForce = */ false) ==
1774 : OGRERR_FAILURE)
1775 : {
1776 0 : if (m_poLayer->GetExtent(&sEnvelope, /* bForce = */ true) ==
1777 : OGRERR_FAILURE)
1778 : {
1779 0 : CPLError(CE_Failure, CPLE_AppDefined,
1780 : "Cannot get layer extent");
1781 0 : return false;
1782 : }
1783 0 : CPLError(CE_Warning, CPLE_AppDefined,
1784 : "Could get layer extent, but using a slower method");
1785 : }
1786 :
1787 20 : const double dfXSize = (sEnvelope.MaxX - sEnvelope.MinX) / dfResX;
1788 20 : if (!(dfXSize >= 0 && dfXSize < INT_MAX))
1789 : {
1790 1 : CPLError(CE_Failure, CPLE_AppDefined,
1791 : "Too small %s, or wrong layer extent", MD_RESX);
1792 1 : return false;
1793 : }
1794 :
1795 19 : const double dfYSize = (sEnvelope.MaxY - sEnvelope.MinY) / dfResY;
1796 19 : if (!(dfYSize >= 0 && dfYSize < INT_MAX))
1797 : {
1798 1 : CPLError(CE_Failure, CPLE_AppDefined,
1799 : "Too small %s, or wrong layer extent", MD_RESY);
1800 1 : return false;
1801 : }
1802 :
1803 18 : m_gt.xorig = sEnvelope.MinX;
1804 18 : m_gt.xscale = dfResX;
1805 18 : m_gt.xrot = 0;
1806 18 : m_gt.yorig = sEnvelope.MaxY;
1807 18 : m_gt.yrot = 0;
1808 18 : m_gt.yscale = -dfResY;
1809 18 : nRasterXSize = static_cast<int>(std::ceil(dfXSize));
1810 18 : nRasterYSize = static_cast<int>(std::ceil(dfYSize));
1811 : }
1812 :
1813 231 : if (nBandCount == 0 && !pszBandCount)
1814 : {
1815 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s metadata item missing",
1816 : MD_BAND_COUNT);
1817 0 : return false;
1818 : }
1819 :
1820 231 : if (!GDALCheckBandCount(nBandCount, false))
1821 1 : return false;
1822 :
1823 230 : if (aeDataTypes.empty() && !pszDataType)
1824 : {
1825 9 : aeDataTypes.resize(nBandCount, GDT_UInt8);
1826 : }
1827 221 : else if (pszDataType)
1828 : {
1829 10 : aeDataTypes.clear();
1830 10 : const CPLStringList aosTokens(CSLTokenizeString2(pszDataType, ", ", 0));
1831 10 : if (aosTokens.size() == 1)
1832 : {
1833 8 : const auto eDataType = GDALGetDataTypeByName(aosTokens[0]);
1834 8 : if (eDataType == GDT_Unknown)
1835 : {
1836 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for %s",
1837 : MD_DATA_TYPE);
1838 1 : return false;
1839 : }
1840 7 : aeDataTypes.resize(nBandCount, eDataType);
1841 : }
1842 2 : else if (aosTokens.size() == nBandCount)
1843 : {
1844 2 : for (int i = 0; i < nBandCount; ++i)
1845 : {
1846 2 : const auto eDataType = GDALGetDataTypeByName(aosTokens[i]);
1847 2 : if (eDataType == GDT_Unknown)
1848 : {
1849 1 : CPLError(CE_Failure, CPLE_AppDefined,
1850 : "Invalid value for %s", MD_DATA_TYPE);
1851 1 : return false;
1852 : }
1853 1 : aeDataTypes.push_back(eDataType);
1854 : }
1855 : }
1856 : else
1857 : {
1858 1 : CPLError(CE_Failure, CPLE_AppDefined,
1859 : "Number of values in %s must be 1 or %s", MD_DATA_TYPE,
1860 : MD_BAND_COUNT);
1861 1 : return false;
1862 : }
1863 : }
1864 :
1865 227 : const char *pszNoData = GetOption(MD_NODATA);
1866 227 : if (pszNoData)
1867 : {
1868 22 : const auto IsValidNoDataStr = [](const char *pszStr)
1869 : {
1870 22 : if (EQUAL(pszStr, "inf") || EQUAL(pszStr, "-inf") ||
1871 18 : EQUAL(pszStr, "nan"))
1872 6 : return true;
1873 16 : const auto eType = CPLGetValueType(pszStr);
1874 16 : return eType == CPL_VALUE_INTEGER || eType == CPL_VALUE_REAL;
1875 : };
1876 :
1877 20 : aNoData.clear();
1878 20 : const CPLStringList aosTokens(CSLTokenizeString2(pszNoData, ", ", 0));
1879 20 : if (aosTokens.size() == 1)
1880 : {
1881 16 : if (!EQUAL(aosTokens[0], "NONE"))
1882 : {
1883 13 : if (!IsValidNoDataStr(aosTokens[0]))
1884 : {
1885 1 : CPLError(CE_Failure, CPLE_AppDefined,
1886 : "Invalid value for %s", MD_NODATA);
1887 1 : return false;
1888 : }
1889 12 : aNoData.resize(nBandCount,
1890 24 : std::pair(true, CPLAtof(aosTokens[0])));
1891 : }
1892 : }
1893 4 : else if (aosTokens.size() == nBandCount)
1894 : {
1895 12 : for (int i = 0; i < nBandCount; ++i)
1896 : {
1897 10 : if (EQUAL(aosTokens[i], "NONE"))
1898 : {
1899 1 : aNoData.emplace_back(false, 0);
1900 : }
1901 9 : else if (IsValidNoDataStr(aosTokens[i]))
1902 : {
1903 8 : aNoData.emplace_back(true, CPLAtof(aosTokens[i]));
1904 : }
1905 : else
1906 : {
1907 1 : CPLError(CE_Failure, CPLE_AppDefined,
1908 : "Invalid value for %s", MD_NODATA);
1909 1 : return false;
1910 : }
1911 : }
1912 : }
1913 : else
1914 : {
1915 1 : CPLError(CE_Failure, CPLE_AppDefined,
1916 : "Number of values in %s must be 1 or %s", MD_NODATA,
1917 : MD_BAND_COUNT);
1918 1 : return false;
1919 : }
1920 : }
1921 :
1922 224 : if (pszColorInterp)
1923 : {
1924 11 : aeColorInterp.clear();
1925 : const CPLStringList aosTokens(
1926 11 : CSLTokenizeString2(pszColorInterp, ", ", 0));
1927 11 : if (aosTokens.size() == 1)
1928 : {
1929 7 : const auto eInterp = GDALGetColorInterpretationByName(aosTokens[0]);
1930 12 : if (eInterp == GCI_Undefined &&
1931 5 : !EQUAL(aosTokens[0],
1932 : GDALGetColorInterpretationName(GCI_Undefined)))
1933 : {
1934 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for %s",
1935 : MD_COLOR_INTERPRETATION);
1936 1 : return false;
1937 : }
1938 6 : aeColorInterp.resize(nBandCount, eInterp);
1939 : }
1940 4 : else if (aosTokens.size() == nBandCount)
1941 : {
1942 11 : for (int i = 0; i < nBandCount; ++i)
1943 : {
1944 : const auto eInterp =
1945 9 : GDALGetColorInterpretationByName(aosTokens[i]);
1946 11 : if (eInterp == GCI_Undefined &&
1947 2 : !EQUAL(aosTokens[i],
1948 : GDALGetColorInterpretationName(GCI_Undefined)))
1949 : {
1950 1 : CPLError(CE_Failure, CPLE_AppDefined,
1951 : "Invalid value for %s", MD_COLOR_INTERPRETATION);
1952 1 : return false;
1953 : }
1954 8 : aeColorInterp.emplace_back(eInterp);
1955 : }
1956 : }
1957 : else
1958 : {
1959 1 : CPLError(CE_Failure, CPLE_AppDefined,
1960 : "Number of values in %s must be 1 or "
1961 : "%s",
1962 : MD_COLOR_INTERPRETATION, MD_BAND_COUNT);
1963 1 : return false;
1964 : }
1965 : }
1966 :
1967 : /* -------------------------------------------------------------------- */
1968 : /* Create bands. */
1969 : /* -------------------------------------------------------------------- */
1970 221 : if (aeDataTypes.size() != static_cast<size_t>(nBandCount))
1971 : {
1972 1 : CPLError(
1973 : CE_Failure, CPLE_AppDefined,
1974 : "Number of data types values found not matching number of bands");
1975 1 : return false;
1976 : }
1977 220 : if (!aNoData.empty() && aNoData.size() != static_cast<size_t>(nBandCount))
1978 : {
1979 1 : CPLError(CE_Failure, CPLE_AppDefined,
1980 : "Number of nodata values found not matching number of bands");
1981 1 : return false;
1982 : }
1983 428 : if (!aeColorInterp.empty() &&
1984 209 : aeColorInterp.size() != static_cast<size_t>(nBandCount))
1985 : {
1986 1 : CPLError(CE_Failure, CPLE_AppDefined,
1987 : "Number of color interpretation values found not matching "
1988 : "number of bands");
1989 1 : return false;
1990 : }
1991 :
1992 218 : int nBlockXSize = 256;
1993 218 : const char *pszBlockXSize = GetOption(MD_BLOCK_X_SIZE);
1994 218 : if (pszBlockXSize)
1995 : {
1996 3 : nBlockXSize = atoi(pszBlockXSize);
1997 3 : if (nBlockXSize <= 0)
1998 : {
1999 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid %s",
2000 : MD_BLOCK_X_SIZE);
2001 1 : return false;
2002 : }
2003 : }
2004 :
2005 217 : int nBlockYSize = 256;
2006 217 : const char *pszBlockYSize = GetOption(MD_BLOCK_Y_SIZE);
2007 217 : if (pszBlockYSize)
2008 : {
2009 3 : nBlockYSize = atoi(pszBlockYSize);
2010 3 : if (nBlockYSize <= 0)
2011 : {
2012 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid %s",
2013 : MD_BLOCK_Y_SIZE);
2014 1 : return false;
2015 : }
2016 : }
2017 :
2018 216 : if (nBlockXSize > INT_MAX / nBlockYSize)
2019 : {
2020 1 : CPLError(CE_Failure, CPLE_AppDefined, "Too big %s * %s",
2021 : MD_BLOCK_X_SIZE, MD_BLOCK_Y_SIZE);
2022 1 : return false;
2023 : }
2024 :
2025 215 : if (dfOvrFactor > 1.0)
2026 : {
2027 4 : m_gt.xscale *= dfOvrFactor;
2028 4 : m_gt.yscale *= dfOvrFactor;
2029 4 : nRasterXSize = static_cast<int>(std::ceil(nRasterXSize / dfOvrFactor));
2030 4 : nRasterYSize = static_cast<int>(std::ceil(nRasterYSize / dfOvrFactor));
2031 : }
2032 :
2033 430 : std::vector<std::string> aosDescriptions;
2034 430 : std::vector<double> adfCenterWavelength;
2035 430 : std::vector<double> adfFullWidthHalfMax;
2036 430 : std::vector<double> adfScale;
2037 430 : std::vector<double> adfOffset;
2038 215 : if (bIsStacGeoParquet && poFeature)
2039 : {
2040 7 : int nBandsIdx = poLayerDefn->GetFieldIndex(
2041 7 : CPLSPrintf("assets.%s.eo:bands", osAssetName.c_str()));
2042 7 : if (nBandsIdx < 0)
2043 2 : nBandsIdx = poLayerDefn->GetFieldIndex("bands");
2044 7 : if (nBandsIdx >= 0 &&
2045 14 : poLayerDefn->GetFieldDefn(nBandsIdx)->GetSubType() == OFSTJSON &&
2046 7 : poFeature->IsFieldSet(nBandsIdx))
2047 : {
2048 7 : const char *pszStr = poFeature->GetFieldAsString(nBandsIdx);
2049 14 : CPLJSONDocument oDoc;
2050 21 : if (oDoc.LoadMemory(pszStr) &&
2051 14 : oDoc.GetRoot().GetType() == CPLJSONObject::Type::Array)
2052 : {
2053 14 : const auto oArray = oDoc.GetRoot().ToArray();
2054 7 : if (oArray.Size() == nBandCount)
2055 : {
2056 7 : int i = 0;
2057 7 : aosDescriptions.resize(nBandCount);
2058 7 : adfCenterWavelength.resize(nBandCount);
2059 7 : adfFullWidthHalfMax.resize(nBandCount);
2060 19 : for (const auto &oObj : oArray)
2061 : {
2062 12 : if (oObj.GetType() == CPLJSONObject::Type::Object)
2063 : {
2064 36 : auto osCommonName = oObj.GetString("common_name");
2065 12 : if (osCommonName.empty())
2066 4 : osCommonName = oObj.GetString("eo:common_name");
2067 : const auto eInterp =
2068 12 : GDALGetColorInterpFromSTACCommonName(
2069 : osCommonName.c_str());
2070 12 : if (eInterp != GCI_Undefined)
2071 11 : aeColorInterp[i] = eInterp;
2072 :
2073 12 : aosDescriptions[i] = oObj.GetString("name");
2074 :
2075 : std::string osDescription =
2076 36 : oObj.GetString("description");
2077 12 : if (!osDescription.empty())
2078 : {
2079 1 : if (aosDescriptions[i].empty())
2080 0 : aosDescriptions[i] =
2081 0 : std::move(osDescription);
2082 : else
2083 1 : aosDescriptions[i]
2084 1 : .append(" (")
2085 1 : .append(osDescription)
2086 1 : .append(")");
2087 : }
2088 :
2089 24 : adfCenterWavelength[i] =
2090 12 : oObj.GetDouble("center_wavelength");
2091 12 : if (adfCenterWavelength[i] == 0)
2092 16 : adfCenterWavelength[i] =
2093 8 : oObj.GetDouble("eo:center_wavelength");
2094 24 : adfFullWidthHalfMax[i] =
2095 12 : oObj.GetDouble("full_width_half_max");
2096 12 : if (adfFullWidthHalfMax[i] == 0)
2097 16 : adfFullWidthHalfMax[i] =
2098 8 : oObj.GetDouble("eo:full_width_half_max");
2099 : }
2100 12 : ++i;
2101 : }
2102 : }
2103 : }
2104 : }
2105 :
2106 7 : int nRasterBandsIdx = poLayerDefn->GetFieldIndex(
2107 7 : CPLSPrintf("assets.%s.raster:bands", osAssetName.c_str()));
2108 7 : if (nRasterBandsIdx < 0)
2109 3 : nRasterBandsIdx = poLayerDefn->GetFieldIndex("bands");
2110 6 : if (nRasterBandsIdx >= 0 &&
2111 6 : poLayerDefn->GetFieldDefn(nRasterBandsIdx)->GetSubType() ==
2112 13 : OFSTJSON &&
2113 6 : poFeature->IsFieldSet(nRasterBandsIdx))
2114 : {
2115 6 : const char *pszStr = poFeature->GetFieldAsString(nRasterBandsIdx);
2116 12 : CPLJSONDocument oDoc;
2117 18 : if (oDoc.LoadMemory(pszStr) &&
2118 12 : oDoc.GetRoot().GetType() == CPLJSONObject::Type::Array)
2119 : {
2120 12 : const auto oArray = oDoc.GetRoot().ToArray();
2121 6 : if (oArray.Size() == nBandCount)
2122 : {
2123 6 : int i = 0;
2124 6 : adfScale.resize(nBandCount,
2125 6 : std::numeric_limits<double>::quiet_NaN());
2126 6 : adfOffset.resize(nBandCount,
2127 6 : std::numeric_limits<double>::quiet_NaN());
2128 14 : for (const auto &oObj : oArray)
2129 : {
2130 8 : if (oObj.GetType() == CPLJSONObject::Type::Object)
2131 : {
2132 8 : adfScale[i] = oObj.GetDouble(
2133 : "scale",
2134 : std::numeric_limits<double>::quiet_NaN());
2135 8 : adfOffset[i] = oObj.GetDouble(
2136 : "offset",
2137 : std::numeric_limits<double>::quiet_NaN());
2138 : }
2139 8 : ++i;
2140 : }
2141 : }
2142 : }
2143 : }
2144 : }
2145 :
2146 215 : GDALTileIndexBand *poFirstBand = nullptr;
2147 512 : for (int i = 0; i < nBandCount; ++i)
2148 : {
2149 297 : GDALDataType eDataType = aeDataTypes[i];
2150 297 : if (!apoXMLNodeBands.empty())
2151 : {
2152 4 : const char *pszVal = CPLGetXMLValue(apoXMLNodeBands[i],
2153 : GTI_XML_BAND_DATATYPE, nullptr);
2154 4 : if (pszVal)
2155 : {
2156 3 : eDataType = GDALGetDataTypeByName(pszVal);
2157 3 : if (eDataType == GDT_Unknown)
2158 0 : return false;
2159 : }
2160 : }
2161 : auto poBandUniquePtr = std::make_unique<GDALTileIndexBand>(
2162 594 : this, i + 1, eDataType, nBlockXSize, nBlockYSize);
2163 297 : auto poBand = poBandUniquePtr.get();
2164 297 : SetBand(i + 1, poBandUniquePtr.release());
2165 297 : if (!poFirstBand)
2166 215 : poFirstBand = poBand;
2167 297 : if (poBand->GetRasterDataType() != poFirstBand->GetRasterDataType())
2168 : {
2169 0 : m_bSameDataType = false;
2170 : }
2171 :
2172 297 : if (!aosDescriptions.empty() && !aosDescriptions[i].empty())
2173 : {
2174 12 : poBand->GDALRasterBand::SetDescription(aosDescriptions[i].c_str());
2175 : }
2176 297 : if (!apoXMLNodeBands.empty())
2177 : {
2178 4 : const char *pszVal = CPLGetXMLValue(
2179 4 : apoXMLNodeBands[i], GTI_XML_BAND_DESCRIPTION, nullptr);
2180 4 : if (pszVal)
2181 : {
2182 2 : poBand->GDALRasterBand::SetDescription(pszVal);
2183 : }
2184 : }
2185 :
2186 297 : if (!aNoData.empty() && aNoData[i].first)
2187 : {
2188 31 : poBand->m_bNoDataValueSet = true;
2189 31 : poBand->m_dfNoDataValue = aNoData[i].second;
2190 : }
2191 297 : if (!apoXMLNodeBands.empty())
2192 : {
2193 4 : const char *pszVal = CPLGetXMLValue(
2194 4 : apoXMLNodeBands[i], GTI_XML_BAND_NODATAVALUE, nullptr);
2195 4 : if (pszVal)
2196 : {
2197 3 : poBand->m_bNoDataValueSet = true;
2198 3 : poBand->m_dfNoDataValue = CPLAtof(pszVal);
2199 : }
2200 : }
2201 593 : if (poBand->m_bNoDataValueSet != poFirstBand->m_bNoDataValueSet ||
2202 296 : !IsSameNaNAware(poBand->m_dfNoDataValue,
2203 : poFirstBand->m_dfNoDataValue))
2204 : {
2205 6 : m_bSameNoData = false;
2206 : }
2207 :
2208 297 : if (!aeColorInterp.empty())
2209 : {
2210 284 : poBand->m_eColorInterp = aeColorInterp[i];
2211 : }
2212 297 : if (!apoXMLNodeBands.empty())
2213 : {
2214 4 : const char *pszVal = CPLGetXMLValue(
2215 4 : apoXMLNodeBands[i], GTI_XML_BAND_COLORINTERP, nullptr);
2216 4 : if (pszVal)
2217 : {
2218 4 : poBand->m_eColorInterp =
2219 4 : GDALGetColorInterpretationByName(pszVal);
2220 : }
2221 : }
2222 :
2223 305 : if (static_cast<int>(adfScale.size()) == nBandCount &&
2224 8 : !std::isnan(adfScale[i]))
2225 : {
2226 4 : poBand->m_dfScale = adfScale[i];
2227 : }
2228 297 : if (const char *pszScale =
2229 297 : GetOption(CPLSPrintf("BAND_%d_%s", i + 1, MD_BAND_SCALE)))
2230 : {
2231 6 : poBand->m_dfScale = CPLAtof(pszScale);
2232 : }
2233 297 : if (!apoXMLNodeBands.empty())
2234 : {
2235 : const char *pszVal =
2236 4 : CPLGetXMLValue(apoXMLNodeBands[i], GTI_XML_BAND_SCALE, nullptr);
2237 4 : if (pszVal)
2238 : {
2239 2 : poBand->m_dfScale = CPLAtof(pszVal);
2240 : }
2241 : }
2242 :
2243 305 : if (static_cast<int>(adfOffset.size()) == nBandCount &&
2244 8 : !std::isnan(adfOffset[i]))
2245 : {
2246 4 : poBand->m_dfOffset = adfOffset[i];
2247 : }
2248 297 : if (const char *pszOffset =
2249 297 : GetOption(CPLSPrintf("BAND_%d_%s", i + 1, MD_BAND_OFFSET)))
2250 : {
2251 6 : poBand->m_dfOffset = CPLAtof(pszOffset);
2252 : }
2253 297 : if (!apoXMLNodeBands.empty())
2254 : {
2255 4 : const char *pszVal = CPLGetXMLValue(apoXMLNodeBands[i],
2256 : GTI_XML_BAND_OFFSET, nullptr);
2257 4 : if (pszVal)
2258 : {
2259 2 : poBand->m_dfOffset = CPLAtof(pszVal);
2260 : }
2261 : }
2262 :
2263 297 : if (const char *pszUnit =
2264 297 : GetOption(CPLSPrintf("BAND_%d_%s", i + 1, MD_BAND_UNITTYPE)))
2265 : {
2266 6 : poBand->m_osUnit = pszUnit;
2267 : }
2268 297 : if (!apoXMLNodeBands.empty())
2269 : {
2270 4 : const char *pszVal = CPLGetXMLValue(apoXMLNodeBands[i],
2271 : GTI_XML_BAND_UNITTYPE, nullptr);
2272 4 : if (pszVal)
2273 : {
2274 2 : poBand->m_osUnit = pszVal;
2275 : }
2276 : }
2277 :
2278 297 : if (!apoXMLNodeBands.empty())
2279 : {
2280 4 : const CPLXMLNode *psBandNode = apoXMLNodeBands[i];
2281 4 : poBand->oMDMD.XMLInit(psBandNode, TRUE);
2282 :
2283 4 : if (const CPLXMLNode *psCategoryNames =
2284 4 : CPLGetXMLNode(psBandNode, GTI_XML_CATEGORYNAMES))
2285 : {
2286 : poBand->m_aosCategoryNames =
2287 2 : VRTParseCategoryNames(psCategoryNames);
2288 : }
2289 :
2290 4 : if (const CPLXMLNode *psColorTable =
2291 4 : CPLGetXMLNode(psBandNode, GTI_XML_COLORTABLE))
2292 : {
2293 2 : poBand->m_poColorTable = VRTParseColorTable(psColorTable);
2294 : }
2295 :
2296 4 : if (const CPLXMLNode *psRAT =
2297 4 : CPLGetXMLNode(psBandNode, GTI_XML_RAT))
2298 : {
2299 : poBand->m_poRAT =
2300 2 : std::make_unique<GDALDefaultRasterAttributeTable>();
2301 2 : poBand->m_poRAT->XMLInit(psRAT, "");
2302 : }
2303 : }
2304 :
2305 309 : if (static_cast<int>(adfCenterWavelength.size()) == nBandCount &&
2306 12 : adfCenterWavelength[i] != 0)
2307 : {
2308 5 : poBand->GDALRasterBand::SetMetadataItem(
2309 : "CENTRAL_WAVELENGTH_UM",
2310 5 : CPLSPrintf("%g", adfCenterWavelength[i]), "IMAGERY");
2311 : }
2312 :
2313 309 : if (static_cast<int>(adfFullWidthHalfMax.size()) == nBandCount &&
2314 12 : adfFullWidthHalfMax[i] != 0)
2315 : {
2316 5 : poBand->GDALRasterBand::SetMetadataItem(
2317 5 : "FWHM_UM", CPLSPrintf("%g", adfFullWidthHalfMax[i]), "IMAGERY");
2318 : }
2319 : }
2320 :
2321 215 : if (nBandCount == 1 && poFirstBand && poSingleColorTable &&
2322 0 : !poFirstBand->m_poColorTable)
2323 0 : poFirstBand->m_poColorTable = std::move(poSingleColorTable);
2324 :
2325 215 : const char *pszMaskBand = GetOption(MD_MASK_BAND);
2326 215 : if (pszMaskBand)
2327 7 : bHasMaskBand = CPLTestBool(pszMaskBand);
2328 215 : if (bHasMaskBand)
2329 : {
2330 8 : m_poMaskBand = std::make_unique<GDALTileIndexBand>(
2331 16 : this, 0, GDT_UInt8, nBlockXSize, nBlockYSize);
2332 : }
2333 :
2334 215 : if (dfOvrFactor == 1.0)
2335 : {
2336 211 : if (psRoot)
2337 : {
2338 84 : for (const CPLXMLNode *psIter = psRoot->psChild; psIter;
2339 66 : psIter = psIter->psNext)
2340 : {
2341 67 : if (psIter->eType == CXT_Element &&
2342 67 : strcmp(psIter->pszValue, GTI_XML_OVERVIEW_ELEMENT) == 0)
2343 : {
2344 9 : const char *pszDataset = CPLGetXMLValue(
2345 : psIter, GTI_XML_OVERVIEW_DATASET, nullptr);
2346 : const char *pszLayer =
2347 9 : CPLGetXMLValue(psIter, GTI_XML_OVERVIEW_LAYER, nullptr);
2348 9 : const char *pszFactor = CPLGetXMLValue(
2349 : psIter, GTI_XML_OVERVIEW_FACTOR, nullptr);
2350 9 : if (!pszDataset && !pszLayer && !pszFactor)
2351 : {
2352 1 : CPLError(
2353 : CE_Failure, CPLE_AppDefined,
2354 : "At least one of %s, %s or %s element "
2355 : "must be present as an %s child",
2356 : GTI_XML_OVERVIEW_DATASET, GTI_XML_OVERVIEW_LAYER,
2357 : GTI_XML_OVERVIEW_FACTOR, GTI_XML_OVERVIEW_ELEMENT);
2358 1 : return false;
2359 : }
2360 : m_aoOverviewDescriptor.emplace_back(
2361 16 : std::string(pszDataset ? pszDataset : ""),
2362 16 : CPLStringList(
2363 : GDALDeserializeOpenOptionsFromXML(psIter)),
2364 16 : std::string(pszLayer ? pszLayer : ""),
2365 24 : pszFactor ? CPLAtof(pszFactor) : 0.0);
2366 : }
2367 : }
2368 : }
2369 : else
2370 : {
2371 193 : for (int iOvr = 0;; ++iOvr)
2372 : {
2373 : const char *pszOvrDSName =
2374 387 : GetOption(CPLSPrintf("OVERVIEW_%d_DATASET", iOvr));
2375 : const char *pszOpenOptions =
2376 387 : GetOption(CPLSPrintf("OVERVIEW_%d_OPEN_OPTIONS", iOvr));
2377 : const char *pszOvrLayer =
2378 387 : GetOption(CPLSPrintf("OVERVIEW_%d_LAYER", iOvr));
2379 : const char *pszOvrFactor =
2380 387 : GetOption(CPLSPrintf("OVERVIEW_%d_FACTOR", iOvr));
2381 387 : if (!pszOvrDSName && !pszOvrLayer && !pszOvrFactor)
2382 : {
2383 : // Before GDAL 3.9.2, we started the iteration at 1.
2384 380 : if (iOvr == 0)
2385 187 : continue;
2386 193 : break;
2387 : }
2388 : m_aoOverviewDescriptor.emplace_back(
2389 14 : std::string(pszOvrDSName ? pszOvrDSName : ""),
2390 14 : pszOpenOptions ? CPLStringList(CSLTokenizeString2(
2391 : pszOpenOptions, ",", 0))
2392 : : CPLStringList(),
2393 14 : std::string(pszOvrLayer ? pszOvrLayer : ""),
2394 21 : pszOvrFactor ? CPLAtof(pszOvrFactor) : 0.0);
2395 194 : }
2396 : }
2397 : }
2398 :
2399 214 : if (psRoot)
2400 : {
2401 19 : oMDMD.XMLInit(psRoot, TRUE);
2402 : }
2403 : else
2404 : {
2405 : // Set on the dataset all metadata items from the index layer which are
2406 : // not "reserved" keywords.
2407 195 : CSLConstList papszLayerMD = m_poLayer->GetMetadata();
2408 542 : for (const auto &[pszKey, pszValue] :
2409 737 : cpl::IterateNameValue(papszLayerMD))
2410 : {
2411 271 : if (STARTS_WITH_CI(pszKey, "OVERVIEW_"))
2412 : {
2413 10 : continue;
2414 : }
2415 :
2416 261 : bool bIsVRTItem = false;
2417 3863 : for (const char *pszTest : apszTIOptions)
2418 : {
2419 3803 : if (EQUAL(pszKey, pszTest))
2420 : {
2421 201 : bIsVRTItem = true;
2422 201 : break;
2423 : }
2424 : }
2425 261 : if (!bIsVRTItem)
2426 : {
2427 60 : if (STARTS_WITH_CI(pszKey, "BAND_"))
2428 : {
2429 52 : const int nBandNr = atoi(pszKey + strlen("BAND_"));
2430 : const char *pszNextUnderscore =
2431 52 : strchr(pszKey + strlen("BAND_"), '_');
2432 52 : if (pszNextUnderscore && nBandNr >= 1 && nBandNr <= nBands)
2433 : {
2434 42 : const char *pszKeyWithoutBand = pszNextUnderscore + 1;
2435 42 : bool bIsReservedBandItem = false;
2436 132 : for (const char *pszItem : apszReservedBandItems)
2437 : {
2438 108 : if (EQUAL(pszKeyWithoutBand, pszItem))
2439 : {
2440 18 : bIsReservedBandItem = true;
2441 18 : break;
2442 : }
2443 : }
2444 42 : if (!bIsReservedBandItem)
2445 : {
2446 24 : GetRasterBand(nBandNr)
2447 24 : ->GDALRasterBand::SetMetadataItem(
2448 : pszKeyWithoutBand, pszValue);
2449 : }
2450 : }
2451 : }
2452 : else
2453 : {
2454 8 : GDALDataset::SetMetadataItem(pszKey, pszValue);
2455 : }
2456 : }
2457 : }
2458 : }
2459 :
2460 214 : const char *pszInterleave = GetOption(MD_INTERLEAVE);
2461 214 : if (pszInterleave)
2462 : {
2463 5 : GDALDataset::SetMetadataItem("INTERLEAVE", pszInterleave,
2464 : "IMAGE_STRUCTURE");
2465 5 : m_bBandInterleave = EQUAL(pszInterleave, "BAND");
2466 : }
2467 209 : else if (nBandCount > 1 && !GetMetadata("IMAGE_STRUCTURE"))
2468 : {
2469 34 : GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
2470 : }
2471 :
2472 : /* -------------------------------------------------------------------- */
2473 : /* Initialize any PAM information. */
2474 : /* -------------------------------------------------------------------- */
2475 214 : SetDescription(poOpenInfo->pszFilename);
2476 214 : TryLoadXML();
2477 :
2478 : /* -------------------------------------------------------------------- */
2479 : /* Check for overviews. */
2480 : /* -------------------------------------------------------------------- */
2481 214 : oOvManager.Initialize(this, poOpenInfo->pszFilename);
2482 :
2483 214 : m_osWarpMemory = CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
2484 214 : "WARPING_MEMORY_SIZE", "");
2485 :
2486 214 : return true;
2487 : }
2488 :
2489 : /************************************************************************/
2490 : /* GetMetadataItem() */
2491 : /************************************************************************/
2492 :
2493 107 : const char *GDALTileIndexDataset::GetMetadataItem(const char *pszName,
2494 : const char *pszDomain)
2495 : {
2496 107 : if (pszName && pszDomain && EQUAL(pszDomain, "__DEBUG__"))
2497 : {
2498 20 : if (EQUAL(pszName, "SCANNED_ONE_FEATURE_AT_OPENING"))
2499 : {
2500 4 : return m_bScannedOneFeatureAtOpening ? "YES" : "NO";
2501 : }
2502 16 : else if (EQUAL(pszName, "NUMBER_OF_CONTRIBUTING_SOURCES"))
2503 : {
2504 5 : return CPLSPrintf("%d", static_cast<int>(m_aoSourceDesc.size()));
2505 : }
2506 11 : else if (EQUAL(pszName, "MULTI_THREADED_RASTERIO_LAST_USED"))
2507 : {
2508 11 : return m_bLastMustUseMultiThreading ? "1" : "0";
2509 : }
2510 : }
2511 87 : return GDALPamDataset::GetMetadataItem(pszName, pszDomain);
2512 : }
2513 :
2514 : /************************************************************************/
2515 : /* TileIndexSupportsEditingLayerMetadata() */
2516 : /************************************************************************/
2517 :
2518 17 : bool GDALTileIndexDataset::TileIndexSupportsEditingLayerMetadata() const
2519 : {
2520 27 : return eAccess == GA_Update && m_poVectorDS->GetDriver() &&
2521 27 : EQUAL(m_poVectorDS->GetDriver()->GetDescription(), "GPKG");
2522 : }
2523 :
2524 : /************************************************************************/
2525 : /* SetMetadataItem() */
2526 : /************************************************************************/
2527 :
2528 3 : CPLErr GDALTileIndexDataset::SetMetadataItem(const char *pszName,
2529 : const char *pszValue,
2530 : const char *pszDomain)
2531 : {
2532 3 : if (m_bXMLUpdatable)
2533 : {
2534 1 : m_bXMLModified = true;
2535 1 : return GDALDataset::SetMetadataItem(pszName, pszValue, pszDomain);
2536 : }
2537 2 : else if (TileIndexSupportsEditingLayerMetadata())
2538 : {
2539 1 : m_poLayer->SetMetadataItem(pszName, pszValue, pszDomain);
2540 1 : return GDALDataset::SetMetadataItem(pszName, pszValue, pszDomain);
2541 : }
2542 : else
2543 : {
2544 1 : return GDALPamDataset::SetMetadataItem(pszName, pszValue, pszDomain);
2545 : }
2546 : }
2547 :
2548 : /************************************************************************/
2549 : /* SetMetadata() */
2550 : /************************************************************************/
2551 :
2552 3 : CPLErr GDALTileIndexDataset::SetMetadata(CSLConstList papszMD,
2553 : const char *pszDomain)
2554 : {
2555 3 : if (m_bXMLUpdatable)
2556 : {
2557 1 : m_bXMLModified = true;
2558 1 : return GDALDataset::SetMetadata(papszMD, pszDomain);
2559 : }
2560 2 : else if (TileIndexSupportsEditingLayerMetadata())
2561 : {
2562 2 : if (!pszDomain || pszDomain[0] == 0)
2563 : {
2564 4 : CPLStringList aosMD(CSLDuplicate(papszMD));
2565 :
2566 : // Reinject dataset reserved items
2567 46 : for (const char *pszItem : apszTIOptions)
2568 : {
2569 44 : if (!aosMD.FetchNameValue(pszItem))
2570 : {
2571 44 : const char *pszValue = m_poLayer->GetMetadataItem(pszItem);
2572 44 : if (pszValue)
2573 : {
2574 2 : aosMD.SetNameValue(pszItem, pszValue);
2575 : }
2576 : }
2577 : }
2578 :
2579 : // Reinject band metadata
2580 2 : CSLConstList papszExistingLayerMD = m_poLayer->GetMetadata();
2581 17 : for (int i = 0; papszExistingLayerMD && papszExistingLayerMD[i];
2582 : ++i)
2583 : {
2584 15 : if (STARTS_WITH_CI(papszExistingLayerMD[i], "BAND_"))
2585 : {
2586 12 : aosMD.AddString(papszExistingLayerMD[i]);
2587 : }
2588 : }
2589 :
2590 4 : m_poLayer->SetMetadata(aosMD.List(), pszDomain);
2591 : }
2592 : else
2593 : {
2594 0 : m_poLayer->SetMetadata(papszMD, pszDomain);
2595 : }
2596 2 : return GDALDataset::SetMetadata(papszMD, pszDomain);
2597 : }
2598 : else
2599 : {
2600 0 : return GDALPamDataset::SetMetadata(papszMD, pszDomain);
2601 : }
2602 : }
2603 :
2604 : /************************************************************************/
2605 : /* GDALTileIndexDatasetIdentify() */
2606 : /************************************************************************/
2607 :
2608 93045 : static int GDALTileIndexDatasetIdentify(GDALOpenInfo *poOpenInfo)
2609 : {
2610 93045 : if (STARTS_WITH(poOpenInfo->pszFilename, GTI_PREFIX))
2611 28 : return true;
2612 :
2613 93017 : if (STARTS_WITH(poOpenInfo->pszFilename, "<GDALTileIndexDataset"))
2614 50 : return true;
2615 :
2616 92967 : if (poOpenInfo->nHeaderBytes >= 100 &&
2617 33247 : STARTS_WITH(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
2618 : "SQLite format 3"))
2619 : {
2620 954 : if (ENDS_WITH_CI(poOpenInfo->pszFilename, ".gti.gpkg"))
2621 : {
2622 : // Most likely handled by GTI driver, but we can't be sure
2623 511 : return GDAL_IDENTIFY_UNKNOWN;
2624 : }
2625 445 : else if (poOpenInfo->IsSingleAllowedDriver("GTI") &&
2626 2 : poOpenInfo->IsExtensionEqualToCI("gpkg"))
2627 : {
2628 2 : return true;
2629 : }
2630 : }
2631 :
2632 92454 : if (poOpenInfo->nHeaderBytes > 0 &&
2633 33352 : (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) != 0)
2634 : {
2635 66704 : if (strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
2636 33340 : "<GDALTileIndexDataset") ||
2637 66692 : ENDS_WITH_CI(poOpenInfo->pszFilename, ".gti.fgb") ||
2638 33340 : ENDS_WITH_CI(poOpenInfo->pszFilename, ".gti.parquet"))
2639 : {
2640 12 : return true;
2641 : }
2642 33340 : else if (poOpenInfo->IsSingleAllowedDriver("GTI") &&
2643 0 : (poOpenInfo->IsExtensionEqualToCI("fgb") ||
2644 0 : poOpenInfo->IsExtensionEqualToCI("parquet")))
2645 : {
2646 0 : return true;
2647 : }
2648 : }
2649 :
2650 92442 : return false;
2651 : }
2652 :
2653 : /************************************************************************/
2654 : /* GDALTileIndexDatasetOpen() */
2655 : /************************************************************************/
2656 :
2657 286 : static GDALDataset *GDALTileIndexDatasetOpen(GDALOpenInfo *poOpenInfo)
2658 : {
2659 286 : if (GDALTileIndexDatasetIdentify(poOpenInfo) == GDAL_IDENTIFY_FALSE)
2660 0 : return nullptr;
2661 572 : auto poDS = std::make_unique<GDALTileIndexDataset>();
2662 286 : if (!poDS->Open(poOpenInfo))
2663 72 : return nullptr;
2664 214 : return poDS.release();
2665 : }
2666 :
2667 : /************************************************************************/
2668 : /* ~GDALTileIndexDataset() */
2669 : /************************************************************************/
2670 :
2671 572 : GDALTileIndexDataset::~GDALTileIndexDataset()
2672 : {
2673 286 : if (m_poVectorDS && m_bIsSQLResultLayer)
2674 4 : m_poVectorDS->ReleaseResultSet(m_poLayer);
2675 :
2676 286 : GDALTileIndexDataset::FlushCache(true);
2677 572 : }
2678 :
2679 : /************************************************************************/
2680 : /* FlushCache() */
2681 : /************************************************************************/
2682 :
2683 287 : CPLErr GDALTileIndexDataset::FlushCache(bool bAtClosing)
2684 : {
2685 287 : CPLErr eErr = CE_None;
2686 287 : if (bAtClosing && m_bXMLModified)
2687 : {
2688 : CPLXMLNode *psRoot =
2689 1 : CPLGetXMLNode(m_psXMLTree.get(), "=GDALTileIndexDataset");
2690 :
2691 : // Suppress existing dataset metadata
2692 : while (true)
2693 : {
2694 1 : CPLXMLNode *psExistingMetadata = CPLGetXMLNode(psRoot, "Metadata");
2695 1 : if (!psExistingMetadata)
2696 1 : break;
2697 0 : CPLRemoveXMLChild(psRoot, psExistingMetadata);
2698 0 : }
2699 :
2700 : // Serialize new dataset metadata
2701 1 : if (CPLXMLNode *psMD = oMDMD.Serialize())
2702 1 : CPLAddXMLChild(psRoot, psMD);
2703 :
2704 : // Update existing band metadata
2705 1 : if (CPLGetXMLNode(psRoot, GTI_XML_BAND_ELEMENT))
2706 : {
2707 0 : for (CPLXMLNode *psIter = psRoot->psChild; psIter;
2708 0 : psIter = psIter->psNext)
2709 : {
2710 0 : if (psIter->eType == CXT_Element &&
2711 0 : strcmp(psIter->pszValue, GTI_XML_BAND_ELEMENT))
2712 : {
2713 : const char *pszBand =
2714 0 : CPLGetXMLValue(psIter, GTI_XML_BAND_NUMBER, nullptr);
2715 0 : if (pszBand)
2716 : {
2717 0 : const int nBand = atoi(pszBand);
2718 0 : if (nBand >= 1 && nBand <= nBands)
2719 : {
2720 : while (true)
2721 : {
2722 : CPLXMLNode *psExistingMetadata =
2723 0 : CPLGetXMLNode(psIter, "Metadata");
2724 0 : if (!psExistingMetadata)
2725 0 : break;
2726 0 : CPLRemoveXMLChild(psIter, psExistingMetadata);
2727 0 : }
2728 :
2729 0 : auto poBand = cpl::down_cast<GDALTileIndexBand *>(
2730 0 : papoBands[nBand - 1]);
2731 0 : if (CPLXMLNode *psMD = poBand->oMDMD.Serialize())
2732 0 : CPLAddXMLChild(psIter, psMD);
2733 : }
2734 : }
2735 : }
2736 : }
2737 : }
2738 : else
2739 : {
2740 : // Create new band objects if they have metadata
2741 2 : std::vector<CPLXMLTreeCloser> aoBandXML;
2742 1 : bool bHasBandMD = false;
2743 2 : for (int i = 1; i <= nBands; ++i)
2744 : {
2745 : auto poBand =
2746 1 : cpl::down_cast<GDALTileIndexBand *>(papoBands[i - 1]);
2747 1 : auto psMD = poBand->oMDMD.Serialize();
2748 1 : if (psMD)
2749 1 : bHasBandMD = true;
2750 1 : aoBandXML.emplace_back(CPLXMLTreeCloser(psMD));
2751 : }
2752 1 : if (bHasBandMD)
2753 : {
2754 2 : for (int i = 1; i <= nBands; ++i)
2755 : {
2756 : auto poBand =
2757 1 : cpl::down_cast<GDALTileIndexBand *>(papoBands[i - 1]);
2758 :
2759 1 : CPLXMLNode *psBand = CPLCreateXMLNode(psRoot, CXT_Element,
2760 : GTI_XML_BAND_ELEMENT);
2761 1 : CPLAddXMLAttributeAndValue(psBand, GTI_XML_BAND_NUMBER,
2762 : CPLSPrintf("%d", i));
2763 1 : CPLAddXMLAttributeAndValue(
2764 : psBand, GTI_XML_BAND_DATATYPE,
2765 : GDALGetDataTypeName(poBand->GetRasterDataType()));
2766 :
2767 1 : const char *pszDescription = poBand->GetDescription();
2768 1 : if (pszDescription && pszDescription[0])
2769 0 : CPLSetXMLValue(psBand, GTI_XML_BAND_DESCRIPTION,
2770 : pszDescription);
2771 :
2772 1 : const auto eColorInterp = poBand->GetColorInterpretation();
2773 1 : if (eColorInterp != GCI_Undefined)
2774 1 : CPLSetXMLValue(
2775 : psBand, GTI_XML_BAND_COLORINTERP,
2776 : GDALGetColorInterpretationName(eColorInterp));
2777 :
2778 1 : if (!std::isnan(poBand->m_dfOffset))
2779 0 : CPLSetXMLValue(psBand, GTI_XML_BAND_OFFSET,
2780 : CPLSPrintf("%.16g", poBand->m_dfOffset));
2781 :
2782 1 : if (!std::isnan(poBand->m_dfScale))
2783 0 : CPLSetXMLValue(psBand, GTI_XML_BAND_SCALE,
2784 : CPLSPrintf("%.16g", poBand->m_dfScale));
2785 :
2786 1 : if (!poBand->m_osUnit.empty())
2787 0 : CPLSetXMLValue(psBand, GTI_XML_BAND_UNITTYPE,
2788 : poBand->m_osUnit.c_str());
2789 :
2790 1 : if (poBand->m_bNoDataValueSet)
2791 : {
2792 0 : CPLSetXMLValue(
2793 : psBand, GTI_XML_BAND_NODATAVALUE,
2794 0 : VRTSerializeNoData(poBand->m_dfNoDataValue,
2795 : poBand->GetRasterDataType(), 18)
2796 : .c_str());
2797 : }
2798 1 : if (aoBandXML[i - 1])
2799 : {
2800 1 : CPLAddXMLChild(psBand, aoBandXML[i - 1].release());
2801 : }
2802 : }
2803 : }
2804 : }
2805 :
2806 1 : if (!CPLSerializeXMLTreeToFile(m_psXMLTree.get(), GetDescription()))
2807 0 : eErr = CE_Failure;
2808 : }
2809 :
2810 : // We also clear the cache of opened sources, in case the user would
2811 : // change the content of a source and would want the GTI dataset to see
2812 : // the refreshed content.
2813 287 : m_oMapSharedSources.clear();
2814 287 : m_dfLastMinXFilter = std::numeric_limits<double>::quiet_NaN();
2815 287 : m_dfLastMinYFilter = std::numeric_limits<double>::quiet_NaN();
2816 287 : m_dfLastMaxXFilter = std::numeric_limits<double>::quiet_NaN();
2817 287 : m_dfLastMaxYFilter = std::numeric_limits<double>::quiet_NaN();
2818 287 : m_anLastBands.clear();
2819 287 : m_aoSourceDesc.clear();
2820 287 : if (GDALPamDataset::FlushCache(bAtClosing) != CE_None)
2821 0 : eErr = CE_Failure;
2822 287 : return eErr;
2823 : }
2824 :
2825 : /************************************************************************/
2826 : /* LoadOverviews() */
2827 : /************************************************************************/
2828 :
2829 44 : void GDALTileIndexDataset::LoadOverviews()
2830 : {
2831 44 : if (m_apoOverviews.empty() && !m_aoOverviewDescriptor.empty())
2832 : {
2833 28 : for (const auto &[osDSName, aosOpenOptions, osLyrName, dfFactor] :
2834 42 : m_aoOverviewDescriptor)
2835 : {
2836 28 : CPLStringList aosNewOpenOptions(aosOpenOptions);
2837 14 : if (dfFactor != 0)
2838 : {
2839 : aosNewOpenOptions.SetNameValue("@FACTOR",
2840 4 : CPLSPrintf("%.17g", dfFactor));
2841 : }
2842 14 : if (!osLyrName.empty())
2843 : {
2844 5 : aosNewOpenOptions.SetNameValue("@LAYER", osLyrName.c_str());
2845 : }
2846 :
2847 : std::unique_ptr<GDALDataset> poOvrDS(GDALDataset::Open(
2848 28 : !osDSName.empty() ? osDSName.c_str() : GetDescription(),
2849 : GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr,
2850 56 : aosNewOpenOptions.List(), nullptr));
2851 :
2852 : const auto IsSmaller =
2853 10 : [](const GDALDataset *a, const GDALDataset *b)
2854 : {
2855 10 : return (a->GetRasterXSize() < b->GetRasterXSize() &&
2856 11 : a->GetRasterYSize() <= b->GetRasterYSize()) ||
2857 1 : (a->GetRasterYSize() < b->GetRasterYSize() &&
2858 10 : a->GetRasterXSize() <= b->GetRasterXSize());
2859 : };
2860 :
2861 32 : if (poOvrDS &&
2862 18 : ((m_apoOverviews.empty() && IsSmaller(poOvrDS.get(), this)) ||
2863 1 : ((!m_apoOverviews.empty() &&
2864 14 : IsSmaller(poOvrDS.get(), m_apoOverviews.back().get())))))
2865 : {
2866 8 : if (poOvrDS->GetRasterCount() == GetRasterCount())
2867 : {
2868 8 : m_apoOverviews.emplace_back(std::move(poOvrDS));
2869 : // Add the overviews of the overview, unless the OVERVIEW_LEVEL
2870 : // option option is specified
2871 8 : if (aosOpenOptions.FetchNameValue("OVERVIEW_LEVEL") ==
2872 : nullptr)
2873 : {
2874 7 : const int nOverviewCount = m_apoOverviews.back()
2875 7 : ->GetRasterBand(1)
2876 7 : ->GetOverviewCount();
2877 8 : for (int i = 0; i < nOverviewCount; ++i)
2878 : {
2879 : aosNewOpenOptions.SetNameValue("OVERVIEW_LEVEL",
2880 1 : CPLSPrintf("%d", i));
2881 : std::unique_ptr<GDALDataset> poOvrOfOvrDS(
2882 : GDALDataset::Open(
2883 2 : !osDSName.empty() ? osDSName.c_str()
2884 0 : : GetDescription(),
2885 : GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
2886 1 : nullptr, aosNewOpenOptions.List(),
2887 3 : nullptr));
2888 2 : if (poOvrOfOvrDS &&
2889 1 : poOvrOfOvrDS->GetRasterCount() ==
2890 3 : GetRasterCount() &&
2891 1 : IsSmaller(poOvrOfOvrDS.get(),
2892 1 : m_apoOverviews.back().get()))
2893 : {
2894 : m_apoOverviews.emplace_back(
2895 1 : std::move(poOvrOfOvrDS));
2896 : }
2897 : }
2898 : }
2899 : }
2900 : else
2901 : {
2902 0 : CPLError(CE_Warning, CPLE_AppDefined,
2903 : "%s has not the same number of bands as %s",
2904 0 : poOvrDS->GetDescription(), GetDescription());
2905 : }
2906 : }
2907 : }
2908 : }
2909 44 : }
2910 :
2911 : /************************************************************************/
2912 : /* GetOverviewCount() */
2913 : /************************************************************************/
2914 :
2915 68 : int GDALTileIndexBand::GetOverviewCount()
2916 : {
2917 68 : const int nPAMOverviews = GDALPamRasterBand::GetOverviewCount();
2918 68 : if (nPAMOverviews)
2919 24 : return nPAMOverviews;
2920 :
2921 44 : m_poDS->LoadOverviews();
2922 44 : return static_cast<int>(m_poDS->m_apoOverviews.size());
2923 : }
2924 :
2925 : /************************************************************************/
2926 : /* GetOverview() */
2927 : /************************************************************************/
2928 :
2929 35 : GDALRasterBand *GDALTileIndexBand::GetOverview(int iOvr)
2930 : {
2931 35 : if (iOvr < 0 || iOvr >= GetOverviewCount())
2932 6 : return nullptr;
2933 :
2934 29 : const int nPAMOverviews = GDALPamRasterBand::GetOverviewCount();
2935 29 : if (nPAMOverviews)
2936 16 : return GDALPamRasterBand::GetOverview(iOvr);
2937 :
2938 13 : if (nBand == 0)
2939 : {
2940 1 : auto poBand = m_poDS->m_apoOverviews[iOvr]->GetRasterBand(1);
2941 1 : if (!poBand)
2942 0 : return nullptr;
2943 1 : return poBand->GetMaskBand();
2944 : }
2945 : else
2946 : {
2947 12 : return m_poDS->m_apoOverviews[iOvr]->GetRasterBand(nBand);
2948 : }
2949 : }
2950 :
2951 : /************************************************************************/
2952 : /* GetGeoTransform() */
2953 : /************************************************************************/
2954 :
2955 23 : CPLErr GDALTileIndexDataset::GetGeoTransform(GDALGeoTransform >) const
2956 : {
2957 23 : gt = m_gt;
2958 23 : return CE_None;
2959 : }
2960 :
2961 : /************************************************************************/
2962 : /* GetSpatialRef() */
2963 : /************************************************************************/
2964 :
2965 18 : const OGRSpatialReference *GDALTileIndexDataset::GetSpatialRef() const
2966 : {
2967 18 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
2968 : }
2969 :
2970 : /************************************************************************/
2971 : /* GDALTileIndexBand() */
2972 : /************************************************************************/
2973 :
2974 305 : GDALTileIndexBand::GDALTileIndexBand(GDALTileIndexDataset *poDSIn, int nBandIn,
2975 : GDALDataType eDT, int nBlockXSizeIn,
2976 305 : int nBlockYSizeIn)
2977 : {
2978 305 : m_poDS = poDSIn;
2979 305 : nBand = nBandIn;
2980 305 : eDataType = eDT;
2981 305 : nRasterXSize = poDSIn->GetRasterXSize();
2982 305 : nRasterYSize = poDSIn->GetRasterYSize();
2983 305 : nBlockXSize = nBlockXSizeIn;
2984 305 : nBlockYSize = nBlockYSizeIn;
2985 305 : }
2986 :
2987 : /************************************************************************/
2988 : /* IReadBlock() */
2989 : /************************************************************************/
2990 :
2991 13 : CPLErr GDALTileIndexBand::IReadBlock(int nBlockXOff, int nBlockYOff,
2992 : void *pImage)
2993 :
2994 : {
2995 13 : const int nPixelSize = GDALGetDataTypeSizeBytes(eDataType);
2996 :
2997 13 : int nReadXSize = nBlockXSize;
2998 13 : int nReadYSize = nBlockYSize;
2999 13 : GetActualBlockSize(nBlockXOff, nBlockYOff, &nReadXSize, &nReadYSize);
3000 :
3001 : GDALRasterIOExtraArg sExtraArg;
3002 13 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
3003 :
3004 26 : return IRasterIO(
3005 13 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize, nReadXSize,
3006 : nReadYSize, pImage, nReadXSize, nReadYSize, eDataType, nPixelSize,
3007 26 : static_cast<GSpacing>(nPixelSize) * nBlockXSize, &sExtraArg);
3008 : }
3009 :
3010 : /************************************************************************/
3011 : /* IRasterIO() */
3012 : /************************************************************************/
3013 :
3014 150 : CPLErr GDALTileIndexBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
3015 : int nXSize, int nYSize, void *pData,
3016 : int nBufXSize, int nBufYSize,
3017 : GDALDataType eBufType, GSpacing nPixelSpace,
3018 : GSpacing nLineSpace,
3019 : GDALRasterIOExtraArg *psExtraArg)
3020 : {
3021 150 : int anBand[] = {nBand};
3022 :
3023 150 : return m_poDS->IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
3024 : nBufXSize, nBufYSize, eBufType, 1, anBand,
3025 300 : nPixelSpace, nLineSpace, 0, psExtraArg);
3026 : }
3027 :
3028 : /************************************************************************/
3029 : /* IGetDataCoverageStatus() */
3030 : /************************************************************************/
3031 :
3032 : #ifndef HAVE_GEOS
3033 : int GDALTileIndexBand::IGetDataCoverageStatus(int /* nXOff */, int /* nYOff */,
3034 : int /* nXSize */,
3035 : int /* nYSize */,
3036 : int /* nMaskFlagStop */,
3037 : double *pdfDataPct)
3038 : {
3039 : if (pdfDataPct != nullptr)
3040 : *pdfDataPct = -1.0;
3041 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
3042 : GDAL_DATA_COVERAGE_STATUS_DATA;
3043 : }
3044 : #else
3045 9 : int GDALTileIndexBand::IGetDataCoverageStatus(int nXOff, int nYOff, int nXSize,
3046 : int nYSize, int nMaskFlagStop,
3047 : double *pdfDataPct)
3048 : {
3049 9 : if (pdfDataPct != nullptr)
3050 9 : *pdfDataPct = -1.0;
3051 :
3052 9 : const double dfMinX = m_poDS->m_gt.xorig + nXOff * m_poDS->m_gt.xscale;
3053 9 : const double dfMaxX = dfMinX + nXSize * m_poDS->m_gt.xscale;
3054 9 : const double dfMaxY = m_poDS->m_gt.yorig + nYOff * m_poDS->m_gt.yscale;
3055 9 : const double dfMinY = dfMaxY + nYSize * m_poDS->m_gt.yscale;
3056 :
3057 9 : OGRLayer *poSQLLayer = nullptr;
3058 9 : if (!m_poDS->m_osSpatialSQL.empty())
3059 : {
3060 : const std::string osSQL =
3061 2 : CPLString(m_poDS->m_osSpatialSQL)
3062 4 : .replaceAll("{XMIN}", CPLSPrintf("%.17g", dfMinX))
3063 4 : .replaceAll("{YMIN}", CPLSPrintf("%.17g", dfMinY))
3064 4 : .replaceAll("{XMAX}", CPLSPrintf("%.17g", dfMaxX))
3065 4 : .replaceAll("{YMAX}", CPLSPrintf("%.17g", dfMaxY));
3066 : poSQLLayer =
3067 2 : m_poDS->m_poVectorDS->ExecuteSQL(osSQL.c_str(), nullptr, nullptr);
3068 2 : if (!poSQLLayer)
3069 1 : return 0;
3070 : }
3071 : else
3072 : {
3073 7 : m_poDS->m_poLayer->SetSpatialFilterRect(dfMinX, dfMinY, dfMaxX, dfMaxY);
3074 7 : m_poDS->m_poLayer->ResetReading();
3075 : }
3076 :
3077 8 : OGRLayer *const poLayer = poSQLLayer ? poSQLLayer : m_poDS->m_poLayer;
3078 :
3079 8 : int nStatus = 0;
3080 :
3081 16 : auto poPolyNonCoveredBySources = std::make_unique<OGRPolygon>();
3082 : {
3083 16 : auto poLR = std::make_unique<OGRLinearRing>();
3084 8 : poLR->addPoint(nXOff, nYOff);
3085 8 : poLR->addPoint(nXOff, nYOff + nYSize);
3086 8 : poLR->addPoint(nXOff + nXSize, nYOff + nYSize);
3087 8 : poLR->addPoint(nXOff + nXSize, nYOff);
3088 8 : poLR->addPoint(nXOff, nYOff);
3089 8 : poPolyNonCoveredBySources->addRingDirectly(poLR.release());
3090 : }
3091 : while (true)
3092 : {
3093 12 : auto poFeature = std::unique_ptr<OGRFeature>(poLayer->GetNextFeature());
3094 12 : if (!poFeature)
3095 2 : break;
3096 10 : if (!poFeature->IsFieldSetAndNotNull(m_poDS->m_nLocationFieldIndex))
3097 : {
3098 0 : continue;
3099 : }
3100 :
3101 10 : const auto poGeom = poFeature->GetGeometryRef();
3102 10 : if (!poGeom || poGeom->IsEmpty())
3103 0 : continue;
3104 :
3105 10 : OGREnvelope sSourceEnvelope;
3106 10 : poGeom->getEnvelope(&sSourceEnvelope);
3107 :
3108 : const double dfDstXOff = std::max<double>(
3109 : nXOff,
3110 10 : (sSourceEnvelope.MinX - m_poDS->m_gt.xorig) / m_poDS->m_gt.xscale);
3111 : const double dfDstXOff2 = std::min<double>(
3112 10 : nXOff + nXSize,
3113 10 : (sSourceEnvelope.MaxX - m_poDS->m_gt.xorig) / m_poDS->m_gt.xscale);
3114 : const double dfDstYOff = std::max<double>(
3115 : nYOff,
3116 10 : (sSourceEnvelope.MaxY - m_poDS->m_gt.yorig) / m_poDS->m_gt.yscale);
3117 : const double dfDstYOff2 = std::min<double>(
3118 10 : nYOff + nYSize,
3119 10 : (sSourceEnvelope.MinY - m_poDS->m_gt.yorig) / m_poDS->m_gt.yscale);
3120 :
3121 : // CPLDebug("GTI", "dfDstXOff=%f, dfDstXOff2=%f, dfDstYOff=%f, dfDstYOff2=%f",
3122 : // dfDstXOff, dfDstXOff2, dfDstYOff, dfDstXOff2);
3123 :
3124 : // Check if the AOI is fully inside the source
3125 10 : if (nXOff >= dfDstXOff && nYOff >= dfDstYOff &&
3126 7 : nXOff + nXSize <= dfDstXOff2 && nYOff + nYSize <= dfDstYOff2)
3127 : {
3128 4 : if (pdfDataPct)
3129 4 : *pdfDataPct = 100.0;
3130 4 : return GDAL_DATA_COVERAGE_STATUS_DATA;
3131 : }
3132 :
3133 : // Check intersection of bounding boxes.
3134 6 : if (dfDstXOff2 > nXOff && dfDstYOff2 > nYOff &&
3135 6 : dfDstXOff < nXOff + nXSize && dfDstYOff < nYOff + nYSize)
3136 : {
3137 6 : nStatus |= GDAL_DATA_COVERAGE_STATUS_DATA;
3138 6 : if (poPolyNonCoveredBySources)
3139 : {
3140 6 : OGRPolygon oPolySource;
3141 6 : auto poLR = std::make_unique<OGRLinearRing>();
3142 6 : poLR->addPoint(dfDstXOff, dfDstYOff);
3143 6 : poLR->addPoint(dfDstXOff, dfDstYOff2);
3144 6 : poLR->addPoint(dfDstXOff2, dfDstYOff2);
3145 6 : poLR->addPoint(dfDstXOff2, dfDstYOff);
3146 6 : poLR->addPoint(dfDstXOff, dfDstYOff);
3147 6 : oPolySource.addRingDirectly(poLR.release());
3148 : auto poRes = std::unique_ptr<OGRGeometry>(
3149 6 : poPolyNonCoveredBySources->Difference(&oPolySource));
3150 6 : if (poRes && poRes->IsEmpty())
3151 : {
3152 2 : if (pdfDataPct)
3153 2 : *pdfDataPct = 100.0;
3154 2 : return GDAL_DATA_COVERAGE_STATUS_DATA;
3155 : }
3156 4 : else if (poRes && poRes->getGeometryType() == wkbPolygon)
3157 : {
3158 4 : poPolyNonCoveredBySources.reset(
3159 : poRes.release()->toPolygon());
3160 : }
3161 : else
3162 : {
3163 0 : poPolyNonCoveredBySources.reset();
3164 : }
3165 : }
3166 : }
3167 4 : if (nMaskFlagStop != 0 && (nStatus & nMaskFlagStop) != 0)
3168 : {
3169 0 : if (poSQLLayer)
3170 0 : m_poDS->ReleaseResultSet(poSQLLayer);
3171 0 : return nStatus;
3172 : }
3173 4 : }
3174 :
3175 2 : if (poSQLLayer)
3176 0 : m_poDS->ReleaseResultSet(poSQLLayer);
3177 :
3178 2 : if (poPolyNonCoveredBySources)
3179 : {
3180 2 : if (!poPolyNonCoveredBySources->IsEmpty())
3181 2 : nStatus |= GDAL_DATA_COVERAGE_STATUS_EMPTY;
3182 2 : if (pdfDataPct)
3183 2 : *pdfDataPct = 100.0 * (1.0 - poPolyNonCoveredBySources->get_Area() /
3184 2 : nXSize / nYSize);
3185 : }
3186 2 : return nStatus;
3187 : }
3188 : #endif // HAVE_GEOS
3189 :
3190 : /************************************************************************/
3191 : /* GetMetadataDomainList() */
3192 : /************************************************************************/
3193 :
3194 1 : char **GDALTileIndexBand::GetMetadataDomainList()
3195 : {
3196 1 : return CSLAddString(GDALRasterBand::GetMetadataDomainList(),
3197 1 : "LocationInfo");
3198 : }
3199 :
3200 : /************************************************************************/
3201 : /* GetMetadataItem() */
3202 : /************************************************************************/
3203 :
3204 44 : const char *GDALTileIndexBand::GetMetadataItem(const char *pszName,
3205 : const char *pszDomain)
3206 :
3207 : {
3208 : /* ==================================================================== */
3209 : /* LocationInfo handling. */
3210 : /* ==================================================================== */
3211 44 : if (pszDomain != nullptr && EQUAL(pszDomain, "LocationInfo") &&
3212 19 : (STARTS_WITH_CI(pszName, "Pixel_") ||
3213 6 : STARTS_WITH_CI(pszName, "GeoPixel_")))
3214 : {
3215 : // What pixel are we aiming at?
3216 18 : int iPixel = 0;
3217 18 : int iLine = 0;
3218 :
3219 18 : if (STARTS_WITH_CI(pszName, "Pixel_"))
3220 : {
3221 13 : pszName += strlen("Pixel_");
3222 13 : iPixel = atoi(pszName);
3223 13 : const char *const pszUnderscore = strchr(pszName, '_');
3224 13 : if (!pszUnderscore)
3225 2 : return nullptr;
3226 11 : iLine = atoi(pszUnderscore + 1);
3227 : }
3228 5 : else if (STARTS_WITH_CI(pszName, "GeoPixel_"))
3229 : {
3230 5 : pszName += strlen("GeoPixel_");
3231 5 : const double dfGeoX = CPLAtof(pszName);
3232 5 : const char *const pszUnderscore = strchr(pszName, '_');
3233 5 : if (!pszUnderscore)
3234 2 : return nullptr;
3235 3 : const double dfGeoY = CPLAtof(pszUnderscore + 1);
3236 :
3237 3 : double adfInvGeoTransform[6] = {0.0};
3238 3 : if (!GDALInvGeoTransform(m_poDS->m_gt.data(), adfInvGeoTransform))
3239 0 : return nullptr;
3240 :
3241 3 : iPixel = static_cast<int>(floor(adfInvGeoTransform[0] +
3242 3 : adfInvGeoTransform[1] * dfGeoX +
3243 3 : adfInvGeoTransform[2] * dfGeoY));
3244 3 : iLine = static_cast<int>(floor(adfInvGeoTransform[3] +
3245 3 : adfInvGeoTransform[4] * dfGeoX +
3246 3 : adfInvGeoTransform[5] * dfGeoY));
3247 : }
3248 : else
3249 : {
3250 0 : return nullptr;
3251 : }
3252 :
3253 23 : if (iPixel < 0 || iLine < 0 || iPixel >= GetXSize() ||
3254 9 : iLine >= GetYSize())
3255 6 : return nullptr;
3256 :
3257 8 : const int anBand[] = {nBand};
3258 8 : if (!m_poDS->CollectSources(iPixel, iLine, 1, 1, 1, anBand,
3259 : /* bMultiThreadAllowed = */ false))
3260 0 : return nullptr;
3261 :
3262 : // Format into XML.
3263 8 : m_osLastLocationInfo = "<LocationInfo>";
3264 :
3265 8 : if (!m_poDS->m_aoSourceDesc.empty())
3266 : {
3267 : const auto AddSource =
3268 6 : [&](const GDALTileIndexDataset::SourceDesc &oSourceDesc)
3269 : {
3270 6 : m_osLastLocationInfo += "<File>";
3271 : char *const pszXMLEscaped =
3272 6 : CPLEscapeString(oSourceDesc.osName.c_str(), -1, CPLES_XML);
3273 6 : m_osLastLocationInfo += pszXMLEscaped;
3274 6 : CPLFree(pszXMLEscaped);
3275 6 : m_osLastLocationInfo += "</File>";
3276 11 : };
3277 :
3278 5 : if (!m_poDS->NeedInitBuffer(1, anBand))
3279 : {
3280 4 : AddSource(m_poDS->m_aoSourceDesc.back());
3281 : }
3282 : else
3283 : {
3284 3 : for (const auto &oSourceDesc : m_poDS->m_aoSourceDesc)
3285 : {
3286 2 : if (oSourceDesc.poDS)
3287 2 : AddSource(oSourceDesc);
3288 : }
3289 : }
3290 : }
3291 :
3292 8 : m_osLastLocationInfo += "</LocationInfo>";
3293 :
3294 8 : return m_osLastLocationInfo.c_str();
3295 : }
3296 :
3297 26 : return GDALPamRasterBand::GetMetadataItem(pszName, pszDomain);
3298 : }
3299 :
3300 : /************************************************************************/
3301 : /* SetMetadataItem() */
3302 : /************************************************************************/
3303 :
3304 13 : CPLErr GDALTileIndexBand::SetMetadataItem(const char *pszName,
3305 : const char *pszValue,
3306 : const char *pszDomain)
3307 : {
3308 13 : if (nBand > 0 && m_poDS->m_bXMLUpdatable)
3309 : {
3310 1 : m_poDS->m_bXMLModified = true;
3311 1 : return GDALRasterBand::SetMetadataItem(pszName, pszValue, pszDomain);
3312 : }
3313 12 : else if (nBand > 0 && m_poDS->TileIndexSupportsEditingLayerMetadata())
3314 : {
3315 6 : m_poDS->m_poLayer->SetMetadataItem(
3316 6 : CPLSPrintf("BAND_%d_%s", nBand, pszName), pszValue, pszDomain);
3317 6 : return GDALRasterBand::SetMetadataItem(pszName, pszValue, pszDomain);
3318 : }
3319 : else
3320 : {
3321 6 : return GDALPamRasterBand::SetMetadataItem(pszName, pszValue, pszDomain);
3322 : }
3323 : }
3324 :
3325 : /************************************************************************/
3326 : /* SetMetadata() */
3327 : /************************************************************************/
3328 :
3329 2 : CPLErr GDALTileIndexBand::SetMetadata(CSLConstList papszMD,
3330 : const char *pszDomain)
3331 : {
3332 2 : if (nBand > 0 && m_poDS->m_bXMLUpdatable)
3333 : {
3334 1 : m_poDS->m_bXMLModified = true;
3335 1 : return GDALRasterBand::SetMetadata(papszMD, pszDomain);
3336 : }
3337 1 : else if (nBand > 0 && m_poDS->TileIndexSupportsEditingLayerMetadata())
3338 : {
3339 2 : CPLStringList aosMD;
3340 :
3341 1 : if (!pszDomain || pszDomain[0] == 0)
3342 : {
3343 : // Reinject dataset metadata
3344 : CSLConstList papszLayerMD =
3345 1 : m_poDS->m_poLayer->GetMetadata(pszDomain);
3346 14 : for (const char *const *papszIter = papszLayerMD;
3347 14 : papszIter && *papszIter; ++papszIter)
3348 : {
3349 13 : if (!STARTS_WITH(*papszIter, "BAND_") ||
3350 12 : STARTS_WITH(*papszIter, MD_BAND_COUNT))
3351 1 : aosMD.AddString(*papszIter);
3352 : }
3353 : }
3354 :
3355 8 : for (int i = 0; papszMD && papszMD[i]; ++i)
3356 : {
3357 7 : aosMD.AddString(CPLSPrintf("BAND_%d_%s", nBand, papszMD[i]));
3358 : }
3359 :
3360 1 : if (!pszDomain || pszDomain[0] == 0)
3361 : {
3362 4 : for (const char *pszItem : apszReservedBandItems)
3363 : {
3364 3 : const char *pszKey = CPLSPrintf("BAND_%d_%s", nBand, pszItem);
3365 3 : if (!aosMD.FetchNameValue(pszKey))
3366 : {
3367 3 : if (const char *pszVal =
3368 3 : m_poDS->m_poLayer->GetMetadataItem(pszKey))
3369 : {
3370 3 : aosMD.SetNameValue(pszKey, pszVal);
3371 : }
3372 : }
3373 : }
3374 : }
3375 :
3376 1 : m_poDS->m_poLayer->SetMetadata(aosMD.List(), pszDomain);
3377 1 : return GDALRasterBand::SetMetadata(papszMD, pszDomain);
3378 : }
3379 : else
3380 : {
3381 0 : return GDALPamRasterBand::SetMetadata(papszMD, pszDomain);
3382 : }
3383 : }
3384 :
3385 : /************************************************************************/
3386 : /* GetSrcDstWin() */
3387 : /************************************************************************/
3388 :
3389 370 : static bool GetSrcDstWin(const GDALGeoTransform &tileGT, int nTileXSize,
3390 : int nTileYSize, const GDALGeoTransform &vrtGT,
3391 : int nVRTXSize, int nVRTYSize, double *pdfSrcXOff,
3392 : double *pdfSrcYOff, double *pdfSrcXSize,
3393 : double *pdfSrcYSize, double *pdfDstXOff,
3394 : double *pdfDstYOff, double *pdfDstXSize,
3395 : double *pdfDstYSize)
3396 : {
3397 370 : const double minX = vrtGT.xorig;
3398 370 : const double we_res = vrtGT.xscale;
3399 370 : const double maxX = minX + nVRTXSize * we_res;
3400 370 : const double maxY = vrtGT.yorig;
3401 370 : const double ns_res = vrtGT.yscale;
3402 370 : const double minY = maxY + nVRTYSize * ns_res;
3403 :
3404 : /* Check that the destination bounding box intersects the source bounding
3405 : * box */
3406 370 : if (tileGT.xorig + nTileXSize * tileGT.xscale <= minX)
3407 0 : return false;
3408 370 : if (tileGT.xorig >= maxX)
3409 1 : return false;
3410 369 : if (tileGT.yorig + nTileYSize * tileGT.yscale >= maxY)
3411 0 : return false;
3412 369 : if (tileGT.yorig <= minY)
3413 0 : return false;
3414 :
3415 369 : if (tileGT.xorig < minX)
3416 : {
3417 1 : *pdfSrcXOff = (minX - tileGT.xorig) / tileGT.xscale;
3418 1 : *pdfDstXOff = 0.0;
3419 : }
3420 : else
3421 : {
3422 368 : *pdfSrcXOff = 0.0;
3423 368 : *pdfDstXOff = ((tileGT.xorig - minX) / we_res);
3424 : }
3425 369 : if (maxY < tileGT.yorig)
3426 : {
3427 1 : *pdfSrcYOff = (tileGT.yorig - maxY) / -tileGT.yscale;
3428 1 : *pdfDstYOff = 0.0;
3429 : }
3430 : else
3431 : {
3432 368 : *pdfSrcYOff = 0.0;
3433 368 : *pdfDstYOff = ((maxY - tileGT.yorig) / -ns_res);
3434 : }
3435 :
3436 369 : *pdfSrcXSize = nTileXSize;
3437 369 : *pdfSrcYSize = nTileYSize;
3438 369 : if (*pdfSrcXOff > 0)
3439 1 : *pdfSrcXSize -= *pdfSrcXOff;
3440 369 : if (*pdfSrcYOff > 0)
3441 1 : *pdfSrcYSize -= *pdfSrcYOff;
3442 :
3443 369 : const double dfSrcToDstXSize = tileGT.xscale / we_res;
3444 369 : *pdfDstXSize = *pdfSrcXSize * dfSrcToDstXSize;
3445 369 : const double dfSrcToDstYSize = tileGT.yscale / ns_res;
3446 369 : *pdfDstYSize = *pdfSrcYSize * dfSrcToDstYSize;
3447 :
3448 369 : if (*pdfDstXOff + *pdfDstXSize > nVRTXSize)
3449 : {
3450 3 : *pdfDstXSize = nVRTXSize - *pdfDstXOff;
3451 3 : *pdfSrcXSize = *pdfDstXSize / dfSrcToDstXSize;
3452 : }
3453 :
3454 369 : if (*pdfDstYOff + *pdfDstYSize > nVRTYSize)
3455 : {
3456 1 : *pdfDstYSize = nVRTYSize - *pdfDstYOff;
3457 1 : *pdfSrcYSize = *pdfDstYSize / dfSrcToDstYSize;
3458 : }
3459 :
3460 738 : return *pdfSrcXSize > 0 && *pdfDstXSize > 0 && *pdfSrcYSize > 0 &&
3461 738 : *pdfDstYSize > 0;
3462 : }
3463 :
3464 : /************************************************************************/
3465 : /* GDALDatasetCastToGTIDataset() */
3466 : /************************************************************************/
3467 :
3468 3 : GDALTileIndexDataset *GDALDatasetCastToGTIDataset(GDALDataset *poDS)
3469 : {
3470 3 : return dynamic_cast<GDALTileIndexDataset *>(poDS);
3471 : }
3472 :
3473 : /************************************************************************/
3474 : /* GTIGetSourcesMoreRecentThan() */
3475 : /************************************************************************/
3476 :
3477 : std::vector<GTISourceDesc>
3478 2 : GTIGetSourcesMoreRecentThan(GDALTileIndexDataset *poDS, int64_t mTime)
3479 : {
3480 2 : return poDS->GetSourcesMoreRecentThan(mTime);
3481 : }
3482 :
3483 : /************************************************************************/
3484 : /* GetSourcesMoreRecentThan() */
3485 : /************************************************************************/
3486 :
3487 : std::vector<GTISourceDesc>
3488 2 : GDALTileIndexDataset::GetSourcesMoreRecentThan(int64_t mTime)
3489 : {
3490 2 : std::vector<GTISourceDesc> oRes;
3491 :
3492 2 : m_poLayer->SetSpatialFilter(nullptr);
3493 6 : for (auto &&poFeature : m_poLayer)
3494 : {
3495 4 : if (!poFeature->IsFieldSetAndNotNull(m_nLocationFieldIndex))
3496 : {
3497 2 : continue;
3498 : }
3499 :
3500 4 : auto poGeom = poFeature->GetGeometryRef();
3501 4 : if (!poGeom || poGeom->IsEmpty())
3502 0 : continue;
3503 :
3504 4 : OGREnvelope sEnvelope;
3505 4 : poGeom->getEnvelope(&sEnvelope);
3506 :
3507 4 : double dfXOff = (sEnvelope.MinX - m_gt.xorig) / m_gt.xscale;
3508 4 : if (dfXOff >= nRasterXSize)
3509 0 : continue;
3510 :
3511 4 : double dfYOff = (sEnvelope.MaxY - m_gt.yorig) / m_gt.yscale;
3512 4 : if (dfYOff >= nRasterYSize)
3513 0 : continue;
3514 :
3515 4 : double dfXSize = (sEnvelope.MaxX - sEnvelope.MinX) / m_gt.xscale;
3516 4 : if (dfXOff < 0)
3517 : {
3518 0 : dfXSize += dfXOff;
3519 0 : dfXOff = 0;
3520 0 : if (dfXSize <= 0)
3521 0 : continue;
3522 : }
3523 :
3524 4 : double dfYSize =
3525 4 : (sEnvelope.MaxY - sEnvelope.MinY) / std::fabs(m_gt.yscale);
3526 4 : if (dfYOff < 0)
3527 : {
3528 0 : dfYSize += dfYOff;
3529 0 : dfYOff = 0;
3530 0 : if (dfYSize <= 0)
3531 0 : continue;
3532 : }
3533 :
3534 : const char *pszTileName =
3535 4 : poFeature->GetFieldAsString(m_nLocationFieldIndex);
3536 : std::string osTileName(GetAbsoluteFileName(
3537 4 : pszTileName, GetDescription(), m_bSTACCollection));
3538 : VSIStatBufL sStatSource;
3539 8 : if (VSIStatL(osTileName.c_str(), &sStatSource) != 0 ||
3540 4 : sStatSource.st_mtime <= mTime)
3541 : {
3542 2 : continue;
3543 : }
3544 :
3545 2 : constexpr double EPS = 1e-8;
3546 4 : GTISourceDesc oSourceDesc;
3547 2 : oSourceDesc.osFilename = std::move(osTileName);
3548 2 : oSourceDesc.nDstXOff = static_cast<int>(dfXOff + EPS);
3549 2 : oSourceDesc.nDstYOff = static_cast<int>(dfYOff + EPS);
3550 2 : oSourceDesc.nDstXSize = static_cast<int>(dfXSize + 0.5);
3551 2 : oSourceDesc.nDstYSize = static_cast<int>(dfYSize + 0.5);
3552 2 : oRes.emplace_back(std::move(oSourceDesc));
3553 : }
3554 :
3555 2 : return oRes;
3556 : }
3557 :
3558 : /************************************************************************/
3559 : /* GetSourceDesc() */
3560 : /************************************************************************/
3561 :
3562 373 : bool GDALTileIndexDataset::GetSourceDesc(const std::string &osTileName,
3563 : SourceDesc &oSourceDesc,
3564 : std::mutex *pMutex, int nBandCount,
3565 : const int *panBandMap)
3566 : {
3567 :
3568 373 : if (pMutex)
3569 138 : pMutex->lock();
3570 373 : std::shared_ptr<SharedDataset> sharedDS;
3571 746 : GTISharedSourceKey key;
3572 373 : key.osTileName = osTileName;
3573 373 : if (m_bBandInterleave)
3574 0 : key.anBands.insert(key.anBands.end(), panBandMap,
3575 9 : panBandMap + nBandCount);
3576 373 : m_oMapSharedSources.tryGet(key, sharedDS);
3577 373 : if (pMutex)
3578 138 : pMutex->unlock();
3579 :
3580 373 : std::shared_ptr<GDALDataset> poTileDS;
3581 373 : GDALDataset *poUnreprojectedDS = nullptr;
3582 373 : bool bBandMapTakenIntoAccount = false;
3583 :
3584 373 : if (sharedDS)
3585 : {
3586 126 : poTileDS = sharedDS->poDS;
3587 126 : poUnreprojectedDS = sharedDS->poUnreprojectedDS;
3588 126 : bBandMapTakenIntoAccount = sharedDS->bBandMapTakenIntoAccount;
3589 : }
3590 : else
3591 : {
3592 494 : poTileDS = std::shared_ptr<GDALDataset>(
3593 : GDALProxyPoolDataset::Create(
3594 : osTileName.c_str(), nullptr, GA_ReadOnly,
3595 : /* bShared = */ true, m_osUniqueHandle.c_str()),
3596 247 : GDALDatasetUniquePtrReleaser());
3597 247 : if (!poTileDS)
3598 : {
3599 3 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot open source %s",
3600 : osTileName.c_str());
3601 3 : return false;
3602 : }
3603 244 : poUnreprojectedDS = poTileDS.get();
3604 244 : if (poTileDS->GetRasterCount() == 0)
3605 : {
3606 0 : CPLError(CE_Failure, CPLE_AppDefined,
3607 : "Source %s has no raster bands", osTileName.c_str());
3608 0 : return false;
3609 : }
3610 :
3611 : // do palette -> RGB(A) expansion if needed
3612 244 : if (!GTIDoPaletteExpansionIfNeeded(poTileDS, nBands))
3613 0 : return false;
3614 :
3615 244 : bool bWarpVRT = false;
3616 244 : bool bExportSRS = false;
3617 244 : bool bAddAlphaToVRT = false;
3618 244 : const OGRSpatialReference *poTileSRS = poTileDS->GetSpatialRef();
3619 244 : GDALGeoTransform tileGT;
3620 426 : if (!m_oSRS.IsEmpty() && poTileSRS != nullptr &&
3621 182 : !m_oSRS.IsSame(poTileSRS))
3622 : {
3623 13 : CPLDebug("GTI",
3624 : "Tile %s has not the same SRS as the VRT. "
3625 : "Proceed to on-the-fly warping",
3626 : osTileName.c_str());
3627 13 : bWarpVRT = true;
3628 13 : bExportSRS = true;
3629 13 : if (poTileDS->GetRasterBand(poTileDS->GetRasterCount())
3630 21 : ->GetColorInterpretation() != GCI_AlphaBand &&
3631 8 : GetRasterBand(nBands)->GetColorInterpretation() ==
3632 : GCI_AlphaBand)
3633 : {
3634 0 : bAddAlphaToVRT = true;
3635 : }
3636 : }
3637 231 : else if (poTileDS->GetGeoTransform(tileGT) == CE_None &&
3638 232 : tileGT.yscale > 0 &&
3639 1 : ((m_oSRS.IsEmpty() && poTileSRS == nullptr) ||
3640 0 : (!m_oSRS.IsEmpty() && poTileSRS && m_oSRS.IsSame(poTileSRS))))
3641 :
3642 : {
3643 1 : CPLDebug("GTI",
3644 : "Tile %s is south-up oriented. "
3645 : "Proceed to on-the-fly warping",
3646 : osTileName.c_str());
3647 1 : bWarpVRT = true;
3648 : }
3649 :
3650 244 : if (bWarpVRT)
3651 : {
3652 14 : CPLStringList aosOptions;
3653 14 : aosOptions.AddString("-of");
3654 14 : aosOptions.AddString("VRT");
3655 :
3656 14 : if ((poTileDS->GetRasterBand(1)->GetColorTable() == nullptr &&
3657 14 : poTileDS->GetRasterBand(1)->GetCategoryNames() == nullptr) ||
3658 0 : m_eResampling == GRIORA_Mode)
3659 : {
3660 14 : aosOptions.AddString("-r");
3661 14 : aosOptions.AddString(m_osResampling.c_str());
3662 : }
3663 :
3664 14 : if (bExportSRS)
3665 : {
3666 13 : if (m_osWKT.empty())
3667 : {
3668 2 : char *pszWKT = nullptr;
3669 2 : const char *const apszWKTOptions[] = {"FORMAT=WKT2_2019",
3670 : nullptr};
3671 2 : m_oSRS.exportToWkt(&pszWKT, apszWKTOptions);
3672 2 : if (pszWKT)
3673 2 : m_osWKT = pszWKT;
3674 2 : CPLFree(pszWKT);
3675 :
3676 2 : if (m_osWKT.empty())
3677 : {
3678 0 : CPLError(CE_Failure, CPLE_AppDefined,
3679 : "Cannot export VRT SRS to WKT2");
3680 0 : return false;
3681 : }
3682 : }
3683 :
3684 13 : aosOptions.AddString("-t_srs");
3685 13 : aosOptions.AddString(m_osWKT.c_str());
3686 : }
3687 :
3688 : // First pass to get the extent of the tile in the
3689 : // target VRT SRS
3690 : GDALWarpAppOptions *psWarpOptions =
3691 14 : GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
3692 14 : GDALDatasetH ahSrcDS[] = {GDALDataset::ToHandle(poTileDS.get())};
3693 14 : int bUsageError = false;
3694 : auto poWarpDS =
3695 : std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(GDALWarp(
3696 14 : "", nullptr, 1, ahSrcDS, psWarpOptions, &bUsageError)));
3697 14 : GDALWarpAppOptionsFree(psWarpOptions);
3698 14 : if (!poWarpDS)
3699 : {
3700 0 : return false;
3701 : }
3702 :
3703 : // Second pass to create a warped source VRT whose
3704 : // extent is aligned on the one of the target GTI
3705 14 : GDALGeoTransform warpDSGT;
3706 14 : const auto eErr = poWarpDS->GetGeoTransform(warpDSGT);
3707 14 : CPL_IGNORE_RET_VAL(eErr);
3708 14 : CPLAssert(eErr == CE_None);
3709 14 : const double dfGTIMinX = m_gt.xorig;
3710 14 : const double dfGTIResX = m_gt.xscale;
3711 14 : const double dfGTIMaxY = m_gt.yorig;
3712 14 : const double dfGTIResYAbs = -m_gt.yscale;
3713 14 : const double dfWarpMinX =
3714 14 : std::floor((warpDSGT.xorig - dfGTIMinX) / dfGTIResX) *
3715 : dfGTIResX +
3716 : dfGTIMinX;
3717 : const double dfWarpMaxX =
3718 28 : std::ceil((warpDSGT.xorig +
3719 14 : warpDSGT.xscale * poWarpDS->GetRasterXSize() -
3720 : dfGTIMinX) /
3721 14 : dfGTIResX) *
3722 : dfGTIResX +
3723 14 : dfGTIMinX;
3724 14 : const double dfWarpMaxY =
3725 : dfGTIMaxY -
3726 14 : std::floor((dfGTIMaxY - warpDSGT.yorig) / dfGTIResYAbs) *
3727 : dfGTIResYAbs;
3728 : const double dfWarpMinY =
3729 : dfGTIMaxY -
3730 14 : std::ceil((dfGTIMaxY -
3731 28 : (warpDSGT.yorig +
3732 14 : warpDSGT.yscale * poWarpDS->GetRasterYSize())) /
3733 14 : dfGTIResYAbs) *
3734 14 : dfGTIResYAbs;
3735 :
3736 14 : aosOptions.AddString("-te");
3737 14 : aosOptions.AddString(CPLSPrintf("%.17g", dfWarpMinX));
3738 14 : aosOptions.AddString(CPLSPrintf("%.17g", dfWarpMinY));
3739 14 : aosOptions.AddString(CPLSPrintf("%.17g", dfWarpMaxX));
3740 14 : aosOptions.AddString(CPLSPrintf("%.17g", dfWarpMaxY));
3741 :
3742 14 : aosOptions.AddString("-tr");
3743 14 : aosOptions.AddString(CPLSPrintf("%.17g", dfGTIResX));
3744 14 : aosOptions.AddString(CPLSPrintf("%.17g", dfGTIResYAbs));
3745 :
3746 14 : if (m_bBandInterleave)
3747 : {
3748 9 : bBandMapTakenIntoAccount = true;
3749 23 : for (int i = 0; i < nBandCount; ++i)
3750 : {
3751 14 : aosOptions.AddString("-b");
3752 14 : aosOptions.AddString(CPLSPrintf("%d", panBandMap[i]));
3753 : }
3754 : }
3755 :
3756 14 : if (bAddAlphaToVRT)
3757 0 : aosOptions.AddString("-dstalpha");
3758 :
3759 14 : if (!m_osWarpMemory.empty())
3760 : {
3761 0 : aosOptions.AddString("-wm");
3762 0 : aosOptions.AddString(m_osWarpMemory.c_str());
3763 : }
3764 :
3765 14 : psWarpOptions = GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
3766 14 : poWarpDS.reset(GDALDataset::FromHandle(GDALWarp(
3767 : "", nullptr, 1, ahSrcDS, psWarpOptions, &bUsageError)));
3768 14 : GDALWarpAppOptionsFree(psWarpOptions);
3769 14 : if (!poWarpDS)
3770 : {
3771 0 : return false;
3772 : }
3773 :
3774 14 : poTileDS = std::move(poWarpDS);
3775 : }
3776 :
3777 244 : sharedDS = std::make_shared<SharedDataset>();
3778 244 : sharedDS->poDS = poTileDS;
3779 244 : sharedDS->poUnreprojectedDS = poUnreprojectedDS;
3780 244 : sharedDS->bBandMapTakenIntoAccount = bBandMapTakenIntoAccount;
3781 :
3782 244 : if (pMutex)
3783 70 : pMutex->lock();
3784 244 : m_oMapSharedSources.insert(key, sharedDS);
3785 244 : if (pMutex)
3786 70 : pMutex->unlock();
3787 : }
3788 :
3789 370 : GDALGeoTransform gtTile;
3790 370 : if (poTileDS->GetGeoTransform(gtTile) != CE_None)
3791 : {
3792 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s lacks geotransform",
3793 : osTileName.c_str());
3794 0 : return false;
3795 : }
3796 :
3797 370 : bool bHasNoData = false;
3798 370 : bool bSameNoData = true;
3799 370 : double dfNoDataValue = 0;
3800 370 : GDALRasterBand *poMaskBand = nullptr;
3801 370 : const int nTileBandCount = poTileDS->GetRasterCount();
3802 1265 : for (int iBand = 0; iBand < nTileBandCount; ++iBand)
3803 : {
3804 895 : auto poTileBand = poTileDS->GetRasterBand(iBand + 1);
3805 895 : int bThisBandHasNoData = false;
3806 : const double dfThisBandNoDataValue =
3807 895 : poTileBand->GetNoDataValue(&bThisBandHasNoData);
3808 895 : if (bThisBandHasNoData)
3809 : {
3810 24 : bHasNoData = true;
3811 24 : dfNoDataValue = dfThisBandNoDataValue;
3812 : }
3813 1420 : if (iBand > 0 &&
3814 525 : (static_cast<int>(bThisBandHasNoData) !=
3815 525 : static_cast<int>(bHasNoData) ||
3816 12 : (bHasNoData &&
3817 12 : !IsSameNaNAware(dfNoDataValue, dfThisBandNoDataValue))))
3818 : {
3819 0 : bSameNoData = false;
3820 : }
3821 :
3822 895 : if (poTileBand->GetMaskFlags() == GMF_PER_DATASET)
3823 2 : poMaskBand = poTileBand->GetMaskBand();
3824 893 : else if (poTileBand->GetColorInterpretation() == GCI_AlphaBand)
3825 36 : poMaskBand = poTileBand;
3826 : }
3827 :
3828 0 : std::unique_ptr<VRTSimpleSource> poSource;
3829 370 : if (!bHasNoData)
3830 : {
3831 358 : poSource = std::make_unique<VRTSimpleSource>();
3832 : }
3833 : else
3834 : {
3835 24 : auto poComplexSource = std::make_unique<VRTComplexSource>();
3836 12 : poComplexSource->SetNoDataValue(dfNoDataValue);
3837 12 : poSource = std::move(poComplexSource);
3838 : }
3839 :
3840 740 : GetSrcDstWin(gtTile, poTileDS->GetRasterXSize(), poTileDS->GetRasterYSize(),
3841 370 : m_gt, GetRasterXSize(), GetRasterYSize(),
3842 370 : &poSource->m_dfSrcXOff, &poSource->m_dfSrcYOff,
3843 370 : &poSource->m_dfSrcXSize, &poSource->m_dfSrcYSize,
3844 370 : &poSource->m_dfDstXOff, &poSource->m_dfDstYOff,
3845 370 : &poSource->m_dfDstXSize, &poSource->m_dfDstYSize);
3846 :
3847 370 : oSourceDesc.osName = osTileName;
3848 370 : oSourceDesc.poDS = std::move(poTileDS);
3849 370 : oSourceDesc.poUnreprojectedDS = poUnreprojectedDS;
3850 370 : oSourceDesc.bBandMapTakenIntoAccount = bBandMapTakenIntoAccount;
3851 370 : oSourceDesc.poSource = std::move(poSource);
3852 370 : oSourceDesc.bHasNoData = bHasNoData;
3853 370 : oSourceDesc.bSameNoData = bSameNoData;
3854 370 : if (bSameNoData)
3855 370 : oSourceDesc.dfSameNoData = dfNoDataValue;
3856 370 : oSourceDesc.poMaskBand = poMaskBand;
3857 370 : return true;
3858 : }
3859 :
3860 : /************************************************************************/
3861 : /* GetNumThreads() */
3862 : /************************************************************************/
3863 :
3864 8 : int GDALTileIndexDataset::GetNumThreads() const
3865 : {
3866 : const char *pszNumThreads =
3867 8 : CSLFetchNameValueDef(GetOpenOptions(), "NUM_THREADS", nullptr);
3868 8 : if (!pszNumThreads)
3869 8 : pszNumThreads = CPLGetConfigOption("GTI_NUM_THREADS", nullptr);
3870 8 : return GDALGetNumThreads(pszNumThreads, GDALGetMaxDatasetPoolSize(),
3871 8 : /* bDefaultAllCPUs = */ true);
3872 : }
3873 :
3874 : /************************************************************************/
3875 : /* CollectSources() */
3876 : /************************************************************************/
3877 :
3878 202 : bool GDALTileIndexDataset::CollectSources(double dfXOff, double dfYOff,
3879 : double dfXSize, double dfYSize,
3880 : int nBandCount, const int *panBandMap,
3881 : bool bMultiThreadAllowed)
3882 : {
3883 202 : const double dfMinX = m_gt.xorig + dfXOff * m_gt.xscale;
3884 202 : const double dfMaxX = dfMinX + dfXSize * m_gt.xscale;
3885 202 : const double dfMaxY = m_gt.yorig + dfYOff * m_gt.yscale;
3886 202 : const double dfMinY = dfMaxY + dfYSize * m_gt.yscale;
3887 :
3888 74 : if (dfMinX == m_dfLastMinXFilter && dfMinY == m_dfLastMinYFilter &&
3889 335 : dfMaxX == m_dfLastMaxXFilter && dfMaxY == m_dfLastMaxYFilter &&
3890 59 : (!m_bBandInterleave ||
3891 6 : (m_anLastBands ==
3892 208 : std::vector<int>(panBandMap, panBandMap + nBandCount))))
3893 : {
3894 54 : return true;
3895 : }
3896 :
3897 148 : m_dfLastMinXFilter = dfMinX;
3898 148 : m_dfLastMinYFilter = dfMinY;
3899 148 : m_dfLastMaxXFilter = dfMaxX;
3900 148 : m_dfLastMaxYFilter = dfMaxY;
3901 148 : if (m_bBandInterleave)
3902 9 : m_anLastBands = std::vector<int>(panBandMap, panBandMap + nBandCount);
3903 148 : m_bLastMustUseMultiThreading = false;
3904 :
3905 148 : OGRLayer *poSQLLayer = nullptr;
3906 148 : if (!m_osSpatialSQL.empty())
3907 : {
3908 : const std::string osSQL =
3909 3 : CPLString(m_osSpatialSQL)
3910 6 : .replaceAll("{XMIN}", CPLSPrintf("%.17g", dfMinX))
3911 6 : .replaceAll("{YMIN}", CPLSPrintf("%.17g", dfMinY))
3912 6 : .replaceAll("{XMAX}", CPLSPrintf("%.17g", dfMaxX))
3913 6 : .replaceAll("{YMAX}", CPLSPrintf("%.17g", dfMaxY));
3914 3 : poSQLLayer = m_poVectorDS->ExecuteSQL(osSQL.c_str(), nullptr, nullptr);
3915 3 : if (!poSQLLayer)
3916 1 : return 0;
3917 : }
3918 : else
3919 : {
3920 145 : m_poLayer->SetSpatialFilterRect(dfMinX, dfMinY, dfMaxX, dfMaxY);
3921 145 : m_poLayer->ResetReading();
3922 : }
3923 :
3924 147 : OGRLayer *const poLayer = poSQLLayer ? poSQLLayer : m_poLayer;
3925 :
3926 147 : m_aoSourceDesc.clear();
3927 : while (true)
3928 : {
3929 487 : auto poFeature = std::unique_ptr<OGRFeature>(poLayer->GetNextFeature());
3930 487 : if (!poFeature)
3931 147 : break;
3932 340 : if (!poFeature->IsFieldSetAndNotNull(m_nLocationFieldIndex))
3933 : {
3934 1 : continue;
3935 : }
3936 :
3937 339 : SourceDesc oSourceDesc;
3938 339 : oSourceDesc.poFeature = std::move(poFeature);
3939 339 : m_aoSourceDesc.emplace_back(std::move(oSourceDesc));
3940 :
3941 339 : if (m_aoSourceDesc.size() > 10 * 1000 * 1000)
3942 : {
3943 : // Safety belt...
3944 0 : CPLError(CE_Failure, CPLE_AppDefined,
3945 : "More than 10 million contributing sources to a "
3946 : "single RasterIO() request is not supported");
3947 0 : return false;
3948 : }
3949 340 : }
3950 :
3951 147 : if (poSQLLayer)
3952 2 : ReleaseResultSet(poSQLLayer);
3953 :
3954 147 : constexpr int MINIMUM_PIXEL_COUNT_FOR_THREADED_IO = 1000 * 1000;
3955 216 : if (bMultiThreadAllowed && m_aoSourceDesc.size() > 1 &&
3956 69 : dfXSize * dfYSize > MINIMUM_PIXEL_COUNT_FOR_THREADED_IO)
3957 : {
3958 8 : if (m_nNumThreads < 0)
3959 8 : m_nNumThreads = GetNumThreads();
3960 8 : bMultiThreadAllowed = m_nNumThreads > 1;
3961 : }
3962 : else
3963 : {
3964 139 : bMultiThreadAllowed = false;
3965 : }
3966 :
3967 147 : if (bMultiThreadAllowed)
3968 : {
3969 : CPLRectObj sGlobalBounds;
3970 5 : sGlobalBounds.minx = dfXOff;
3971 5 : sGlobalBounds.miny = dfYOff;
3972 5 : sGlobalBounds.maxx = dfXOff + dfXSize;
3973 5 : sGlobalBounds.maxy = dfYOff + dfYSize;
3974 5 : CPLQuadTree *hQuadTree = CPLQuadTreeCreate(&sGlobalBounds, nullptr);
3975 :
3976 5 : bool bCompatibleOfMultiThread = true;
3977 5 : std::set<std::string> oSetTileNames;
3978 77 : for (const auto &oSourceDesc : m_aoSourceDesc)
3979 : {
3980 : const char *pszTileName =
3981 73 : oSourceDesc.poFeature->GetFieldAsString(m_nLocationFieldIndex);
3982 73 : if (oSetTileNames.find(pszTileName) != oSetTileNames.end())
3983 : {
3984 0 : bCompatibleOfMultiThread = false;
3985 1 : break;
3986 : }
3987 73 : oSetTileNames.insert(pszTileName);
3988 :
3989 73 : const auto poGeom = oSourceDesc.poFeature->GetGeometryRef();
3990 73 : if (!poGeom || poGeom->IsEmpty())
3991 0 : continue;
3992 :
3993 73 : OGREnvelope sEnvelope;
3994 73 : poGeom->getEnvelope(&sEnvelope);
3995 :
3996 : CPLRectObj sSourceBounds;
3997 73 : sSourceBounds.minx = (sEnvelope.MinX - m_gt.xorig) / m_gt.xscale;
3998 73 : sSourceBounds.maxx = (sEnvelope.MaxX - m_gt.xorig) / m_gt.xscale;
3999 : // Yes use of MaxY to compute miny is intended given that MaxY is
4000 : // in georeferenced space whereas miny is in pixel space.
4001 73 : sSourceBounds.miny = (sEnvelope.MaxY - m_gt.yorig) / m_gt.yscale;
4002 : // Same here for maxy vs Miny
4003 73 : sSourceBounds.maxy = (sEnvelope.MinY - m_gt.yorig) / m_gt.yscale;
4004 :
4005 : // Clamp to global bounds and some epsilon to avoid adjacent tiles
4006 : // to be considered as overlapping
4007 73 : constexpr double EPSILON = 0.1;
4008 73 : sSourceBounds.minx =
4009 73 : std::max(sGlobalBounds.minx, sSourceBounds.minx) + EPSILON;
4010 73 : sSourceBounds.maxx =
4011 73 : std::min(sGlobalBounds.maxx, sSourceBounds.maxx) - EPSILON;
4012 73 : sSourceBounds.miny =
4013 73 : std::max(sGlobalBounds.miny, sSourceBounds.miny) + EPSILON;
4014 73 : sSourceBounds.maxy =
4015 73 : std::min(sGlobalBounds.maxy, sSourceBounds.maxy) - EPSILON;
4016 :
4017 : // Check that the new source doesn't overlap an existing one.
4018 73 : if (CPLQuadTreeHasMatch(hQuadTree, &sSourceBounds))
4019 : {
4020 1 : bCompatibleOfMultiThread = false;
4021 1 : break;
4022 : }
4023 :
4024 72 : CPLQuadTreeInsertWithBounds(
4025 : hQuadTree,
4026 : const_cast<void *>(static_cast<const void *>(&oSourceDesc)),
4027 : &sSourceBounds);
4028 : }
4029 :
4030 5 : CPLQuadTreeDestroy(hQuadTree);
4031 :
4032 5 : if (bCompatibleOfMultiThread)
4033 : {
4034 4 : m_bLastMustUseMultiThreading = true;
4035 4 : return true;
4036 : }
4037 : }
4038 :
4039 143 : if (m_aoSourceDesc.size() > 1)
4040 : {
4041 66 : SortSourceDesc();
4042 : }
4043 :
4044 : // Try to find the last (most priority) fully opaque source covering
4045 : // the whole AOI. We only need to start rendering from it.
4046 143 : size_t i = m_aoSourceDesc.size();
4047 284 : while (i > 0)
4048 : {
4049 235 : --i;
4050 235 : auto &poFeature = m_aoSourceDesc[i].poFeature;
4051 : const char *pszTileName =
4052 235 : poFeature->GetFieldAsString(m_nLocationFieldIndex);
4053 : const std::string osTileName(GetAbsoluteFileName(
4054 235 : pszTileName, GetDescription(), m_bSTACCollection));
4055 :
4056 235 : SourceDesc oSourceDesc;
4057 235 : if (!GetSourceDesc(osTileName, oSourceDesc, nullptr, nBandCount,
4058 : panBandMap))
4059 2 : return false;
4060 :
4061 : // Check consistency of bounding box in tile index vs actual
4062 : // extent of the tile.
4063 233 : GDALGeoTransform tileGT;
4064 233 : if (oSourceDesc.poDS->GetGeoTransform(tileGT) == CE_None &&
4065 233 : tileGT.xrot == 0 && tileGT.yrot == 0)
4066 : {
4067 233 : OGREnvelope sActualTileExtent;
4068 233 : sActualTileExtent.MinX = tileGT.xorig;
4069 233 : sActualTileExtent.MaxX =
4070 466 : sActualTileExtent.MinX +
4071 233 : oSourceDesc.poDS->GetRasterXSize() * tileGT.xscale;
4072 233 : sActualTileExtent.MaxY = tileGT.yorig;
4073 233 : sActualTileExtent.MinY =
4074 466 : sActualTileExtent.MaxY +
4075 233 : oSourceDesc.poDS->GetRasterYSize() * tileGT.yscale;
4076 233 : const auto poGeom = poFeature->GetGeometryRef();
4077 233 : if (poGeom && !poGeom->IsEmpty())
4078 : {
4079 233 : OGREnvelope sGeomTileExtent;
4080 233 : poGeom->getEnvelope(&sGeomTileExtent);
4081 233 : sGeomTileExtent.MinX -= m_gt.xscale;
4082 233 : sGeomTileExtent.MaxX += m_gt.xscale;
4083 233 : sGeomTileExtent.MinY -= std::fabs(m_gt.yscale);
4084 233 : sGeomTileExtent.MaxY += std::fabs(m_gt.yscale);
4085 233 : if (!sGeomTileExtent.Contains(sActualTileExtent))
4086 : {
4087 3 : if (!sGeomTileExtent.Intersects(sActualTileExtent))
4088 : {
4089 1 : CPLError(CE_Warning, CPLE_AppDefined,
4090 : "Tile index is out of sync with actual "
4091 : "extent of %s. Bounding box from tile index "
4092 : "is (%.15g, %.15g, %.15g, %.15g) does not "
4093 : "intersect at all bounding box from tile "
4094 : "(%.15g, %.15g, %.15g, %.15g)",
4095 : osTileName.c_str(), sGeomTileExtent.MinX,
4096 : sGeomTileExtent.MinY, sGeomTileExtent.MaxX,
4097 : sGeomTileExtent.MaxY, sActualTileExtent.MinX,
4098 : sActualTileExtent.MinY, sActualTileExtent.MaxX,
4099 : sActualTileExtent.MaxY);
4100 2 : continue;
4101 : }
4102 :
4103 : // The above test assumes, in the case of reprojection, that
4104 : // the reprojected geometry in the index is computed the
4105 : // same way as we do here, that is using GDALWarp()
4106 : // Which wasn't the case for example before GDAL 3.12.3 when
4107 : // using "gdal raster index", which uses a simple 4-corner
4108 : // reprojection logic. So also test using that method,
4109 : // before emitting any warning.
4110 2 : if (oSourceDesc.poUnreprojectedDS != oSourceDesc.poDS.get())
4111 : {
4112 : const int nXSize =
4113 1 : oSourceDesc.poUnreprojectedDS->GetRasterXSize();
4114 : const int nYSize =
4115 1 : oSourceDesc.poUnreprojectedDS->GetRasterYSize();
4116 1 : GDALGeoTransform gt;
4117 : const auto poSrcSRS =
4118 1 : oSourceDesc.poUnreprojectedDS->GetSpatialRef();
4119 2 : if (poSrcSRS && !m_oSRS.IsEmpty() &&
4120 1 : oSourceDesc.poUnreprojectedDS->GetGeoTransform(
4121 1 : gt) == CE_None)
4122 : {
4123 1 : double adfX[4] = {0.0, 0.0, 0.0, 0.0};
4124 1 : double adfY[4] = {0.0, 0.0, 0.0, 0.0};
4125 1 : adfX[0] = gt.xorig + 0 * gt.xscale + 0 * gt.xrot;
4126 1 : adfY[0] = gt.yorig + 0 * gt.yrot + 0 * gt.yscale;
4127 :
4128 1 : adfX[1] =
4129 1 : gt.xorig + nXSize * gt.xscale + 0 * gt.xrot;
4130 1 : adfY[1] =
4131 1 : gt.yorig + nXSize * gt.yrot + 0 * gt.yscale;
4132 :
4133 1 : adfX[2] = gt.xorig + nXSize * gt.xscale +
4134 1 : nYSize * gt.xrot;
4135 1 : adfY[2] = gt.yorig + nXSize * gt.yrot +
4136 1 : nYSize * gt.yscale;
4137 :
4138 1 : adfX[3] =
4139 1 : gt.xorig + 0 * gt.xscale + nYSize * gt.xrot;
4140 1 : adfY[3] =
4141 1 : gt.yorig + 0 * gt.yrot + nYSize * gt.yscale;
4142 :
4143 : auto poCT =
4144 : std::unique_ptr<OGRCoordinateTransformation>(
4145 : OGRCreateCoordinateTransformation(poSrcSRS,
4146 1 : &m_oSRS));
4147 1 : if (poCT && poCT->Transform(4, adfX, adfY, nullptr))
4148 : {
4149 1 : OGREnvelope sActualTileExtent2;
4150 0 : sActualTileExtent2.MinX = std::min(
4151 1 : {adfX[0], adfX[1], adfX[2], adfX[3]});
4152 0 : sActualTileExtent2.MinY = std::min(
4153 1 : {adfY[0], adfY[1], adfY[2], adfY[3]});
4154 0 : sActualTileExtent2.MaxX = std::max(
4155 1 : {adfX[0], adfX[1], adfX[2], adfX[3]});
4156 0 : sActualTileExtent2.MaxY = std::max(
4157 1 : {adfY[0], adfY[1], adfY[2], adfY[3]});
4158 1 : if (sGeomTileExtent.Contains(
4159 1 : sActualTileExtent2))
4160 : {
4161 1 : continue;
4162 : }
4163 : }
4164 : }
4165 : }
4166 :
4167 1 : CPLError(CE_Warning, CPLE_AppDefined,
4168 : "Tile index is out of sync with actual extent "
4169 : "of %s. Bounding box from tile index is "
4170 : "(%.15g, %.15g, %.15g, %.15g) does not fully "
4171 : "contain bounding box from tile "
4172 : "(%.15g, %.15g, %.15g, %.15g)",
4173 : osTileName.c_str(), sGeomTileExtent.MinX,
4174 : sGeomTileExtent.MinY, sGeomTileExtent.MaxX,
4175 : sGeomTileExtent.MaxY, sActualTileExtent.MinX,
4176 : sActualTileExtent.MinY, sActualTileExtent.MaxX,
4177 : sActualTileExtent.MaxY);
4178 : }
4179 : }
4180 : }
4181 :
4182 231 : const auto &poSource = oSourceDesc.poSource;
4183 231 : if (dfXOff >= poSource->m_dfDstXOff + poSource->m_dfDstXSize ||
4184 231 : dfYOff >= poSource->m_dfDstYOff + poSource->m_dfDstYSize ||
4185 690 : poSource->m_dfDstXOff >= dfXOff + dfXSize ||
4186 228 : poSource->m_dfDstYOff >= dfYOff + dfYSize)
4187 : {
4188 : // Can happen as some spatial filters select slightly more features
4189 : // than strictly needed.
4190 3 : continue;
4191 : }
4192 :
4193 : const bool bCoversWholeAOI =
4194 228 : (poSource->m_dfDstXOff <= dfXOff &&
4195 151 : poSource->m_dfDstYOff <= dfYOff &&
4196 149 : poSource->m_dfDstXOff + poSource->m_dfDstXSize >=
4197 515 : dfXOff + dfXSize &&
4198 136 : poSource->m_dfDstYOff + poSource->m_dfDstYSize >=
4199 136 : dfYOff + dfYSize);
4200 228 : oSourceDesc.bCoversWholeAOI = bCoversWholeAOI;
4201 :
4202 228 : m_aoSourceDesc[i] = std::move(oSourceDesc);
4203 :
4204 228 : if (m_aoSourceDesc[i].bCoversWholeAOI &&
4205 228 : !m_aoSourceDesc[i].bHasNoData && !m_aoSourceDesc[i].poMaskBand)
4206 : {
4207 92 : break;
4208 : }
4209 : }
4210 :
4211 141 : if (i > 0)
4212 : {
4213 : // Remove sources that will not be rendered
4214 32 : m_aoSourceDesc.erase(m_aoSourceDesc.begin(),
4215 64 : m_aoSourceDesc.begin() + i);
4216 : }
4217 :
4218 : // Remove elements that have no dataset
4219 0 : m_aoSourceDesc.erase(std::remove_if(m_aoSourceDesc.begin(),
4220 : m_aoSourceDesc.end(),
4221 233 : [](const SourceDesc &desc)
4222 374 : { return desc.poDS == nullptr; }),
4223 282 : m_aoSourceDesc.end());
4224 :
4225 141 : return true;
4226 : }
4227 :
4228 : /************************************************************************/
4229 : /* SortSourceDesc() */
4230 : /************************************************************************/
4231 :
4232 66 : void GDALTileIndexDataset::SortSourceDesc()
4233 : {
4234 66 : const auto eFieldType = m_nSortFieldIndex >= 0
4235 66 : ? m_poLayer->GetLayerDefn()
4236 47 : ->GetFieldDefn(m_nSortFieldIndex)
4237 47 : ->GetType()
4238 66 : : OFTMaxType;
4239 66 : std::sort(
4240 : m_aoSourceDesc.begin(), m_aoSourceDesc.end(),
4241 1828 : [this, eFieldType](const SourceDesc &a, const SourceDesc &b)
4242 : {
4243 419 : const auto &poFeatureA = (m_bSortFieldAsc ? a : b).poFeature;
4244 419 : const auto &poFeatureB = (m_bSortFieldAsc ? b : a).poFeature;
4245 918 : if (m_nSortFieldIndex >= 0 &&
4246 499 : poFeatureA->IsFieldSetAndNotNull(m_nSortFieldIndex) &&
4247 80 : poFeatureB->IsFieldSetAndNotNull(m_nSortFieldIndex))
4248 : {
4249 80 : if (eFieldType == OFTString)
4250 : {
4251 : const int nCmp =
4252 5 : strcmp(poFeatureA->GetFieldAsString(m_nSortFieldIndex),
4253 : poFeatureB->GetFieldAsString(m_nSortFieldIndex));
4254 5 : if (nCmp < 0)
4255 1 : return true;
4256 4 : if (nCmp > 0)
4257 2 : return false;
4258 : }
4259 75 : else if (eFieldType == OFTInteger || eFieldType == OFTInteger64)
4260 : {
4261 : const auto nA =
4262 45 : poFeatureA->GetFieldAsInteger64(m_nSortFieldIndex);
4263 : const auto nB =
4264 45 : poFeatureB->GetFieldAsInteger64(m_nSortFieldIndex);
4265 45 : if (nA < nB)
4266 3 : return true;
4267 42 : if (nA > nB)
4268 42 : return false;
4269 : }
4270 30 : else if (eFieldType == OFTReal)
4271 : {
4272 : const auto dfA =
4273 3 : poFeatureA->GetFieldAsDouble(m_nSortFieldIndex);
4274 : const auto dfB =
4275 3 : poFeatureB->GetFieldAsDouble(m_nSortFieldIndex);
4276 3 : if (dfA < dfB)
4277 1 : return true;
4278 2 : if (dfA > dfB)
4279 2 : return false;
4280 : }
4281 27 : else if (eFieldType == OFTDate || eFieldType == OFTDateTime)
4282 : {
4283 : const auto poFieldA =
4284 27 : poFeatureA->GetRawFieldRef(m_nSortFieldIndex);
4285 : const auto poFieldB =
4286 27 : poFeatureB->GetRawFieldRef(m_nSortFieldIndex);
4287 :
4288 : #define COMPARE_DATE_COMPONENT(comp) \
4289 : do \
4290 : { \
4291 : if (poFieldA->Date.comp < poFieldB->Date.comp) \
4292 : return true; \
4293 : if (poFieldA->Date.comp > poFieldB->Date.comp) \
4294 : return false; \
4295 : } while (0)
4296 :
4297 27 : COMPARE_DATE_COMPONENT(Year);
4298 21 : COMPARE_DATE_COMPONENT(Month);
4299 15 : COMPARE_DATE_COMPONENT(Day);
4300 9 : COMPARE_DATE_COMPONENT(Hour);
4301 8 : COMPARE_DATE_COMPONENT(Minute);
4302 7 : COMPARE_DATE_COMPONENT(Second);
4303 : }
4304 : else
4305 : {
4306 0 : CPLAssert(false);
4307 : }
4308 : }
4309 347 : return poFeatureA->GetFID() < poFeatureB->GetFID();
4310 : });
4311 66 : }
4312 :
4313 : /************************************************************************/
4314 : /* CompositeSrcWithMaskIntoDest() */
4315 : /************************************************************************/
4316 :
4317 : static void
4318 66 : CompositeSrcWithMaskIntoDest(const int nOutXSize, const int nOutYSize,
4319 : const GDALDataType eBufType,
4320 : const int nBufTypeSize, const GSpacing nPixelSpace,
4321 : const GSpacing nLineSpace, const GByte *pabySrc,
4322 : const GByte *const pabyMask, GByte *const pabyDest)
4323 : {
4324 66 : size_t iMaskIdx = 0;
4325 66 : if (eBufType == GDT_UInt8)
4326 : {
4327 : // Optimization for byte case
4328 136 : for (int iY = 0; iY < nOutYSize; iY++)
4329 : {
4330 86 : GByte *pabyDestLine =
4331 86 : pabyDest + static_cast<GPtrDiff_t>(iY * nLineSpace);
4332 86 : int iX = 0;
4333 : #ifdef USE_SSE2_OPTIM
4334 86 : if (nPixelSpace == 1)
4335 : {
4336 : // SSE2 version up to 6 times faster than portable version
4337 86 : const auto xmm_zero = _mm_setzero_si128();
4338 86 : constexpr int SIZEOF_REG = static_cast<int>(sizeof(xmm_zero));
4339 110 : for (; iX + SIZEOF_REG <= nOutXSize; iX += SIZEOF_REG)
4340 : {
4341 48 : auto xmm_mask = _mm_loadu_si128(
4342 : reinterpret_cast<__m128i const *>(pabyMask + iMaskIdx));
4343 24 : const auto xmm_src = _mm_loadu_si128(
4344 : reinterpret_cast<__m128i const *>(pabySrc));
4345 24 : auto xmm_dst = _mm_loadu_si128(
4346 : reinterpret_cast<__m128i const *>(pabyDestLine));
4347 : #ifdef USE_SSE41_OPTIM
4348 : xmm_dst = _mm_blendv_epi8(xmm_dst, xmm_src, xmm_mask);
4349 : #else
4350 : // mask[i] = 0 becomes 255, and mask[i] != 0 becomes 0
4351 24 : xmm_mask = _mm_cmpeq_epi8(xmm_mask, xmm_zero);
4352 : // dst_data[i] = (mask[i] & dst_data[i]) |
4353 : // (~mask[i] & src_data[i])
4354 : // That is:
4355 : // dst_data[i] = dst_data[i] when mask[i] = 255
4356 : // dst_data[i] = src_data[i] when mask[i] = 0
4357 72 : xmm_dst = _mm_or_si128(_mm_and_si128(xmm_mask, xmm_dst),
4358 : _mm_andnot_si128(xmm_mask, xmm_src));
4359 : #endif
4360 : _mm_storeu_si128(reinterpret_cast<__m128i *>(pabyDestLine),
4361 : xmm_dst);
4362 24 : pabyDestLine += SIZEOF_REG;
4363 24 : pabySrc += SIZEOF_REG;
4364 24 : iMaskIdx += SIZEOF_REG;
4365 : }
4366 : }
4367 : #endif
4368 342 : for (; iX < nOutXSize; iX++)
4369 : {
4370 256 : if (pabyMask[iMaskIdx])
4371 : {
4372 218 : *pabyDestLine = *pabySrc;
4373 : }
4374 256 : pabyDestLine += static_cast<GPtrDiff_t>(nPixelSpace);
4375 256 : pabySrc++;
4376 256 : iMaskIdx++;
4377 : }
4378 : }
4379 : }
4380 : else
4381 : {
4382 38 : for (int iY = 0; iY < nOutYSize; iY++)
4383 : {
4384 22 : GByte *pabyDestLine =
4385 22 : pabyDest + static_cast<GPtrDiff_t>(iY * nLineSpace);
4386 54 : for (int iX = 0; iX < nOutXSize; iX++)
4387 : {
4388 32 : if (pabyMask[iMaskIdx])
4389 : {
4390 16 : memcpy(pabyDestLine, pabySrc, nBufTypeSize);
4391 : }
4392 32 : pabyDestLine += static_cast<GPtrDiff_t>(nPixelSpace);
4393 32 : pabySrc += nBufTypeSize;
4394 32 : iMaskIdx++;
4395 : }
4396 : }
4397 : }
4398 66 : }
4399 :
4400 : /************************************************************************/
4401 : /* NeedInitBuffer() */
4402 : /************************************************************************/
4403 :
4404 : // Must be called after CollectSources()
4405 190 : bool GDALTileIndexDataset::NeedInitBuffer(int nBandCount,
4406 : const int *panBandMap) const
4407 : {
4408 190 : bool bNeedInitBuffer = true;
4409 : // If the last source (that is the most priority one) covers at least
4410 : // the window of interest and is fully opaque, then we don't need to
4411 : // initialize the buffer, and can directly render that source.
4412 190 : int bHasNoData = false;
4413 376 : if (!m_aoSourceDesc.empty() && m_aoSourceDesc.back().bCoversWholeAOI &&
4414 173 : (!m_aoSourceDesc.back().bHasNoData ||
4415 : // Also, if there's a single source and that the VRT bands and the
4416 : // source bands have the same nodata value, we can skip initialization.
4417 12 : (m_aoSourceDesc.size() == 1 && m_aoSourceDesc.back().bSameNoData &&
4418 10 : m_bSameNoData && m_bSameDataType &&
4419 5 : IsSameNaNAware(papoBands[0]->GetNoDataValue(&bHasNoData),
4420 5 : m_aoSourceDesc.back().dfSameNoData) &&
4421 381 : bHasNoData)) &&
4422 204 : (!m_aoSourceDesc.back().poMaskBand ||
4423 : // Also, if there's a single source that has a mask band, and the VRT
4424 : // bands have no-nodata or a 0-nodata value, we can skip
4425 : // initialization.
4426 51 : (m_aoSourceDesc.size() == 1 && m_bSameDataType &&
4427 11 : !(nBandCount == 1 && panBandMap[0] == 0) && m_bSameNoData &&
4428 11 : papoBands[0]->GetNoDataValue(&bHasNoData) == 0)))
4429 : {
4430 135 : bNeedInitBuffer = false;
4431 : }
4432 190 : return bNeedInitBuffer;
4433 : }
4434 :
4435 : /************************************************************************/
4436 : /* InitBuffer() */
4437 : /************************************************************************/
4438 :
4439 60 : void GDALTileIndexDataset::InitBuffer(void *pData, int nBufXSize, int nBufYSize,
4440 : GDALDataType eBufType, int nBandCount,
4441 : const int *panBandMap,
4442 : GSpacing nPixelSpace, GSpacing nLineSpace,
4443 : GSpacing nBandSpace) const
4444 : {
4445 60 : const int nBufTypeSize = GDALGetDataTypeSizeBytes(eBufType);
4446 60 : if (m_bSameNoData && nBandCount > 1 &&
4447 18 : ((nPixelSpace == nBufTypeSize &&
4448 18 : nLineSpace == nBufXSize * nPixelSpace &&
4449 18 : nBandSpace == nBufYSize * nLineSpace) ||
4450 0 : (nBandSpace == nBufTypeSize &&
4451 0 : nPixelSpace == nBandCount * nBandSpace &&
4452 0 : nLineSpace == nBufXSize * nPixelSpace)))
4453 : {
4454 18 : const int nBandNr = panBandMap[0];
4455 : auto poVRTBand =
4456 : nBandNr == 0
4457 18 : ? m_poMaskBand.get()
4458 18 : : cpl::down_cast<GDALTileIndexBand *>(papoBands[nBandNr - 1]);
4459 18 : CPLAssert(poVRTBand);
4460 18 : const double dfNoData = poVRTBand->m_dfNoDataValue;
4461 18 : if (dfNoData == 0.0)
4462 : {
4463 16 : memset(pData, 0,
4464 16 : static_cast<size_t>(nBufXSize) * nBufYSize * nBandCount *
4465 16 : nBufTypeSize);
4466 : }
4467 : else
4468 : {
4469 2 : GDALCopyWords64(
4470 : &dfNoData, GDT_Float64, 0, pData, eBufType, nBufTypeSize,
4471 2 : static_cast<size_t>(nBufXSize) * nBufYSize * nBandCount);
4472 18 : }
4473 : }
4474 : else
4475 : {
4476 85 : for (int i = 0; i < nBandCount; ++i)
4477 : {
4478 43 : const int nBandNr = panBandMap[i];
4479 43 : auto poVRTBand = nBandNr == 0 ? m_poMaskBand.get()
4480 41 : : cpl::down_cast<GDALTileIndexBand *>(
4481 41 : papoBands[nBandNr - 1]);
4482 43 : GByte *pabyBandData = static_cast<GByte *>(pData) + i * nBandSpace;
4483 43 : if (nPixelSpace == nBufTypeSize &&
4484 43 : poVRTBand->m_dfNoDataValue == 0.0)
4485 : {
4486 37 : if (nLineSpace == nBufXSize * nPixelSpace)
4487 : {
4488 37 : memset(pabyBandData, 0,
4489 37 : static_cast<size_t>(nBufYSize * nLineSpace));
4490 : }
4491 : else
4492 : {
4493 0 : for (int iLine = 0; iLine < nBufYSize; iLine++)
4494 : {
4495 0 : memset(static_cast<GByte *>(pabyBandData) +
4496 0 : static_cast<GIntBig>(iLine) * nLineSpace,
4497 0 : 0, static_cast<size_t>(nBufXSize * nPixelSpace));
4498 : }
4499 37 : }
4500 : }
4501 : else
4502 : {
4503 6 : double dfWriteValue = poVRTBand->m_dfNoDataValue;
4504 :
4505 1014 : for (int iLine = 0; iLine < nBufYSize; iLine++)
4506 : {
4507 1008 : GDALCopyWords(&dfWriteValue, GDT_Float64, 0,
4508 1008 : static_cast<GByte *>(pabyBandData) +
4509 1008 : static_cast<GIntBig>(nLineSpace) * iLine,
4510 : eBufType, static_cast<int>(nPixelSpace),
4511 : nBufXSize);
4512 : }
4513 : }
4514 : }
4515 : }
4516 60 : }
4517 :
4518 : /************************************************************************/
4519 : /* RenderSource() */
4520 : /************************************************************************/
4521 :
4522 492 : CPLErr GDALTileIndexDataset::RenderSource(
4523 : const SourceDesc &oSourceDesc, bool bNeedInitBuffer, int nBandNrMax,
4524 : int nXOff, int nYOff, int nXSize, int nYSize, double dfXOff, double dfYOff,
4525 : double dfXSize, double dfYSize, int nBufXSize, int nBufYSize, void *pData,
4526 : GDALDataType eBufType, int nBandCount, BANDMAP_TYPE panBandMapIn,
4527 : GSpacing nPixelSpace, GSpacing nLineSpace, GSpacing nBandSpace,
4528 : GDALRasterIOExtraArg *psExtraArg,
4529 : VRTSource::WorkingState &oWorkingState) const
4530 : {
4531 492 : auto &poTileDS = oSourceDesc.poDS;
4532 492 : auto &poSource = oSourceDesc.poSource;
4533 492 : auto poComplexSource = dynamic_cast<VRTComplexSource *>(poSource.get());
4534 492 : CPLErr eErr = CE_None;
4535 :
4536 872 : const auto GetBandFromBandMap = [&oSourceDesc, panBandMapIn](int iBand)
4537 : {
4538 441 : return oSourceDesc.bBandMapTakenIntoAccount ? iBand + 1
4539 441 : : panBandMapIn[iBand];
4540 492 : };
4541 :
4542 984 : std::vector<int> anSerial;
4543 492 : if (oSourceDesc.bBandMapTakenIntoAccount)
4544 : {
4545 25 : for (int i = 0; i < nBandCount; ++i)
4546 15 : anSerial.push_back(i + 1);
4547 : }
4548 : BANDMAP_TYPE panBandMapForRasterIO =
4549 492 : oSourceDesc.bBandMapTakenIntoAccount ? anSerial.data() : panBandMapIn;
4550 :
4551 492 : if (poTileDS->GetRasterCount() + 1 == nBandNrMax &&
4552 496 : papoBands[nBandNrMax - 1]->GetColorInterpretation() == GCI_AlphaBand &&
4553 4 : papoBands[nBandNrMax - 1]->GetRasterDataType() == GDT_UInt8)
4554 : {
4555 : GDALRasterIOExtraArg sExtraArg;
4556 4 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
4557 4 : if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
4558 : {
4559 : // cppcheck-suppress redundantAssignment
4560 0 : sExtraArg.eResampleAlg = psExtraArg->eResampleAlg;
4561 : }
4562 : else
4563 : {
4564 : // cppcheck-suppress redundantAssignment
4565 4 : sExtraArg.eResampleAlg = m_eResampling;
4566 : }
4567 :
4568 : // Special case when there's typically a mix of RGB and RGBA source
4569 : // datasets and we read a RGB one.
4570 14 : for (int iBand = 0; iBand < nBandCount && eErr == CE_None; ++iBand)
4571 : {
4572 10 : const int nBandNr = GetBandFromBandMap(iBand);
4573 10 : if (nBandNr == nBandNrMax)
4574 : {
4575 : // The window we will actually request from the source raster band.
4576 4 : double dfReqXOff = 0.0;
4577 4 : double dfReqYOff = 0.0;
4578 4 : double dfReqXSize = 0.0;
4579 4 : double dfReqYSize = 0.0;
4580 4 : int nReqXOff = 0;
4581 4 : int nReqYOff = 0;
4582 4 : int nReqXSize = 0;
4583 4 : int nReqYSize = 0;
4584 :
4585 : // The window we will actual set _within_ the pData buffer.
4586 4 : int nOutXOff = 0;
4587 4 : int nOutYOff = 0;
4588 4 : int nOutXSize = 0;
4589 4 : int nOutYSize = 0;
4590 :
4591 4 : bool bError = false;
4592 :
4593 4 : auto poTileBand = poTileDS->GetRasterBand(1);
4594 4 : poSource->SetRasterBand(poTileBand, false);
4595 4 : if (poSource->GetSrcDstWindow(
4596 : dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
4597 : sExtraArg.eResampleAlg, &dfReqXOff, &dfReqYOff,
4598 : &dfReqXSize, &dfReqYSize, &nReqXOff, &nReqYOff,
4599 : &nReqXSize, &nReqYSize, &nOutXOff, &nOutYOff,
4600 4 : &nOutXSize, &nOutYSize, bError))
4601 : {
4602 4 : GByte *pabyOut =
4603 : static_cast<GByte *>(pData) +
4604 4 : static_cast<GPtrDiff_t>(iBand * nBandSpace +
4605 4 : nOutXOff * nPixelSpace +
4606 4 : nOutYOff * nLineSpace);
4607 :
4608 4 : constexpr GByte n255 = 255;
4609 8 : for (int iY = 0; iY < nOutYSize; iY++)
4610 : {
4611 4 : GDALCopyWords(
4612 : &n255, GDT_UInt8, 0,
4613 4 : pabyOut + static_cast<GPtrDiff_t>(iY * nLineSpace),
4614 : eBufType, static_cast<int>(nPixelSpace), nOutXSize);
4615 : }
4616 : }
4617 : }
4618 : else
4619 : {
4620 6 : auto poTileBand = poTileDS->GetRasterBand(nBandNr);
4621 6 : if (poComplexSource)
4622 : {
4623 0 : int bHasNoData = false;
4624 : const double dfNoDataValue =
4625 0 : poTileBand->GetNoDataValue(&bHasNoData);
4626 0 : poComplexSource->SetNoDataValue(
4627 0 : bHasNoData ? dfNoDataValue : VRT_NODATA_UNSET);
4628 : }
4629 6 : poSource->SetRasterBand(poTileBand, false);
4630 :
4631 6 : GByte *pabyBandData =
4632 6 : static_cast<GByte *>(pData) + iBand * nBandSpace;
4633 12 : eErr = poSource->RasterIO(
4634 : poTileBand->GetRasterDataType(), nXOff, nYOff, nXSize,
4635 : nYSize, pabyBandData, nBufXSize, nBufYSize, eBufType,
4636 6 : nPixelSpace, nLineSpace, &sExtraArg, oWorkingState);
4637 : }
4638 : }
4639 4 : return eErr;
4640 : }
4641 966 : else if (!oSourceDesc.bBandMapTakenIntoAccount &&
4642 478 : poTileDS->GetRasterCount() < nBandNrMax)
4643 : {
4644 2 : CPLError(CE_Failure, CPLE_AppDefined, "%s has not enough bands.",
4645 : oSourceDesc.osName.c_str());
4646 2 : return CE_Failure;
4647 : }
4648 :
4649 486 : if ((oSourceDesc.poMaskBand && bNeedInitBuffer) || nBandNrMax == 0)
4650 : {
4651 : // The window we will actually request from the source raster band.
4652 55 : double dfReqXOff = 0.0;
4653 55 : double dfReqYOff = 0.0;
4654 55 : double dfReqXSize = 0.0;
4655 55 : double dfReqYSize = 0.0;
4656 55 : int nReqXOff = 0;
4657 55 : int nReqYOff = 0;
4658 55 : int nReqXSize = 0;
4659 55 : int nReqYSize = 0;
4660 :
4661 : // The window we will actual set _within_ the pData buffer.
4662 55 : int nOutXOff = 0;
4663 55 : int nOutYOff = 0;
4664 55 : int nOutXSize = 0;
4665 55 : int nOutYSize = 0;
4666 :
4667 55 : bool bError = false;
4668 :
4669 55 : auto poFirstTileBand = poTileDS->GetRasterBand(1);
4670 55 : poSource->SetRasterBand(poFirstTileBand, false);
4671 :
4672 : GDALRasterIOExtraArg sExtraArg;
4673 55 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
4674 55 : CPL_IGNORE_RET_VAL(sExtraArg.bFloatingPointWindowValidity);
4675 55 : CPL_IGNORE_RET_VAL(sExtraArg.eResampleAlg);
4676 55 : if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
4677 : {
4678 0 : sExtraArg.eResampleAlg = psExtraArg->eResampleAlg;
4679 : }
4680 : else
4681 : {
4682 55 : sExtraArg.eResampleAlg = m_eResampling;
4683 : }
4684 :
4685 55 : if (poSource->GetSrcDstWindow(
4686 : dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
4687 : sExtraArg.eResampleAlg, &dfReqXOff, &dfReqYOff, &dfReqXSize,
4688 : &dfReqYSize, &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
4689 55 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError))
4690 : {
4691 55 : int iMaskBandIdx = -1;
4692 55 : if (eBufType == GDT_UInt8 && nBandNrMax == 0)
4693 : {
4694 : // when called from m_poMaskBand
4695 4 : iMaskBandIdx = 0;
4696 : }
4697 51 : else if (oSourceDesc.poMaskBand)
4698 : {
4699 : // If we request a Byte buffer and the mask band is actually
4700 : // one of the queried bands of this request, we can save
4701 : // requesting it separately.
4702 51 : const int nMaskBandNr = oSourceDesc.poMaskBand->GetBand();
4703 39 : if (eBufType == GDT_UInt8 && nMaskBandNr >= 1 &&
4704 129 : nMaskBandNr <= poTileDS->GetRasterCount() &&
4705 39 : poTileDS->GetRasterBand(nMaskBandNr) ==
4706 39 : oSourceDesc.poMaskBand)
4707 : {
4708 61 : for (int iBand = 0; iBand < nBandCount; ++iBand)
4709 : {
4710 44 : if (panBandMapIn[iBand] == nMaskBandNr)
4711 : {
4712 20 : iMaskBandIdx = iBand;
4713 20 : break;
4714 : }
4715 : }
4716 : }
4717 : }
4718 :
4719 55 : sExtraArg.bFloatingPointWindowValidity = TRUE;
4720 55 : sExtraArg.dfXOff = dfReqXOff;
4721 55 : sExtraArg.dfYOff = dfReqYOff;
4722 55 : sExtraArg.dfXSize = dfReqXSize;
4723 55 : sExtraArg.dfYSize = dfReqYSize;
4724 :
4725 76 : if (iMaskBandIdx < 0 && oSourceDesc.abyMask.empty() &&
4726 21 : oSourceDesc.poMaskBand)
4727 : {
4728 : // Fetch the mask band
4729 : try
4730 : {
4731 21 : oSourceDesc.abyMask.resize(static_cast<size_t>(nOutXSize) *
4732 21 : nOutYSize);
4733 : }
4734 0 : catch (const std::bad_alloc &)
4735 : {
4736 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
4737 : "Cannot allocate working buffer for mask");
4738 0 : return CE_Failure;
4739 : }
4740 :
4741 21 : if (oSourceDesc.poMaskBand->RasterIO(
4742 : GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
4743 21 : oSourceDesc.abyMask.data(), nOutXSize, nOutYSize,
4744 21 : GDT_UInt8, 0, 0, &sExtraArg) != CE_None)
4745 : {
4746 0 : oSourceDesc.abyMask.clear();
4747 0 : return CE_Failure;
4748 : }
4749 : }
4750 :
4751 : // Allocate a temporary contiguous buffer to receive pixel data
4752 55 : const int nBufTypeSize = GDALGetDataTypeSizeBytes(eBufType);
4753 55 : const size_t nWorkBufferBandSize =
4754 55 : static_cast<size_t>(nOutXSize) * nOutYSize * nBufTypeSize;
4755 55 : std::vector<GByte> abyWorkBuffer;
4756 : try
4757 : {
4758 55 : abyWorkBuffer.resize(nBandCount * nWorkBufferBandSize);
4759 : }
4760 0 : catch (const std::bad_alloc &)
4761 : {
4762 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
4763 : "Cannot allocate working buffer");
4764 0 : return CE_Failure;
4765 : }
4766 :
4767 : const GByte *const pabyMask =
4768 : iMaskBandIdx >= 0
4769 24 : ? abyWorkBuffer.data() + iMaskBandIdx * nWorkBufferBandSize
4770 79 : : oSourceDesc.abyMask.data();
4771 :
4772 55 : if (nBandNrMax == 0)
4773 : {
4774 : // Special case when called from m_poMaskBand
4775 12 : if (poTileDS->GetRasterBand(1)->GetMaskBand()->RasterIO(
4776 : GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
4777 6 : abyWorkBuffer.data(), nOutXSize, nOutYSize, eBufType, 0,
4778 6 : 0, &sExtraArg) != CE_None)
4779 : {
4780 0 : return CE_Failure;
4781 : }
4782 : }
4783 98 : else if (poTileDS->RasterIO(GF_Read, nReqXOff, nReqYOff, nReqXSize,
4784 49 : nReqYSize, abyWorkBuffer.data(),
4785 : nOutXSize, nOutYSize, eBufType,
4786 : nBandCount, panBandMapForRasterIO, 0, 0,
4787 49 : 0, &sExtraArg) != CE_None)
4788 : {
4789 0 : return CE_Failure;
4790 : }
4791 :
4792 : // Compose the temporary contiguous buffer into the target
4793 : // buffer, taking into account the mask
4794 55 : GByte *pabyOut = static_cast<GByte *>(pData) +
4795 55 : static_cast<GPtrDiff_t>(nOutXOff * nPixelSpace +
4796 55 : nOutYOff * nLineSpace);
4797 :
4798 121 : for (int iBand = 0; iBand < nBandCount && eErr == CE_None; ++iBand)
4799 : {
4800 66 : GByte *pabyDestBand =
4801 66 : pabyOut + static_cast<GPtrDiff_t>(iBand * nBandSpace);
4802 : const GByte *pabySrc =
4803 66 : abyWorkBuffer.data() + iBand * nWorkBufferBandSize;
4804 :
4805 66 : CompositeSrcWithMaskIntoDest(
4806 : nOutXSize, nOutYSize, eBufType, nBufTypeSize, nPixelSpace,
4807 : nLineSpace, pabySrc, pabyMask, pabyDestBand);
4808 : }
4809 55 : }
4810 : }
4811 431 : else if (m_bSameDataType && !bNeedInitBuffer && oSourceDesc.bHasNoData)
4812 : {
4813 : // We create a non-VRTComplexSource SimpleSource copy of poSource
4814 : // to be able to call DatasetRasterIO()
4815 4 : VRTSimpleSource oSimpleSource(poSource.get(), 1.0, 1.0);
4816 :
4817 : GDALRasterIOExtraArg sExtraArg;
4818 4 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
4819 4 : if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
4820 : {
4821 : // cppcheck-suppress redundantAssignment
4822 0 : sExtraArg.eResampleAlg = psExtraArg->eResampleAlg;
4823 : }
4824 : else
4825 : {
4826 : // cppcheck-suppress redundantAssignment
4827 4 : sExtraArg.eResampleAlg = m_eResampling;
4828 : }
4829 :
4830 4 : auto poTileBand = poTileDS->GetRasterBand(GetBandFromBandMap(0));
4831 4 : oSimpleSource.SetRasterBand(poTileBand, false);
4832 4 : eErr = oSimpleSource.DatasetRasterIO(
4833 4 : papoBands[0]->GetRasterDataType(), nXOff, nYOff, nXSize, nYSize,
4834 : pData, nBufXSize, nBufYSize, eBufType, nBandCount,
4835 : panBandMapForRasterIO, nPixelSpace, nLineSpace, nBandSpace,
4836 4 : &sExtraArg);
4837 : }
4838 427 : else if (m_bSameDataType && !poComplexSource)
4839 : {
4840 418 : auto poTileBand = poTileDS->GetRasterBand(GetBandFromBandMap(0));
4841 418 : poSource->SetRasterBand(poTileBand, false);
4842 :
4843 : GDALRasterIOExtraArg sExtraArg;
4844 418 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
4845 418 : if (poTileBand->GetColorTable())
4846 : {
4847 : // cppcheck-suppress redundantAssignment
4848 0 : sExtraArg.eResampleAlg = GRIORA_NearestNeighbour;
4849 : }
4850 418 : else if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
4851 : {
4852 : // cppcheck-suppress redundantAssignment
4853 0 : sExtraArg.eResampleAlg = psExtraArg->eResampleAlg;
4854 : }
4855 : else
4856 : {
4857 : // cppcheck-suppress redundantAssignment
4858 418 : sExtraArg.eResampleAlg = m_eResampling;
4859 : }
4860 :
4861 836 : eErr = poSource->DatasetRasterIO(
4862 418 : papoBands[0]->GetRasterDataType(), nXOff, nYOff, nXSize, nYSize,
4863 : pData, nBufXSize, nBufYSize, eBufType, nBandCount,
4864 : panBandMapForRasterIO, nPixelSpace, nLineSpace, nBandSpace,
4865 418 : &sExtraArg);
4866 : }
4867 : else
4868 : {
4869 18 : for (int i = 0; i < nBandCount && eErr == CE_None; ++i)
4870 : {
4871 9 : const int nBandNr = GetBandFromBandMap(i);
4872 9 : GByte *pabyBandData = static_cast<GByte *>(pData) + i * nBandSpace;
4873 9 : auto poTileBand = poTileDS->GetRasterBand(nBandNr);
4874 9 : if (poComplexSource)
4875 : {
4876 9 : int bHasNoData = false;
4877 : const double dfNoDataValue =
4878 9 : poTileBand->GetNoDataValue(&bHasNoData);
4879 9 : poComplexSource->SetNoDataValue(bHasNoData ? dfNoDataValue
4880 9 : : VRT_NODATA_UNSET);
4881 : }
4882 9 : poSource->SetRasterBand(poTileBand, false);
4883 :
4884 : GDALRasterIOExtraArg sExtraArg;
4885 9 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
4886 9 : if (poTileBand->GetColorTable())
4887 : {
4888 : // cppcheck-suppress redundantAssignment
4889 0 : sExtraArg.eResampleAlg = GRIORA_NearestNeighbour;
4890 : }
4891 9 : else if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
4892 : {
4893 : // cppcheck-suppress redundantAssignment
4894 0 : sExtraArg.eResampleAlg = psExtraArg->eResampleAlg;
4895 : }
4896 : else
4897 : {
4898 : // cppcheck-suppress redundantAssignment
4899 9 : sExtraArg.eResampleAlg = m_eResampling;
4900 : }
4901 :
4902 18 : eErr = poSource->RasterIO(
4903 9 : papoBands[nBandNr - 1]->GetRasterDataType(), nXOff, nYOff,
4904 : nXSize, nYSize, pabyBandData, nBufXSize, nBufYSize, eBufType,
4905 9 : nPixelSpace, nLineSpace, &sExtraArg, oWorkingState);
4906 : }
4907 : }
4908 486 : return eErr;
4909 : }
4910 :
4911 : /************************************************************************/
4912 : /* IRasterIO() */
4913 : /************************************************************************/
4914 :
4915 196 : CPLErr GDALTileIndexDataset::IRasterIO(
4916 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
4917 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
4918 : int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
4919 : GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg)
4920 : {
4921 196 : if (eRWFlag != GF_Read)
4922 0 : return CE_Failure;
4923 :
4924 196 : if (nBufXSize < nXSize && nBufYSize < nYSize && AreOverviewsEnabled())
4925 : {
4926 2 : int bTried = FALSE;
4927 2 : const CPLErr eErr = TryOverviewRasterIO(
4928 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
4929 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
4930 : nBandSpace, psExtraArg, &bTried);
4931 2 : if (bTried)
4932 2 : return eErr;
4933 : }
4934 :
4935 194 : double dfXOff = nXOff;
4936 194 : double dfYOff = nYOff;
4937 194 : double dfXSize = nXSize;
4938 194 : double dfYSize = nYSize;
4939 194 : if (psExtraArg->bFloatingPointWindowValidity)
4940 : {
4941 6 : dfXOff = psExtraArg->dfXOff;
4942 6 : dfYOff = psExtraArg->dfYOff;
4943 6 : dfXSize = psExtraArg->dfXSize;
4944 6 : dfYSize = psExtraArg->dfYSize;
4945 : }
4946 :
4947 194 : if (!CollectSources(dfXOff, dfYOff, dfXSize, dfYSize, nBandCount,
4948 : panBandMap,
4949 : /* bMultiThreadAllowed = */ true))
4950 : {
4951 3 : return CE_Failure;
4952 : }
4953 :
4954 : // We might be called with nBandCount == 1 && panBandMap[0] == 0
4955 : // to mean m_poMaskBand
4956 191 : int nBandNrMax = 0;
4957 436 : for (int i = 0; i < nBandCount; ++i)
4958 : {
4959 245 : const int nBandNr = panBandMap[i];
4960 245 : nBandNrMax = std::max(nBandNrMax, nBandNr);
4961 : }
4962 :
4963 : const bool bNeedInitBuffer =
4964 191 : m_bLastMustUseMultiThreading || NeedInitBuffer(nBandCount, panBandMap);
4965 :
4966 191 : if (!bNeedInitBuffer)
4967 : {
4968 131 : return RenderSource(
4969 131 : m_aoSourceDesc.back(), bNeedInitBuffer, nBandNrMax, nXOff, nYOff,
4970 : nXSize, nYSize, dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize,
4971 : nBufYSize, pData, eBufType, nBandCount, panBandMap, nPixelSpace,
4972 262 : nLineSpace, nBandSpace, psExtraArg, m_oWorkingState);
4973 : }
4974 : else
4975 : {
4976 60 : InitBuffer(pData, nBufXSize, nBufYSize, eBufType, nBandCount,
4977 : panBandMap, nPixelSpace, nLineSpace, nBandSpace);
4978 :
4979 60 : if (m_bLastMustUseMultiThreading)
4980 : {
4981 12 : CPLErrorAccumulator oErrorAccumulator;
4982 6 : std::atomic<bool> bSuccess = true;
4983 : const int nContributingSources =
4984 6 : static_cast<int>(m_aoSourceDesc.size());
4985 6 : CPLWorkerThreadPool *psThreadPool = GDALGetGlobalThreadPool(
4986 6 : std::min(nContributingSources, m_nNumThreads));
4987 : const int nThreads =
4988 6 : std::min(nContributingSources, psThreadPool->GetThreadCount());
4989 6 : CPLDebugOnly("GTI",
4990 : "IRasterIO(): use optimized "
4991 : "multi-threaded code path. "
4992 : "Using %d threads",
4993 : nThreads);
4994 :
4995 : {
4996 12 : std::lock_guard oLock(m_oQueueWorkingStates.oMutex);
4997 6 : if (m_oQueueWorkingStates.oStates.size() <
4998 6 : static_cast<size_t>(nThreads))
4999 : {
5000 4 : m_oQueueWorkingStates.oStates.resize(nThreads);
5001 : }
5002 22 : for (int i = 0; i < nThreads; ++i)
5003 : {
5004 16 : if (!m_oQueueWorkingStates.oStates[i])
5005 10 : m_oQueueWorkingStates.oStates[i] =
5006 20 : std::make_unique<VRTSource::WorkingState>();
5007 : }
5008 : }
5009 :
5010 6 : auto oQueue = psThreadPool->CreateJobQueue();
5011 6 : std::atomic<int> nCompletedJobs = 0;
5012 144 : for (auto &oSourceDesc : m_aoSourceDesc)
5013 : {
5014 138 : auto psJob = new RasterIOJob();
5015 138 : psJob->bSTACCollection = m_bSTACCollection;
5016 138 : psJob->poDS = this;
5017 138 : psJob->pbSuccess = &bSuccess;
5018 138 : psJob->poErrorAccumulator = &oErrorAccumulator;
5019 138 : psJob->pnCompletedJobs = &nCompletedJobs;
5020 138 : psJob->poQueueWorkingStates = &m_oQueueWorkingStates;
5021 138 : psJob->nBandNrMax = nBandNrMax;
5022 138 : psJob->nXOff = nXOff;
5023 138 : psJob->nYOff = nYOff;
5024 138 : psJob->nXSize = nXSize;
5025 138 : psJob->nYSize = nYSize;
5026 138 : psJob->pData = pData;
5027 138 : psJob->nBufXSize = nBufXSize;
5028 138 : psJob->nBufYSize = nBufYSize;
5029 138 : psJob->eBufType = eBufType;
5030 138 : psJob->nBandCount = nBandCount;
5031 138 : psJob->panBandMap = panBandMap;
5032 138 : psJob->nPixelSpace = nPixelSpace;
5033 138 : psJob->nLineSpace = nLineSpace;
5034 138 : psJob->nBandSpace = nBandSpace;
5035 138 : psJob->psExtraArg = psExtraArg;
5036 :
5037 : psJob->osTileName = oSourceDesc.poFeature->GetFieldAsString(
5038 138 : m_nLocationFieldIndex);
5039 :
5040 138 : if (!oQueue->SubmitJob(RasterIOJob::Func, psJob))
5041 : {
5042 0 : delete psJob;
5043 0 : bSuccess = false;
5044 0 : break;
5045 : }
5046 : }
5047 :
5048 58 : while (oQueue->WaitEvent())
5049 : {
5050 : // Quite rough progress callback. We could do better by counting
5051 : // the number of contributing pixels.
5052 52 : if (psExtraArg->pfnProgress)
5053 : {
5054 100 : psExtraArg->pfnProgress(double(nCompletedJobs.load()) /
5055 : nContributingSources,
5056 : "", psExtraArg->pProgressData);
5057 : }
5058 : }
5059 :
5060 6 : oErrorAccumulator.ReplayErrors();
5061 :
5062 6 : if (bSuccess && psExtraArg->pfnProgress)
5063 : {
5064 4 : psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
5065 : }
5066 :
5067 6 : return bSuccess ? CE_None : CE_Failure;
5068 : }
5069 : else
5070 : {
5071 : // Now render from bottom of the stack to top.
5072 278 : for (auto &oSourceDesc : m_aoSourceDesc)
5073 : {
5074 448 : if (oSourceDesc.poDS &&
5075 224 : RenderSource(oSourceDesc, bNeedInitBuffer, nBandNrMax,
5076 : nXOff, nYOff, nXSize, nYSize, dfXOff, dfYOff,
5077 : dfXSize, dfYSize, nBufXSize, nBufYSize, pData,
5078 : eBufType, nBandCount, panBandMap, nPixelSpace,
5079 : nLineSpace, nBandSpace, psExtraArg,
5080 448 : m_oWorkingState) != CE_None)
5081 0 : return CE_Failure;
5082 : }
5083 :
5084 54 : if (psExtraArg->pfnProgress)
5085 : {
5086 4 : psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
5087 : }
5088 :
5089 54 : return CE_None;
5090 : }
5091 : }
5092 : }
5093 :
5094 : /************************************************************************/
5095 : /* GDALTileIndexDataset::RasterIOJob::Func() */
5096 : /************************************************************************/
5097 :
5098 138 : void GDALTileIndexDataset::RasterIOJob::Func(void *pData)
5099 : {
5100 : auto psJob =
5101 276 : std::unique_ptr<RasterIOJob>(static_cast<RasterIOJob *>(pData));
5102 138 : if (*psJob->pbSuccess)
5103 : {
5104 : const std::string osTileName(GetAbsoluteFileName(
5105 414 : psJob->osTileName.c_str(), psJob->poDS->GetDescription(),
5106 276 : psJob->bSTACCollection));
5107 :
5108 276 : SourceDesc oSourceDesc;
5109 :
5110 276 : auto oAccumulator = psJob->poErrorAccumulator->InstallForCurrentScope();
5111 138 : CPL_IGNORE_RET_VAL(oAccumulator);
5112 :
5113 : const bool bCanOpenSource =
5114 138 : psJob->poDS->GetSourceDesc(osTileName, oSourceDesc,
5115 138 : &psJob->poQueueWorkingStates->oMutex,
5116 413 : psJob->nBandCount, psJob->panBandMap) &&
5117 137 : oSourceDesc.poDS;
5118 :
5119 138 : if (!bCanOpenSource)
5120 : {
5121 1 : *psJob->pbSuccess = false;
5122 : }
5123 : else
5124 : {
5125 137 : GDALRasterIOExtraArg sArg = *(psJob->psExtraArg);
5126 137 : sArg.pfnProgress = nullptr;
5127 137 : sArg.pProgressData = nullptr;
5128 :
5129 137 : std::unique_ptr<VRTSource::WorkingState> poWorkingState;
5130 : {
5131 274 : std::lock_guard oLock(psJob->poQueueWorkingStates->oMutex);
5132 : poWorkingState =
5133 137 : std::move(psJob->poQueueWorkingStates->oStates.back());
5134 137 : psJob->poQueueWorkingStates->oStates.pop_back();
5135 137 : CPLAssert(poWorkingState.get());
5136 : }
5137 :
5138 137 : double dfXOff = psJob->nXOff;
5139 137 : double dfYOff = psJob->nYOff;
5140 137 : double dfXSize = psJob->nXSize;
5141 137 : double dfYSize = psJob->nYSize;
5142 137 : if (psJob->psExtraArg->bFloatingPointWindowValidity)
5143 : {
5144 0 : dfXOff = psJob->psExtraArg->dfXOff;
5145 0 : dfYOff = psJob->psExtraArg->dfYOff;
5146 0 : dfXSize = psJob->psExtraArg->dfXSize;
5147 0 : dfYSize = psJob->psExtraArg->dfYSize;
5148 : }
5149 :
5150 : const bool bRenderOK =
5151 274 : psJob->poDS->RenderSource(
5152 137 : oSourceDesc, /*bNeedInitBuffer = */ true, psJob->nBandNrMax,
5153 137 : psJob->nXOff, psJob->nYOff, psJob->nXSize, psJob->nYSize,
5154 137 : dfXOff, dfYOff, dfXSize, dfYSize, psJob->nBufXSize,
5155 137 : psJob->nBufYSize, psJob->pData, psJob->eBufType,
5156 137 : psJob->nBandCount, psJob->panBandMap, psJob->nPixelSpace,
5157 137 : psJob->nLineSpace, psJob->nBandSpace, &sArg,
5158 137 : *(poWorkingState.get())) == CE_None;
5159 :
5160 137 : if (!bRenderOK)
5161 : {
5162 1 : *psJob->pbSuccess = false;
5163 : }
5164 :
5165 : {
5166 274 : std::lock_guard oLock(psJob->poQueueWorkingStates->oMutex);
5167 274 : psJob->poQueueWorkingStates->oStates.push_back(
5168 137 : std::move(poWorkingState));
5169 : }
5170 : }
5171 : }
5172 :
5173 138 : ++(*psJob->pnCompletedJobs);
5174 138 : }
5175 :
5176 : #ifdef GDAL_ENABLE_ALGORITHMS
5177 :
5178 : /************************************************************************/
5179 : /* GDALGTICreateAlgorithm */
5180 : /************************************************************************/
5181 :
5182 : class GDALGTICreateAlgorithm final : public GDALRasterIndexAlgorithm
5183 : {
5184 : public:
5185 : static constexpr const char *NAME = "create";
5186 : static constexpr const char *DESCRIPTION =
5187 : "Create an index of raster datasets compatible of the GDAL Tile Index "
5188 : "(GTI) driver.";
5189 : static constexpr const char *HELP_URL =
5190 : "/programs/gdal_driver_gti_create.html";
5191 :
5192 : GDALGTICreateAlgorithm();
5193 :
5194 : protected:
5195 : bool AddExtraOptions(CPLStringList &aosOptions) override;
5196 :
5197 : private:
5198 : std::string m_xmlFilename{};
5199 : std::vector<double> m_resolution{};
5200 : std::vector<double> m_bbox{};
5201 : std::string m_dataType{};
5202 : int m_bandCount = 0;
5203 : std::vector<double> m_nodata{};
5204 : std::vector<std::string> m_colorInterpretation{};
5205 : bool m_mask = false;
5206 : std::vector<std::string> m_fetchedMetadata{};
5207 : };
5208 :
5209 : /************************************************************************/
5210 : /* GDALGTICreateAlgorithm::GDALGTICreateAlgorithm() */
5211 : /************************************************************************/
5212 :
5213 146 : GDALGTICreateAlgorithm::GDALGTICreateAlgorithm()
5214 146 : : GDALRasterIndexAlgorithm(NAME, DESCRIPTION, HELP_URL)
5215 : {
5216 146 : AddProgressArg();
5217 146 : AddInputDatasetArg(&m_inputDatasets, GDAL_OF_RASTER)
5218 146 : .SetAutoOpenDataset(false);
5219 146 : GDALVectorOutputAbstractAlgorithm::AddAllOutputArgs();
5220 :
5221 146 : AddCommonOptions();
5222 :
5223 : AddArg("xml-filename", 0,
5224 : _("Filename of the XML Virtual Tile Index file to generate, that "
5225 : "can be used as an input for the GDAL GTI / Virtual Raster Tile "
5226 : "Index driver"),
5227 292 : &m_xmlFilename)
5228 146 : .SetMinCharCount(1);
5229 :
5230 : AddArg("resolution", 0,
5231 : _("Resolution (in destination CRS units) of the virtual mosaic"),
5232 292 : &m_resolution)
5233 146 : .SetMinCount(2)
5234 146 : .SetMaxCount(2)
5235 146 : .SetMinValueExcluded(0)
5236 146 : .SetRepeatedArgAllowed(false)
5237 146 : .SetDisplayHintAboutRepetition(false)
5238 146 : .SetMetaVar("<xres>,<yres>");
5239 :
5240 : AddBBOXArg(
5241 : &m_bbox,
5242 146 : _("Bounding box (in destination CRS units) of the virtual mosaic"));
5243 146 : AddOutputDataTypeArg(&m_dataType, _("Datatype of the virtual mosaic"));
5244 : AddArg("band-count", 0, _("Number of bands of the virtual mosaic"),
5245 292 : &m_bandCount)
5246 146 : .SetMinValueIncluded(1);
5247 : AddArg("nodata", 0, _("Nodata value(s) of the bands of the virtual mosaic"),
5248 146 : &m_nodata);
5249 : AddArg("color-interpretation", 0,
5250 : _("Color interpretation(s) of the bands of the virtual mosaic"),
5251 292 : &m_colorInterpretation)
5252 146 : .SetChoices("red", "green", "blue", "alpha", "gray", "undefined");
5253 : AddArg("mask", 0, _("Defines that the virtual mosaic has a mask band"),
5254 146 : &m_mask);
5255 : AddArg("fetch-metadata", 0,
5256 : _("Fetch a metadata item from source rasters and write it as a "
5257 : "field in the index."),
5258 292 : &m_fetchedMetadata)
5259 292 : .SetMetaVar("<gdal-metadata-name>,<field-name>,<field-type>")
5260 146 : .SetPackedValuesAllowed(false)
5261 : .AddValidationAction(
5262 6 : [this]()
5263 : {
5264 6 : for (const std::string &s : m_fetchedMetadata)
5265 : {
5266 : const CPLStringList aosTokens(
5267 4 : CSLTokenizeString2(s.c_str(), ",", 0));
5268 4 : if (aosTokens.size() != 3)
5269 : {
5270 1 : ReportError(
5271 : CE_Failure, CPLE_IllegalArg,
5272 : "'%s' is not of the form "
5273 : "<gdal-metadata-name>,<field-name>,<field-type>",
5274 : s.c_str());
5275 1 : return false;
5276 : }
5277 3 : bool ok = false;
5278 18 : for (const char *type : {"String", "Integer", "Integer64",
5279 21 : "Real", "Date", "DateTime"})
5280 : {
5281 18 : if (EQUAL(aosTokens[2], type))
5282 2 : ok = true;
5283 : }
5284 3 : if (!ok)
5285 : {
5286 1 : ReportError(CE_Failure, CPLE_IllegalArg,
5287 : "'%s' has an invalid field type '%s'. It "
5288 : "should be one of 'String', 'Integer', "
5289 : "'Integer64', 'Real', 'Date', 'DateTime'.",
5290 : s.c_str(), aosTokens[2]);
5291 1 : return false;
5292 : }
5293 : }
5294 2 : return true;
5295 146 : });
5296 146 : }
5297 :
5298 : /************************************************************************/
5299 : /* GDALGTICreateAlgorithm::AddExtraOptions() */
5300 : /************************************************************************/
5301 :
5302 6 : bool GDALGTICreateAlgorithm::AddExtraOptions(CPLStringList &aosOptions)
5303 : {
5304 6 : if (!m_xmlFilename.empty())
5305 : {
5306 1 : aosOptions.push_back("-gti_filename");
5307 1 : aosOptions.push_back(m_xmlFilename);
5308 : }
5309 6 : if (!m_resolution.empty())
5310 : {
5311 1 : aosOptions.push_back("-tr");
5312 1 : aosOptions.push_back(CPLSPrintf("%.17g", m_resolution[0]));
5313 1 : aosOptions.push_back(CPLSPrintf("%.17g", m_resolution[1]));
5314 : }
5315 6 : if (!m_bbox.empty())
5316 : {
5317 1 : aosOptions.push_back("-te");
5318 1 : aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[0]));
5319 1 : aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[1]));
5320 1 : aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[2]));
5321 1 : aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[3]));
5322 : }
5323 6 : if (!m_dataType.empty())
5324 : {
5325 1 : aosOptions.push_back("-ot");
5326 1 : aosOptions.push_back(m_dataType);
5327 : }
5328 6 : if (m_bandCount > 0)
5329 : {
5330 3 : aosOptions.push_back("-bandcount");
5331 3 : aosOptions.push_back(CPLSPrintf("%d", m_bandCount));
5332 :
5333 5 : if (!m_nodata.empty() && m_nodata.size() != 1 &&
5334 2 : static_cast<int>(m_nodata.size()) != m_bandCount)
5335 : {
5336 1 : ReportError(CE_Failure, CPLE_IllegalArg,
5337 : "%d nodata values whereas one or %d were expected",
5338 1 : static_cast<int>(m_nodata.size()), m_bandCount);
5339 1 : return false;
5340 : }
5341 :
5342 4 : if (!m_colorInterpretation.empty() &&
5343 4 : m_colorInterpretation.size() != 1 &&
5344 2 : static_cast<int>(m_colorInterpretation.size()) != m_bandCount)
5345 : {
5346 1 : ReportError(
5347 : CE_Failure, CPLE_IllegalArg,
5348 : "%d color interpretations whereas one or %d were expected",
5349 1 : static_cast<int>(m_colorInterpretation.size()), m_bandCount);
5350 1 : return false;
5351 : }
5352 : }
5353 4 : if (!m_nodata.empty())
5354 : {
5355 2 : std::string val;
5356 3 : for (double v : m_nodata)
5357 : {
5358 2 : if (!val.empty())
5359 1 : val += ',';
5360 2 : val += CPLSPrintf("%.17g", v);
5361 : }
5362 1 : aosOptions.push_back("-nodata");
5363 1 : aosOptions.push_back(val);
5364 : }
5365 4 : if (!m_colorInterpretation.empty())
5366 : {
5367 2 : std::string val;
5368 3 : for (const std::string &s : m_colorInterpretation)
5369 : {
5370 2 : if (!val.empty())
5371 1 : val += ',';
5372 2 : val += s;
5373 : }
5374 1 : aosOptions.push_back("-colorinterp");
5375 1 : aosOptions.push_back(val);
5376 : }
5377 4 : if (m_mask)
5378 1 : aosOptions.push_back("-mask");
5379 5 : for (const std::string &s : m_fetchedMetadata)
5380 : {
5381 1 : aosOptions.push_back("-fetch_md");
5382 2 : const CPLStringList aosTokens(CSLTokenizeString2(s.c_str(), ",", 0));
5383 4 : for (const char *token : aosTokens)
5384 : {
5385 3 : aosOptions.push_back(token);
5386 : }
5387 : }
5388 4 : return true;
5389 : }
5390 :
5391 : /************************************************************************/
5392 : /* GDALTileIndexInstantiateAlgorithm() */
5393 : /************************************************************************/
5394 :
5395 : static GDALAlgorithm *
5396 146 : GDALTileIndexInstantiateAlgorithm(const std::vector<std::string> &aosPath)
5397 : {
5398 146 : if (aosPath.size() == 1 && aosPath[0] == "create")
5399 : {
5400 146 : return std::make_unique<GDALGTICreateAlgorithm>().release();
5401 : }
5402 : else
5403 : {
5404 0 : return nullptr;
5405 : }
5406 : }
5407 :
5408 : #endif
5409 :
5410 : /************************************************************************/
5411 : /* GDALRegister_GTI() */
5412 : /************************************************************************/
5413 :
5414 2058 : void GDALRegister_GTI()
5415 : {
5416 2058 : if (GDALGetDriverByName("GTI") != nullptr)
5417 263 : return;
5418 :
5419 3590 : auto poDriver = std::make_unique<GDALDriver>();
5420 :
5421 1795 : poDriver->SetDescription("GTI");
5422 1795 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
5423 1795 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "GDAL Raster Tile Index");
5424 1795 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "gti.gpkg gti.fgb gti");
5425 1795 : poDriver->SetMetadataItem(GDAL_DMD_CONNECTION_PREFIX, GTI_PREFIX);
5426 1795 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/gti.html");
5427 :
5428 1795 : poDriver->pfnOpen = GDALTileIndexDatasetOpen;
5429 1795 : poDriver->pfnIdentify = GDALTileIndexDatasetIdentify;
5430 :
5431 1795 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
5432 :
5433 1795 : poDriver->SetMetadataItem(
5434 : GDAL_DMD_OPENOPTIONLIST,
5435 : "<OpenOptionList>"
5436 : " <Option name='LAYER' type='string'/>"
5437 : " <Option name='SQL' type='string'/>"
5438 : " <Option name='SPATIAL_SQL' type='string'/>"
5439 : " <Option name='LOCATION_FIELD' type='string'/>"
5440 : " <Option name='SORT_FIELD' type='string'/>"
5441 : " <Option name='SORT_FIELD_ASC' type='boolean'/>"
5442 : " <Option name='FILTER' type='string'/>"
5443 : " <Option name='SRS' type='string'/>"
5444 : " <Option name='RESX' type='float'/>"
5445 : " <Option name='RESY' type='float'/>"
5446 : " <Option name='MINX' type='float'/>"
5447 : " <Option name='MINY' type='float'/>"
5448 : " <Option name='MAXX' type='float'/>"
5449 : " <Option name='MAXY' type='float'/>"
5450 : " <Option name='NUM_THREADS' type='string' description="
5451 : "'Number of worker threads for reading. Can be set to ALL_CPUS' "
5452 : "default='ALL_CPUS'/>"
5453 : " <Option name='WARPING_MEMORY_SIZE' type='string' description="
5454 : "'Set the amount of memory that the warp API is allowed to use for "
5455 : "caching' default='64MB'/>"
5456 1795 : "</OpenOptionList>");
5457 :
5458 : #ifdef GDAL_ENABLE_ALGORITHMS
5459 3590 : poDriver->DeclareAlgorithm({"create"});
5460 1795 : poDriver->pfnInstantiateAlgorithm = GDALTileIndexInstantiateAlgorithm;
5461 : #endif
5462 :
5463 : #ifdef BUILT_AS_PLUGIN
5464 : // Used by gdaladdo and test_gdaladdo.py
5465 : poDriver->SetMetadataItem("IS_PLUGIN", "YES");
5466 : #endif
5467 :
5468 1795 : GetGDALDriverManager()->RegisterDriver(poDriver.release());
5469 : }
|
