GDAL VRT - How to create RGB composite from two input GeoTIFF files (input1, input2, input1/input2)

Question

I'm trying to create an RGB composite from two input GeoTIFF files. The first two bands in the VRT are straightforward and look like the following:

<VRTDataset rasterXSize="10980" rasterYSize="10980">
  <SRS dataAxisToSRSAxisMapping="1,2">...</SRS>
  <GeoTransform>...</GeoTransform>
  <VRTRasterBand dataType="Float32" band="1">
    <NoDataValue>nan</NoDataValue>
    <ColorInterp>Red</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">input1.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="10980" RasterYSize="10980" DataType="Float32" BlockXSize="512" BlockYSize="512" />
      <SrcRect xOff="0" yOff="0" xSize="10980" ySize="10980" />
      <DstRect xOff="0" yOff="0" xSize="10980" ySize="10980" />
      <NODATA>nan</NODATA>
    </ComplexSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="2">
    <NoDataValue>nan</NoDataValue>
    <ColorInterp>Green</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">input2.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="10980" RasterYSize="10980" DataType="Float32" BlockXSize="512" BlockYSize="512" />
      <SrcRect xOff="0" yOff="0" xSize="10980" ySize="10980" />
      <DstRect xOff="0" yOff="0" xSize="10980" ySize="10980" />
      <NODATA>nan</NODATA>
    </ComplexSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="3">
    <NoDataValue>nan</NoDataValue>
    <ColorInterp>Blue</ColorInterp>
    ...?

  <OverviewList resampling="nearest">2 4 8 16 32</OverviewList>
</VRTDataset>

I'm not sure how to best implement the third VRTRasterband, though. As mentioned in the title, this third band is supposed to be a ratio of the first two bands, so I wanted to use one of the existing pixel functions. However, as far as I can see there is no pixel function that divides two input sources. I thought about a possible work-around by using the mul pixel function ("multiply 2 or more raster bands") in combination with the inv pixel function ("inverse (1./x)") to create the ratio.

band3 = mul(input1, inv(input2))

Is it somehow possible to combine pixel functions in such a way in one VRTRasterband? I wasn't able to find anything related in the documentation or in Github issues. Does anyone have another solution that I might have missed in the documentation?

It's also important to note, that I need to solve this without writing a new pixel function and without the option to create a pixel function in Python!

Answer 1

In a future version of GDAL, the div function should be available as a default pixel function as there's work going on to add new and improve existing functions.

In the meantime, make a separate VRT for the inv input2 then use input1 and inverse_input2.vrt as the bands for the mul.

To have it all in one VRT file, you can nest VRT XML in the <SourceFilename etc...> element using a CDATA section which is content that is to be interpreted as text data, not as XML.

E.g. <![CDATA[vrt xml string]]>.

This trick works as GDAL understands how to open a Dataset from a string of VRT XML as well as a file.

E.g. (tested and works with dummy data, obviously you'll need to rewrite it to suit your data, rows, cols, CRS etc...):

<VRTDataset rasterXSize="256" rasterYSize="256">
    <VRTRasterBand dataType="Float32" band="1">
        <Description>input 1</Description>
        <SimpleSource>
            <SourceFilename relativeToVRT="1">input1.tif</SourceFilename>
            <SourceBand>1</SourceBand>
            <SourceProperties RasterXSize="256" RasterYSize="256" DataType="Float32" />
            <SrcRect xOff="0" yOff="0" xSize="256" ySize="256" />
            <DstRect xOff="0" yOff="0" xSize="256" ySize="256" />
        </SimpleSource>
    </VRTRasterBand>
    <VRTRasterBand dataType="Float32" band="2">
        <Description>input 2</Description>
        <SimpleSource>
            <SourceFilename relativeToVRT="1">input2.tif</SourceFilename>
            <SourceBand>1</SourceBand>
            <SourceProperties RasterXSize="256" RasterYSize="256" DataType="Float32" />
            <SrcRect xOff="0" yOff="0" xSize="256" ySize="256" />
            <DstRect xOff="0" yOff="0" xSize="256" ySize="256" />
        </SimpleSource>
    </VRTRasterBand>
    <VRTRasterBand dataType="Float32" band="3" subClass="VRTDerivedRasterBand">
        <Description>multiply input 1 and inverse of input 2</Description>
        <PixelFunctionType>mul</PixelFunctionType>
        <SimpleSource>
            <SourceFilename relativeToVRT="1">input1.tif</SourceFilename>
            <SourceBand>1</SourceBand>
            <SrcRect xOff="0" yOff="0" xSize="256" ySize="256"/>
            <DstRect xOff="0" yOff="0" xSize="256" ySize="256"/>
        </SimpleSource>
        <SimpleSource>
            <SourceFilename relativeToVRT="1"><![CDATA[
                <VRTDataset rasterXSize="256" rasterYSize="256">
                    <VRTRasterBand dataType="Float32" band="1" subClass="VRTDerivedRasterBand">
                        <Description>inverse</Description>
                        <PixelFunctionType>inv</PixelFunctionType>
                        <SimpleSource>
                          <SourceFilename relativeToVRT="1">input2.tif</SourceFilename>
                          <SourceBand>1</SourceBand>
                          <SourceProperties RasterXSize="256" RasterYSize="256" DataType="Float32" />
                          <SrcRect xOff="0" yOff="0" xSize="256" ySize="256" />
                          <DstRect xOff="0" yOff="0" xSize="256" ySize="256" />
                        </SimpleSource>
                    </VRTRasterBand>
                </VRTDataset>
            ]]></SourceFilename>
            <SourceBand>1</SourceBand>
            <SrcRect xOff="0" yOff="0" xSize="256" ySize="256"/>
            <DstRect xOff="0" yOff="0" xSize="256" ySize="256"/>
        </SimpleSource>
    </VRTRasterBand>
</VRTDataset>