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