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 51040 : int STACITDataset::Identify(GDALOpenInfo *poOpenInfo)
94 : {
95 51040 : if (STARTS_WITH(poOpenInfo->pszFilename, "STACIT:"))
96 : {
97 12 : return true;
98 : }
99 :
100 51028 : const bool bIsSingleDriver = poOpenInfo->IsSingleAllowedDriver("STACIT");
101 51023 : if (bIsSingleDriver && (STARTS_WITH(poOpenInfo->pszFilename, "http://") ||
102 4 : STARTS_WITH(poOpenInfo->pszFilename, "https://")))
103 : {
104 1 : return true;
105 : }
106 :
107 51022 : if (poOpenInfo->nHeaderBytes == 0)
108 : {
109 48070 : return false;
110 : }
111 :
112 8796 : for (int i = 0; i < 2; i++)
113 : {
114 : // TryToIngest() may reallocate pabyHeader, so do not move this
115 : // before the loop.
116 5874 : const char *pszHeader =
117 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
118 6270 : while (*pszHeader != 0 &&
119 4512 : std::isspace(static_cast<unsigned char>(*pszHeader)))
120 396 : ++pszHeader;
121 5874 : if (bIsSingleDriver)
122 : {
123 4 : return pszHeader[0] == '{';
124 : }
125 :
126 5870 : if (strstr(pszHeader, "\"stac_version\"") != nullptr &&
127 26 : strstr(pszHeader, "\"proj:transform\"") != nullptr)
128 : {
129 26 : return true;
130 : }
131 :
132 5844 : if (i == 0)
133 : {
134 : // Should be enough for a STACIT .json file
135 2922 : poOpenInfo->TryToIngest(32768);
136 : }
137 : }
138 :
139 2922 : 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 : if (dfXRes == 0 || dfYRes == 0)
462 : {
463 0 : CPLError(CE_Failure, CPLE_AppDefined,
464 : "Invalid computed dataset dimensions");
465 0 : return false;
466 : }
467 17 : double dfXSize = std::round((dfXMax - dfXMin) / dfXRes);
468 17 : double dfYSize = std::round((dfYMax - dfYMin) / dfYRes);
469 17 : if (dfXSize <= 0 || dfYSize <= 0 || dfXSize > INT_MAX || dfYSize > INT_MAX)
470 : {
471 0 : CPLError(CE_Failure, CPLE_AppDefined,
472 : "Invalid computed dataset dimensions");
473 0 : return false;
474 : }
475 17 : nRasterXSize = static_cast<int>(dfXSize);
476 17 : nRasterYSize = static_cast<int>(dfYSize);
477 :
478 : // Set geotransform
479 : double adfGeoTransform[6];
480 17 : adfGeoTransform[0] = dfXMin;
481 17 : adfGeoTransform[1] = dfXRes;
482 17 : adfGeoTransform[2] = 0;
483 17 : adfGeoTransform[3] = dfYMax;
484 17 : adfGeoTransform[4] = 0;
485 17 : adfGeoTransform[5] = -dfYRes;
486 17 : SetGeoTransform(adfGeoTransform);
487 :
488 : // Set SRS
489 34 : OGRSpatialReference oSRS;
490 17 : if (oSRS.SetFromUserInput(
491 : assetByProj.osProjUserString.c_str(),
492 17 : OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get()) ==
493 : OGRERR_NONE)
494 : {
495 17 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
496 17 : SetSpatialRef(&oSRS);
497 : }
498 :
499 : // Open of the items to find the number of bands, their data type
500 : // and nodata.
501 : const auto BuildVSICurlFilename =
502 49 : [&osSTACITFilename, &collection](const std::string &osFilename)
503 : {
504 49 : std::string osRet;
505 49 : if (STARTS_WITH(osFilename.c_str(), "http"))
506 : {
507 0 : if (STARTS_WITH(osSTACITFilename.c_str(),
508 : "https://planetarycomputer.microsoft.com/api/"))
509 : {
510 0 : osRet = "/vsicurl?pc_url_signing=yes&pc_collection=";
511 0 : osRet += collection.osName;
512 0 : osRet += "&url=";
513 : char *pszEncoded =
514 0 : CPLEscapeString(osFilename.c_str(), -1, CPLES_URL);
515 0 : CPLString osEncoded(pszEncoded);
516 0 : CPLFree(pszEncoded);
517 : // something is confused if the whole URL appears as
518 : // /vsicurl...blabla_without_slash.tif" so transform back %2F to
519 : // /
520 0 : osEncoded.replaceAll("%2F", '/');
521 0 : osRet += osEncoded;
522 : }
523 : else
524 : {
525 0 : osRet = "/vsicurl/";
526 0 : osRet += osFilename;
527 : }
528 : }
529 49 : else if (STARTS_WITH(osFilename.c_str(), "file://"))
530 : {
531 3 : osRet = osFilename.substr(strlen("file://"));
532 : }
533 46 : else if (STARTS_WITH(osFilename.c_str(), "s3://"))
534 : {
535 0 : osRet = "/vsis3/";
536 0 : osRet += osFilename.substr(strlen("s3://"));
537 : }
538 :
539 : else
540 : {
541 46 : osRet = osFilename;
542 : }
543 49 : return osRet;
544 17 : };
545 :
546 34 : const auto osFirstItemName(BuildVSICurlFilename(items.front().osFilename));
547 : auto poItemDS = std::unique_ptr<GDALDataset>(
548 34 : GDALDataset::Open(osFirstItemName.c_str()));
549 17 : if (!poItemDS)
550 : {
551 0 : CPLError(CE_Failure, CPLE_AppDefined,
552 : "Cannot open %s to retrieve band characteristics",
553 : osFirstItemName.c_str());
554 0 : return false;
555 : }
556 :
557 : // Sort by ascending datetime
558 17 : std::sort(items.begin(), items.end(),
559 26 : [](const AssetItem &a, const AssetItem &b)
560 : {
561 26 : if (!a.osDateTime.empty() && !b.osDateTime.empty())
562 26 : return a.osDateTime < b.osDateTime;
563 0 : return &a < &b;
564 : });
565 :
566 : // Create VRT bands and add sources
567 17 : bool bAtLeastOneBandHasNoData = false;
568 34 : for (int i = 0; i < poItemDS->GetRasterCount(); i++)
569 : {
570 17 : auto poItemBand = poItemDS->GetRasterBand(i + 1);
571 17 : AddBand(poItemBand->GetRasterDataType(), nullptr);
572 : auto poVRTBand =
573 17 : cpl::down_cast<VRTSourcedRasterBand *>(GetRasterBand(i + 1));
574 17 : int bHasNoData = FALSE;
575 17 : const double dfNoData = poItemBand->GetNoDataValue(&bHasNoData);
576 17 : if (bHasNoData)
577 : {
578 4 : bAtLeastOneBandHasNoData = true;
579 4 : poVRTBand->SetNoDataValue(dfNoData);
580 : }
581 :
582 17 : const auto eInterp = poItemBand->GetColorInterpretation();
583 17 : if (eInterp != GCI_Undefined)
584 17 : poVRTBand->SetColorInterpretation(eInterp);
585 :
586 : // Set band properties
587 34 : if (asset.eoBands.IsValid() &&
588 17 : asset.eoBands.Size() == poItemDS->GetRasterCount())
589 : {
590 34 : const auto &eoBand = asset.eoBands[i];
591 51 : const auto osBandName = eoBand["name"].ToString();
592 17 : if (!osBandName.empty())
593 17 : poVRTBand->SetDescription(osBandName.c_str());
594 :
595 51 : const auto osCommonName = eoBand["common_name"].ToString();
596 17 : if (!osCommonName.empty())
597 : {
598 : const auto eInterpFromCommonName =
599 13 : GDALGetColorInterpFromSTACCommonName(osCommonName.c_str());
600 13 : if (eInterpFromCommonName != GCI_Undefined)
601 13 : poVRTBand->SetColorInterpretation(eInterpFromCommonName);
602 : }
603 :
604 73 : for (const auto &eoBandChild : eoBand.GetChildren())
605 : {
606 112 : const auto osChildName = eoBandChild.GetName();
607 56 : if (osChildName != "name" && osChildName != "common_name")
608 : {
609 26 : poVRTBand->SetMetadataItem(osChildName.c_str(),
610 52 : eoBandChild.ToString().c_str());
611 : }
612 : }
613 : }
614 :
615 : // Add items as VRT sources
616 49 : for (const auto &assetItem : items)
617 : {
618 64 : const auto osItemName(BuildVSICurlFilename(assetItem.osFilename));
619 32 : const double dfDstXOff = (assetItem.dfXMin - dfXMin) / dfXRes;
620 32 : const double dfDstXSize =
621 32 : (assetItem.dfXMax - assetItem.dfXMin) / dfXRes;
622 32 : const double dfDstYOff = (dfYMax - assetItem.dfYMax) / dfYRes;
623 32 : const double dfDstYSize =
624 32 : (assetItem.dfYMax - assetItem.dfYMin) / dfYRes;
625 32 : if (!bHasNoData)
626 : {
627 24 : poVRTBand->AddSimpleSource(osItemName.c_str(), i + 1, 0, 0,
628 24 : assetItem.nXSize, assetItem.nYSize,
629 : dfDstXOff, dfDstYOff, dfDstXSize,
630 : dfDstYSize);
631 : }
632 : else
633 : {
634 8 : poVRTBand->AddComplexSource(osItemName.c_str(), i + 1, 0, 0,
635 8 : assetItem.nXSize, assetItem.nYSize,
636 : dfDstXOff, dfDstYOff, dfDstXSize,
637 : dfDstYSize,
638 : 0.0, // offset
639 : 1.0, // scale
640 : dfNoData);
641 : }
642 : }
643 :
644 : const char *pszOverlapStrategy =
645 17 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
646 : "OVERLAP_STRATEGY", "REMOVE_IF_NO_NODATA");
647 17 : if ((EQUAL(pszOverlapStrategy, "REMOVE_IF_NO_NODATA") &&
648 13 : !bAtLeastOneBandHasNoData) ||
649 6 : EQUAL(pszOverlapStrategy, "USE_MOST_RECENT"))
650 : {
651 13 : const char *const apszOptions[] = {
652 : "EMIT_ERROR_IF_GEOS_NOT_AVAILABLE=NO", nullptr};
653 13 : poVRTBand->RemoveCoveredSources(apszOptions);
654 : }
655 : }
656 17 : return true;
657 : }
658 :
659 : /************************************************************************/
660 : /* SetSubdatasets() */
661 : /************************************************************************/
662 :
663 1 : void STACITDataset::SetSubdatasets(
664 : const std::string &osFilename,
665 : const std::map<std::string, Collection> &oMapCollection)
666 : {
667 2 : CPLStringList aosSubdatasets;
668 1 : int nCount = 1;
669 3 : for (const auto &collectionKV : oMapCollection)
670 : {
671 5 : for (const auto &assetKV : collectionKV.second.assets)
672 : {
673 6 : std::string osCollectionAssetArg;
674 3 : if (oMapCollection.size() > 1)
675 : osCollectionAssetArg +=
676 3 : "collection=" + collectionKV.first + ",";
677 3 : osCollectionAssetArg += "asset=" + assetKV.first;
678 :
679 6 : std::string osCollectionAssetText;
680 3 : if (oMapCollection.size() > 1)
681 : osCollectionAssetText +=
682 3 : "Collection " + collectionKV.first + ", ";
683 3 : osCollectionAssetText += "Asset " + assetKV.first;
684 :
685 3 : if (assetKV.second.assets.size() == 1)
686 : {
687 : aosSubdatasets.AddString(CPLSPrintf(
688 : "SUBDATASET_%d_NAME=STACIT:\"%s\":%s", nCount,
689 2 : osFilename.c_str(), osCollectionAssetArg.c_str()));
690 : aosSubdatasets.AddString(CPLSPrintf(
691 : "SUBDATASET_%d_DESC=%s of %s", nCount,
692 2 : osCollectionAssetText.c_str(), osFilename.c_str()));
693 2 : nCount++;
694 : }
695 : else
696 : {
697 3 : for (const auto &assetsByProjKV : assetKV.second.assets)
698 : {
699 : aosSubdatasets.AddString(CPLSPrintf(
700 : "SUBDATASET_%d_NAME=STACIT:\"%s\":%s,crs=%s", nCount,
701 : osFilename.c_str(), osCollectionAssetArg.c_str(),
702 2 : SanitizeCRSValue(assetsByProjKV.first).c_str()));
703 : aosSubdatasets.AddString(CPLSPrintf(
704 : "SUBDATASET_%d_DESC=%s of %s in CRS %s", nCount,
705 : osCollectionAssetText.c_str(), osFilename.c_str(),
706 2 : assetsByProjKV.first.c_str()));
707 2 : nCount++;
708 : }
709 : }
710 : }
711 : }
712 1 : GDALDataset::SetMetadata(aosSubdatasets.List(), "SUBDATASETS");
713 1 : }
714 :
715 : /************************************************************************/
716 : /* Open() */
717 : /************************************************************************/
718 :
719 20 : bool STACITDataset::Open(GDALOpenInfo *poOpenInfo)
720 : {
721 40 : CPLString osFilename(poOpenInfo->pszFilename);
722 : std::string osFilteredCollection =
723 40 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "COLLECTION", "");
724 : std::string osFilteredAsset =
725 40 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "ASSET", "");
726 : std::string osFilteredCRS = SanitizeCRSValue(
727 60 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "CRS", ""));
728 20 : if (STARTS_WITH(poOpenInfo->pszFilename, "STACIT:"))
729 : {
730 : const CPLStringList aosTokens(CSLTokenizeString2(
731 6 : poOpenInfo->pszFilename, ":", CSLT_HONOURSTRINGS));
732 6 : if (aosTokens.size() != 2 && aosTokens.size() != 3)
733 0 : return false;
734 6 : osFilename = aosTokens[1];
735 6 : if (aosTokens.size() >= 3)
736 : {
737 : const CPLStringList aosFilters(
738 12 : CSLTokenizeString2(aosTokens[2], ",", 0));
739 : osFilteredCollection = aosFilters.FetchNameValueDef(
740 6 : "collection", osFilteredCollection.c_str());
741 : osFilteredAsset =
742 6 : aosFilters.FetchNameValueDef("asset", osFilteredAsset.c_str());
743 : osFilteredCRS =
744 6 : aosFilters.FetchNameValueDef("crs", osFilteredCRS.c_str());
745 : }
746 : }
747 :
748 40 : std::map<std::string, Collection> oMapCollection;
749 20 : GIntBig nItemIter = 0;
750 20 : GIntBig nMaxItems = CPLAtoGIntBig(CSLFetchNameValueDef(
751 20 : poOpenInfo->papszOpenOptions, "MAX_ITEMS", "1000"));
752 :
753 : const bool bMaxItemsSpecified =
754 20 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, "MAX_ITEMS") != nullptr;
755 20 : if (!bMaxItemsSpecified)
756 : {
757 : // If the URL includes a limit parameter, and it's larger than our
758 : // default MAX_ITEMS value, then increase the later to the former.
759 19 : const auto nPos = osFilename.ifind("&limit=");
760 19 : if (nPos != std::string::npos)
761 : {
762 0 : const auto nLimit = CPLAtoGIntBig(
763 0 : osFilename.substr(nPos + strlen("&limit=")).c_str());
764 0 : nMaxItems = std::max(nMaxItems, nLimit);
765 : }
766 : }
767 :
768 40 : auto osCurFilename = osFilename;
769 40 : std::string osMethod = "GET";
770 40 : CPLJSONObject oHeaders;
771 40 : CPLJSONObject oBody;
772 20 : bool bMerge = false;
773 20 : int nLoops = 0;
774 3 : do
775 : {
776 23 : ++nLoops;
777 23 : if (nMaxItems > 0 && nLoops > nMaxItems)
778 : {
779 20 : break;
780 : }
781 :
782 23 : CPLJSONDocument oDoc;
783 :
784 45 : if (STARTS_WITH(osCurFilename, "http://") ||
785 22 : STARTS_WITH(osCurFilename, "https://"))
786 : {
787 : // Cf // Cf https://github.com/radiantearth/stac-api-spec/tree/release/v1.0.0/item-search#pagination
788 1 : CPLStringList aosOptions;
789 2 : if (oBody.IsValid() &&
790 1 : oBody.GetType() == CPLJSONObject::Type::Object)
791 : {
792 1 : if (bMerge)
793 0 : CPLDebug("STACIT",
794 : "Ignoring 'merge' attribute from next link");
795 : const std::string osPostContent =
796 2 : oBody.Format(CPLJSONObject::PrettyFormat::Pretty);
797 1 : aosOptions.SetNameValue("POSTFIELDS", osPostContent.c_str());
798 : }
799 1 : aosOptions.SetNameValue("CUSTOMREQUEST", osMethod.c_str());
800 1 : CPLString osHeaders;
801 4 : if (!oHeaders.IsValid() ||
802 2 : oHeaders.GetType() != CPLJSONObject::Type::Object ||
803 2 : oHeaders["Content-Type"].ToString().empty())
804 : {
805 1 : osHeaders = "Content-Type: application/json";
806 : }
807 2 : if (oHeaders.IsValid() &&
808 1 : oHeaders.GetType() == CPLJSONObject::Type::Object)
809 : {
810 2 : for (const auto &obj : oHeaders.GetChildren())
811 : {
812 1 : osHeaders += "\r\n";
813 1 : osHeaders += obj.GetName();
814 1 : osHeaders += ": ";
815 1 : osHeaders += obj.ToString();
816 : }
817 : }
818 1 : aosOptions.SetNameValue("HEADERS", osHeaders.c_str());
819 : CPLHTTPResult *psResult =
820 1 : CPLHTTPFetch(osCurFilename.c_str(), aosOptions.List());
821 1 : if (!psResult)
822 0 : return false;
823 1 : if (!psResult->pabyData)
824 : {
825 0 : CPLHTTPDestroyResult(psResult);
826 0 : return false;
827 : }
828 2 : const bool bOK = oDoc.LoadMemory(
829 1 : reinterpret_cast<const char *>(psResult->pabyData));
830 : // CPLDebug("STACIT", "Response: %s", reinterpret_cast<const char*>(psResult->pabyData));
831 1 : CPLHTTPDestroyResult(psResult);
832 1 : if (!bOK)
833 0 : return false;
834 : }
835 : else
836 : {
837 22 : if (!oDoc.Load(osCurFilename))
838 0 : return false;
839 : }
840 23 : const auto oRoot = oDoc.GetRoot();
841 46 : auto oFeatures = oRoot.GetArray("features");
842 23 : if (!oFeatures.IsValid())
843 : {
844 1 : if (oRoot.GetString("type") == "Feature")
845 : {
846 1 : oFeatures = CPLJSONArray();
847 1 : oFeatures.Add(oRoot);
848 : }
849 : else
850 : {
851 0 : CPLError(CE_Failure, CPLE_AppDefined, "Missing features");
852 0 : return false;
853 : }
854 : }
855 72 : for (const auto &oFeature : oFeatures)
856 : {
857 49 : nItemIter++;
858 49 : if (nMaxItems > 0 && nItemIter > nMaxItems)
859 : {
860 0 : break;
861 : }
862 :
863 98 : auto oStacExtensions = oFeature.GetArray("stac_extensions");
864 49 : if (!oStacExtensions.IsValid())
865 : {
866 0 : CPLDebug("STACIT",
867 : "Skipping Feature that lacks stac_extensions");
868 0 : continue;
869 : }
870 49 : bool bProjExtensionFound = false;
871 147 : for (const auto &oStacExtension : oStacExtensions)
872 : {
873 245 : if (oStacExtension.ToString() == "proj" ||
874 147 : oStacExtension.ToString().find(
875 : "https://stac-extensions.github.io/projection/") == 0)
876 : {
877 49 : bProjExtensionFound = true;
878 49 : break;
879 : }
880 : }
881 49 : if (!bProjExtensionFound)
882 : {
883 0 : CPLDebug(
884 : "STACIT",
885 : "Skipping Feature that lacks the 'proj' STAC extension");
886 0 : continue;
887 : }
888 :
889 98 : auto jAssets = oFeature["assets"];
890 98 : if (!jAssets.IsValid() ||
891 49 : jAssets.GetType() != CPLJSONObject::Type::Object)
892 : {
893 0 : CPLError(CE_Warning, CPLE_AppDefined,
894 : "Missing assets on a Feature");
895 0 : continue;
896 : }
897 :
898 98 : auto oProperties = oFeature["properties"];
899 98 : if (!oProperties.IsValid() ||
900 49 : oProperties.GetType() != CPLJSONObject::Type::Object)
901 : {
902 0 : CPLError(CE_Warning, CPLE_AppDefined,
903 : "Missing properties on a Feature");
904 0 : continue;
905 : }
906 :
907 98 : const auto osCollection = oFeature["collection"].ToString();
908 64 : if (!osFilteredCollection.empty() &&
909 15 : osFilteredCollection != osCollection)
910 8 : continue;
911 :
912 92 : for (const auto &jAsset : jAssets.GetChildren())
913 : {
914 51 : const auto osAssetName = jAsset.GetName();
915 51 : if (!osFilteredAsset.empty() && osFilteredAsset != osAssetName)
916 8 : continue;
917 :
918 43 : ParseAsset(jAsset, oProperties, osCollection, osFilteredCRS,
919 : oMapCollection);
920 : }
921 : }
922 23 : if (nMaxItems > 0 && nItemIter >= nMaxItems)
923 : {
924 1 : if (!bMaxItemsSpecified)
925 : {
926 0 : CPLError(CE_Warning, CPLE_AppDefined,
927 : "Maximum number of items (" CPL_FRMT_GIB
928 : ") allowed to be retrieved has been hit",
929 : nMaxItems);
930 : }
931 : else
932 : {
933 1 : CPLDebug("STACIT",
934 : "Maximum number of items (" CPL_FRMT_GIB
935 : ") allowed to be retrieved has been hit",
936 : nMaxItems);
937 : }
938 1 : break;
939 : }
940 :
941 : // Follow next link
942 : // Cf https://github.com/radiantearth/stac-api-spec/tree/release/v1.0.0/item-search#pagination
943 44 : const auto oLinks = oRoot.GetArray("links");
944 22 : if (!oLinks.IsValid())
945 19 : break;
946 6 : std::string osNewFilename;
947 6 : for (const auto &oLink : oLinks)
948 : {
949 6 : const auto osType = oLink["type"].ToString();
950 6 : if (oLink["rel"].ToString() == "next" &&
951 3 : (osType.empty() || osType == "application/geo+json"))
952 : {
953 3 : osMethod = oLink.GetString("method", "GET");
954 3 : osNewFilename = oLink["href"].ToString();
955 3 : oHeaders = oLink["headers"];
956 3 : oBody = oLink["body"];
957 3 : bMerge = oLink.GetBool("merge", false);
958 3 : if (osType == "application/geo+json")
959 : {
960 0 : break;
961 : }
962 : }
963 : }
964 6 : if (!osNewFilename.empty() &&
965 3 : (osNewFilename != osCurFilename ||
966 0 : (oBody.IsValid() &&
967 0 : oBody.GetType() == CPLJSONObject::Type::Object)))
968 : {
969 3 : osCurFilename = osNewFilename;
970 : }
971 : else
972 0 : osCurFilename.clear();
973 3 : } while (!osCurFilename.empty());
974 :
975 20 : if (oMapCollection.empty())
976 : {
977 2 : CPLError(CE_Failure, CPLE_AppDefined, "No compatible asset found");
978 2 : return false;
979 : }
980 :
981 35 : if (oMapCollection.size() > 1 ||
982 35 : oMapCollection.begin()->second.assets.size() > 1 ||
983 35 : oMapCollection.begin()->second.assets.begin()->second.assets.size() > 1)
984 : {
985 : // If there's more than one asset type or more than one SRS, expose
986 : // subdatasets.
987 1 : SetSubdatasets(osFilename, oMapCollection);
988 1 : return true;
989 : }
990 : else
991 : {
992 17 : return SetupDataset(poOpenInfo, osFilename, oMapCollection);
993 : }
994 : }
995 :
996 : /************************************************************************/
997 : /* OpenStatic() */
998 : /************************************************************************/
999 :
1000 20 : GDALDataset *STACITDataset::OpenStatic(GDALOpenInfo *poOpenInfo)
1001 : {
1002 20 : if (!Identify(poOpenInfo))
1003 0 : return nullptr;
1004 40 : auto poDS = std::make_unique<STACITDataset>();
1005 20 : if (!poDS->Open(poOpenInfo))
1006 2 : return nullptr;
1007 18 : return poDS.release();
1008 : }
1009 :
1010 : /************************************************************************/
1011 : /* GDALRegister_STACIT() */
1012 : /************************************************************************/
1013 :
1014 1682 : void GDALRegister_STACIT()
1015 :
1016 : {
1017 1682 : if (GDALGetDriverByName("STACIT") != nullptr)
1018 301 : return;
1019 :
1020 1381 : GDALDriver *poDriver = new GDALDriver();
1021 :
1022 1381 : poDriver->SetDescription("STACIT");
1023 1381 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1024 1381 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1025 1381 : "Spatio-Temporal Asset Catalog Items");
1026 1381 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/stacit.html");
1027 1381 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1028 1381 : poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");
1029 1381 : poDriver->SetMetadataItem(
1030 : GDAL_DMD_OPENOPTIONLIST,
1031 : "<OpenOptionList>"
1032 : " <Option name='MAX_ITEMS' type='int' default='1000' "
1033 : "description='Maximum number of items fetched. 0=unlimited'/>"
1034 : " <Option name='COLLECTION' type='string' "
1035 : "description='Name of collection to filter items'/>"
1036 : " <Option name='ASSET' type='string' "
1037 : "description='Name of asset to filter items'/>"
1038 : " <Option name='CRS' type='string' "
1039 : "description='Name of CRS to filter items'/>"
1040 : " <Option name='RESOLUTION' type='string-select' default='AVERAGE' "
1041 : "description='Strategy to use to determine dataset resolution'>"
1042 : " <Value>AVERAGE</Value>"
1043 : " <Value>HIGHEST</Value>"
1044 : " <Value>LOWEST</Value>"
1045 : " </Option>"
1046 : " <Option name='OVERLAP_STRATEGY' type='string-select' "
1047 : "default='REMOVE_IF_NO_NODATA' "
1048 : "description='Strategy to use when some sources are fully "
1049 : "covered by others'>"
1050 : " <Value>REMOVE_IF_NO_NODATA</Value>"
1051 : " <Value>USE_ALL</Value>"
1052 : " <Value>USE_MOST_RECENT</Value>"
1053 : " </Option>"
1054 1381 : "</OpenOptionList>");
1055 :
1056 1381 : poDriver->pfnOpen = STACITDataset::OpenStatic;
1057 1381 : poDriver->pfnIdentify = STACITDataset::Identify;
1058 :
1059 1381 : GetGDALDriverManager()->RegisterDriver(poDriver);
1060 : }
|
