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