GDAL VRT – Creating RGB Composite from Two GeoTIFF Files

gdalvrt

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!

Best Answer

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>
Related Question