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);
}

Compile the function

gcc -fPIC -c TestFunction.c
gcc -shared TestFunction.o -o TestFunction.so

Create a directory for the GDAL plugins (if not already existing), and copy the compiled shared library for the pixel function there:

mkdir /usr/local/share/gdal/gdalplugins
cp TestFunction.so /usr/local/share/gdal/gdalplugins/gdal_TestFunction.so

Reference the GDAL plugin directory for use in MapServer (with Apache)

Add a line to your Apache config file httpd.conf or apache2.conf like

   SetEnv GDAL_DRIVER_PATH /usr/local/share/gdal/gdalplugins

Restart the Apache server and the VRT images with the pixel functions should be read without problems from MapServer. At least for simple raster algebra I have not noticed performance degradation due to the derived bands.

Important to note that the pixel function works also for the image overviews (at least using the standard internal and external .ovr files, ERDAS .aux types have not been tested with this).

Use of VRT derived bands with PHP MapScript

The above mentioned