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