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