VRT: add a <UseMaskBand>true</UseMaskBand> element of <ComplexSource>…

… to allow compositing of overlapping sources (fixes OSGeo#1148)

autotest/gcore/vrt_read.py

@@ -1529,3 +1529,118 @@ def test_vrt_implicit_ovr_with_hidenodatavalue():
     assert got_data[32*64+32] == 10
 
     gdal.Unlink('/vsimem/src.tif')
+
+
+def test_vrt_usemaskband():
+
+    ds = gdal.GetDriverByName('GTiff').Create('/vsimem/src1.tif', 3, 1)
+    ds.GetRasterBand(1).Fill(255)
+    ds.CreateMaskBand(0)
+    ds.GetRasterBand(1).GetMaskBand().WriteRaster(0, 0, 1, 1, b'\xff')
+    ds = None
+
+    ds = gdal.GetDriverByName('GTiff').Create('/vsimem/src2.tif', 3, 1)
+    ds.GetRasterBand(1).Fill(127)
+    ds.CreateMaskBand(0)
+    ds.GetRasterBand(1).GetMaskBand().WriteRaster(1, 0, 1, 1, b'\xff')
+    ds = None
+
+    vrt_text = """<VRTDataset rasterXSize="3" rasterYSize="1">
+  <VRTRasterBand dataType="Byte" band="1">
+    <ComplexSource>
+      <SourceFilename>/vsimem/src1.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SrcRect xOff="0" yOff="0" xSize="3" ySize="1" />
+      <DstRect xOff="0" yOff="0" xSize="3" ySize="1" />
+      <UseMaskBand>true</UseMaskBand>
+    </ComplexSource>
+    <ComplexSource>
+      <SourceFilename>/vsimem/src2.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SrcRect xOff="0" yOff="0" xSize="3" ySize="1" />
+      <DstRect xOff="0" yOff="0" xSize="3" ySize="1" />
+      <UseMaskBand>true</UseMaskBand>
+    </ComplexSource>
+  </VRTRasterBand>
+  <MaskBand>
+    <VRTRasterBand dataType="Byte">
+        <ComplexSource>
+            <SourceFilename>/vsimem/src1.tif</SourceFilename>
+            <SourceBand>mask,1</SourceBand>
+            <SrcRect xOff="0" yOff="0" xSize="3" ySize="1" />
+            <DstRect xOff="0" yOff="0" xSize="3" ySize="1" />
+            <UseMaskBand>true</UseMaskBand>
+        </ComplexSource>
+        <ComplexSource>
+            <SourceFilename>/vsimem/src2.tif</SourceFilename>
+            <SourceBand>mask,1</SourceBand>
+            <SrcRect xOff="0" yOff="0" xSize="3" ySize="1" />
+            <DstRect xOff="0" yOff="0" xSize="3" ySize="1" />
+            <UseMaskBand>true</UseMaskBand>
+        </ComplexSource>
+    </VRTRasterBand>
+  </MaskBand>
+</VRTDataset>"""
+    ds = gdal.Open(vrt_text)
+    assert struct.unpack('B' * 3, ds.ReadRaster()) == (255, 127, 0)
+    assert struct.unpack('B' * 3, ds.GetRasterBand(1).GetMaskBand().ReadRaster()) == (255, 255, 0)
+
+    gdal.GetDriverByName('GTiff').Delete('/vsimem/src1.tif')
+    gdal.GetDriverByName('GTiff').Delete('/vsimem/src2.tif')
+
+
+def test_vrt_usemaskband_alpha():
+
+    ds = gdal.GetDriverByName('GTiff').Create('/vsimem/src1.tif', 3, 1, 2)
+    ds.GetRasterBand(1).Fill(255)
+    ds.GetRasterBand(1).GetMaskBand().WriteRaster(0, 0, 1, 1, b'\xff')
+    ds.GetRasterBand(2).SetColorInterpretation(gdal.GCI_AlphaBand)
+    ds.GetRasterBand(2).WriteRaster(0, 0, 1, 1, b'\xff')
+
+    ds = gdal.GetDriverByName('GTiff').Create('/vsimem/src2.tif', 3, 1, 2)
+    ds.GetRasterBand(1).Fill(127)
+    ds.GetRasterBand(2).SetColorInterpretation(gdal.GCI_AlphaBand)
+    ds.GetRasterBand(2).WriteRaster(1, 0, 1, 1, b'\xff')
+    ds = None
+
+    vrt_text = """<VRTDataset rasterXSize="3" rasterYSize="1">
+  <VRTRasterBand dataType="Byte" band="1">
+    <ComplexSource>
+      <SourceFilename>/vsimem/src1.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SrcRect xOff="0" yOff="0" xSize="3" ySize="1" />
+      <DstRect xOff="0" yOff="0" xSize="3" ySize="1" />
+      <UseMaskBand>true</UseMaskBand>
+    </ComplexSource>
+    <ComplexSource>
+      <SourceFilename>/vsimem/src2.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SrcRect xOff="0" yOff="0" xSize="3" ySize="1" />
+      <DstRect xOff="0" yOff="0" xSize="3" ySize="1" />
+      <UseMaskBand>true</UseMaskBand>
+    </ComplexSource>
+  </VRTRasterBand>
+  <VRTRasterBand dataType="Byte" band="2">
+    <ColorInterp>Alpha</ColorInterp>
+    <ComplexSource>
+        <SourceFilename>/vsimem/src1.tif</SourceFilename>
+        <SourceBand>2</SourceBand>
+        <SrcRect xOff="0" yOff="0" xSize="3" ySize="1" />
+        <DstRect xOff="0" yOff="0" xSize="3" ySize="1" />
+        <UseMaskBand>true</UseMaskBand>
+    </ComplexSource>
+    <ComplexSource>
+        <SourceFilename>/vsimem/src2.tif</SourceFilename>
+        <SourceBand>2</SourceBand>
+        <SrcRect xOff="0" yOff="0" xSize="3" ySize="1" />
+        <DstRect xOff="0" yOff="0" xSize="3" ySize="1" />
+        <UseMaskBand>true</UseMaskBand>
+    </ComplexSource>
+  </VRTRasterBand>
+</VRTDataset>"""
+    ds = gdal.Open(vrt_text)
+    assert struct.unpack('B' * 3, ds.GetRasterBand(1).ReadRaster()) == (255, 127, 0)
+    assert struct.unpack('B' * 3, ds.GetRasterBand(2).ReadRaster()) == (255, 255, 0)
+
+    gdal.GetDriverByName('GTiff').Delete('/vsimem/src1.tif')
+    gdal.GetDriverByName('GTiff').Delete('/vsimem/src2.tif')

autotest/utilities/test_gdalbuildvrt_lib.py

@@ -277,3 +277,43 @@ def test_gdalbuildvrt_lib_separate_nodata_4():
 
     assert b'<NoDataValue>' not in data
     assert b'<NODATA>' not in data
+
+
+###############################################################################
+def test_gdalbuildvrt_lib_usemaskband_on_mask_band():
+
+    src1_ds = gdal.GetDriverByName('MEM').Create('src1', 3, 1)
+    src1_ds.SetGeoTransform([2,0.001,0,49,0,-0.001])
+    src1_ds.GetRasterBand(1).Fill(255)
+    src1_ds.CreateMaskBand(0)
+    src1_ds.GetRasterBand(1).GetMaskBand().WriteRaster(0, 0, 1, 1, b'\xff')
+
+    src2_ds = gdal.GetDriverByName('MEM').Create('src2', 3, 1)
+    src2_ds.SetGeoTransform([2,0.001,0,49,0,-0.001])
+    src2_ds.GetRasterBand(1).Fill(127)
+    src2_ds.CreateMaskBand(0)
+    src2_ds.GetRasterBand(1).GetMaskBand().WriteRaster(1, 0, 1, 1, b'\xff')
+
+    ds = gdal.BuildVRT('', [src1_ds, src2_ds])
+    assert struct.unpack('B' * 3, ds.ReadRaster()) == (255, 127, 0)
+    assert struct.unpack('B' * 3, ds.GetRasterBand(1).GetMaskBand().ReadRaster()) == (255, 255, 0)
+
+
+###############################################################################
+def test_gdalbuildvrt_lib_usemaskband_on_alpha_band():
+
+    src1_ds = gdal.GetDriverByName('MEM').Create('src1', 3, 1, 2)
+    src1_ds.SetGeoTransform([2,0.001,0,49,0,-0.001])
+    src1_ds.GetRasterBand(1).Fill(255)
+    src1_ds.GetRasterBand(2).SetColorInterpretation(gdal.GCI_AlphaBand)
+    src1_ds.GetRasterBand(2).WriteRaster(0, 0, 1, 1, b'\xff')
+
+    src2_ds = gdal.GetDriverByName('MEM').Create('src2', 3, 1, 2)
+    src2_ds.SetGeoTransform([2,0.001,0,49,0,-0.001])
+    src2_ds.GetRasterBand(1).Fill(127)
+    src2_ds.GetRasterBand(2).SetColorInterpretation(gdal.GCI_AlphaBand)
+    src2_ds.GetRasterBand(2).WriteRaster(1, 0, 1, 1, b'\xff')
+
+    ds = gdal.BuildVRT('', [src1_ds, src2_ds])
+    assert struct.unpack('B' * 3, ds.GetRasterBand(1).ReadRaster()) == (255, 127, 0)
+    assert struct.unpack('B' * 3, ds.GetRasterBand(2).ReadRaster()) == (255, 255, 0)

gdal/apps/gdalbuildvrt_bin.cpp

@@ -51,6 +51,7 @@ static void Usage(const char* pszErrorMsg)
             "                    [-allow_projection_difference] [-q]\n"
             "                    [-addalpha] [-hidenodata]\n"
             "                    [-srcnodata \"value [value...]\"] [-vrtnodata \"value [value...]\"] \n"
+            "                    [-ignore_srcmaskband]\n"
             "                    [-a_srs srs_def]\n"
             "                    [-r {nearest,bilinear,cubic,cubicspline,lanczos,average,mode}]\n"
             "                    [-oo NAME=VALUE]*\n"

gdal/apps/gdalbuildvrt_lib.cpp

@@ -90,6 +90,7 @@ struct DatasetProperty
     std::vector<bool>   abHasOffset{};
     std::vector<double> adfOffset{};
     std::vector<bool>   abHasScale{};
+    std::vector<bool>   abHasMaskBand{};
     std::vector<double> adfScale{};
     int    bHasDatasetMask = 0;
     int    nMaskBlockXSize = 0;
@@ -240,6 +241,7 @@ class VRTBuilder
     char               *pszOutputSRS = nullptr;
     char               *pszResampling = nullptr;
     char              **papszOpenOptions = nullptr;
+    bool                bUseSrcMaskBand = true;
 
     /* Internal variables */
     char               *pszProjectionRef = nullptr;
@@ -278,6 +280,7 @@ class VRTBuilder
                            int bSeparate, int bAllowProjectionDifference,
                            int bAddAlpha, int bHideNoData, int nSubdataset,
                            const char* pszSrcNoData, const char* pszVRTNoData,
+                           bool bUseSrcMaskBand,
                            const char* pszOutputSRS,
                            const char* pszResampling,
                            const char* const* papszOpenOptionsIn );
@@ -303,6 +306,7 @@ VRTBuilder::VRTBuilder(const char* pszOutputFilenameIn,
                        int bSeparateIn, int bAllowProjectionDifferenceIn,
                        int bAddAlphaIn, int bHideNoDataIn, int nSubdatasetIn,
                        const char* pszSrcNoDataIn, const char* pszVRTNoDataIn,
+                       bool bUseSrcMaskBandIn,
                        const char* pszOutputSRSIn,
                        const char* pszResamplingIn,
                        const char* const * papszOpenOptionsIn )
@@ -359,6 +363,7 @@ VRTBuilder::VRTBuilder(const char* pszOutputFilenameIn,
     pszVRTNoData = (pszVRTNoDataIn) ? CPLStrdup(pszVRTNoDataIn) : nullptr;
     pszOutputSRS = (pszOutputSRSIn) ? CPLStrdup(pszOutputSRSIn) : nullptr;
     pszResampling = (pszResamplingIn) ? CPLStrdup(pszResamplingIn) : nullptr;
+    bUseSrcMaskBand = bUseSrcMaskBandIn;
 }
 
 /************************************************************************/
@@ -606,6 +611,8 @@ int VRTBuilder::AnalyseRaster( GDALDatasetH hDS, DatasetProperty* psDatasetPrope
     psDatasetProperties->adfScale.resize(_nBands);
     psDatasetProperties->abHasScale.resize(_nBands);
 
+    psDatasetProperties->abHasMaskBand.resize(_nBands);
+
     psDatasetProperties->bHasDatasetMask = poFirstBand->GetMaskFlags() == GMF_PER_DATASET;
     if (psDatasetProperties->bHasDatasetMask)
         bHasDatasetMask = TRUE;
@@ -668,6 +675,11 @@ int VRTBuilder::AnalyseRaster( GDALDatasetH hDS, DatasetProperty* psDatasetPrope
         psDatasetProperties->adfScale[j] = poBand->GetScale(&bHasScale);
         psDatasetProperties->abHasScale[j] = bHasScale != 0 &&
                             psDatasetProperties->adfScale[j] != 1.0;
+
+        const int nMaskFlags = poBand->GetMaskFlags();
+        psDatasetProperties->abHasMaskBand[j] =
+                    (nMaskFlags != GMF_ALL_VALID && nMaskFlags != GMF_NODATA) ||
+                    poBand->GetColorInterpretation() == GCI_AlphaBand;
     }
 
     if (bFirst)
@@ -1004,6 +1016,12 @@ void VRTBuilder::CreateVRTSeparate(VRTDatasetH hVRTDS)
                 poSimpleSource->SetNoDataValue( psDatasetProperties->adfNoDataValues[0] );
             }
         }
+        else if( bUseSrcMaskBand && psDatasetProperties->abHasMaskBand[0] )
+        {
+            auto poSource = new VRTComplexSource();
+            poSource->SetUseMaskBand(true);
+            poSimpleSource = poSource;
+        }
         else
             poSimpleSource = new VRTSimpleSource();
 
@@ -1173,6 +1191,12 @@ void VRTBuilder::CreateVRTNonSeparate(VRTDatasetH hVRTDS)
                 poSimpleSource = new VRTComplexSource();
                 poSimpleSource->SetNoDataValue( psDatasetProperties->adfNoDataValues[nSelBand] );
             }
+            else if( bUseSrcMaskBand && psDatasetProperties->abHasMaskBand[nSelBand] )
+            {
+                auto poSource = new VRTComplexSource();
+                poSource->SetUseMaskBand(true);
+                poSimpleSource = poSource;
+            }
             else
                 poSimpleSource = new VRTSimpleSource();
             if( pszResampling )
@@ -1204,18 +1228,28 @@ void VRTBuilder::CreateVRTNonSeparate(VRTDatasetH hVRTDS)
         }
         else if (bHasDatasetMask)
         {
-            VRTSimpleSource* poSimpleSource = new VRTSimpleSource();
+            VRTSimpleSource* poSource;
+            if( bUseSrcMaskBand )
+            {
+                auto poComplexSource = new VRTComplexSource();
+                poComplexSource->SetUseMaskBand(true);
+                poSource = poComplexSource;
+            }
+            else
+            {
+                poSource =  new VRTSimpleSource();
+            }
             if( pszResampling )
-                poSimpleSource->SetResampling(pszResampling);
-            poMaskVRTBand->ConfigureSource( poSimpleSource,
+                poSource->SetResampling(pszResampling);
+            poMaskVRTBand->ConfigureSource( poSource,
                                             static_cast<GDALRasterBand*>(GDALGetRasterBand(hSourceDS, 1)),
                                             TRUE,
                                             dfSrcXOff, dfSrcYOff,
                                             dfSrcXSize, dfSrcYSize,
                                             dfDstXOff, dfDstYOff,
                                             dfDstXSize, dfDstYSize );
 
-            poMaskVRTBand->AddSource( poSimpleSource );
+            poMaskVRTBand->AddSource( poSource );
         }
 
         if( bDropRef )
@@ -1583,6 +1617,7 @@ struct GDALBuildVRTOptions
     int nMaxBandNo;
     char* pszResampling;
     char** papszOpenOptions;
+    bool bUseSrcMaskBand;
 
     /*! allow or suppress progress monitor and other non-error output */
     int bQuiet;
@@ -1727,6 +1762,7 @@ GDALDatasetH GDALBuildVRT( const char *pszDest,
                         psOptions->bSeparate, psOptions->bAllowProjectionDifference,
                         psOptions->bAddAlpha, psOptions->bHideNoData, psOptions->nSubdataset,
                         psOptions->pszSrcNoData, psOptions->pszVRTNoData,
+                        psOptions->bUseSrcMaskBand,
                         psOptions->pszOutputSRS, psOptions->pszResampling,
                         psOptions->papszOpenOptions);
 
@@ -1795,6 +1831,7 @@ GDALBuildVRTOptions *GDALBuildVRTOptionsNew(char** papszArgv,
     psOptions->bQuiet = TRUE;
     psOptions->pfnProgress = GDALDummyProgress;
     psOptions->pProgressData = nullptr;
+    psOptions->bUseSrcMaskBand = true;
 
 /* -------------------------------------------------------------------- */
 /*      Parse arguments.                                                */
@@ -1971,6 +2008,10 @@ GDALBuildVRTOptions *GDALBuildVRTOptionsNew(char** papszArgv,
                     CSLAddString( psOptions->papszOpenOptions,
                                   papszArgv[++iArg] );
         }
+        else if( EQUAL(papszArgv[iArg],"-ignore_srcmaskband") )
+        {
+            psOptions->bUseSrcMaskBand = false;
+        }
         else if( papszArgv[iArg][0] == '-' )
         {
             CPLError(CE_Failure, CPLE_NotSupported,

gdal/data/gdalvrt.xsd

@@ -389,7 +389,8 @@
                 <xs:element name="SrcMax" type="xs:double"/>
                 <xs:element name="DstMin" type="xs:double"/>
                 <xs:element name="DstMax" type="xs:double"/>
-                <xs:element name="NODATA" type="DoubleOrNanType"/>
+                <xs:element name="NODATA" type="DoubleOrNanType"/> <!-- NODATA and UseMaskBand are mutually exclusive -->
+                <xs:element name="UseMaskBand" type="xs:boolean"/> <!-- NODATA and UseMaskBand are mutually exclusive -->
                 <xs:element name="LUT" type="xs:string"/>
             </xs:choice>
         </xs:sequence>

gdal/doc/source/drivers/raster/vrt.rst

@@ -162,7 +162,7 @@ The allowed subelements for VRTRasterBand are :
 
   <ColorInterp>Gray</ColorInterp>:
 
-- **NoDataValue**: If this element exists a raster band has a nodata value associated with, of the value given as data in the element.
+- **NoDataValue**: If this element exists a raster band has a nodata value associated with, of the value given as data in the element. This must not be confused with the NODATA element of a VRTComplexSource element.
 
 .. code-block:: xml
 
@@ -271,7 +271,7 @@ The allowed subelements for VRTRasterBand are :
 
 - **AveragedSource**: The AveragedSource is derived from the SimpleSource and shares the same properties except that it uses an averaging resampling instead of a nearest neighbour algorithm as in SimpleSource, when the size of the destination rectangle is not the same as the size of the source rectangle. Note: a more general mechanism to specify resampling algorithms can be used. See above paragraph about the 'resampling' attribute.
 
-- **ComplexSource**: The ComplexSource_ is derived from the SimpleSource (so it shares the SourceFilename, SourceBand, SrcRect and DestRect elements), but it provides support to rescale and offset the range of the source values. Certain regions of the source can be masked by specifying the NODATA value.
+- **ComplexSource**: The ComplexSource_ is derived from the SimpleSource (so it shares the SourceFilename, SourceBand, SrcRect and DestRect elements), but it provides support to rescale and offset the range of the source values. Certain regions of the source can be masked by specifying the NODATA value, or starting with GDAL 3.3, with the <UseMaskBand>true</UseMaskBand> element.
 
 - **KernelFilteredSource**: The KernelFilteredSource_ is a pixel source derived from the Simple Source (so it shares the SourceFilename, SourceBand, SrcRect and DestRect elements, but it also passes the data through a simple filtering kernel specified with the Kernel element.
 
@@ -375,7 +375,10 @@ the blue band or 4 for the alpha band.
 When transforming the source values the operations are executed
 in the following order:
 
-- Nodata masking
+- Masking, if the NODATA element is set or, starting with GDAL 3.3,
+  if the UseMaskBand is set to true and the source band has a mask band.
+  Note that this is binary masking only, so no alpha blending is done if the
+  mask band is actually an alpha band with non-0 or non-255 values.
 - Color table expansion
 - For linear scaling, applying the scale ratio, then scale offset
 - For non-linear scaling, apply (DstMax-DstMin) * pow( (SrcValue-SrcMin) / (SrcMax-SrcMin), Exponent) + DstMin
@@ -390,7 +393,7 @@ in the following order:
       <ScaleRatio>1</ScaleRatio>
       <ColorTableComponent>1</ColorTableComponent>
       <LUT>0:0,2345.12:64,56789.5:128,2364753.02:255</LUT>
-      <NODATA>0</NODATA>
+      <NODATA>0</NODATA>  <!-- if the mask is a mask or alpha band, use <UseMaskBand>true</UseMaskBand> -->
       <SrcRect xOff="0" yOff="0" xSize="512" ySize="512"/>
       <DstRect xOff="0" yOff="0" xSize="512" ySize="512"/>
     </ComplexSource>