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