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