Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: STACIT (Spatio-Temporal Asset Catalog ITems) driver
5 : * Author: Even Rouault, <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2021, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_json.h"
14 : #include "cpl_http.h"
15 : #include "vrtdataset.h"
16 : #include "ogr_spatialref.h"
17 :
18 : #include <algorithm>
19 : #include <limits>
20 : #include <map>
21 : #include <string>
22 :
23 : namespace
24 : {
25 :
26 : struct AssetItem
27 : {
28 : std::string osFilename{};
29 : std::string osDateTime{};
30 : int nXSize = 0;
31 : int nYSize = 0;
32 : double dfXMin = 0;
33 : double dfYMin = 0;
34 : double dfXMax = 0;
35 : double dfYMax = 0;
36 : };
37 :
38 : struct AssetSetByProjection
39 : {
40 : std::string osProjUserString{};
41 : std::vector<AssetItem> assets{};
42 : };
43 :
44 : struct Asset
45 : {
46 : std::string osName{};
47 : CPLJSONArray eoBands{};
48 : std::map<std::string, AssetSetByProjection> assets{};
49 : };
50 :
51 : struct Collection
52 : {
53 : std::string osName{};
54 : std::map<std::string, Asset> assets{};
55 : };
56 : } // namespace
57 :
58 : /************************************************************************/
59 : /* STACITDataset */
60 : /************************************************************************/
61 :
62 : class STACITDataset final : public VRTDataset
63 : {
64 : bool Open(GDALOpenInfo *poOpenInfo);
65 : bool SetupDataset(GDALOpenInfo *poOpenInfo,
66 : const std::string &osSTACITFilename,
67 : std::map<std::string, Collection> &oMapCollection);
68 : void
69 : SetSubdatasets(const std::string &osFilename,
70 : const std::map<std::string, Collection> &oMapCollection);
71 :
72 : public:
73 : STACITDataset();
74 :
75 : static int Identify(GDALOpenInfo *poOpenInfo);
76 : static GDALDataset *OpenStatic(GDALOpenInfo *poOpenInfo);
77 : };
78 :
79 : /************************************************************************/
80 : /* STACITDataset() */
81 : /************************************************************************/
82 :
83 20 : STACITDataset::STACITDataset() : VRTDataset(0, 0)
84 : {
85 20 : poDriver = nullptr; // cancel what the VRTDataset did
86 20 : SetWritable(false);
87 20 : }
88 :
89 : /************************************************************************/
90 : /* Identify() */
91 : /************************************************************************/
92 :
93 50475 : int STACITDataset::Identify(GDALOpenInfo *poOpenInfo)
94 : {
95 50475 : if (STARTS_WITH(poOpenInfo->pszFilename, "STACIT:"))
96 : {
97 12 : return true;
98 : }
99 :
100 50463 : const bool bIsSingleDriver = poOpenInfo->IsSingleAllowedDriver("STACIT");
101 50455 : if (bIsSingleDriver && (STARTS_WITH(poOpenInfo->pszFilename, "http://") ||
102 4 : STARTS_WITH(poOpenInfo->pszFilename, "https://")))
103 : {
104 1 : return true;
105 : }
106 :
107 50454 : if (poOpenInfo->nHeaderBytes == 0)
108 : {
109 47595 : return false;
110 : }
111 :
112 8514 : for (int i = 0; i < 2; i++)
113 : {
114 : // TryToIngest() may reallocate pabyHeader, so do not move this
115 : // before the loop.
116 5686 : const char *pszHeader =
117 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
118 6002 : while (*pszHeader != 0 &&
119 4306 : std::isspace(static_cast<unsigned char>(*pszHeader)))
120 316 : ++pszHeader;
121 5686 : if (bIsSingleDriver)
122 : {
123 4 : return pszHeader[0] == '{';
124 : }
125 :
126 5682 : if (strstr(pszHeader, "\"stac_version\"") != nullptr &&
127 26 : strstr(pszHeader, "\"proj:transform\"") != nullptr)
128 : {
129 26 : return true;
130 : }
131 :
132 5656 : if (i == 0)
133 : {
134 : // Should be enough for a STACIT .json file
135 2828 : poOpenInfo->TryToIngest(32768);
136 : }
137 : }
138 :
139 2828 : return false;
140 : }
141 :
142 : /************************************************************************/
143 : /* SanitizeCRSValue() */
144 : /************************************************************************/
145 :
146 26 : static std::string SanitizeCRSValue(const std::string &v)
147 : {
148 26 : std::string ret;
149 26 : bool lastWasAlphaNum = true;
150 86 : for (char ch : v)
151 : {
152 60 : if (!isalnum(static_cast<unsigned char>(ch)))
153 : {
154 6 : if (lastWasAlphaNum)
155 6 : ret += '_';
156 6 : lastWasAlphaNum = false;
157 : }
158 : else
159 : {
160 54 : ret += ch;
161 54 : lastWasAlphaNum = true;
162 : }
163 : }
164 26 : if (!ret.empty() && ret.back() == '_')
165 0 : ret.pop_back();
166 26 : return ret;
167 : }
168 :
169 : /************************************************************************/
170 : /* ParseAsset() */
171 : /************************************************************************/
172 :
173 43 : static void ParseAsset(const CPLJSONObject &jAsset,
174 : const CPLJSONObject &oProperties,
175 : const std::string &osCollection,
176 : const std::string &osFilteredCRS,
177 : std::map<std::string, Collection> &oMapCollection)
178 : {
179 : // Skip assets that are obviously not images
180 86 : const auto osType = jAsset["type"].ToString();
181 81 : if (osType == "application/json" || osType == "application/xml" ||
182 38 : osType == "text/plain")
183 : {
184 5 : return;
185 : }
186 :
187 : // Skip assets whose role is obviously non georeferenced
188 76 : const auto oRoles = jAsset.GetArray("roles");
189 38 : if (oRoles.IsValid())
190 : {
191 76 : for (const auto &oRole : oRoles)
192 : {
193 76 : const auto osRole = oRole.ToString();
194 76 : if (osRole == "thumbnail" || osRole == "info" ||
195 38 : osRole == "metadata")
196 : {
197 0 : return;
198 : }
199 : }
200 : }
201 :
202 38 : const auto osAssetName = jAsset.GetName();
203 :
204 76 : const auto osHref = jAsset["href"].ToString();
205 38 : if (osHref.empty())
206 : {
207 0 : CPLError(CE_Warning, CPLE_AppDefined, "Missing href on asset %s",
208 : osAssetName.c_str());
209 0 : return;
210 : }
211 :
212 : const auto GetAssetOrFeatureProperty =
213 217 : [&oProperties, &jAsset](const char *pszName)
214 : {
215 438 : auto obj = jAsset[pszName];
216 146 : if (obj.IsValid())
217 75 : return obj;
218 71 : return oProperties[pszName];
219 38 : };
220 :
221 38 : auto oProjEPSG = GetAssetOrFeatureProperty("proj:epsg");
222 38 : std::string osProjUserString;
223 38 : if (oProjEPSG.IsValid() && oProjEPSG.GetType() != CPLJSONObject::Type::Null)
224 : {
225 38 : osProjUserString = "EPSG:" + oProjEPSG.ToString();
226 : }
227 : else
228 : {
229 0 : auto oProjWKT2 = GetAssetOrFeatureProperty("proj:wkt2");
230 0 : if (oProjWKT2.IsValid() &&
231 0 : oProjWKT2.GetType() == CPLJSONObject::Type::String)
232 : {
233 0 : osProjUserString = oProjWKT2.ToString();
234 : }
235 : else
236 : {
237 0 : auto oProjPROJJSON = GetAssetOrFeatureProperty("proj:projjson");
238 0 : if (oProjPROJJSON.IsValid() &&
239 0 : oProjPROJJSON.GetType() == CPLJSONObject::Type::Object)
240 : {
241 0 : osProjUserString = oProjPROJJSON.ToString();
242 : }
243 : else
244 : {
245 0 : CPLDebug("STACIT",
246 : "Skipping asset %s that lacks a valid CRS member",
247 : osAssetName.c_str());
248 0 : return;
249 : }
250 : }
251 : }
252 :
253 42 : if (!osFilteredCRS.empty() &&
254 42 : osFilteredCRS != SanitizeCRSValue(osProjUserString))
255 : {
256 2 : return;
257 : }
258 :
259 36 : AssetItem item;
260 36 : item.osFilename = osHref;
261 36 : item.osDateTime = oProperties["datetime"].ToString();
262 :
263 : // Figure out item bounds and width/height
264 36 : auto oProjBBOX = GetAssetOrFeatureProperty("proj:bbox").ToArray();
265 36 : auto oProjShape = GetAssetOrFeatureProperty("proj:shape").ToArray();
266 36 : auto oProjTransform = GetAssetOrFeatureProperty("proj:transform").ToArray();
267 36 : const bool bIsBBOXValid = oProjBBOX.IsValid() && oProjBBOX.Size() == 4;
268 36 : const bool bIsShapeValid = oProjShape.IsValid() && oProjShape.Size() == 2;
269 : const bool bIsTransformValid =
270 72 : oProjTransform.IsValid() &&
271 36 : (oProjTransform.Size() == 6 || oProjTransform.Size() == 9);
272 36 : std::vector<double> bbox;
273 36 : if (bIsBBOXValid)
274 : {
275 125 : for (const auto &oItem : oProjBBOX)
276 100 : bbox.push_back(oItem.ToDouble());
277 25 : CPLAssert(bbox.size() == 4);
278 : }
279 36 : std::vector<int> shape;
280 36 : if (bIsShapeValid)
281 : {
282 33 : for (const auto &oItem : oProjShape)
283 22 : shape.push_back(oItem.ToInteger());
284 11 : CPLAssert(shape.size() == 2);
285 : }
286 36 : std::vector<double> transform;
287 36 : if (bIsTransformValid)
288 : {
289 345 : for (const auto &oItem : oProjTransform)
290 309 : transform.push_back(oItem.ToDouble());
291 36 : CPLAssert(transform.size() == 6 || transform.size() == 9);
292 : #if defined(__GNUC__)
293 : #pragma GCC diagnostic push
294 : #pragma GCC diagnostic ignored "-Wnull-dereference"
295 : #endif
296 72 : if (transform[0] <= 0 || transform[1] != 0 || transform[3] != 0 ||
297 139 : transform[4] >= 0 ||
298 36 : (transform.size() == 9 &&
299 31 : (transform[6] != 0 || transform[7] != 0 || transform[8] != 1)))
300 : {
301 0 : CPLError(
302 : CE_Warning, CPLE_AppDefined,
303 : "Skipping asset %s because its proj:transform is "
304 : "not of the form [xres,0,xoffset,0,yres<0,yoffset[,0,0,1]]",
305 : osAssetName.c_str());
306 0 : return;
307 : }
308 : #if defined(__GNUC__)
309 : #pragma GCC diagnostic pop
310 : #endif
311 : }
312 :
313 36 : if (bIsBBOXValid && bIsShapeValid)
314 : {
315 0 : item.nXSize = shape[1];
316 0 : item.nYSize = shape[0];
317 0 : item.dfXMin = bbox[0];
318 0 : item.dfYMin = bbox[1];
319 0 : item.dfXMax = bbox[2];
320 0 : item.dfYMax = bbox[3];
321 : }
322 36 : else if (bIsBBOXValid && bIsTransformValid)
323 : {
324 25 : item.dfXMin = bbox[0];
325 25 : item.dfYMin = bbox[1];
326 25 : item.dfXMax = bbox[2];
327 25 : item.dfYMax = bbox[3];
328 25 : if (item.dfXMin != transform[2] || item.dfYMax != transform[5])
329 : {
330 0 : CPLError(CE_Warning, CPLE_AppDefined,
331 : "Skipping asset %s because the origin of "
332 : "proj:transform and proj:bbox are not consistent",
333 : osAssetName.c_str());
334 0 : return;
335 : }
336 25 : double dfXSize = (item.dfXMax - item.dfXMin) / transform[0];
337 25 : double dfYSize = (item.dfYMax - item.dfYMin) / -transform[4];
338 25 : if (!(dfXSize < INT_MAX && dfYSize < INT_MAX))
339 0 : return;
340 25 : item.nXSize = static_cast<int>(dfXSize);
341 25 : item.nYSize = static_cast<int>(dfYSize);
342 : }
343 11 : else if (bIsShapeValid && bIsTransformValid)
344 : {
345 11 : item.nXSize = shape[1];
346 11 : item.nYSize = shape[0];
347 11 : item.dfXMin = transform[2];
348 11 : item.dfYMax = transform[5];
349 11 : item.dfXMax = item.dfXMin + item.nXSize * transform[0];
350 11 : item.dfYMin = item.dfYMax + item.nYSize * transform[4];
351 : }
352 : else
353 : {
354 0 : CPLDebug("STACIT",
355 : "Skipping asset %s that lacks at least 2 members among "
356 : "proj:bbox, proj:shape and proj:transform",
357 : osAssetName.c_str());
358 0 : return;
359 : }
360 :
361 36 : if (item.nXSize <= 0 || item.nYSize <= 0)
362 : {
363 0 : CPLError(CE_Warning, CPLE_AppDefined,
364 : "Skipping asset %s because the size is invalid",
365 : osAssetName.c_str());
366 0 : return;
367 : }
368 :
369 : // Create/fetch collection
370 36 : if (oMapCollection.find(osCollection) == oMapCollection.end())
371 : {
372 38 : Collection collection;
373 19 : collection.osName = osCollection;
374 19 : oMapCollection[osCollection] = std::move(collection);
375 : }
376 36 : auto &collection = oMapCollection[osCollection];
377 :
378 : // Create/fetch asset in collection
379 36 : if (collection.assets.find(osAssetName) == collection.assets.end())
380 : {
381 40 : Asset asset;
382 20 : asset.osName = osAssetName;
383 20 : asset.eoBands = jAsset.GetArray("eo:bands");
384 :
385 20 : collection.assets[osAssetName] = std::move(asset);
386 : }
387 36 : auto &asset = collection.assets[osAssetName];
388 :
389 : // Create/fetch projection in asset
390 36 : if (asset.assets.find(osProjUserString) == asset.assets.end())
391 : {
392 42 : AssetSetByProjection assetByProj;
393 21 : assetByProj.osProjUserString = osProjUserString;
394 21 : asset.assets[osProjUserString] = std::move(assetByProj);
395 : }
396 36 : auto &assets = asset.assets[osProjUserString];
397 :
398 : // Add item
399 36 : assets.assets.emplace_back(item);
400 : }
401 :
402 : /************************************************************************/
403 : /* SetupDataset() */
404 : /************************************************************************/
405 :
406 17 : bool STACITDataset::SetupDataset(
407 : GDALOpenInfo *poOpenInfo, const std::string &osSTACITFilename,
408 : std::map<std::string, Collection> &oMapCollection)
409 : {
410 17 : auto &collection = oMapCollection.begin()->second;
411 17 : auto &asset = collection.assets.begin()->second;
412 17 : auto &assetByProj = asset.assets.begin()->second;
413 17 : auto &items = assetByProj.assets;
414 :
415 : // Compute global bounds and resolution
416 17 : double dfXMin = std::numeric_limits<double>::max();
417 17 : double dfYMin = std::numeric_limits<double>::max();
418 17 : double dfXMax = -std::numeric_limits<double>::max();
419 17 : double dfYMax = -std::numeric_limits<double>::max();
420 17 : double dfXRes = 0;
421 17 : double dfYRes = 0;
422 34 : const char *pszResolution = CSLFetchNameValueDef(
423 17 : poOpenInfo->papszOpenOptions, "RESOLUTION", "AVERAGE");
424 49 : for (const auto &assetItem : items)
425 : {
426 32 : dfXMin = std::min(dfXMin, assetItem.dfXMin);
427 32 : dfYMin = std::min(dfYMin, assetItem.dfYMin);
428 32 : dfXMax = std::max(dfXMax, assetItem.dfXMax);
429 32 : dfYMax = std::max(dfYMax, assetItem.dfYMax);
430 32 : const double dfThisXRes =
431 32 : (assetItem.dfXMax - assetItem.dfXMin) / assetItem.nXSize;
432 32 : const double dfThisYRes =
433 32 : (assetItem.dfYMax - assetItem.dfYMin) / assetItem.nYSize;
434 : #ifdef DEBUG_VERBOSE
435 : CPLDebug("STACIT", "%s -> resx=%f resy=%f",
436 : assetItem.osFilename.c_str(), dfThisXRes, dfThisYRes);
437 : #endif
438 32 : if (dfXRes != 0 && EQUAL(pszResolution, "HIGHEST"))
439 : {
440 0 : dfXRes = std::min(dfXRes, dfThisXRes);
441 0 : dfYRes = std::min(dfYRes, dfThisYRes);
442 : }
443 32 : else if (dfXRes != 0 && EQUAL(pszResolution, "LOWEST"))
444 : {
445 0 : dfXRes = std::max(dfXRes, dfThisXRes);
446 0 : dfYRes = std::max(dfYRes, dfThisYRes);
447 : }
448 : else
449 : {
450 32 : dfXRes += dfThisXRes;
451 32 : dfYRes += dfThisYRes;
452 : }
453 : }
454 17 : if (EQUAL(pszResolution, "AVERAGE"))
455 : {
456 17 : dfXRes /= static_cast<int>(items.size());
457 17 : dfYRes /= static_cast<int>(items.size());
458 : }
459 :
460 : // Set raster size
461 17 : double dfXSize = std::round((dfXMax - dfXMin) / dfXRes);
462 17 : double dfYSize = std::round((dfYMax - dfYMin) / dfYRes);
463 17 : if (dfXSize <= 0 || dfYSize <= 0 || dfXSize > INT_MAX || dfYSize > INT_MAX)
464 : {
465 0 : CPLError(CE_Failure, CPLE_AppDefined,
466 : "Invalid computed dataset dimensions");
467 0 : return false;
468 : }
469 17 : nRasterXSize = static_cast<int>(dfXSize);
470 17 : nRasterYSize = static_cast<int>(dfYSize);
471 :
472 : // Set geotransform
473 : double adfGeoTransform[6];
474 17 : adfGeoTransform[0] = dfXMin;
475 17 : adfGeoTransform[1] = dfXRes;
476 17 : adfGeoTransform[2] = 0;
477 17 : adfGeoTransform[3] = dfYMax;
478 17 : adfGeoTransform[4] = 0;
479 17 : adfGeoTransform[5] = -dfYRes;
480 17 : SetGeoTransform(adfGeoTransform);
481 :
482 : // Set SRS
483 34 : OGRSpatialReference oSRS;
484 17 : if (oSRS.SetFromUserInput(
485 : assetByProj.osProjUserString.c_str(),
486 17 : OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get()) ==
487 : OGRERR_NONE)
488 : {
489 17 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
490 17 : SetSpatialRef(&oSRS);
491 : }
492 :
493 : // Open of the items to find the number of bands, their data type
494 : // and nodata.
495 : const auto BuildVSICurlFilename =
496 49 : [&osSTACITFilename, &collection](const std::string &osFilename)
497 : {
498 49 : std::string osRet;
499 49 : if (STARTS_WITH(osFilename.c_str(), "http"))
500 : {
501 0 : if (STARTS_WITH(osSTACITFilename.c_str(),
502 : "https://planetarycomputer.microsoft.com/api/"))
503 : {
504 0 : osRet = "/vsicurl?pc_url_signing=yes&pc_collection=";
505 0 : osRet += collection.osName;
506 0 : osRet += "&url=";
507 : char *pszEncoded =
508 0 : CPLEscapeString(osFilename.c_str(), -1, CPLES_URL);
509 0 : CPLString osEncoded(pszEncoded);
510 0 : CPLFree(pszEncoded);
511 : // something is confused if the whole URL appears as
512 : // /vsicurl...blabla_without_slash.tif" so transform back %2F to
513 : // /
514 0 : osEncoded.replaceAll("%2F", '/');
515 0 : osRet += osEncoded;
516 : }
517 : else
518 : {
519 0 : osRet = "/vsicurl/";
520 0 : osRet += osFilename;
521 : }
522 : }
523 49 : else if (STARTS_WITH(osFilename.c_str(), "file://"))
524 : {
525 3 : osRet = osFilename.substr(strlen("file://"));
526 : }
527 46 : else if (STARTS_WITH(osFilename.c_str(), "s3://"))
528 : {
529 0 : osRet = "/vsis3/";
530 0 : osRet += osFilename.substr(strlen("s3://"));
531 : }
532 :
533 : else
534 : {
535 46 : osRet = osFilename;
536 : }
537 49 : return osRet;
538 17 : };
539 :
540 34 : const auto osFirstItemName(BuildVSICurlFilename(items.front().osFilename));
541 : auto poItemDS = std::unique_ptr<GDALDataset>(
542 34 : GDALDataset::Open(osFirstItemName.c_str()));
543 17 : if (!poItemDS)
544 : {
545 0 : CPLError(CE_Failure, CPLE_AppDefined,
546 : "Cannot open %s to retrieve band characteristics",
547 : osFirstItemName.c_str());
548 0 : return false;
549 : }
550 :
551 : // Sort by ascending datetime
552 17 : std::sort(items.begin(), items.end(),
553 26 : [](const AssetItem &a, const AssetItem &b)
554 : {
555 26 : if (!a.osDateTime.empty() && !b.osDateTime.empty())
556 26 : return a.osDateTime < b.osDateTime;
557 0 : return &a < &b;
558 : });
559 :
560 : // Create VRT bands and add sources
561 17 : bool bAtLeastOneBandHasNoData = false;
562 34 : for (int i = 0; i < poItemDS->GetRasterCount(); i++)
563 : {
564 17 : auto poItemBand = poItemDS->GetRasterBand(i + 1);
565 17 : AddBand(poItemBand->GetRasterDataType(), nullptr);
566 : auto poVRTBand =
567 17 : cpl::down_cast<VRTSourcedRasterBand *>(GetRasterBand(i + 1));
568 17 : int bHasNoData = FALSE;
569 17 : const double dfNoData = poItemBand->GetNoDataValue(&bHasNoData);
570 17 : if (bHasNoData)
571 : {
572 4 : bAtLeastOneBandHasNoData = true;
573 4 : poVRTBand->SetNoDataValue(dfNoData);
574 : }
575 :
576 17 : const auto eInterp = poItemBand->GetColorInterpretation();
577 17 : if (eInterp != GCI_Undefined)
578 17 : poVRTBand->SetColorInterpretation(eInterp);
579 :
580 : // Set band properties
581 34 : if (asset.eoBands.IsValid() &&
582 17 : asset.eoBands.Size() == poItemDS->GetRasterCount())
583 : {
584 34 : const auto &eoBand = asset.eoBands[i];
585 51 : const auto osBandName = eoBand["name"].ToString();
586 17 : if (!osBandName.empty())
587 17 : poVRTBand->SetDescription(osBandName.c_str());
588 :
589 51 : const auto osCommonName = eoBand["common_name"].ToString();
590 17 : if (!osCommonName.empty())
591 : {
592 : const auto eInterpFromCommonName =
593 13 : GDALGetColorInterpFromSTACCommonName(osCommonName.c_str());
594 13 : if (eInterpFromCommonName != GCI_Undefined)
595 13 : poVRTBand->SetColorInterpretation(eInterpFromCommonName);
596 : }
597 :
598 73 : for (const auto &eoBandChild : eoBand.GetChildren())
599 : {
600 112 : const auto osChildName = eoBandChild.GetName();
601 56 : if (osChildName != "name" && osChildName != "common_name")
602 : {
603 26 : poVRTBand->SetMetadataItem(osChildName.c_str(),
604 52 : eoBandChild.ToString().c_str());
605 : }
606 : }
607 : }
608 :
609 : // Add items as VRT sources
610 49 : for (const auto &assetItem : items)
611 : {
612 64 : const auto osItemName(BuildVSICurlFilename(assetItem.osFilename));
613 32 : const double dfDstXOff = (assetItem.dfXMin - dfXMin) / dfXRes;
614 32 : const double dfDstXSize =
615 32 : (assetItem.dfXMax - assetItem.dfXMin) / dfXRes;
616 32 : const double dfDstYOff = (dfYMax - assetItem.dfYMax) / dfYRes;
617 32 : const double dfDstYSize =
618 32 : (assetItem.dfYMax - assetItem.dfYMin) / dfYRes;
619 32 : if (!bHasNoData)
620 : {
621 24 : poVRTBand->AddSimpleSource(osItemName.c_str(), i + 1, 0, 0,
622 24 : assetItem.nXSize, assetItem.nYSize,
623 : dfDstXOff, dfDstYOff, dfDstXSize,
624 : dfDstYSize);
625 : }
626 : else
627 : {
628 8 : poVRTBand->AddComplexSource(osItemName.c_str(), i + 1, 0, 0,
629 8 : assetItem.nXSize, assetItem.nYSize,
630 : dfDstXOff, dfDstYOff, dfDstXSize,
631 : dfDstYSize,
632 : 0.0, // offset
633 : 1.0, // scale
634 : dfNoData);
635 : }
636 : }
637 :
638 : const char *pszOverlapStrategy =
639 17 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
640 : "OVERLAP_STRATEGY", "REMOVE_IF_NO_NODATA");
641 17 : if ((EQUAL(pszOverlapStrategy, "REMOVE_IF_NO_NODATA") &&
642 13 : !bAtLeastOneBandHasNoData) ||
643 6 : EQUAL(pszOverlapStrategy, "USE_MOST_RECENT"))
644 : {
645 13 : const char *const apszOptions[] = {
646 : "EMIT_ERROR_IF_GEOS_NOT_AVAILABLE=NO", nullptr};
647 13 : poVRTBand->RemoveCoveredSources(apszOptions);
648 : }
649 : }
650 17 : return true;
651 : }
652 :
653 : /************************************************************************/
654 : /* SetSubdatasets() */
655 : /************************************************************************/
656 :
657 1 : void STACITDataset::SetSubdatasets(
658 : const std::string &osFilename,
659 : const std::map<std::string, Collection> &oMapCollection)
660 : {
661 2 : CPLStringList aosSubdatasets;
662 1 : int nCount = 1;
663 3 : for (const auto &collectionKV : oMapCollection)
664 : {
665 5 : for (const auto &assetKV : collectionKV.second.assets)
666 : {
667 6 : std::string osCollectionAssetArg;
668 3 : if (oMapCollection.size() > 1)
669 : osCollectionAssetArg +=
670 3 : "collection=" + collectionKV.first + ",";
671 3 : osCollectionAssetArg += "asset=" + assetKV.first;
672 :
673 6 : std::string osCollectionAssetText;
674 3 : if (oMapCollection.size() > 1)
675 : osCollectionAssetText +=
676 3 : "Collection " + collectionKV.first + ", ";
677 3 : osCollectionAssetText += "Asset " + assetKV.first;
678 :
679 3 : if (assetKV.second.assets.size() == 1)
680 : {
681 : aosSubdatasets.AddString(CPLSPrintf(
682 : "SUBDATASET_%d_NAME=STACIT:\"%s\":%s", nCount,
683 2 : osFilename.c_str(), osCollectionAssetArg.c_str()));
684 : aosSubdatasets.AddString(CPLSPrintf(
685 : "SUBDATASET_%d_DESC=%s of %s", nCount,
686 2 : osCollectionAssetText.c_str(), osFilename.c_str()));
687 2 : nCount++;
688 : }
689 : else
690 : {
691 3 : for (const auto &assetsByProjKV : assetKV.second.assets)
692 : {
693 : aosSubdatasets.AddString(CPLSPrintf(
694 : "SUBDATASET_%d_NAME=STACIT:\"%s\":%s,crs=%s", nCount,
695 : osFilename.c_str(), osCollectionAssetArg.c_str(),
696 2 : SanitizeCRSValue(assetsByProjKV.first).c_str()));
697 : aosSubdatasets.AddString(CPLSPrintf(
698 : "SUBDATASET_%d_DESC=%s of %s in CRS %s", nCount,
699 : osCollectionAssetText.c_str(), osFilename.c_str(),
700 2 : assetsByProjKV.first.c_str()));
701 2 : nCount++;
702 : }
703 : }
704 : }
705 : }
706 1 : GDALDataset::SetMetadata(aosSubdatasets.List(), "SUBDATASETS");
707 1 : }
708 :
709 : /************************************************************************/
710 : /* Open() */
711 : /************************************************************************/
712 :
713 20 : bool STACITDataset::Open(GDALOpenInfo *poOpenInfo)
714 : {
715 40 : CPLString osFilename(poOpenInfo->pszFilename);
716 : std::string osFilteredCollection =
717 40 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "COLLECTION", "");
718 : std::string osFilteredAsset =
719 40 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "ASSET", "");
720 : std::string osFilteredCRS = SanitizeCRSValue(
721 60 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "CRS", ""));
722 20 : if (STARTS_WITH(poOpenInfo->pszFilename, "STACIT:"))
723 : {
724 : const CPLStringList aosTokens(CSLTokenizeString2(
725 6 : poOpenInfo->pszFilename, ":", CSLT_HONOURSTRINGS));
726 6 : if (aosTokens.size() != 2 && aosTokens.size() != 3)
727 0 : return false;
728 6 : osFilename = aosTokens[1];
729 6 : if (aosTokens.size() >= 3)
730 : {
731 : const CPLStringList aosFilters(
732 12 : CSLTokenizeString2(aosTokens[2], ",", 0));
733 : osFilteredCollection = aosFilters.FetchNameValueDef(
734 6 : "collection", osFilteredCollection.c_str());
735 : osFilteredAsset =
736 6 : aosFilters.FetchNameValueDef("asset", osFilteredAsset.c_str());
737 : osFilteredCRS =
738 6 : aosFilters.FetchNameValueDef("crs", osFilteredCRS.c_str());
739 : }
740 : }
741 :
742 40 : std::map<std::string, Collection> oMapCollection;
743 20 : GIntBig nItemIter = 0;
744 20 : GIntBig nMaxItems = CPLAtoGIntBig(CSLFetchNameValueDef(
745 20 : poOpenInfo->papszOpenOptions, "MAX_ITEMS", "1000"));
746 :
747 : const bool bMaxItemsSpecified =
748 20 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MAX_ITEMS") != nullptr;
749 20 : if (!bMaxItemsSpecified)
750 : {
751 : // If the URL includes a limit parameter, and it's larger than our
752 : // default MAX_ITEMS value, then increase the later to the former.
753 19 : const auto nPos = osFilename.ifind("&limit=");
754 19 : if (nPos != std::string::npos)
755 : {
756 0 : const auto nLimit = CPLAtoGIntBig(
757 0 : osFilename.substr(nPos + strlen("&limit=")).c_str());
758 0 : nMaxItems = std::max(nMaxItems, nLimit);
759 : }
760 : }
761 :
762 40 : auto osCurFilename = osFilename;
763 40 : std::string osMethod = "GET";
764 40 : CPLJSONObject oHeaders;
765 40 : CPLJSONObject oBody;
766 20 : bool bMerge = false;
767 20 : int nLoops = 0;
768 3 : do
769 : {
770 23 : ++nLoops;
771 23 : if (nMaxItems > 0 && nLoops > nMaxItems)
772 : {
773 20 : break;
774 : }
775 :
776 23 : CPLJSONDocument oDoc;
777 :
778 45 : if (STARTS_WITH(osCurFilename, "http://") ||
779 22 : STARTS_WITH(osCurFilename, "https://"))
780 : {
781 : // Cf // Cf https://github.com/radiantearth/stac-api-spec/tree/release/v1.0.0/item-search#pagination
782 1 : CPLStringList aosOptions;
783 2 : if (oBody.IsValid() &&
784 1 : oBody.GetType() == CPLJSONObject::Type::Object)
785 : {
786 1 : if (bMerge)
787 0 : CPLDebug("STACIT",
788 : "Ignoring 'merge' attribute from next link");
789 : const std::string osPostContent =
790 2 : oBody.Format(CPLJSONObject::PrettyFormat::Pretty);
791 1 : aosOptions.SetNameValue("POSTFIELDS", osPostContent.c_str());
792 : }
793 1 : aosOptions.SetNameValue("CUSTOMREQUEST", osMethod.c_str());
794 1 : CPLString osHeaders;
795 4 : if (!oHeaders.IsValid() ||
796 2 : oHeaders.GetType() != CPLJSONObject::Type::Object ||
797 2 : oHeaders["Content-Type"].ToString().empty())
798 : {
799 1 : osHeaders = "Content-Type: application/json";
800 : }
801 2 : if (oHeaders.IsValid() &&
802 1 : oHeaders.GetType() == CPLJSONObject::Type::Object)
803 : {
804 2 : for (const auto &obj : oHeaders.GetChildren())
805 : {
806 1 : osHeaders += "\r\n";
807 1 : osHeaders += obj.GetName();
808 1 : osHeaders += ": ";
809 1 : osHeaders += obj.ToString();
810 : }
811 : }
812 1 : aosOptions.SetNameValue("HEADERS", osHeaders.c_str());
813 : CPLHTTPResult *psResult =
814 1 : CPLHTTPFetch(osCurFilename.c_str(), aosOptions.List());
815 1 : if (!psResult)
816 0 : return false;
817 1 : if (!psResult->pabyData)
818 : {
819 0 : CPLHTTPDestroyResult(psResult);
820 0 : return false;
821 : }
822 2 : const bool bOK = oDoc.LoadMemory(
823 1 : reinterpret_cast<const char *>(psResult->pabyData));
824 : // CPLDebug("STACIT", "Response: %s", reinterpret_cast<const char*>(psResult->pabyData));
825 1 : CPLHTTPDestroyResult(psResult);
826 1 : if (!bOK)
827 0 : return false;
828 : }
829 : else
830 : {
831 22 : if (!oDoc.Load(osCurFilename))
832 0 : return false;
833 : }
834 23 : const auto oRoot = oDoc.GetRoot();
835 46 : auto oFeatures = oRoot.GetArray("features");
836 23 : if (!oFeatures.IsValid())
837 : {
838 1 : if (oRoot.GetString("type") == "Feature")
839 : {
840 1 : oFeatures = CPLJSONArray();
841 1 : oFeatures.Add(oRoot);
842 : }
843 : else
844 : {
845 0 : CPLError(CE_Failure, CPLE_AppDefined, "Missing features");
846 0 : return false;
847 : }
848 : }
849 72 : for (const auto &oFeature : oFeatures)
850 : {
851 49 : nItemIter++;
852 49 : if (nMaxItems > 0 && nItemIter > nMaxItems)
853 : {
854 0 : break;
855 : }
856 :
857 98 : auto oStacExtensions = oFeature.GetArray("stac_extensions");
858 49 : if (!oStacExtensions.IsValid())
859 : {
860 0 : CPLDebug("STACIT",
861 : "Skipping Feature that lacks stac_extensions");
862 0 : continue;
863 : }
864 49 : bool bProjExtensionFound = false;
865 147 : for (const auto &oStacExtension : oStacExtensions)
866 : {
867 245 : if (oStacExtension.ToString() == "proj" ||
868 147 : oStacExtension.ToString().find(
869 : "https://stac-extensions.github.io/projection/") == 0)
870 : {
871 49 : bProjExtensionFound = true;
872 49 : break;
873 : }
874 : }
875 49 : if (!bProjExtensionFound)
876 : {
877 0 : CPLDebug(
878 : "STACIT",
879 : "Skipping Feature that lacks the 'proj' STAC extension");
880 0 : continue;
881 : }
882 :
883 98 : auto jAssets = oFeature["assets"];
884 98 : if (!jAssets.IsValid() ||
885 49 : jAssets.GetType() != CPLJSONObject::Type::Object)
886 : {
887 0 : CPLError(CE_Warning, CPLE_AppDefined,
888 : "Missing assets on a Feature");
889 0 : continue;
890 : }
891 :
892 98 : auto oProperties = oFeature["properties"];
893 98 : if (!oProperties.IsValid() ||
894 49 : oProperties.GetType() != CPLJSONObject::Type::Object)
895 : {
896 0 : CPLError(CE_Warning, CPLE_AppDefined,
897 : "Missing properties on a Feature");
898 0 : continue;
899 : }
900 :
901 98 : const auto osCollection = oFeature["collection"].ToString();
902 64 : if (!osFilteredCollection.empty() &&
903 15 : osFilteredCollection != osCollection)
904 8 : continue;
905 :
906 92 : for (const auto &jAsset : jAssets.GetChildren())
907 : {
908 51 : const auto osAssetName = jAsset.GetName();
909 51 : if (!osFilteredAsset.empty() && osFilteredAsset != osAssetName)
910 8 : continue;
911 :
912 43 : ParseAsset(jAsset, oProperties, osCollection, osFilteredCRS,
913 : oMapCollection);
914 : }
915 : }
916 23 : if (nMaxItems > 0 && nItemIter >= nMaxItems)
917 : {
918 1 : if (!bMaxItemsSpecified)
919 : {
920 0 : CPLError(CE_Warning, CPLE_AppDefined,
921 : "Maximum number of items (" CPL_FRMT_GIB
922 : ") allowed to be retrieved has been hit",
923 : nMaxItems);
924 : }
925 : else
926 : {
927 1 : CPLDebug("STACIT",
928 : "Maximum number of items (" CPL_FRMT_GIB
929 : ") allowed to be retrieved has been hit",
930 : nMaxItems);
931 : }
932 1 : break;
933 : }
934 :
935 : // Follow next link
936 : // Cf https://github.com/radiantearth/stac-api-spec/tree/release/v1.0.0/item-search#pagination
937 44 : const auto oLinks = oRoot.GetArray("links");
938 22 : if (!oLinks.IsValid())
939 19 : break;
940 6 : std::string osNewFilename;
941 6 : for (const auto &oLink : oLinks)
942 : {
943 6 : const auto osType = oLink["type"].ToString();
944 6 : if (oLink["rel"].ToString() == "next" &&
945 3 : (osType.empty() || osType == "application/geo+json"))
946 : {
947 3 : osMethod = oLink.GetString("method", "GET");
948 3 : osNewFilename = oLink["href"].ToString();
949 3 : oHeaders = oLink["headers"];
950 3 : oBody = oLink["body"];
951 3 : bMerge = oLink.GetBool("merge", false);
952 3 : if (osType == "application/geo+json")
953 : {
954 0 : break;
955 : }
956 : }
957 : }
958 6 : if (!osNewFilename.empty() &&
959 3 : (osNewFilename != osCurFilename ||
960 0 : (oBody.IsValid() &&
961 0 : oBody.GetType() == CPLJSONObject::Type::Object)))
962 : {
963 3 : osCurFilename = osNewFilename;
964 : }
965 : else
966 0 : osCurFilename.clear();
967 3 : } while (!osCurFilename.empty());
968 :
969 20 : if (oMapCollection.empty())
970 : {
971 2 : CPLError(CE_Failure, CPLE_AppDefined, "No compatible asset found");
972 2 : return false;
973 : }
974 :
975 35 : if (oMapCollection.size() > 1 ||
976 35 : oMapCollection.begin()->second.assets.size() > 1 ||
977 35 : oMapCollection.begin()->second.assets.begin()->second.assets.size() > 1)
978 : {
979 : // If there's more than one asset type or more than one SRS, expose
980 : // subdatasets.
981 1 : SetSubdatasets(osFilename, oMapCollection);
982 1 : return true;
983 : }
984 : else
985 : {
986 17 : return SetupDataset(poOpenInfo, osFilename, oMapCollection);
987 : }
988 : }
989 :
990 : /************************************************************************/
991 : /* OpenStatic() */
992 : /************************************************************************/
993 :
994 20 : GDALDataset *STACITDataset::OpenStatic(GDALOpenInfo *poOpenInfo)
995 : {
996 20 : if (!Identify(poOpenInfo))
997 0 : return nullptr;
998 40 : auto poDS = std::make_unique<STACITDataset>();
999 20 : if (!poDS->Open(poOpenInfo))
1000 2 : return nullptr;
1001 18 : return poDS.release();
1002 : }
1003 :
1004 : /************************************************************************/
1005 : /* GDALRegister_STACIT() */
1006 : /************************************************************************/
1007 :
1008 1595 : void GDALRegister_STACIT()
1009 :
1010 : {
1011 1595 : if (GDALGetDriverByName("STACIT") != nullptr)
1012 302 : return;
1013 :
1014 1293 : GDALDriver *poDriver = new GDALDriver();
1015 :
1016 1293 : poDriver->SetDescription("STACIT");
1017 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1018 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1019 1293 : "Spatio-Temporal Asset Catalog Items");
1020 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/stacit.html");
1021 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1022 1293 : poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");
1023 1293 : poDriver->SetMetadataItem(
1024 : GDAL_DMD_OPENOPTIONLIST,
1025 : "<OpenOptionList>"
1026 : " <Option name='MAX_ITEMS' type='int' default='1000' "
1027 : "description='Maximum number of items fetched. 0=unlimited'/>"
1028 : " <Option name='COLLECTION' type='string' "
1029 : "description='Name of collection to filter items'/>"
1030 : " <Option name='ASSET' type='string' "
1031 : "description='Name of asset to filter items'/>"
1032 : " <Option name='CRS' type='string' "
1033 : "description='Name of CRS to filter items'/>"
1034 : " <Option name='RESOLUTION' type='string-select' default='AVERAGE' "
1035 : "description='Strategy to use to determine dataset resolution'>"
1036 : " <Value>AVERAGE</Value>"
1037 : " <Value>HIGHEST</Value>"
1038 : " <Value>LOWEST</Value>"
1039 : " </Option>"
1040 : " <Option name='OVERLAP_STRATEGY' type='string-select' "
1041 : "default='REMOVE_IF_NO_NODATA' "
1042 : "description='Strategy to use when some sources are fully "
1043 : "covered by others'>"
1044 : " <Value>REMOVE_IF_NO_NODATA</Value>"
1045 : " <Value>USE_ALL</Value>"
1046 : " <Value>USE_MOST_RECENT</Value>"
1047 : " </Option>"
1048 1293 : "</OpenOptionList>");
1049 :
1050 1293 : poDriver->pfnOpen = STACITDataset::OpenStatic;
1051 1293 : poDriver->pfnIdentify = STACITDataset::Identify;
1052 :
1053 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1054 : }
|
