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