STACIT: add a OVERLAP_STRATEGY opening option to determine how to …

autotest/gdrivers/data/stacit/overlapping_sources_with_nodata.json

@@ -0,0 +1,103 @@
+{
+  "type": "FeatureCollection",
+  "stac_version": "1.0.0-beta.2",
+  "stac_extensions": [],
+  "features": [
+    {
+      "type": "Feature",
+      "stac_version": "1.0.0-beta.2",
+      "stac_extensions": [
+        "eo",
+        "proj"
+      ],
+      "id": "byte",
+      "geometry": null,
+      "properties": {
+        "datetime": "2021-07-25T00:00:00Z",
+        "proj:epsg": 26711,
+      },
+      "collection": "my_collection",
+      "assets": {
+        "B01": {
+          "title": "Band 1 (coastal)",
+          "type": "image/tiff; application=geotiff; profile=cloud-optimized",
+          "roles": [
+            "data"
+          ],
+          "eo:bands": [
+            {
+              "name": "B01",
+              "common_name": "coastal",
+              "center_wavelength": 0.4439,
+              "full_width_half_max": 0.027
+            }
+          ],
+          "href": "data/byte_nodata_0.tif",
+          "proj:bbox": [
+            440720.000, 3750120.000,
+            441920.000, 3751320.000
+          ],
+          "proj:transform": [
+            60,
+            0,
+            440720,
+            0,
+            -60,
+            3751320,
+            0,
+            0,
+            1
+          ]
+        }
+      }
+    },
+    {
+      "type": "Feature",
+      "stac_version": "1.0.0-beta.2",
+      "stac_extensions": [
+        "eo",
+        "proj"
+      ],
+      "id": "under",
+      "geometry": null,
+      "properties": {
+        "datetime": "2021-07-19T10:57:30Z",
+        "proj:epsg": 26711,
+      },
+      "collection": "my_collection",
+      "assets": {
+        "B01": {
+          "title": "Band 1 (coastal)",
+          "type": "image/tiff; application=geotiff; profile=cloud-optimized",
+          "roles": [
+            "data"
+          ],
+          "eo:bands": [
+            {
+              "name": "B01",
+              "common_name": "coastal",
+              "center_wavelength": 0.4439,
+              "full_width_half_max": 0.027
+            }
+          ],
+          "href": "data/byte.tif",
+          "proj:bbox": [
+            440720.000, 3750120.000,
+            441920.000, 3751320.000
+          ],
+          "proj:transform": [
+            60,
+            0,
+            440720,
+            0,
+            -60,
+            3751320,
+            0,
+            0,
+            1
+          ]
+        }
+      }
+    }
+  ]
+}

autotest/gdrivers/stacit.py

@@ -156,7 +156,7 @@ def test_stacit_overlapping_sources():
 
     # Check that the source covered by another one is not listed
     vrt = ds.GetMetadata("xml:VRT")[0]
-    placement_vrt = """
+    only_one_simple_source = """
     <ColorInterp>Gray</ColorInterp>
     <SimpleSource>
       <SourceFilename relativeToVRT="0">data/byte.tif</SourceFilename>
@@ -166,4 +166,78 @@ def test_stacit_overlapping_sources():
     </SimpleSource>
   </VRTRasterBand>"""
     # print(vrt)
-    assert placement_vrt in vrt
+    assert only_one_simple_source in vrt
+
+    ds = gdal.OpenEx(
+        "data/stacit/overlapping_sources.json",
+        open_options=["OVERLAP_STRATEGY=REMOVE_IF_NO_NODATA"],
+    )
+    assert ds is not None
+    vrt = ds.GetMetadata("xml:VRT")[0]
+    assert only_one_simple_source in vrt
+
+    ds = gdal.OpenEx(
+        "data/stacit/overlapping_sources.json",
+        open_options=["OVERLAP_STRATEGY=USE_MOST_RECENT"],
+    )
+    assert ds is not None
+    vrt = ds.GetMetadata("xml:VRT")[0]
+    assert only_one_simple_source in vrt
+
+    ds = gdal.OpenEx(
+        "data/stacit/overlapping_sources.json",
+        open_options=["OVERLAP_STRATEGY=USE_ALL"],
+    )
+    assert ds is not None
+    assert len(ds.GetFileList()) == 4
+    vrt = ds.GetMetadata("xml:VRT")[0]
+
+
+@pytest.mark.require_geos
+def test_stacit_overlapping_sources_with_nodata():
+
+    ds = gdal.Open("data/stacit/overlapping_sources_with_nodata.json")
+    assert ds is not None
+    assert len(ds.GetFileList()) == 3
+    vrt = ds.GetMetadata("xml:VRT")[0]
+    # print(vrt)
+    two_sources = """<ComplexSource>
+      <SourceFilename relativeToVRT="0">data/byte.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SrcRect xOff="0" yOff="0" xSize="20" ySize="20" />
+      <DstRect xOff="0" yOff="0" xSize="20" ySize="20" />
+      <NODATA>0</NODATA>
+    </ComplexSource>
+    <ComplexSource>
+      <SourceFilename relativeToVRT="0">data/byte_nodata_0.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SrcRect xOff="0" yOff="0" xSize="20" ySize="20" />
+      <DstRect xOff="0" yOff="0" xSize="20" ySize="20" />
+      <NODATA>0</NODATA>
+    </ComplexSource>"""
+    assert two_sources in vrt
+
+    ds = gdal.OpenEx(
+        "data/stacit/overlapping_sources_with_nodata.json",
+        open_options=["OVERLAP_STRATEGY=REMOVE_IF_NO_NODATA"],
+    )
+    assert ds is not None
+    vrt = ds.GetMetadata("xml:VRT")[0]
+    assert len(ds.GetFileList()) == 3
+    assert two_sources in vrt
+
+    ds = gdal.OpenEx(
+        "data/stacit/overlapping_sources_with_nodata.json",
+        open_options=["OVERLAP_STRATEGY=USE_MOST_RECENT"],
+    )
+    assert ds is not None
+    assert len(ds.GetFileList()) == 2
+
+    ds = gdal.OpenEx(
+        "data/stacit/overlapping_sources_with_nodata.json",
+        open_options=["OVERLAP_STRATEGY=USE_ALL"],
+    )
+    assert ds is not None
+    vrt = ds.GetMetadata("xml:VRT")[0]
+    assert len(ds.GetFileList()) == 3
+    assert two_sources in vrt

doc/source/drivers/raster/stacit.rst

@@ -18,10 +18,6 @@ Thus, translating it into VRT will result in a VRT file that directly references
 
 Note that `STAC API ItemCollections <https://github.com/radiantearth/stac-api-spec/blob/main/fragments/itemcollection/README.md>`_ are not the same as  `STAC Collections <https://github.com/radiantearth/stac-spec/tree/master/collection-spec>`_. STAC API ItemCollections are GeoJSON FeatureCollections enhanced with STAC entities.
 
-Note that when the ItemCollections contains overlapping items, and that some items
-are fully covered by other items that are more recent, the STACIT virtual mosaic will
-not list those fully covered items not participating to the pixel values of the mosaic.
-
 Open syntax
 -----------
 
@@ -67,6 +63,28 @@ The following open options are supported:
 
       Strategy to use to determine dataset resolution.
 
+-  .. oo:: OVERLAP_STRATEGY
+      :choices: REMOVE_IF_NO_NODATA, USE_ALL, USE_MOST_RECENT
+      :default: REMOVE_IF_NO_NODATA
+      :since: 3.9.1
+
+      Strategy to use when the ItemCollections contains overlapping items, and
+      that some items are fully covered by other items that are more recent.
+
+      Starting with GDAL 3.9.1, the ``REMOVE_IF_NO_NODATA`` strategy is applied
+      by default. The STACIT virtual mosaic will omit fully covered items,
+      only if no band declares a nodata value.
+      (Note that the determination whether a band has a nodata value of not is
+      done by opening one of the items, and assuming it is representative of
+      the characteristics of the others in the collection).
+
+      This strategy can be forced in all cases by selecting the ``USE_MOST_RECENT``
+      strategy (this was the strategy applied prior to 3.9.1)
+
+      The ``USE_ALL`` strategy always causes all items to be listed in the virtual
+      mosaic, with the most recent ones being rendered on top of the less recent ones.
+
+
 Subdatasets
 -----------
 

frmts/stacit/stacitdataset.cpp

@@ -550,6 +550,7 @@ bool STACITDataset::SetupDataset(
               });
 
     // Create VRT bands and add sources
+    bool bAtLeastOneBandHasNoData = false;
     for (int i = 0; i < poItemDS->GetRasterCount(); i++)
     {
         auto poItemBand = poItemDS->GetRasterBand(i + 1);
@@ -559,7 +560,10 @@ bool STACITDataset::SetupDataset(
         int bHasNoData = FALSE;
         const double dfNoData = poItemBand->GetNoDataValue(&bHasNoData);
         if (bHasNoData)
+        {
+            bAtLeastOneBandHasNoData = true;
             poVRTBand->SetNoDataValue(dfNoData);
+        }
 
         const auto eInterp = poItemBand->GetColorInterpretation();
         if (eInterp != GCI_Undefined)
@@ -627,9 +631,17 @@ bool STACITDataset::SetupDataset(
             }
         }
 
-        const char *apszOptions[] = {"EMIT_ERROR_IF_GEOS_NOT_AVAILABLE=NO",
-                                     nullptr};
-        poVRTBand->RemoveCoveredSources(apszOptions);
+        const char *pszOverlapStrategy =
+            CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
+                                 "OVERLAP_STRATEGY", "REMOVE_IF_NO_NODATA");
+        if ((EQUAL(pszOverlapStrategy, "REMOVE_IF_NO_NODATA") &&
+             !bAtLeastOneBandHasNoData) ||
+            EQUAL(pszOverlapStrategy, "USE_MOST_RECENT"))
+        {
+            const char *const apszOptions[] = {
+                "EMIT_ERROR_IF_GEOS_NOT_AVAILABLE=NO", nullptr};
+            poVRTBand->RemoveCoveredSources(apszOptions);
+        }
     }
     return true;
 }
@@ -940,6 +952,14 @@ void GDALRegister_STACIT()
         "       <Value>HIGHEST</Value>"
         "       <Value>LOWEST</Value>"
         "   </Option>"
+        "   <Option name='OVERLAP_STRATEGY' type='string-select' "
+        "default='REMOVE_IF_NO_NODATA' "
+        "description='Strategy to use when some sources are fully "
+        "covered by others'>"
+        "       <Value>REMOVE_IF_NO_NODATA</Value>"
+        "       <Value>USE_ALL</Value>"
+        "       <Value>USE_MOST_RECENT</Value>"
+        "   </Option>"
         "</OpenOptionList>");
 
     poDriver->pfnOpen = STACITDataset::OpenStatic;