Using GDAL VRT images with derived bands in MapServer

The GDAL virtual image format (VRT) allows the creation of derived bands that derive their pixel information from their source bands. A general overview is given at the ​GDAL web site. This how-to describes the steps to make this type of raster format usable in MapServer, based on the sample code of the GDAL site. The emphasis lies on the creation and registration of the required pixel function used for the derived band.

Definition of the VRT image

The VRT band definition starts with

<VRTDataset rasterXSize="1000" rasterYSize="1000">
    <VRTRasterBand dataType="Float32" band="1" subClass="VRTDerivedRasterBand">
        <Description>Magnitude</Description>
        <PixelFunctionType>TestFunction</PixelFunctionType>
        <SourceTransferType>CFloat64</SourceTransferType>
    ...

The important part is the refernce of the pixel function, here name TestFunction

Creation of the pixel function

Create a file TestFunction.c with the following contents

#include "gdal.h"

CPLErr TestFunction(void **papoSources, int nSources, void *pData,
                    int nXSize, int nYSize,
                    GDALDataType eSrcType, GDALDataType eBufType,
                    int nPixelSpace, int nLineSpace)
{
    int ii, iLine, iCol;
    double pix_val;
    double x0, x3, x4, x8;

    // ---- Init ----
    if (nSources != 4) return CE_Failure;

    // ---- Set pixels ----
    for( iLine = 0; iLine < nYSize; iLine++ )
    {
        for( iCol = 0; iCol < nXSize; iCol++ )
        {
            ii = iLine * nXSize + iCol;
            /* Source raster pixels may be obtained with SRCVAL macro */
            x0 = SRCVAL(papoSources[0], eSrcType, ii);
            x3 = SRCVAL(papoSources[1], eSrcType, ii);
            x4 = SRCVAL(papoSources[2], eSrcType, ii);
            x8 = SRCVAL(papoSources[3], eSrcType, ii);

            pix_val = sqrt((x3*x3+x4*x4)/(x0*x8));
            
            GDALCopyWords(&pix_val, GDT_Float64, 0,
                          ((GByte *)pData) + nLineSpace * iLine + iCol * nPixelSpace,
                          eBufType, nPixelSpace, 1);
        }
    }

    // ---- Return success ----
    return CE_None;
}

void GDALRegister_TestFunction() {
     GDALAddDerivedBandPixelFunc("TestFunction", TestFunction);
}