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 50501 : int STACTADataset::Identify(GDALOpenInfo *poOpenInfo)
636 : {
637 50501 : if (STARTS_WITH(poOpenInfo->pszFilename, "STACTA:"))
638 : {
639 6 : return true;
640 : }
641 :
642 50495 : const bool bIsSingleDriver = poOpenInfo->IsSingleAllowedDriver("STACTA");
643 50481 : if (bIsSingleDriver && (STARTS_WITH(poOpenInfo->pszFilename, "http://") ||
644 4 : STARTS_WITH(poOpenInfo->pszFilename, "https://")))
645 : {
646 1 : return true;
647 : }
648 :
649 50491 : if (
650 : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
651 100967 : (!bIsSingleDriver &&
652 50913 : !EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "json")) ||
653 : #endif
654 433 : poOpenInfo->nHeaderBytes == 0)
655 : {
656 50339 : return false;
657 : }
658 :
659 357 : for (int i = 0; i < 2; i++)
660 : {
661 : // TryToIngest() may reallocate pabyHeader, so do not move this
662 : // before the loop.
663 252 : const char *pszHeader =
664 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
665 252 : while (*pszHeader != 0 &&
666 252 : std::isspace(static_cast<unsigned char>(*pszHeader)))
667 0 : ++pszHeader;
668 252 : if (bIsSingleDriver)
669 : {
670 4 : return pszHeader[0] == '{';
671 : }
672 :
673 248 : if (strstr(pszHeader, "\"stac_extensions\"") != nullptr &&
674 64 : (strstr(pszHeader, "\"tiled-assets\"") != nullptr ||
675 46 : strstr(pszHeader,
676 : "https://stac-extensions.github.io/tiled-assets/") !=
677 : nullptr))
678 : {
679 38 : return true;
680 : }
681 :
682 210 : if (i == 0)
683 : {
684 : // Should be enough for a STACTA .json file
685 105 : poOpenInfo->TryToIngest(32768);
686 : }
687 : }
688 :
689 105 : return false;
690 : }
691 :
692 : /************************************************************************/
693 : /* Open() */
694 : /************************************************************************/
695 :
696 23 : bool STACTADataset::Open(GDALOpenInfo *poOpenInfo)
697 : {
698 46 : CPLString osFilename(poOpenInfo->pszFilename);
699 46 : CPLString osAssetName;
700 46 : CPLString osTMS;
701 23 : if (STARTS_WITH(poOpenInfo->pszFilename, "STACTA:"))
702 : {
703 : const CPLStringList aosTokens(CSLTokenizeString2(
704 3 : poOpenInfo->pszFilename, ":", CSLT_HONOURSTRINGS));
705 4 : if (aosTokens.size() != 2 && aosTokens.size() != 3 &&
706 1 : aosTokens.size() != 4)
707 0 : return false;
708 3 : osFilename = aosTokens[1];
709 3 : if (aosTokens.size() >= 3)
710 3 : osAssetName = aosTokens[2];
711 3 : if (aosTokens.size() == 4)
712 1 : osTMS = aosTokens[3];
713 : }
714 :
715 46 : CPLJSONDocument oDoc;
716 46 : if (STARTS_WITH(osFilename, "http://") ||
717 23 : STARTS_WITH(osFilename, "https://"))
718 : {
719 0 : if (!oDoc.LoadUrl(osFilename, nullptr))
720 0 : return false;
721 : }
722 : else
723 : {
724 23 : if (!oDoc.Load(osFilename))
725 0 : return false;
726 : }
727 46 : const auto oRoot = oDoc.GetRoot();
728 69 : const auto oProperties = oRoot["properties"];
729 46 : if (!oProperties.IsValid() ||
730 23 : oProperties.GetType() != CPLJSONObject::Type::Object)
731 : {
732 0 : CPLError(CE_Failure, CPLE_AppDefined, "Missing properties");
733 0 : return false;
734 : }
735 :
736 69 : const auto oAssetTemplates = oRoot["asset_templates"];
737 46 : if (!oAssetTemplates.IsValid() ||
738 23 : oAssetTemplates.GetType() != CPLJSONObject::Type::Object)
739 : {
740 0 : CPLError(CE_Failure, CPLE_AppDefined, "Missing asset_templates");
741 0 : return false;
742 : }
743 :
744 46 : const auto aoAssetTemplates = oAssetTemplates.GetChildren();
745 23 : if (aoAssetTemplates.size() == 0)
746 : {
747 0 : CPLError(CE_Failure, CPLE_AppDefined, "Empty asset_templates");
748 0 : return false;
749 : }
750 :
751 69 : const auto oTMSs = oProperties.GetObj("tiles:tile_matrix_sets");
752 23 : if (!oTMSs.IsValid() || oTMSs.GetType() != CPLJSONObject::Type::Object)
753 : {
754 0 : CPLError(CE_Failure, CPLE_AppDefined,
755 : "Missing properties[\"tiles:tile_matrix_sets\"]");
756 0 : return false;
757 : }
758 46 : const auto aoTMSs = oTMSs.GetChildren();
759 23 : if (aoTMSs.empty())
760 : {
761 0 : CPLError(CE_Failure, CPLE_AppDefined,
762 : "Empty properties[\"tiles:tile_matrix_sets\"]");
763 0 : return false;
764 : }
765 :
766 48 : if ((aoAssetTemplates.size() >= 2 || aoTMSs.size() >= 2) &&
767 48 : osAssetName.empty() && osTMS.empty())
768 : {
769 2 : int nSDSCount = 0;
770 5 : for (const auto &oAssetTemplate : aoAssetTemplates)
771 : {
772 6 : const CPLString osAssetNameSubDS = oAssetTemplate.GetName();
773 3 : const char *pszAssetNameSubDS = osAssetNameSubDS.c_str();
774 3 : if (aoTMSs.size() >= 2)
775 : {
776 3 : for (const auto &oTMS : aoTMSs)
777 : {
778 2 : const CPLString osTMSSubDS = oTMS.GetName();
779 2 : const char *pszTMSSubDS = osTMSSubDS.c_str();
780 2 : GDALDataset::SetMetadataItem(
781 : CPLSPrintf("SUBDATASET_%d_NAME", nSDSCount + 1),
782 : CPLSPrintf("STACTA:\"%s\":%s:%s", osFilename.c_str(),
783 : pszAssetNameSubDS, pszTMSSubDS),
784 : "SUBDATASETS");
785 2 : GDALDataset::SetMetadataItem(
786 : CPLSPrintf("SUBDATASET_%d_DESC", nSDSCount + 1),
787 : CPLSPrintf("Asset %s, tile matrix set %s",
788 : pszAssetNameSubDS, pszTMSSubDS),
789 : "SUBDATASETS");
790 2 : nSDSCount++;
791 : }
792 : }
793 : else
794 : {
795 2 : GDALDataset::SetMetadataItem(
796 : CPLSPrintf("SUBDATASET_%d_NAME", nSDSCount + 1),
797 : CPLSPrintf("STACTA:\"%s\":%s", osFilename.c_str(),
798 : pszAssetNameSubDS),
799 : "SUBDATASETS");
800 2 : GDALDataset::SetMetadataItem(
801 : CPLSPrintf("SUBDATASET_%d_DESC", nSDSCount + 1),
802 : CPLSPrintf("Asset %s", pszAssetNameSubDS), "SUBDATASETS");
803 2 : nSDSCount++;
804 : }
805 : }
806 2 : return true;
807 : }
808 :
809 21 : if (osAssetName.empty())
810 : {
811 18 : osAssetName = aoAssetTemplates[0].GetName();
812 : }
813 42 : const auto oAssetTemplate = oAssetTemplates.GetObj(osAssetName);
814 42 : if (!oAssetTemplate.IsValid() ||
815 21 : oAssetTemplate.GetType() != CPLJSONObject::Type::Object)
816 : {
817 0 : CPLError(CE_Failure, CPLE_AppDefined,
818 : "Cannot find asset_templates[\"%s\"]", osAssetName.c_str());
819 0 : return false;
820 : }
821 :
822 21 : if (osTMS.empty())
823 : {
824 20 : osTMS = aoTMSs[0].GetName();
825 : }
826 42 : const auto oTMS = oTMSs.GetObj(osTMS);
827 21 : if (!oTMS.IsValid() || oTMS.GetType() != CPLJSONObject::Type::Object)
828 : {
829 0 : CPLError(CE_Failure, CPLE_AppDefined,
830 : "Cannot find properties[\"tiles:tile_matrix_sets\"][\"%s\"]",
831 : osTMS.c_str());
832 0 : return false;
833 : }
834 :
835 : auto poTMS = gdal::TileMatrixSet::parse(
836 42 : oTMS.Format(CPLJSONObject::PrettyFormat::Plain).c_str());
837 21 : if (poTMS == nullptr)
838 0 : return false;
839 :
840 63 : CPLString osURLTemplate = oAssetTemplate.GetString("href");
841 21 : if (osURLTemplate.empty())
842 : {
843 0 : CPLError(CE_Failure, CPLE_AppDefined,
844 : "Cannot find asset_templates[\"%s\"][\"href\"]",
845 : osAssetName.c_str());
846 : }
847 21 : osURLTemplate.replaceAll("{TileMatrixSet}", osTMS);
848 21 : if (STARTS_WITH(osURLTemplate, "file://"))
849 : {
850 0 : osURLTemplate = osURLTemplate.substr(strlen("file://"));
851 : }
852 21 : else if (STARTS_WITH(osURLTemplate, "s3://"))
853 : {
854 0 : osURLTemplate = "/vsis3/" + osURLTemplate.substr(strlen("s3://"));
855 : }
856 :
857 41 : if (!STARTS_WITH(osURLTemplate, "http://") &&
858 20 : !STARTS_WITH(osURLTemplate, "https://"))
859 : {
860 20 : if (STARTS_WITH(osURLTemplate, "./"))
861 20 : osURLTemplate = osURLTemplate.substr(2);
862 : osURLTemplate = CPLProjectRelativeFilename(CPLGetDirname(osFilename),
863 20 : osURLTemplate);
864 : }
865 :
866 : // Parse optional tile matrix set limits
867 42 : std::map<CPLString, Limits> oMapLimits;
868 63 : const auto oTMLinks = oProperties.GetObj("tiles:tile_matrix_links");
869 21 : if (oTMLinks.IsValid())
870 : {
871 18 : if (oTMLinks.GetType() != CPLJSONObject::Type::Object)
872 : {
873 0 : CPLError(
874 : CE_Failure, CPLE_AppDefined,
875 : "Invalid type for properties[\"tiles:tile_matrix_links\"]");
876 0 : return false;
877 : }
878 :
879 54 : auto oLimits = oTMLinks[osTMS]["limits"];
880 36 : if (oLimits.IsValid() &&
881 18 : oLimits.GetType() == CPLJSONObject::Type::Object)
882 : {
883 72 : for (const auto &oLimit : oLimits.GetChildren())
884 : {
885 54 : Limits limits;
886 54 : limits.min_tile_col = oLimit.GetInteger("min_tile_col");
887 54 : limits.max_tile_col = oLimit.GetInteger("max_tile_col");
888 54 : limits.min_tile_row = oLimit.GetInteger("min_tile_row");
889 54 : limits.max_tile_row = oLimit.GetInteger("max_tile_row");
890 54 : oMapLimits[oLimit.GetName()] = limits;
891 : }
892 : }
893 : }
894 21 : const auto &tmsList = poTMS->tileMatrixList();
895 21 : if (tmsList.empty())
896 0 : return false;
897 :
898 42 : m_bSkipMissingMetaTile = CPLTestBool(CSLFetchNameValueDef(
899 21 : poOpenInfo->papszOpenOptions, "SKIP_MISSING_METATILE",
900 : CPLGetConfigOption("GDAL_STACTA_SKIP_MISSING_METATILE", "NO")));
901 :
902 : // STAC 1.1 uses bands instead of eo:bands and raster:bands
903 63 : const auto oBands = oAssetTemplate.GetArray("bands");
904 :
905 : // Check if there are both eo:bands and raster:bands extension
906 : // If so, we don't need to fetch a prototype metatile to derive the
907 : // information we need (number of bands, data type and nodata value)
908 : const auto oEoBands =
909 62 : oBands.IsValid() ? oBands : oAssetTemplate.GetArray("eo:bands");
910 : const auto oRasterBands =
911 62 : oBands.IsValid() ? oBands : oAssetTemplate.GetArray("raster:bands");
912 :
913 42 : std::vector<GDALDataType> aeDT;
914 42 : std::vector<double> adfNoData;
915 42 : std::vector<bool> abSetNoData;
916 21 : int nExpectedBandCount = 0;
917 21 : if (oRasterBands.IsValid())
918 : {
919 9 : if (oEoBands.IsValid() && oEoBands.Size() != oRasterBands.Size())
920 : {
921 1 : CPLError(CE_Warning, CPLE_AppDefined,
922 : "Number of bands in eo:bands and raster:bands is not "
923 : "identical. Ignoring the later");
924 : }
925 : else
926 : {
927 8 : nExpectedBandCount = oRasterBands.Size();
928 :
929 : const struct
930 : {
931 : const char *pszStacDataType;
932 : GDALDataType eGDALDataType;
933 8 : } aDataTypeMapping[] = {
934 : {"int8", GDT_Int8},
935 : {"int16", GDT_Int16},
936 : {"int32", GDT_Int32},
937 : {"int64", GDT_Int64},
938 : {"uint8", GDT_Byte},
939 : {"uint16", GDT_UInt16},
940 : {"uint32", GDT_UInt32},
941 : {"uint64", GDT_UInt64},
942 : // float16: 16-bit float; unhandled
943 : {"float32", GDT_Float32},
944 : {"float64", GDT_Float64},
945 : {"cint16", GDT_CInt16},
946 : {"cint32", GDT_CInt32},
947 : {"cfloat32", GDT_CFloat32},
948 : {"cfloat64", GDT_CFloat64},
949 : };
950 :
951 38 : for (int i = 0; i < nExpectedBandCount; ++i)
952 : {
953 33 : if (oRasterBands[i].GetType() != CPLJSONObject::Type::Object)
954 : {
955 1 : CPLError(CE_Failure, CPLE_AppDefined,
956 : "Wrong raster:bands[%d]", i);
957 3 : return false;
958 : }
959 : const std::string osDataType =
960 64 : oRasterBands[i].GetString("data_type");
961 32 : GDALDataType eDT = GDT_Unknown;
962 280 : for (const auto &oTuple : aDataTypeMapping)
963 : {
964 278 : if (osDataType == oTuple.pszStacDataType)
965 : {
966 30 : eDT = oTuple.eGDALDataType;
967 30 : break;
968 : }
969 : }
970 32 : if (eDT == GDT_Unknown)
971 : {
972 2 : CPLError(CE_Failure, CPLE_AppDefined,
973 : "Wrong raster:bands[%d].data_type = %s", i,
974 : osDataType.c_str());
975 2 : return false;
976 : }
977 30 : aeDT.push_back(eDT);
978 :
979 90 : const auto oNoData = oRasterBands[i].GetObj("nodata");
980 30 : if (oNoData.GetType() == CPLJSONObject::Type::String)
981 : {
982 48 : const std::string osNoData = oNoData.ToString();
983 16 : if (osNoData == "inf")
984 : {
985 5 : abSetNoData.push_back(true);
986 5 : adfNoData.push_back(
987 5 : std::numeric_limits<double>::infinity());
988 : }
989 11 : else if (osNoData == "-inf")
990 : {
991 5 : abSetNoData.push_back(true);
992 5 : adfNoData.push_back(
993 5 : -std::numeric_limits<double>::infinity());
994 : }
995 6 : else if (osNoData == "nan")
996 : {
997 5 : abSetNoData.push_back(true);
998 5 : adfNoData.push_back(
999 5 : std::numeric_limits<double>::quiet_NaN());
1000 : }
1001 : else
1002 : {
1003 1 : CPLError(CE_Warning, CPLE_AppDefined,
1004 : "Invalid raster:bands[%d].nodata = %s", i,
1005 : osNoData.c_str());
1006 1 : abSetNoData.push_back(false);
1007 1 : adfNoData.push_back(
1008 1 : std::numeric_limits<double>::quiet_NaN());
1009 : }
1010 : }
1011 14 : else if (oNoData.GetType() == CPLJSONObject::Type::Integer ||
1012 25 : oNoData.GetType() == CPLJSONObject::Type::Long ||
1013 11 : oNoData.GetType() == CPLJSONObject::Type::Double)
1014 : {
1015 8 : abSetNoData.push_back(true);
1016 8 : adfNoData.push_back(oNoData.ToDouble());
1017 : }
1018 6 : else if (!oNoData.IsValid())
1019 : {
1020 5 : abSetNoData.push_back(false);
1021 5 : adfNoData.push_back(
1022 5 : std::numeric_limits<double>::quiet_NaN());
1023 : }
1024 : else
1025 : {
1026 1 : CPLError(CE_Warning, CPLE_AppDefined,
1027 : "Invalid raster:bands[%d].nodata", i);
1028 1 : abSetNoData.push_back(false);
1029 1 : adfNoData.push_back(
1030 1 : std::numeric_limits<double>::quiet_NaN());
1031 : }
1032 : }
1033 :
1034 5 : CPLAssert(aeDT.size() == abSetNoData.size());
1035 5 : CPLAssert(adfNoData.size() == abSetNoData.size());
1036 : }
1037 : }
1038 :
1039 18 : std::unique_ptr<GDALDataset> poProtoDS;
1040 18 : if (aeDT.empty())
1041 : {
1042 13 : for (int i = 0; i < static_cast<int>(tmsList.size()); i++)
1043 : {
1044 : // Open a metatile to get mostly its band data type
1045 13 : int nProtoTileCol = 0;
1046 13 : int nProtoTileRow = 0;
1047 13 : auto oIterLimit = oMapLimits.find(tmsList[i].mId);
1048 13 : if (oIterLimit != oMapLimits.end())
1049 : {
1050 10 : nProtoTileCol = oIterLimit->second.min_tile_col;
1051 10 : nProtoTileRow = oIterLimit->second.min_tile_row;
1052 : }
1053 : const CPLString osURL =
1054 13 : CPLString(osURLTemplate)
1055 26 : .replaceAll("{TileMatrix}", tmsList[i].mId)
1056 26 : .replaceAll("{TileRow}", CPLSPrintf("%d", nProtoTileRow))
1057 26 : .replaceAll("{TileCol}", CPLSPrintf("%d", nProtoTileCol));
1058 38 : CPLString osProtoDSName = (STARTS_WITH(osURL, "http://") ||
1059 12 : STARTS_WITH(osURL, "https://"))
1060 26 : ? CPLString("/vsicurl/" + osURL)
1061 26 : : osURL;
1062 : CPLConfigOptionSetter oSetter("GDAL_DISABLE_READDIR_ON_OPEN",
1063 : "EMPTY_DIR",
1064 13 : /* bSetOnlyIfUndefined = */ true);
1065 13 : if (m_bSkipMissingMetaTile)
1066 2 : CPLPushErrorHandler(CPLQuietErrorHandler);
1067 13 : poProtoDS.reset(GDALDataset::Open(osProtoDSName.c_str()));
1068 13 : if (m_bSkipMissingMetaTile)
1069 2 : CPLPopErrorHandler();
1070 13 : if (poProtoDS != nullptr)
1071 : {
1072 11 : break;
1073 : }
1074 2 : if (!m_bSkipMissingMetaTile)
1075 : {
1076 2 : CPLError(CE_Failure, CPLE_OpenFailed, "Cannot open %s",
1077 : osURL.c_str());
1078 2 : return false;
1079 : }
1080 : }
1081 11 : if (poProtoDS == nullptr)
1082 : {
1083 0 : CPLError(CE_Failure, CPLE_AppDefined,
1084 : "Cannot find prototype dataset");
1085 0 : return false;
1086 : }
1087 11 : nExpectedBandCount = poProtoDS->GetRasterCount();
1088 : }
1089 :
1090 : // Iterate over tile matrices to create corresponding STACTARawDataset
1091 : // objects
1092 64 : for (int i = static_cast<int>(tmsList.size() - 1); i >= 0; i--)
1093 : {
1094 48 : const auto &oTM = tmsList[i];
1095 48 : int nMatrixWidth = oTM.mMatrixWidth;
1096 48 : int nMatrixHeight = oTM.mMatrixHeight;
1097 48 : auto oIterLimit = oMapLimits.find(tmsList[i].mId);
1098 48 : if (oIterLimit != oMapLimits.end())
1099 : {
1100 39 : nMatrixWidth = oIterLimit->second.max_tile_col -
1101 39 : oIterLimit->second.min_tile_col + 1;
1102 39 : nMatrixHeight = oIterLimit->second.max_tile_row -
1103 39 : oIterLimit->second.min_tile_row + 1;
1104 : }
1105 48 : if (nMatrixWidth <= 0 || oTM.mTileWidth > INT_MAX / nMatrixWidth ||
1106 48 : nMatrixHeight <= 0 || oTM.mTileHeight > INT_MAX / nMatrixHeight)
1107 : {
1108 0 : continue;
1109 : }
1110 48 : auto poRawDS = std::make_unique<STACTARawDataset>();
1111 96 : if (!poRawDS->InitRaster(poProtoDS.get(), aeDT, abSetNoData, adfNoData,
1112 48 : poTMS.get(), tmsList[i].mId, oTM, oMapLimits))
1113 : {
1114 0 : return false;
1115 : }
1116 48 : poRawDS->m_osURLTemplate = osURLTemplate;
1117 48 : poRawDS->m_osURLTemplate.replaceAll("{TileMatrix}", tmsList[i].mId);
1118 48 : poRawDS->m_poMasterDS = this;
1119 :
1120 48 : if (m_poDS == nullptr)
1121 : {
1122 16 : nRasterXSize = poRawDS->GetRasterXSize();
1123 16 : nRasterYSize = poRawDS->GetRasterYSize();
1124 16 : m_oSRS = poRawDS->m_oSRS;
1125 16 : memcpy(&m_adfGeoTransform[0], &poRawDS->m_adfGeoTransform[0],
1126 : 6 * sizeof(double));
1127 16 : m_poDS = std::move(poRawDS);
1128 : }
1129 : else
1130 : {
1131 32 : const double dfMinX = m_adfGeoTransform[0];
1132 : const double dfMaxX =
1133 32 : m_adfGeoTransform[0] + GetRasterXSize() * m_adfGeoTransform[1];
1134 32 : const double dfMaxY = m_adfGeoTransform[3];
1135 : const double dfMinY =
1136 32 : m_adfGeoTransform[3] + GetRasterYSize() * m_adfGeoTransform[5];
1137 :
1138 32 : const double dfOvrMinX = poRawDS->m_adfGeoTransform[0];
1139 : const double dfOvrMaxX =
1140 32 : poRawDS->m_adfGeoTransform[0] +
1141 32 : poRawDS->GetRasterXSize() * poRawDS->m_adfGeoTransform[1];
1142 32 : const double dfOvrMaxY = poRawDS->m_adfGeoTransform[3];
1143 : const double dfOvrMinY =
1144 32 : poRawDS->m_adfGeoTransform[3] +
1145 32 : poRawDS->GetRasterYSize() * poRawDS->m_adfGeoTransform[5];
1146 :
1147 32 : if (fabs(dfMinX - dfOvrMinX) < 1e-10 * fabs(dfMinX) &&
1148 30 : fabs(dfMinY - dfOvrMinY) < 1e-10 * fabs(dfMinY) &&
1149 30 : fabs(dfMaxX - dfOvrMaxX) < 1e-10 * fabs(dfMaxX) &&
1150 30 : fabs(dfMaxY - dfOvrMaxY) < 1e-10 * fabs(dfMaxY))
1151 : {
1152 30 : m_apoOverviewDS.emplace_back(std::move(poRawDS));
1153 : }
1154 : else
1155 : {
1156 : // If this zoom level doesn't share the same origin and extent
1157 : // as the most resoluted one, then subset it
1158 2 : CPLStringList aosOptions;
1159 2 : aosOptions.AddString("-of");
1160 2 : aosOptions.AddString("VRT");
1161 2 : aosOptions.AddString("-projwin");
1162 2 : aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
1163 2 : aosOptions.AddString(CPLSPrintf("%.17g", dfMaxY));
1164 2 : aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
1165 2 : aosOptions.AddString(CPLSPrintf("%.17g", dfMinY));
1166 : auto psOptions =
1167 2 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
1168 : auto hDS =
1169 2 : GDALTranslate("", GDALDataset::ToHandle(poRawDS.get()),
1170 : psOptions, nullptr);
1171 2 : GDALTranslateOptionsFree(psOptions);
1172 2 : if (hDS == nullptr)
1173 0 : continue;
1174 2 : m_apoIntermediaryDS.emplace_back(std::move(poRawDS));
1175 2 : m_apoOverviewDS.emplace_back(GDALDataset::FromHandle(hDS));
1176 : }
1177 : }
1178 : }
1179 16 : if (m_poDS == nullptr)
1180 : {
1181 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find valid tile matrix");
1182 0 : return false;
1183 : }
1184 :
1185 : // Create main bands
1186 79 : for (int i = 0; i < m_poDS->GetRasterCount(); i++)
1187 : {
1188 63 : auto poSrcBand = m_poDS->GetRasterBand(i + 1);
1189 63 : auto poBand = new STACTARasterBand(this, i + 1, poSrcBand);
1190 63 : if (oEoBands.IsValid() && oEoBands.Size() == nExpectedBandCount)
1191 : {
1192 : // Set band metadata
1193 48 : if (oEoBands[i].GetType() == CPLJSONObject::Type::Object)
1194 : {
1195 136 : for (const auto &oItem : oEoBands[i].GetChildren())
1196 : {
1197 88 : if (oBands.IsValid())
1198 : {
1199 : // STAC 1.1
1200 22 : if (STARTS_WITH(oItem.GetName().c_str(), "eo:"))
1201 : {
1202 2 : poBand->GDALRasterBand::SetMetadataItem(
1203 2 : oItem.GetName().c_str() + strlen("eo:"),
1204 2 : oItem.ToString().c_str());
1205 : }
1206 36 : else if (oItem.GetName() != "data_type" &&
1207 46 : oItem.GetName() != "nodata" &&
1208 40 : oItem.GetName() != "unit" &&
1209 38 : oItem.GetName() != "raster:scale" &&
1210 65 : oItem.GetName() != "raster:offset" &&
1211 28 : oItem.GetName() != "raster:bits_per_sample")
1212 : {
1213 12 : poBand->GDALRasterBand::SetMetadataItem(
1214 12 : oItem.GetName().c_str(),
1215 12 : oItem.ToString().c_str());
1216 : }
1217 : }
1218 : else
1219 : {
1220 : // STAC 1.0
1221 132 : poBand->GDALRasterBand::SetMetadataItem(
1222 198 : oItem.GetName().c_str(), oItem.ToString().c_str());
1223 : }
1224 : }
1225 : }
1226 : }
1227 219 : if (oRasterBands.IsValid() &&
1228 93 : oRasterBands.Size() == nExpectedBandCount &&
1229 93 : oRasterBands[i].GetType() == CPLJSONObject::Type::Object)
1230 : {
1231 30 : poBand->m_osUnit = oRasterBands[i].GetString("unit");
1232 60 : const double dfScale = oRasterBands[i].GetDouble(
1233 30 : oBands.IsValid() ? "raster:scale" : "scale");
1234 30 : if (dfScale != 0)
1235 5 : poBand->m_dfScale = dfScale;
1236 60 : poBand->m_dfOffset = oRasterBands[i].GetDouble(
1237 30 : oBands.IsValid() ? "raster:offset" : "offset");
1238 60 : const int nBitsPerSample = oRasterBands[i].GetInteger(
1239 30 : oBands.IsValid() ? "raster:bits_per_sample"
1240 : : "bits_per_sample");
1241 5 : if (((nBitsPerSample >= 1 && nBitsPerSample <= 7) &&
1242 35 : poBand->GetRasterDataType() == GDT_Byte) ||
1243 0 : ((nBitsPerSample >= 9 && nBitsPerSample <= 15) &&
1244 0 : poBand->GetRasterDataType() == GDT_UInt16))
1245 : {
1246 5 : poBand->GDALRasterBand::SetMetadataItem(
1247 : "NBITS", CPLSPrintf("%d", nBitsPerSample),
1248 : "IMAGE_STRUCTURE");
1249 : }
1250 : }
1251 63 : SetBand(i + 1, poBand);
1252 : }
1253 :
1254 : // Set dataset metadata
1255 109 : for (const auto &oItem : oProperties.GetChildren())
1256 : {
1257 186 : const auto osName = oItem.GetName();
1258 173 : if (osName != "tiles:tile_matrix_links" &&
1259 80 : osName != "tiles:tile_matrix_sets")
1260 : {
1261 64 : GDALDataset::SetMetadataItem(osName.c_str(),
1262 128 : oItem.ToString().c_str());
1263 : }
1264 : }
1265 :
1266 16 : if (poProtoDS)
1267 : {
1268 : const char *pszInterleave =
1269 11 : poProtoDS->GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE");
1270 11 : GDALDataset::SetMetadataItem("INTERLEAVE",
1271 : pszInterleave ? pszInterleave : "PIXEL",
1272 : "IMAGE_STRUCTURE");
1273 : }
1274 : else
1275 : {
1276 : // A bit bold to assume that, but that should be a reasonable
1277 : // setting
1278 5 : GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
1279 : }
1280 :
1281 16 : m_bDownloadWholeMetaTile = CPLTestBool(CSLFetchNameValueDef(
1282 16 : poOpenInfo->papszOpenOptions, "WHOLE_METATILE", "NO"));
1283 :
1284 16 : return true;
1285 : }
1286 :
1287 : /************************************************************************/
1288 : /* ~STACTADataset() */
1289 : /************************************************************************/
1290 :
1291 46 : STACTADataset::~STACTADataset()
1292 : {
1293 23 : m_poDS.reset();
1294 23 : m_apoOverviewDS.clear();
1295 23 : m_apoIntermediaryDS.clear();
1296 46 : }
1297 :
1298 : /************************************************************************/
1299 : /* FlushCache() */
1300 : /************************************************************************/
1301 :
1302 0 : CPLErr STACTADataset::FlushCache(bool bAtClosing)
1303 : {
1304 0 : m_oCacheTileDS.clear();
1305 0 : return GDALDataset::FlushCache(bAtClosing);
1306 : }
1307 :
1308 : /************************************************************************/
1309 : /* InitRaster() */
1310 : /************************************************************************/
1311 :
1312 48 : bool STACTARawDataset::InitRaster(GDALDataset *poProtoDS,
1313 : const std::vector<GDALDataType> &aeDT,
1314 : const std::vector<bool> &abSetNoData,
1315 : const std::vector<double> &adfNoData,
1316 : const gdal::TileMatrixSet *poTMS,
1317 : const std::string &osTMId,
1318 : const gdal::TileMatrixSet::TileMatrix &oTM,
1319 : const std::map<CPLString, Limits> &oMapLimits)
1320 : {
1321 48 : int nMatrixWidth = oTM.mMatrixWidth;
1322 48 : int nMatrixHeight = oTM.mMatrixHeight;
1323 48 : auto oIterLimit = oMapLimits.find(osTMId);
1324 48 : if (oIterLimit != oMapLimits.end())
1325 : {
1326 39 : m_nMinMetaTileCol = oIterLimit->second.min_tile_col;
1327 39 : m_nMinMetaTileRow = oIterLimit->second.min_tile_row;
1328 39 : nMatrixWidth = oIterLimit->second.max_tile_col - m_nMinMetaTileCol + 1;
1329 39 : nMatrixHeight = oIterLimit->second.max_tile_row - m_nMinMetaTileRow + 1;
1330 : }
1331 48 : m_nMetaTileWidth = oTM.mTileWidth;
1332 48 : m_nMetaTileHeight = oTM.mTileHeight;
1333 48 : nRasterXSize = nMatrixWidth * m_nMetaTileWidth;
1334 48 : nRasterYSize = nMatrixHeight * m_nMetaTileHeight;
1335 :
1336 48 : if (poProtoDS)
1337 : {
1338 132 : for (int i = 0; i < poProtoDS->GetRasterCount(); i++)
1339 : {
1340 99 : auto poProtoBand = poProtoDS->GetRasterBand(i + 1);
1341 99 : auto poBand = new STACTARawRasterBand(this, i + 1, poProtoBand);
1342 99 : SetBand(i + 1, poBand);
1343 : }
1344 : }
1345 : else
1346 : {
1347 105 : for (int i = 0; i < static_cast<int>(aeDT.size()); i++)
1348 : {
1349 90 : auto poBand = new STACTARawRasterBand(this, i + 1, aeDT[i],
1350 90 : abSetNoData[i], adfNoData[i]);
1351 90 : SetBand(i + 1, poBand);
1352 : }
1353 : }
1354 :
1355 96 : CPLString osCRS = poTMS->crs().c_str();
1356 48 : if (osCRS == "http://www.opengis.net/def/crs/OGC/1.3/CRS84")
1357 48 : osCRS = "EPSG:4326";
1358 48 : if (m_oSRS.SetFromUserInput(osCRS) != OGRERR_NONE)
1359 : {
1360 0 : return false;
1361 : }
1362 48 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1363 96 : m_adfGeoTransform[0] =
1364 48 : oTM.mTopLeftX + m_nMinMetaTileCol * m_nMetaTileWidth * oTM.mResX;
1365 48 : m_adfGeoTransform[1] = oTM.mResX;
1366 96 : m_adfGeoTransform[3] =
1367 48 : oTM.mTopLeftY - m_nMinMetaTileRow * m_nMetaTileHeight * oTM.mResY;
1368 48 : m_adfGeoTransform[5] = -oTM.mResY;
1369 48 : SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
1370 :
1371 48 : return true;
1372 : }
1373 :
1374 : /************************************************************************/
1375 : /* GetSpatialRef () */
1376 : /************************************************************************/
1377 :
1378 3 : const OGRSpatialReference *STACTADataset::GetSpatialRef() const
1379 : {
1380 3 : return nBands == 0 ? nullptr : &m_oSRS;
1381 : }
1382 :
1383 : /************************************************************************/
1384 : /* GetGeoTransform() */
1385 : /************************************************************************/
1386 :
1387 3 : CPLErr STACTADataset::GetGeoTransform(double *padfGeoTransform)
1388 : {
1389 3 : memcpy(padfGeoTransform, &m_adfGeoTransform[0], 6 * sizeof(double));
1390 3 : return nBands == 0 ? CE_Failure : CE_None;
1391 : }
1392 :
1393 : /************************************************************************/
1394 : /* OpenStatic() */
1395 : /************************************************************************/
1396 :
1397 23 : GDALDataset *STACTADataset::OpenStatic(GDALOpenInfo *poOpenInfo)
1398 : {
1399 23 : if (!Identify(poOpenInfo))
1400 0 : return nullptr;
1401 46 : auto poDS = std::make_unique<STACTADataset>();
1402 23 : if (!poDS->Open(poOpenInfo))
1403 5 : return nullptr;
1404 18 : return poDS.release();
1405 : }
1406 :
1407 : /************************************************************************/
1408 : /* GDALRegister_STACTA() */
1409 : /************************************************************************/
1410 :
1411 1595 : void GDALRegister_STACTA()
1412 :
1413 : {
1414 1595 : if (GDALGetDriverByName("STACTA") != nullptr)
1415 302 : return;
1416 :
1417 1293 : GDALDriver *poDriver = new GDALDriver();
1418 :
1419 1293 : poDriver->SetDescription("STACTA");
1420 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1421 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1422 1293 : "Spatio-Temporal Asset Catalog Tiled Assets");
1423 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/stacta.html");
1424 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "json");
1425 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1426 1293 : poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");
1427 1293 : poDriver->SetMetadataItem(
1428 : GDAL_DMD_OPENOPTIONLIST,
1429 : "<OpenOptionList>"
1430 : " <Option name='WHOLE_METATILE' type='boolean' "
1431 : "description='Whether to download whole metatiles'/>"
1432 : " <Option name='SKIP_MISSING_METATILE' type='boolean' "
1433 : "description='Whether to gracefully skip missing metatiles'/>"
1434 1293 : "</OpenOptionList>");
1435 :
1436 1293 : poDriver->pfnOpen = STACTADataset::OpenStatic;
1437 1293 : poDriver->pfnIdentify = STACTADataset::Identify;
1438 :
1439 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1440 : }
|
