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