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