Using GDAL VRT images with derived bands in MapServer

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

Definition of the VRT image

The VRT band definition starts with

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

The important part is the reference of the pixel function, here named 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;
}

/* Register the function */
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

More information for GDAL_DRIVER_PATH ​see here.

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 experienced noticeable 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 overviews or external .ovr files, ERDAS .aux types have not been tested with this).

Use of VRT derived bands with PHP MapScript

The above mentioned configuration works well using MapServer CGI. If you want to use them via PHP MapScript the GDAL plugin directory is not identified. In order to get this working add the respective GDAL configuration variable in the PHP MapScript code that calls the $map->draw() function or that calls the WxS/OWS functions that returns the WxS/OWS output, like

putenv("GDAL_DRIVER_PATH=/usr/local/share/gdal/gdalplugins");
$_ENV['GDAL_DRIVER_PATH'] = "/usr/local/share/gdal/gdalplugins";

(one of them might not be necessary...)

Restart Apache and everything should work.

For other MapScript versions this additional setting of the GDAL environment variable might be necessary as well.