Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: STACTA (Spatio-Temporal Asset Catalog Tiled Assets) driver
5 : * Author: Even Rouault, <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2020, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * Permission is hereby granted, free of charge, to any person obtaining a
11 : * copy of this software and associated documentation files (the "Software"),
12 : * to deal in the Software without restriction, including without limitation
13 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 : * and/or sell copies of the Software, and to permit persons to whom the
15 : * Software is furnished to do so, subject to the following conditions:
16 : *
17 : * The above copyright notice and this permission notice shall be included
18 : * in all copies or substantial portions of the Software.
19 : *
20 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 : * DEALINGS IN THE SOFTWARE.
27 : ****************************************************************************/
28 :
29 : #include "cpl_json.h"
30 : #include "cpl_mem_cache.h"
31 : #include "cpl_string.h"
32 : #include "gdal_pam.h"
33 : #include "gdal_utils.h"
34 : #include "memdataset.h"
35 : #include "tilematrixset.hpp"
36 : #include "stactadataset.h"
37 :
38 : #include <algorithm>
39 : #include <array>
40 : #include <limits>
41 : #include <map>
42 : #include <memory>
43 : #include <vector>
44 :
45 : extern "C" void GDALRegister_STACTA();
46 :
47 : // Implements a driver for
48 : // https://github.com/stac-extensions/tiled-assets
49 :
50 : /************************************************************************/
51 : /* STACTARasterBand() */
52 : /************************************************************************/
53 :
54 57 : STACTARasterBand::STACTARasterBand(STACTADataset *poDSIn, int nBandIn,
55 57 : GDALRasterBand *poProtoBand)
56 57 : : m_eColorInterp(poProtoBand->GetColorInterpretation())
57 : {
58 57 : poDS = poDSIn;
59 57 : nBand = nBandIn;
60 57 : eDataType = poProtoBand->GetRasterDataType();
61 57 : poProtoBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
62 57 : nRasterXSize = poDSIn->GetRasterXSize();
63 57 : nRasterYSize = poDSIn->GetRasterYSize();
64 57 : m_dfNoData = poProtoBand->GetNoDataValue(&m_bHasNoDataValue);
65 57 : }
66 :
67 : /************************************************************************/
68 : /* IReadBlock() */
69 : /************************************************************************/
70 :
71 0 : CPLErr STACTARasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
72 : void *pImage)
73 : {
74 0 : auto poGDS = cpl::down_cast<STACTADataset *>(poDS);
75 0 : return poGDS->m_poDS->GetRasterBand(nBand)->ReadBlock(nBlockXOff,
76 0 : nBlockYOff, pImage);
77 : }
78 :
79 : /************************************************************************/
80 : /* IRasterIO() */
81 : /************************************************************************/
82 :
83 4 : CPLErr STACTARasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
84 : int nXSize, int nYSize, void *pData,
85 : int nBufXSize, int nBufYSize,
86 : GDALDataType eBufType, GSpacing nPixelSpace,
87 : GSpacing nLineSpace,
88 : GDALRasterIOExtraArg *psExtraArg)
89 : {
90 4 : auto poGDS = cpl::down_cast<STACTADataset *>(poDS);
91 4 : if ((nBufXSize < nXSize || nBufYSize < nYSize) &&
92 8 : poGDS->m_apoOverviewDS.size() >= 1 && eRWFlag == GF_Read)
93 : {
94 : int bTried;
95 1 : CPLErr eErr = TryOverviewRasterIO(
96 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
97 : eBufType, nPixelSpace, nLineSpace, psExtraArg, &bTried);
98 1 : if (bTried)
99 1 : return eErr;
100 : }
101 :
102 3 : return poGDS->m_poDS->GetRasterBand(nBand)->RasterIO(
103 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
104 3 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
105 : }
106 :
107 : /************************************************************************/
108 : /* IRasterIO() */
109 : /************************************************************************/
110 :
111 12 : CPLErr STACTADataset::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
112 : int nXSize, int nYSize, void *pData,
113 : int nBufXSize, int nBufYSize,
114 : GDALDataType eBufType, int nBandCount,
115 : int *panBandMap, GSpacing nPixelSpace,
116 : GSpacing nLineSpace, GSpacing nBandSpace,
117 : GDALRasterIOExtraArg *psExtraArg)
118 : {
119 12 : if ((nBufXSize < nXSize || nBufYSize < nYSize) &&
120 24 : m_apoOverviewDS.size() >= 1 && eRWFlag == GF_Read)
121 : {
122 : int bTried;
123 6 : CPLErr eErr = TryOverviewRasterIO(
124 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
125 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
126 : nBandSpace, psExtraArg, &bTried);
127 6 : if (bTried)
128 2 : return eErr;
129 : }
130 :
131 10 : return m_poDS->RasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
132 : nBufXSize, nBufYSize, eBufType, nBandCount,
133 : panBandMap, nPixelSpace, nLineSpace, nBandSpace,
134 10 : psExtraArg);
135 : }
136 :
137 : /************************************************************************/
138 : /* GetOverviewCount() */
139 : /************************************************************************/
140 :
141 31 : int STACTARasterBand::GetOverviewCount()
142 : {
143 31 : STACTADataset *poGDS = cpl::down_cast<STACTADataset *>(poDS);
144 31 : return static_cast<int>(poGDS->m_apoOverviewDS.size());
145 : }
146 :
147 : /************************************************************************/
148 : /* GetOverview() */
149 : /************************************************************************/
150 :
151 21 : GDALRasterBand *STACTARasterBand::GetOverview(int nIdx)
152 : {
153 21 : STACTADataset *poGDS = cpl::down_cast<STACTADataset *>(poDS);
154 21 : if (nIdx < 0 || nIdx >= GetOverviewCount())
155 0 : return nullptr;
156 21 : return poGDS->m_apoOverviewDS[nIdx]->GetRasterBand(nBand);
157 : }
158 :
159 : /************************************************************************/
160 : /* GetNoDataValue() */
161 : /************************************************************************/
162 :
163 20 : double STACTARasterBand::GetNoDataValue(int *pbHasNoData)
164 : {
165 20 : if (pbHasNoData)
166 20 : *pbHasNoData = m_bHasNoDataValue;
167 20 : return m_dfNoData;
168 : }
169 :
170 : /************************************************************************/
171 : /* STACTARawRasterBand() */
172 : /************************************************************************/
173 :
174 81 : STACTARawRasterBand::STACTARawRasterBand(STACTARawDataset *poDSIn, int nBandIn,
175 81 : GDALRasterBand *poProtoBand)
176 81 : : m_eColorInterp(poProtoBand->GetColorInterpretation())
177 : {
178 81 : poDS = poDSIn;
179 81 : nBand = nBandIn;
180 81 : eDataType = poProtoBand->GetRasterDataType();
181 81 : nBlockXSize = 256;
182 81 : nBlockYSize = 256;
183 : int nProtoBlockXSize;
184 : int nProtoBlockYSize;
185 : // Use tile block size if it divides the metatile dimension.
186 81 : poProtoBand->GetBlockSize(&nProtoBlockXSize, &nProtoBlockYSize);
187 81 : if ((poDSIn->m_nMetaTileWidth % nProtoBlockXSize) == 0 &&
188 81 : (poDSIn->m_nMetaTileHeight % nProtoBlockYSize) == 0)
189 : {
190 81 : nBlockXSize = nProtoBlockXSize;
191 81 : nBlockYSize = nProtoBlockYSize;
192 : }
193 81 : nRasterXSize = poDSIn->GetRasterXSize();
194 81 : nRasterYSize = poDSIn->GetRasterYSize();
195 81 : m_dfNoData = poProtoBand->GetNoDataValue(&m_bHasNoDataValue);
196 81 : }
197 :
198 : /************************************************************************/
199 : /* STACTARawRasterBand() */
200 : /************************************************************************/
201 :
202 90 : STACTARawRasterBand::STACTARawRasterBand(STACTARawDataset *poDSIn, int nBandIn,
203 : GDALDataType eDT, bool bSetNoData,
204 90 : double dfNoData)
205 : {
206 90 : poDS = poDSIn;
207 90 : nBand = nBandIn;
208 90 : eDataType = eDT;
209 90 : nBlockXSize = 256;
210 90 : nBlockYSize = 256;
211 90 : nRasterXSize = poDSIn->GetRasterXSize();
212 90 : nRasterYSize = poDSIn->GetRasterYSize();
213 90 : m_bHasNoDataValue = bSetNoData;
214 90 : m_dfNoData = dfNoData;
215 90 : }
216 :
217 : /************************************************************************/
218 : /* GetNoDataValue() */
219 : /************************************************************************/
220 :
221 81 : double STACTARawRasterBand::GetNoDataValue(int *pbHasNoData)
222 : {
223 81 : if (pbHasNoData)
224 75 : *pbHasNoData = m_bHasNoDataValue;
225 81 : return m_dfNoData;
226 : }
227 :
228 : /************************************************************************/
229 : /* IReadBlock() */
230 : /************************************************************************/
231 :
232 0 : CPLErr STACTARawRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
233 : void *pImage)
234 : {
235 0 : const int nXOff = nBlockXOff * nBlockXSize;
236 0 : const int nYOff = nBlockYOff * nBlockYSize;
237 0 : const int nXSize = std::min(nBlockXSize, nRasterXSize - nXOff);
238 0 : const int nYSize = std::min(nBlockYSize, nRasterYSize - nYOff);
239 : GDALRasterIOExtraArg sExtraArgs;
240 0 : INIT_RASTERIO_EXTRA_ARG(sExtraArgs);
241 0 : const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
242 0 : return IRasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize, pImage, nBlockXSize,
243 : nBlockYSize, eDataType, nDTSize,
244 0 : static_cast<GSpacing>(nDTSize) * nBlockXSize, &sExtraArgs);
245 : }
246 :
247 : /************************************************************************/
248 : /* IRasterIO() */
249 : /************************************************************************/
250 :
251 7 : CPLErr STACTARawRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
252 : int nXSize, int nYSize, void *pData,
253 : int nBufXSize, int nBufYSize,
254 : GDALDataType eBufType,
255 : GSpacing nPixelSpace, GSpacing nLineSpace,
256 : GDALRasterIOExtraArg *psExtraArg)
257 : {
258 7 : CPLDebugOnly("STACTA", "Band %d RasterIO: %d,%d,%d,%d->%d,%d", nBand, nXOff,
259 : nYOff, nXSize, nYSize, nBufXSize, nBufYSize);
260 7 : auto poGDS = cpl::down_cast<STACTARawDataset *>(poDS);
261 :
262 7 : const int nKernelRadius = 3; // up to 3 for Lanczos
263 : const int nRadiusX =
264 7 : nKernelRadius * static_cast<int>(std::ceil(nXSize / nBufXSize));
265 : const int nRadiusY =
266 7 : nKernelRadius * static_cast<int>(std::ceil(nYSize / nBufYSize));
267 7 : const int nXOffMod = std::max(0, nXOff - nRadiusX);
268 7 : const int nYOffMod = std::max(0, nYOff - nRadiusY);
269 7 : const int nXSizeMod = static_cast<int>(std::min(
270 14 : nXOff + nXSize + static_cast<GIntBig>(nRadiusX),
271 7 : static_cast<GIntBig>(nRasterXSize))) -
272 7 : nXOffMod;
273 7 : const int nYSizeMod = static_cast<int>(std::min(
274 14 : nYOff + nYSize + static_cast<GIntBig>(nRadiusY),
275 7 : static_cast<GIntBig>(nRasterYSize))) -
276 7 : nYOffMod;
277 :
278 7 : const bool bRequestFitsInSingleMetaTile =
279 7 : nXOffMod / poGDS->m_nMetaTileWidth ==
280 11 : (nXOffMod + nXSizeMod - 1) / poGDS->m_nMetaTileWidth &&
281 4 : nYOffMod / poGDS->m_nMetaTileHeight ==
282 4 : (nYOffMod + nYSizeMod - 1) / poGDS->m_nMetaTileHeight;
283 :
284 7 : if (eRWFlag != GF_Read || ((nXSize != nBufXSize || nYSize != nBufYSize) &&
285 4 : !bRequestFitsInSingleMetaTile))
286 : {
287 0 : if (!(eRWFlag == GF_Read && nXSizeMod <= 4096 && nYSizeMod <= 4096))
288 : {
289 : // If not reading at nominal resolution, fallback to default block
290 : // reading
291 0 : return GDALRasterBand::IRasterIO(
292 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
293 0 : nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg);
294 : }
295 : }
296 :
297 : // Use optimized dataset level RasterIO()
298 7 : return poGDS->IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
299 : nBufXSize, nBufYSize, eBufType, 1, &nBand,
300 7 : nPixelSpace, nLineSpace, 0, psExtraArg);
301 : }
302 :
303 : /************************************************************************/
304 : /* IRasterIO() */
305 : /************************************************************************/
306 :
307 20 : CPLErr STACTARawDataset::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
308 : int nXSize, int nYSize, void *pData,
309 : int nBufXSize, int nBufYSize,
310 : GDALDataType eBufType, int nBandCount,
311 : int *panBandMap, GSpacing nPixelSpace,
312 : GSpacing nLineSpace, GSpacing nBandSpace,
313 : GDALRasterIOExtraArg *psExtraArg)
314 : {
315 20 : CPLDebugOnly("STACTA", "Dataset RasterIO: %d,%d,%d,%d->%d,%d", nXOff, nYOff,
316 : nXSize, nYSize, nBufXSize, nBufYSize);
317 20 : const int nMinBlockX = nXOff / m_nMetaTileWidth;
318 20 : const int nMaxBlockX = (nXOff + nXSize - 1) / m_nMetaTileWidth;
319 20 : const int nMinBlockY = nYOff / m_nMetaTileHeight;
320 20 : const int nMaxBlockY = (nYOff + nYSize - 1) / m_nMetaTileHeight;
321 :
322 20 : const int nKernelRadius = 3; // up to 3 for Lanczos
323 : const int nRadiusX =
324 20 : nKernelRadius * static_cast<int>(std::ceil(nXSize / nBufXSize));
325 : const int nRadiusY =
326 20 : nKernelRadius * static_cast<int>(std::ceil(nYSize / nBufYSize));
327 20 : const int nXOffMod = std::max(0, nXOff - nRadiusX);
328 20 : const int nYOffMod = std::max(0, nYOff - nRadiusY);
329 20 : const int nXSizeMod = static_cast<int>(std::min(
330 40 : nXOff + nXSize + static_cast<GIntBig>(nRadiusX),
331 20 : static_cast<GIntBig>(nRasterXSize))) -
332 20 : nXOffMod;
333 20 : const int nYSizeMod = static_cast<int>(std::min(
334 40 : nYOff + nYSize + static_cast<GIntBig>(nRadiusY),
335 20 : static_cast<GIntBig>(nRasterYSize))) -
336 20 : nYOffMod;
337 :
338 20 : const bool bRequestFitsInSingleMetaTile =
339 20 : nXOffMod / m_nMetaTileWidth ==
340 27 : (nXOffMod + nXSizeMod - 1) / m_nMetaTileWidth &&
341 7 : nYOffMod / m_nMetaTileHeight ==
342 7 : (nYOffMod + nYSizeMod - 1) / m_nMetaTileHeight;
343 20 : const auto eBandDT = GetRasterBand(1)->GetRasterDataType();
344 20 : const int nDTSize = GDALGetDataTypeSizeBytes(eBandDT);
345 :
346 20 : if (eRWFlag != GF_Read || ((nXSize != nBufXSize || nYSize != nBufYSize) &&
347 9 : !bRequestFitsInSingleMetaTile))
348 : {
349 2 : if (eRWFlag == GF_Read && nXSizeMod <= 4096 && nYSizeMod <= 4096 &&
350 : nBandCount <= 10)
351 : {
352 : // If extracting from a small enough window, do a RasterIO()
353 : // at full resolution into a MEM dataset, and then proceeding to
354 : // resampling on it. This will avoid to fallback on block based
355 : // approach.
356 : GDALRasterIOExtraArg sExtraArgs;
357 2 : INIT_RASTERIO_EXTRA_ARG(sExtraArgs);
358 2 : const size_t nXSizeModeMulYSizeModMulDTSize =
359 2 : static_cast<size_t>(nXSizeMod) * nYSizeMod * nDTSize;
360 : std::vector<GByte> abyBuf(nXSizeModeMulYSizeModMulDTSize *
361 4 : nBandCount);
362 2 : if (IRasterIO(GF_Read, nXOffMod, nYOffMod, nXSizeMod, nYSizeMod,
363 2 : &abyBuf[0], nXSizeMod, nYSizeMod, eBandDT, nBandCount,
364 : panBandMap, nDTSize,
365 2 : static_cast<GSpacing>(nDTSize) * nXSizeMod,
366 2 : static_cast<GSpacing>(nDTSize) * nXSizeMod *
367 2 : nYSizeMod,
368 2 : &sExtraArgs) != CE_None)
369 : {
370 0 : return CE_Failure;
371 : }
372 :
373 : auto poMEMDS = std::unique_ptr<MEMDataset>(MEMDataset::Create(
374 4 : "", nXSizeMod, nYSizeMod, 0, eBandDT, nullptr));
375 8 : for (int i = 0; i < nBandCount; i++)
376 : {
377 6 : auto hBand = MEMCreateRasterBandEx(
378 6 : poMEMDS.get(), i + 1,
379 6 : &abyBuf[0] + i * nXSizeModeMulYSizeModMulDTSize, eBandDT, 0,
380 : 0, false);
381 6 : poMEMDS->AddMEMBand(hBand);
382 : }
383 :
384 2 : sExtraArgs.eResampleAlg = psExtraArg->eResampleAlg;
385 2 : if (psExtraArg->bFloatingPointWindowValidity)
386 : {
387 0 : sExtraArgs.bFloatingPointWindowValidity = true;
388 0 : sExtraArgs.dfXOff = psExtraArg->dfXOff - nXOffMod;
389 0 : sExtraArgs.dfYOff = psExtraArg->dfYOff - nYOffMod;
390 0 : sExtraArgs.dfXSize = psExtraArg->dfXSize;
391 0 : sExtraArgs.dfYSize = psExtraArg->dfYSize;
392 : }
393 2 : return poMEMDS->RasterIO(
394 : GF_Read, nXOff - nXOffMod, nYOff - nYOffMod, nXSize, nYSize,
395 : pData, nBufXSize, nBufYSize, eBufType, nBandCount, nullptr,
396 2 : nPixelSpace, nLineSpace, nBandSpace, &sExtraArgs);
397 : }
398 :
399 : // If not reading at nominal resolution, fallback to default block
400 : // reading
401 0 : return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
402 : pData, nBufXSize, nBufYSize, eBufType,
403 : nBandCount, panBandMap, nPixelSpace,
404 0 : nLineSpace, nBandSpace, psExtraArg);
405 : }
406 :
407 18 : int nBufYOff = 0;
408 :
409 : // If the (uncompressed) size of a metatile is small enough, then download
410 : // it entirely to minimize the number of network requests
411 18 : const bool bDownloadWholeMetaTile =
412 33 : m_poMasterDS->m_bDownloadWholeMetaTile ||
413 15 : (static_cast<GIntBig>(m_nMetaTileWidth) * m_nMetaTileHeight * nBands *
414 15 : nDTSize <
415 : 128 * 1024);
416 :
417 : // Split the request on each metatile that it intersects
418 35 : for (int iY = nMinBlockY; iY <= nMaxBlockY; iY++)
419 : {
420 18 : const int nTileYOff = std::max(0, nYOff - iY * m_nMetaTileHeight);
421 : const int nTileYSize =
422 18 : std::min((iY + 1) * m_nMetaTileHeight, nYOff + nYSize) -
423 18 : std::max(nYOff, iY * m_nMetaTileHeight);
424 :
425 18 : int nBufXOff = 0;
426 46 : for (int iX = nMinBlockX; iX <= nMaxBlockX; iX++)
427 : {
428 29 : CPLString osURL(m_osURLTemplate);
429 : osURL.replaceAll("{TileRow}",
430 29 : CPLSPrintf("%d", iY + m_nMinMetaTileRow));
431 : osURL.replaceAll("{TileCol}",
432 29 : CPLSPrintf("%d", iX + m_nMinMetaTileCol));
433 :
434 29 : const int nTileXOff = std::max(0, nXOff - iX * m_nMetaTileWidth);
435 : const int nTileXSize =
436 29 : std::min((iX + 1) * m_nMetaTileWidth, nXOff + nXSize) -
437 29 : std::max(nXOff, iX * m_nMetaTileWidth);
438 :
439 29 : const int nBufXSizeEffective =
440 29 : bRequestFitsInSingleMetaTile ? nBufXSize : nTileXSize;
441 29 : const int nBufYSizeEffective =
442 29 : bRequestFitsInSingleMetaTile ? nBufYSize : nTileYSize;
443 :
444 29 : bool bMissingTile = false;
445 : do
446 : {
447 : std::unique_ptr<GDALDataset> *ppoTileDS =
448 29 : m_poMasterDS->m_oCacheTileDS.getPtr(osURL);
449 29 : if (ppoTileDS == nullptr)
450 : {
451 :
452 : // Avoid probing side car files
453 : CPLConfigOptionSetter oSetter(
454 : "GDAL_DISABLE_READDIR_ON_OPEN", "EMPTY_DIR",
455 15 : /* bSetOnlyIfUndefined = */ true);
456 :
457 15 : CPLStringList aosAllowedDrivers;
458 15 : aosAllowedDrivers.AddString("GTiff");
459 15 : aosAllowedDrivers.AddString("PNG");
460 15 : aosAllowedDrivers.AddString("JPEG");
461 15 : aosAllowedDrivers.AddString("JP2KAK");
462 15 : aosAllowedDrivers.AddString("JP2ECW");
463 15 : aosAllowedDrivers.AddString("JP2MrSID");
464 15 : aosAllowedDrivers.AddString("JP2OpenJPEG");
465 0 : std::unique_ptr<GDALDataset> poTileDS;
466 15 : if (bDownloadWholeMetaTile && !VSIIsLocal(osURL.c_str()))
467 : {
468 0 : if (m_poMasterDS->m_bSkipMissingMetaTile)
469 0 : CPLPushErrorHandler(CPLQuietErrorHandler);
470 0 : VSILFILE *fp = VSIFOpenL(osURL, "rb");
471 0 : if (m_poMasterDS->m_bSkipMissingMetaTile)
472 0 : CPLPopErrorHandler();
473 0 : if (fp == nullptr)
474 : {
475 0 : if (m_poMasterDS->m_bSkipMissingMetaTile)
476 : {
477 0 : m_poMasterDS->m_oCacheTileDS.insert(osURL,
478 0 : nullptr);
479 0 : bMissingTile = true;
480 0 : break;
481 : }
482 0 : CPLError(CE_Failure, CPLE_OpenFailed,
483 : "Cannot open %s", osURL.c_str());
484 0 : return CE_Failure;
485 : }
486 0 : GByte *pabyBuf = nullptr;
487 0 : vsi_l_offset nSize = 0;
488 0 : if (!VSIIngestFile(fp, nullptr, &pabyBuf, &nSize, -1))
489 : {
490 0 : VSIFCloseL(fp);
491 0 : return CE_Failure;
492 : }
493 0 : VSIFCloseL(fp);
494 0 : const CPLString osMEMFilename("/vsimem/stacta/" +
495 0 : osURL);
496 0 : VSIFCloseL(VSIFileFromMemBuffer(osMEMFilename, pabyBuf,
497 : nSize, TRUE));
498 0 : poTileDS = std::unique_ptr<GDALDataset>(
499 : GDALDataset::Open(osMEMFilename,
500 : GDAL_OF_INTERNAL | GDAL_OF_RASTER,
501 0 : aosAllowedDrivers.List()));
502 0 : if (poTileDS)
503 0 : poTileDS->MarkSuppressOnClose();
504 : else
505 0 : VSIUnlink(osMEMFilename);
506 : }
507 29 : else if (bDownloadWholeMetaTile ||
508 14 : (!STARTS_WITH(osURL, "http://") &&
509 12 : !STARTS_WITH(osURL, "https://")))
510 : {
511 13 : aosAllowedDrivers.AddString("HTTP");
512 13 : if (m_poMasterDS->m_bSkipMissingMetaTile)
513 4 : CPLPushErrorHandler(CPLQuietErrorHandler);
514 : poTileDS =
515 26 : std::unique_ptr<GDALDataset>(GDALDataset::Open(
516 : osURL, GDAL_OF_INTERNAL | GDAL_OF_RASTER,
517 26 : aosAllowedDrivers.List()));
518 13 : if (m_poMasterDS->m_bSkipMissingMetaTile)
519 4 : CPLPopErrorHandler();
520 : }
521 : else
522 : {
523 2 : if (m_poMasterDS->m_bSkipMissingMetaTile)
524 0 : CPLPushErrorHandler(CPLQuietErrorHandler);
525 4 : poTileDS = std::unique_ptr<GDALDataset>(
526 4 : GDALDataset::Open(("/vsicurl/" + osURL).c_str(),
527 : GDAL_OF_INTERNAL | GDAL_OF_RASTER,
528 4 : aosAllowedDrivers.List()));
529 2 : if (m_poMasterDS->m_bSkipMissingMetaTile)
530 0 : CPLPopErrorHandler();
531 : }
532 15 : if (poTileDS == nullptr)
533 : {
534 3 : if (m_poMasterDS->m_bSkipMissingMetaTile)
535 : {
536 2 : m_poMasterDS->m_oCacheTileDS.insert(
537 2 : osURL, std::move(poTileDS));
538 2 : bMissingTile = true;
539 2 : break;
540 : }
541 1 : CPLError(CE_Failure, CPLE_OpenFailed, "Cannot open %s",
542 : osURL.c_str());
543 1 : return CE_Failure;
544 : }
545 12 : ppoTileDS = &m_poMasterDS->m_oCacheTileDS.insert(
546 12 : osURL, std::move(poTileDS));
547 : }
548 26 : std::unique_ptr<GDALDataset> &poTileDS = *ppoTileDS;
549 26 : if (poTileDS == nullptr)
550 : {
551 0 : bMissingTile = true;
552 0 : break;
553 : }
554 :
555 : GDALRasterIOExtraArg sExtraArgs;
556 26 : INIT_RASTERIO_EXTRA_ARG(sExtraArgs);
557 26 : if (bRequestFitsInSingleMetaTile)
558 : {
559 7 : sExtraArgs.eResampleAlg = psExtraArg->eResampleAlg;
560 7 : if (psExtraArg->bFloatingPointWindowValidity)
561 : {
562 5 : sExtraArgs.bFloatingPointWindowValidity = true;
563 5 : sExtraArgs.dfXOff =
564 5 : psExtraArg->dfXOff - iX * m_nMetaTileWidth;
565 5 : sExtraArgs.dfYOff =
566 5 : psExtraArg->dfYOff - iY * m_nMetaTileHeight;
567 5 : sExtraArgs.dfXSize = psExtraArg->dfXSize;
568 5 : sExtraArgs.dfYSize = psExtraArg->dfYSize;
569 : }
570 : }
571 26 : CPLDebugOnly("STACTA", "Reading %d,%d,%d,%d in %s", nTileXOff,
572 : nTileYOff, nTileXSize, nTileYSize, osURL.c_str());
573 26 : if (poTileDS->RasterIO(
574 : GF_Read, nTileXOff, nTileYOff, nTileXSize, nTileYSize,
575 26 : static_cast<GByte *>(pData) + nBufXOff * nPixelSpace +
576 26 : nBufYOff * nLineSpace,
577 : nBufXSizeEffective, nBufYSizeEffective, eBufType,
578 : nBandCount, panBandMap, nPixelSpace, nLineSpace,
579 26 : nBandSpace, &sExtraArgs) != CE_None)
580 : {
581 0 : return CE_Failure;
582 : }
583 : } while (false);
584 :
585 28 : if (bMissingTile)
586 : {
587 2 : CPLDebugOnly("STACTA", "Missing metatile %s", osURL.c_str());
588 8 : for (int iBand = 0; iBand < nBandCount; iBand++)
589 : {
590 6 : int bHasNoData = FALSE;
591 6 : double dfNodata = GetRasterBand(panBandMap[iBand])
592 6 : ->GetNoDataValue(&bHasNoData);
593 6 : if (!bHasNoData)
594 0 : dfNodata = 0;
595 6150 : for (int nYBufOff = 0; nYBufOff < nBufYSizeEffective;
596 : nYBufOff++)
597 : {
598 6144 : GByte *pabyDest = static_cast<GByte *>(pData) +
599 6144 : iBand * nBandSpace +
600 6144 : nBufXOff * nPixelSpace +
601 6144 : (nBufYOff + nYBufOff) * nLineSpace;
602 6144 : GDALCopyWords(&dfNodata, GDT_Float64, 0, pabyDest,
603 : eBufType, static_cast<int>(nPixelSpace),
604 : nBufXSizeEffective);
605 : }
606 : }
607 : }
608 :
609 28 : if (iX == nMinBlockX)
610 : {
611 36 : nBufXOff = m_nMetaTileWidth -
612 18 : std::max(0, nXOff - nMinBlockX * m_nMetaTileWidth);
613 : }
614 : else
615 : {
616 10 : nBufXOff += m_nMetaTileWidth;
617 : }
618 : }
619 :
620 17 : if (iY == nMinBlockY)
621 : {
622 34 : nBufYOff = m_nMetaTileHeight -
623 17 : std::max(0, nYOff - nMinBlockY * m_nMetaTileHeight);
624 : }
625 : else
626 : {
627 0 : nBufYOff += m_nMetaTileHeight;
628 : }
629 : }
630 :
631 17 : return CE_None;
632 : }
633 :
634 : /************************************************************************/
635 : /* GetGeoTransform() */
636 : /************************************************************************/
637 :
638 4 : CPLErr STACTARawDataset::GetGeoTransform(double *padfGeoTransform)
639 : {
640 4 : memcpy(padfGeoTransform, &m_adfGeoTransform[0], 6 * sizeof(double));
641 4 : return CE_None;
642 : }
643 :
644 : /************************************************************************/
645 : /* Identify() */
646 : /************************************************************************/
647 :
648 47854 : int STACTADataset::Identify(GDALOpenInfo *poOpenInfo)
649 : {
650 47854 : if (STARTS_WITH(poOpenInfo->pszFilename, "STACTA:"))
651 : {
652 6 : return true;
653 : }
654 :
655 47843 : if (
656 : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
657 48245 : !EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "json") ||
658 : #endif
659 397 : poOpenInfo->nHeaderBytes == 0)
660 : {
661 47724 : return false;
662 : }
663 :
664 291 : for (int i = 0; i < 2; i++)
665 : {
666 : // TryToIngest() may reallocate pabyHeader, so do not move this
667 : // before the loop.
668 206 : const char *pszHeader =
669 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
670 206 : if (strstr(pszHeader, "\"stac_extensions\"") != nullptr &&
671 44 : (strstr(pszHeader, "\"tiled-assets\"") != nullptr ||
672 28 : strstr(pszHeader,
673 : "https://stac-extensions.github.io/tiled-assets/") !=
674 : nullptr))
675 : {
676 36 : return true;
677 : }
678 :
679 170 : if (i == 0)
680 : {
681 : // Should be enough for a STACTA .json file
682 85 : poOpenInfo->TryToIngest(32768);
683 : }
684 : }
685 :
686 85 : return false;
687 : }
688 :
689 : /************************************************************************/
690 : /* Open() */
691 : /************************************************************************/
692 :
693 21 : bool STACTADataset::Open(GDALOpenInfo *poOpenInfo)
694 : {
695 42 : CPLString osFilename(poOpenInfo->pszFilename);
696 42 : CPLString osAssetName;
697 42 : CPLString osTMS;
698 21 : if (STARTS_WITH(poOpenInfo->pszFilename, "STACTA:"))
699 : {
700 : const CPLStringList aosTokens(CSLTokenizeString2(
701 3 : poOpenInfo->pszFilename, ":", CSLT_HONOURSTRINGS));
702 4 : if (aosTokens.size() != 2 && aosTokens.size() != 3 &&
703 1 : aosTokens.size() != 4)
704 0 : return false;
705 3 : osFilename = aosTokens[1];
706 3 : if (aosTokens.size() >= 3)
707 3 : osAssetName = aosTokens[2];
708 3 : if (aosTokens.size() == 4)
709 1 : osTMS = aosTokens[3];
710 : }
711 :
712 42 : CPLJSONDocument oDoc;
713 42 : if (STARTS_WITH(osFilename, "http://") ||
714 21 : STARTS_WITH(osFilename, "https://"))
715 : {
716 0 : if (!oDoc.LoadUrl(osFilename, nullptr))
717 0 : return false;
718 : }
719 : else
720 : {
721 21 : if (!oDoc.Load(osFilename))
722 0 : return false;
723 : }
724 42 : const auto oRoot = oDoc.GetRoot();
725 63 : const auto oProperties = oRoot["properties"];
726 42 : if (!oProperties.IsValid() ||
727 21 : oProperties.GetType() != CPLJSONObject::Type::Object)
728 : {
729 0 : CPLError(CE_Failure, CPLE_AppDefined, "Missing properties");
730 0 : return false;
731 : }
732 :
733 63 : const auto oAssetTemplates = oRoot["asset_templates"];
734 42 : if (!oAssetTemplates.IsValid() ||
735 21 : oAssetTemplates.GetType() != CPLJSONObject::Type::Object)
736 : {
737 0 : CPLError(CE_Failure, CPLE_AppDefined, "Missing asset_templates");
738 0 : return false;
739 : }
740 :
741 42 : const auto aoAssetTemplates = oAssetTemplates.GetChildren();
742 21 : if (aoAssetTemplates.size() == 0)
743 : {
744 0 : CPLError(CE_Failure, CPLE_AppDefined, "Empty asset_templates");
745 0 : return false;
746 : }
747 :
748 63 : const auto oTMSs = oProperties.GetObj("tiles:tile_matrix_sets");
749 21 : if (!oTMSs.IsValid() || oTMSs.GetType() != CPLJSONObject::Type::Object)
750 : {
751 0 : CPLError(CE_Failure, CPLE_AppDefined,
752 : "Missing properties[\"tiles:tile_matrix_sets\"]");
753 0 : return false;
754 : }
755 42 : const auto aoTMSs = oTMSs.GetChildren();
756 21 : if (aoTMSs.empty())
757 : {
758 0 : CPLError(CE_Failure, CPLE_AppDefined,
759 : "Empty properties[\"tiles:tile_matrix_sets\"]");
760 0 : return false;
761 : }
762 :
763 44 : if ((aoAssetTemplates.size() >= 2 || aoTMSs.size() >= 2) &&
764 44 : osAssetName.empty() && osTMS.empty())
765 : {
766 2 : int nSDSCount = 0;
767 5 : for (const auto &oAssetTemplate : aoAssetTemplates)
768 : {
769 6 : const CPLString osAssetNameSubDS = oAssetTemplate.GetName();
770 3 : const char *pszAssetNameSubDS = osAssetNameSubDS.c_str();
771 3 : if (aoTMSs.size() >= 2)
772 : {
773 3 : for (const auto &oTMS : aoTMSs)
774 : {
775 2 : const CPLString osTMSSubDS = oTMS.GetName();
776 2 : const char *pszTMSSubDS = osTMSSubDS.c_str();
777 2 : GDALDataset::SetMetadataItem(
778 : CPLSPrintf("SUBDATASET_%d_NAME", nSDSCount + 1),
779 : CPLSPrintf("STACTA:\"%s\":%s:%s", osFilename.c_str(),
780 : pszAssetNameSubDS, pszTMSSubDS),
781 : "SUBDATASETS");
782 2 : GDALDataset::SetMetadataItem(
783 : CPLSPrintf("SUBDATASET_%d_DESC", nSDSCount + 1),
784 : CPLSPrintf("Asset %s, tile matrix set %s",
785 : pszAssetNameSubDS, pszTMSSubDS),
786 : "SUBDATASETS");
787 2 : nSDSCount++;
788 : }
789 : }
790 : else
791 : {
792 2 : GDALDataset::SetMetadataItem(
793 : CPLSPrintf("SUBDATASET_%d_NAME", nSDSCount + 1),
794 : CPLSPrintf("STACTA:\"%s\":%s", osFilename.c_str(),
795 : pszAssetNameSubDS),
796 : "SUBDATASETS");
797 2 : GDALDataset::SetMetadataItem(
798 : CPLSPrintf("SUBDATASET_%d_DESC", nSDSCount + 1),
799 : CPLSPrintf("Asset %s", pszAssetNameSubDS), "SUBDATASETS");
800 2 : nSDSCount++;
801 : }
802 : }
803 2 : return true;
804 : }
805 :
806 19 : if (osAssetName.empty())
807 : {
808 16 : osAssetName = aoAssetTemplates[0].GetName();
809 : }
810 38 : const auto oAssetTemplate = oAssetTemplates.GetObj(osAssetName);
811 38 : if (!oAssetTemplate.IsValid() ||
812 19 : oAssetTemplate.GetType() != CPLJSONObject::Type::Object)
813 : {
814 0 : CPLError(CE_Failure, CPLE_AppDefined,
815 : "Cannot find asset_templates[\"%s\"]", osAssetName.c_str());
816 0 : return false;
817 : }
818 :
819 19 : if (osTMS.empty())
820 : {
821 18 : osTMS = aoTMSs[0].GetName();
822 : }
823 38 : const auto oTMS = oTMSs.GetObj(osTMS);
824 19 : if (!oTMS.IsValid() || oTMS.GetType() != CPLJSONObject::Type::Object)
825 : {
826 0 : CPLError(CE_Failure, CPLE_AppDefined,
827 : "Cannot find properties[\"tiles:tile_matrix_sets\"][\"%s\"]",
828 : osTMS.c_str());
829 0 : return false;
830 : }
831 :
832 : auto poTMS = gdal::TileMatrixSet::parse(
833 38 : oTMS.Format(CPLJSONObject::PrettyFormat::Plain).c_str());
834 19 : if (poTMS == nullptr)
835 0 : return false;
836 :
837 57 : CPLString osURLTemplate = oAssetTemplate.GetString("href");
838 19 : if (osURLTemplate.empty())
839 : {
840 0 : CPLError(CE_Failure, CPLE_AppDefined,
841 : "Cannot find asset_templates[\"%s\"][\"href\"]",
842 : osAssetName.c_str());
843 : }
844 19 : osURLTemplate.replaceAll("{TileMatrixSet}", osTMS);
845 19 : if (STARTS_WITH(osURLTemplate, "file://"))
846 : {
847 0 : osURLTemplate = osURLTemplate.substr(strlen("file://"));
848 : }
849 19 : else if (STARTS_WITH(osURLTemplate, "s3://"))
850 : {
851 0 : osURLTemplate = "/vsis3/" + osURLTemplate.substr(strlen("s3://"));
852 : }
853 :
854 37 : if (!STARTS_WITH(osURLTemplate, "http://") &&
855 18 : !STARTS_WITH(osURLTemplate, "https://"))
856 : {
857 18 : if (STARTS_WITH(osURLTemplate, "./"))
858 18 : osURLTemplate = osURLTemplate.substr(2);
859 : osURLTemplate = CPLProjectRelativeFilename(CPLGetDirname(osFilename),
860 18 : osURLTemplate);
861 : }
862 :
863 : // Parse optional tile matrix set limits
864 38 : std::map<CPLString, Limits> oMapLimits;
865 57 : const auto oTMLinks = oProperties.GetObj("tiles:tile_matrix_links");
866 19 : if (oTMLinks.IsValid())
867 : {
868 16 : if (oTMLinks.GetType() != CPLJSONObject::Type::Object)
869 : {
870 0 : CPLError(
871 : CE_Failure, CPLE_AppDefined,
872 : "Invalid type for properties[\"tiles:tile_matrix_links\"]");
873 0 : return false;
874 : }
875 :
876 48 : auto oLimits = oTMLinks[osTMS]["limits"];
877 32 : if (oLimits.IsValid() &&
878 16 : oLimits.GetType() == CPLJSONObject::Type::Object)
879 : {
880 64 : for (const auto &oLimit : oLimits.GetChildren())
881 : {
882 48 : Limits limits;
883 48 : limits.min_tile_col = oLimit.GetInteger("min_tile_col");
884 48 : limits.max_tile_col = oLimit.GetInteger("max_tile_col");
885 48 : limits.min_tile_row = oLimit.GetInteger("min_tile_row");
886 48 : limits.max_tile_row = oLimit.GetInteger("max_tile_row");
887 48 : oMapLimits[oLimit.GetName()] = limits;
888 : }
889 : }
890 : }
891 19 : const auto &tmsList = poTMS->tileMatrixList();
892 19 : if (tmsList.empty())
893 0 : return false;
894 :
895 38 : m_bSkipMissingMetaTile = CPLTestBool(CSLFetchNameValueDef(
896 19 : poOpenInfo->papszOpenOptions, "SKIP_MISSING_METATILE",
897 : CPLGetConfigOption("GDAL_STACTA_SKIP_MISSING_METATILE", "NO")));
898 :
899 : // STAC 1.1 uses bands instead of eo:bands and raster:bands
900 57 : const auto oBands = oAssetTemplate.GetArray("bands");
901 :
902 : // Check if there are both eo:bands and raster:bands extension
903 : // If so, we don't need to fetch a prototype metatile to derive the
904 : // information we need (number of bands, data type and nodata value)
905 : const auto oEoBands =
906 56 : oBands.IsValid() ? oBands : oAssetTemplate.GetArray("eo:bands");
907 : const auto oRasterBands =
908 56 : oBands.IsValid() ? oBands : oAssetTemplate.GetArray("raster:bands");
909 :
910 38 : std::vector<GDALDataType> aeDT;
911 38 : std::vector<double> adfNoData;
912 38 : std::vector<bool> abSetNoData;
913 19 : int nExpectedBandCount = 0;
914 19 : if (oRasterBands.IsValid())
915 : {
916 9 : if (oEoBands.IsValid() && oEoBands.Size() != oRasterBands.Size())
917 : {
918 1 : CPLError(CE_Warning, CPLE_AppDefined,
919 : "Number of bands in eo:bands and raster:bands is not "
920 : "identical. Ignoring the later");
921 : }
922 : else
923 : {
924 8 : nExpectedBandCount = oRasterBands.Size();
925 :
926 : const struct
927 : {
928 : const char *pszStacDataType;
929 : GDALDataType eGDALDataType;
930 8 : } aDataTypeMapping[] = {
931 : {"int8", GDT_Int8},
932 : {"int16", GDT_Int16},
933 : {"int32", GDT_Int32},
934 : {"int64", GDT_Int64},
935 : {"uint8", GDT_Byte},
936 : {"uint16", GDT_UInt16},
937 : {"uint32", GDT_UInt32},
938 : {"uint64", GDT_UInt64},
939 : // float16: 16-bit float; unhandled
940 : {"float32", GDT_Float32},
941 : {"float64", GDT_Float64},
942 : {"cint16", GDT_CInt16},
943 : {"cint32", GDT_CInt32},
944 : {"cfloat32", GDT_CFloat32},
945 : {"cfloat64", GDT_CFloat64},
946 : };
947 :
948 38 : for (int i = 0; i < nExpectedBandCount; ++i)
949 : {
950 33 : if (oRasterBands[i].GetType() != CPLJSONObject::Type::Object)
951 : {
952 1 : CPLError(CE_Failure, CPLE_AppDefined,
953 : "Wrong raster:bands[%d]", i);
954 3 : return false;
955 : }
956 : const std::string osDataType =
957 64 : oRasterBands[i].GetString("data_type");
958 32 : GDALDataType eDT = GDT_Unknown;
959 280 : for (const auto &oTuple : aDataTypeMapping)
960 : {
961 278 : if (osDataType == oTuple.pszStacDataType)
962 : {
963 30 : eDT = oTuple.eGDALDataType;
964 30 : break;
965 : }
966 : }
967 32 : if (eDT == GDT_Unknown)
968 : {
969 2 : CPLError(CE_Failure, CPLE_AppDefined,
970 : "Wrong raster:bands[%d].data_type = %s", i,
971 : osDataType.c_str());
972 2 : return false;
973 : }
974 30 : aeDT.push_back(eDT);
975 :
976 90 : const auto oNoData = oRasterBands[i].GetObj("nodata");
977 30 : if (oNoData.GetType() == CPLJSONObject::Type::String)
978 : {
979 48 : const std::string osNoData = oNoData.ToString();
980 16 : if (osNoData == "inf")
981 : {
982 5 : abSetNoData.push_back(true);
983 5 : adfNoData.push_back(
984 5 : std::numeric_limits<double>::infinity());
985 : }
986 11 : else if (osNoData == "-inf")
987 : {
988 5 : abSetNoData.push_back(true);
989 5 : adfNoData.push_back(
990 5 : -std::numeric_limits<double>::infinity());
991 : }
992 6 : else if (osNoData == "nan")
993 : {
994 5 : abSetNoData.push_back(true);
995 5 : adfNoData.push_back(
996 5 : std::numeric_limits<double>::quiet_NaN());
997 : }
998 : else
999 : {
1000 1 : CPLError(CE_Warning, CPLE_AppDefined,
1001 : "Invalid raster:bands[%d].nodata = %s", i,
1002 : osNoData.c_str());
1003 1 : abSetNoData.push_back(false);
1004 1 : adfNoData.push_back(
1005 1 : std::numeric_limits<double>::quiet_NaN());
1006 : }
1007 : }
1008 14 : else if (oNoData.GetType() == CPLJSONObject::Type::Integer ||
1009 25 : oNoData.GetType() == CPLJSONObject::Type::Long ||
1010 11 : oNoData.GetType() == CPLJSONObject::Type::Double)
1011 : {
1012 8 : abSetNoData.push_back(true);
1013 8 : adfNoData.push_back(oNoData.ToDouble());
1014 : }
1015 6 : else if (!oNoData.IsValid())
1016 : {
1017 5 : abSetNoData.push_back(false);
1018 5 : adfNoData.push_back(
1019 5 : std::numeric_limits<double>::quiet_NaN());
1020 : }
1021 : else
1022 : {
1023 1 : CPLError(CE_Warning, CPLE_AppDefined,
1024 : "Invalid raster:bands[%d].nodata", i);
1025 1 : abSetNoData.push_back(false);
1026 1 : adfNoData.push_back(
1027 1 : std::numeric_limits<double>::quiet_NaN());
1028 : }
1029 : }
1030 :
1031 5 : CPLAssert(aeDT.size() == abSetNoData.size());
1032 5 : CPLAssert(adfNoData.size() == abSetNoData.size());
1033 : }
1034 : }
1035 :
1036 16 : std::unique_ptr<GDALDataset> poProtoDS;
1037 16 : if (aeDT.empty())
1038 : {
1039 11 : for (int i = 0; i < static_cast<int>(tmsList.size()); i++)
1040 : {
1041 : // Open a metatile to get mostly its band data type
1042 11 : int nProtoTileCol = 0;
1043 11 : int nProtoTileRow = 0;
1044 11 : auto oIterLimit = oMapLimits.find(tmsList[i].mId);
1045 11 : if (oIterLimit != oMapLimits.end())
1046 : {
1047 8 : nProtoTileCol = oIterLimit->second.min_tile_col;
1048 8 : nProtoTileRow = oIterLimit->second.min_tile_row;
1049 : }
1050 : const CPLString osURL =
1051 11 : CPLString(osURLTemplate)
1052 22 : .replaceAll("{TileMatrix}", tmsList[i].mId)
1053 22 : .replaceAll("{TileRow}", CPLSPrintf("%d", nProtoTileRow))
1054 22 : .replaceAll("{TileCol}", CPLSPrintf("%d", nProtoTileCol));
1055 32 : CPLString osProtoDSName = (STARTS_WITH(osURL, "http://") ||
1056 10 : STARTS_WITH(osURL, "https://"))
1057 22 : ? CPLString("/vsicurl/" + osURL)
1058 22 : : osURL;
1059 : CPLConfigOptionSetter oSetter("GDAL_DISABLE_READDIR_ON_OPEN",
1060 : "EMPTY_DIR",
1061 11 : /* bSetOnlyIfUndefined = */ true);
1062 11 : if (m_bSkipMissingMetaTile)
1063 2 : CPLPushErrorHandler(CPLQuietErrorHandler);
1064 11 : poProtoDS.reset(GDALDataset::Open(osProtoDSName.c_str()));
1065 11 : if (m_bSkipMissingMetaTile)
1066 2 : CPLPopErrorHandler();
1067 11 : if (poProtoDS != nullptr)
1068 : {
1069 9 : break;
1070 : }
1071 2 : if (!m_bSkipMissingMetaTile)
1072 : {
1073 2 : CPLError(CE_Failure, CPLE_OpenFailed, "Cannot open %s",
1074 : osURL.c_str());
1075 2 : return false;
1076 : }
1077 : }
1078 9 : if (poProtoDS == nullptr)
1079 : {
1080 0 : CPLError(CE_Failure, CPLE_AppDefined,
1081 : "Cannot find prototype dataset");
1082 0 : return false;
1083 : }
1084 9 : nExpectedBandCount = poProtoDS->GetRasterCount();
1085 : }
1086 :
1087 : // Iterate over tile matrices to create corresponding STACTARawDataset
1088 : // objects
1089 56 : for (int i = static_cast<int>(tmsList.size() - 1); i >= 0; i--)
1090 : {
1091 42 : const auto &oTM = tmsList[i];
1092 42 : int nMatrixWidth = oTM.mMatrixWidth;
1093 42 : int nMatrixHeight = oTM.mMatrixHeight;
1094 42 : auto oIterLimit = oMapLimits.find(tmsList[i].mId);
1095 42 : if (oIterLimit != oMapLimits.end())
1096 : {
1097 33 : nMatrixWidth = oIterLimit->second.max_tile_col -
1098 33 : oIterLimit->second.min_tile_col + 1;
1099 33 : nMatrixHeight = oIterLimit->second.max_tile_row -
1100 33 : oIterLimit->second.min_tile_row + 1;
1101 : }
1102 42 : if (nMatrixWidth <= 0 || oTM.mTileWidth > INT_MAX / nMatrixWidth ||
1103 42 : nMatrixHeight <= 0 || oTM.mTileHeight > INT_MAX / nMatrixHeight)
1104 : {
1105 0 : continue;
1106 : }
1107 42 : auto poRawDS = std::make_unique<STACTARawDataset>();
1108 84 : if (!poRawDS->InitRaster(poProtoDS.get(), aeDT, abSetNoData, adfNoData,
1109 42 : poTMS.get(), tmsList[i].mId, oTM, oMapLimits))
1110 : {
1111 0 : return false;
1112 : }
1113 42 : poRawDS->m_osURLTemplate = osURLTemplate;
1114 42 : poRawDS->m_osURLTemplate.replaceAll("{TileMatrix}", tmsList[i].mId);
1115 42 : poRawDS->m_poMasterDS = this;
1116 :
1117 42 : if (m_poDS == nullptr)
1118 : {
1119 14 : nRasterXSize = poRawDS->GetRasterXSize();
1120 14 : nRasterYSize = poRawDS->GetRasterYSize();
1121 14 : m_oSRS = poRawDS->m_oSRS;
1122 14 : memcpy(&m_adfGeoTransform[0], &poRawDS->m_adfGeoTransform[0],
1123 : 6 * sizeof(double));
1124 14 : m_poDS = std::move(poRawDS);
1125 : }
1126 : else
1127 : {
1128 28 : const double dfMinX = m_adfGeoTransform[0];
1129 : const double dfMaxX =
1130 28 : m_adfGeoTransform[0] + GetRasterXSize() * m_adfGeoTransform[1];
1131 28 : const double dfMaxY = m_adfGeoTransform[3];
1132 : const double dfMinY =
1133 28 : m_adfGeoTransform[3] + GetRasterYSize() * m_adfGeoTransform[5];
1134 :
1135 28 : const double dfOvrMinX = poRawDS->m_adfGeoTransform[0];
1136 : const double dfOvrMaxX =
1137 28 : poRawDS->m_adfGeoTransform[0] +
1138 28 : poRawDS->GetRasterXSize() * poRawDS->m_adfGeoTransform[1];
1139 28 : const double dfOvrMaxY = poRawDS->m_adfGeoTransform[3];
1140 : const double dfOvrMinY =
1141 28 : poRawDS->m_adfGeoTransform[3] +
1142 28 : poRawDS->GetRasterYSize() * poRawDS->m_adfGeoTransform[5];
1143 :
1144 28 : if (fabs(dfMinX - dfOvrMinX) < 1e-10 * fabs(dfMinX) &&
1145 26 : fabs(dfMinY - dfOvrMinY) < 1e-10 * fabs(dfMinY) &&
1146 26 : fabs(dfMaxX - dfOvrMaxX) < 1e-10 * fabs(dfMaxX) &&
1147 26 : fabs(dfMaxY - dfOvrMaxY) < 1e-10 * fabs(dfMaxY))
1148 : {
1149 26 : m_apoOverviewDS.emplace_back(std::move(poRawDS));
1150 : }
1151 : else
1152 : {
1153 : // If this zoom level doesn't share the same origin and extent
1154 : // as the most resoluted one, then subset it
1155 2 : CPLStringList aosOptions;
1156 2 : aosOptions.AddString("-of");
1157 2 : aosOptions.AddString("VRT");
1158 2 : aosOptions.AddString("-projwin");
1159 2 : aosOptions.AddString(CPLSPrintf("%.18g", dfMinX));
1160 2 : aosOptions.AddString(CPLSPrintf("%.18g", dfMaxY));
1161 2 : aosOptions.AddString(CPLSPrintf("%.18g", dfMaxX));
1162 2 : aosOptions.AddString(CPLSPrintf("%.18g", dfMinY));
1163 : auto psOptions =
1164 2 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
1165 : auto hDS =
1166 2 : GDALTranslate("", GDALDataset::ToHandle(poRawDS.get()),
1167 : psOptions, nullptr);
1168 2 : GDALTranslateOptionsFree(psOptions);
1169 2 : if (hDS == nullptr)
1170 0 : continue;
1171 2 : m_apoIntermediaryDS.emplace_back(std::move(poRawDS));
1172 2 : m_apoOverviewDS.emplace_back(GDALDataset::FromHandle(hDS));
1173 : }
1174 : }
1175 : }
1176 14 : if (m_poDS == nullptr)
1177 : {
1178 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find valid tile matrix");
1179 0 : return false;
1180 : }
1181 :
1182 : // Create main bands
1183 71 : for (int i = 0; i < m_poDS->GetRasterCount(); i++)
1184 : {
1185 57 : auto poSrcBand = m_poDS->GetRasterBand(i + 1);
1186 57 : auto poBand = new STACTARasterBand(this, i + 1, poSrcBand);
1187 57 : if (oEoBands.IsValid() && oEoBands.Size() == nExpectedBandCount)
1188 : {
1189 : // Set band metadata
1190 42 : if (oEoBands[i].GetType() == CPLJSONObject::Type::Object)
1191 : {
1192 118 : for (const auto &oItem : oEoBands[i].GetChildren())
1193 : {
1194 76 : if (oBands.IsValid())
1195 : {
1196 : // STAC 1.1
1197 22 : if (STARTS_WITH(oItem.GetName().c_str(), "eo:"))
1198 : {
1199 2 : poBand->GDALRasterBand::SetMetadataItem(
1200 2 : oItem.GetName().c_str() + strlen("eo:"),
1201 2 : oItem.ToString().c_str());
1202 : }
1203 36 : else if (oItem.GetName() != "data_type" &&
1204 46 : oItem.GetName() != "nodata" &&
1205 40 : oItem.GetName() != "unit" &&
1206 38 : oItem.GetName() != "raster:scale" &&
1207 65 : oItem.GetName() != "raster:offset" &&
1208 28 : oItem.GetName() != "raster:bits_per_sample")
1209 : {
1210 12 : poBand->GDALRasterBand::SetMetadataItem(
1211 12 : oItem.GetName().c_str(),
1212 12 : oItem.ToString().c_str());
1213 : }
1214 : }
1215 : else
1216 : {
1217 : // STAC 1.0
1218 108 : poBand->GDALRasterBand::SetMetadataItem(
1219 162 : oItem.GetName().c_str(), oItem.ToString().c_str());
1220 : }
1221 : }
1222 : }
1223 : }
1224 201 : if (oRasterBands.IsValid() &&
1225 87 : oRasterBands.Size() == nExpectedBandCount &&
1226 87 : oRasterBands[i].GetType() == CPLJSONObject::Type::Object)
1227 : {
1228 30 : poBand->m_osUnit = oRasterBands[i].GetString("unit");
1229 60 : const double dfScale = oRasterBands[i].GetDouble(
1230 30 : oBands.IsValid() ? "raster:scale" : "scale");
1231 30 : if (dfScale != 0)
1232 5 : poBand->m_dfScale = dfScale;
1233 60 : poBand->m_dfOffset = oRasterBands[i].GetDouble(
1234 30 : oBands.IsValid() ? "raster:offset" : "offset");
1235 60 : const int nBitsPerSample = oRasterBands[i].GetInteger(
1236 30 : oBands.IsValid() ? "raster:bits_per_sample"
1237 : : "bits_per_sample");
1238 5 : if (((nBitsPerSample >= 1 && nBitsPerSample <= 7) &&
1239 35 : poBand->GetRasterDataType() == GDT_Byte) ||
1240 0 : ((nBitsPerSample >= 9 && nBitsPerSample <= 15) &&
1241 0 : poBand->GetRasterDataType() == GDT_UInt16))
1242 : {
1243 5 : poBand->GDALRasterBand::SetMetadataItem(
1244 : "NBITS", CPLSPrintf("%d", nBitsPerSample),
1245 : "IMAGE_STRUCTURE");
1246 : }
1247 : }
1248 57 : SetBand(i + 1, poBand);
1249 : }
1250 :
1251 : // Set dataset metadata
1252 95 : for (const auto &oItem : oProperties.GetChildren())
1253 : {
1254 162 : const auto osName = oItem.GetName();
1255 151 : if (osName != "tiles:tile_matrix_links" &&
1256 70 : osName != "tiles:tile_matrix_sets")
1257 : {
1258 56 : GDALDataset::SetMetadataItem(osName.c_str(),
1259 112 : oItem.ToString().c_str());
1260 : }
1261 : }
1262 :
1263 14 : if (poProtoDS)
1264 : {
1265 : const char *pszInterleave =
1266 9 : poProtoDS->GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE");
1267 9 : GDALDataset::SetMetadataItem("INTERLEAVE",
1268 : pszInterleave ? pszInterleave : "PIXEL",
1269 : "IMAGE_STRUCTURE");
1270 : }
1271 : else
1272 : {
1273 : // A bit bold to assume that, but that should be a reasonable
1274 : // setting
1275 5 : GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
1276 : }
1277 :
1278 14 : m_bDownloadWholeMetaTile = CPLTestBool(CSLFetchNameValueDef(
1279 14 : poOpenInfo->papszOpenOptions, "WHOLE_METATILE", "NO"));
1280 :
1281 14 : return true;
1282 : }
1283 :
1284 : /************************************************************************/
1285 : /* ~STACTADataset() */
1286 : /************************************************************************/
1287 :
1288 42 : STACTADataset::~STACTADataset()
1289 : {
1290 21 : m_poDS.reset();
1291 21 : m_apoOverviewDS.clear();
1292 21 : m_apoIntermediaryDS.clear();
1293 42 : }
1294 :
1295 : /************************************************************************/
1296 : /* FlushCache() */
1297 : /************************************************************************/
1298 :
1299 0 : CPLErr STACTADataset::FlushCache(bool bAtClosing)
1300 : {
1301 0 : m_oCacheTileDS.clear();
1302 0 : return GDALDataset::FlushCache(bAtClosing);
1303 : }
1304 :
1305 : /************************************************************************/
1306 : /* InitRaster() */
1307 : /************************************************************************/
1308 :
1309 42 : bool STACTARawDataset::InitRaster(GDALDataset *poProtoDS,
1310 : const std::vector<GDALDataType> &aeDT,
1311 : const std::vector<bool> &abSetNoData,
1312 : const std::vector<double> &adfNoData,
1313 : const gdal::TileMatrixSet *poTMS,
1314 : const std::string &osTMId,
1315 : const gdal::TileMatrixSet::TileMatrix &oTM,
1316 : const std::map<CPLString, Limits> &oMapLimits)
1317 : {
1318 42 : int nMatrixWidth = oTM.mMatrixWidth;
1319 42 : int nMatrixHeight = oTM.mMatrixHeight;
1320 42 : auto oIterLimit = oMapLimits.find(osTMId);
1321 42 : if (oIterLimit != oMapLimits.end())
1322 : {
1323 33 : m_nMinMetaTileCol = oIterLimit->second.min_tile_col;
1324 33 : m_nMinMetaTileRow = oIterLimit->second.min_tile_row;
1325 33 : nMatrixWidth = oIterLimit->second.max_tile_col - m_nMinMetaTileCol + 1;
1326 33 : nMatrixHeight = oIterLimit->second.max_tile_row - m_nMinMetaTileRow + 1;
1327 : }
1328 42 : m_nMetaTileWidth = oTM.mTileWidth;
1329 42 : m_nMetaTileHeight = oTM.mTileHeight;
1330 42 : nRasterXSize = nMatrixWidth * m_nMetaTileWidth;
1331 42 : nRasterYSize = nMatrixHeight * m_nMetaTileHeight;
1332 :
1333 42 : if (poProtoDS)
1334 : {
1335 108 : for (int i = 0; i < poProtoDS->GetRasterCount(); i++)
1336 : {
1337 81 : auto poProtoBand = poProtoDS->GetRasterBand(i + 1);
1338 81 : auto poBand = new STACTARawRasterBand(this, i + 1, poProtoBand);
1339 81 : SetBand(i + 1, poBand);
1340 : }
1341 : }
1342 : else
1343 : {
1344 105 : for (int i = 0; i < static_cast<int>(aeDT.size()); i++)
1345 : {
1346 90 : auto poBand = new STACTARawRasterBand(this, i + 1, aeDT[i],
1347 90 : abSetNoData[i], adfNoData[i]);
1348 90 : SetBand(i + 1, poBand);
1349 : }
1350 : }
1351 :
1352 84 : CPLString osCRS = poTMS->crs().c_str();
1353 42 : if (osCRS == "http://www.opengis.net/def/crs/OGC/1.3/CRS84")
1354 42 : osCRS = "EPSG:4326";
1355 42 : if (m_oSRS.SetFromUserInput(osCRS) != OGRERR_NONE)
1356 : {
1357 0 : return false;
1358 : }
1359 42 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1360 84 : m_adfGeoTransform[0] =
1361 42 : oTM.mTopLeftX + m_nMinMetaTileCol * m_nMetaTileWidth * oTM.mResX;
1362 42 : m_adfGeoTransform[1] = oTM.mResX;
1363 84 : m_adfGeoTransform[3] =
1364 42 : oTM.mTopLeftY - m_nMinMetaTileRow * m_nMetaTileHeight * oTM.mResY;
1365 42 : m_adfGeoTransform[5] = -oTM.mResY;
1366 42 : SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
1367 :
1368 42 : return true;
1369 : }
1370 :
1371 : /************************************************************************/
1372 : /* GetSpatialRef () */
1373 : /************************************************************************/
1374 :
1375 3 : const OGRSpatialReference *STACTADataset::GetSpatialRef() const
1376 : {
1377 3 : return nBands == 0 ? nullptr : &m_oSRS;
1378 : }
1379 :
1380 : /************************************************************************/
1381 : /* GetGeoTransform() */
1382 : /************************************************************************/
1383 :
1384 3 : CPLErr STACTADataset::GetGeoTransform(double *padfGeoTransform)
1385 : {
1386 3 : memcpy(padfGeoTransform, &m_adfGeoTransform[0], 6 * sizeof(double));
1387 3 : return nBands == 0 ? CE_Failure : CE_None;
1388 : }
1389 :
1390 : /************************************************************************/
1391 : /* OpenStatic() */
1392 : /************************************************************************/
1393 :
1394 21 : GDALDataset *STACTADataset::OpenStatic(GDALOpenInfo *poOpenInfo)
1395 : {
1396 21 : if (!Identify(poOpenInfo))
1397 0 : return nullptr;
1398 42 : auto poDS = std::make_unique<STACTADataset>();
1399 21 : if (!poDS->Open(poOpenInfo))
1400 5 : return nullptr;
1401 16 : return poDS.release();
1402 : }
1403 :
1404 : /************************************************************************/
1405 : /* GDALRegister_STACTA() */
1406 : /************************************************************************/
1407 :
1408 1522 : void GDALRegister_STACTA()
1409 :
1410 : {
1411 1522 : if (GDALGetDriverByName("STACTA") != nullptr)
1412 301 : return;
1413 :
1414 1221 : GDALDriver *poDriver = new GDALDriver();
1415 :
1416 1221 : poDriver->SetDescription("STACTA");
1417 1221 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1418 1221 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1419 1221 : "Spatio-Temporal Asset Catalog Tiled Assets");
1420 1221 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/stacta.html");
1421 1221 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "json");
1422 1221 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1423 1221 : poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");
1424 1221 : poDriver->SetMetadataItem(
1425 : GDAL_DMD_OPENOPTIONLIST,
1426 : "<OpenOptionList>"
1427 : " <Option name='WHOLE_METATILE' type='boolean' "
1428 : "description='Whether to download whole metatiles'/>"
1429 : " <Option name='SKIP_MISSING_METATILE' type='boolean' "
1430 : "description='Whether to gracefully skip missing metatiles'/>"
1431 1221 : "</OpenOptionList>");
1432 :
1433 1221 : poDriver->pfnOpen = STACTADataset::OpenStatic;
1434 1221 : poDriver->pfnIdentify = STACTADataset::Identify;
1435 :
1436 1221 : GetGDALDriverManager()->RegisterDriver(poDriver);
1437 : }
|
