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