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