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